The collision and continued convergence of the Indian and Eurasian plates since ca. 55 Ma have resulted in substantial deformation and uplift within the interior of the Tibetan Plateau, together with the outward expansion of the plateau along its margins (e.g., Molnar and Tapponnier, 1975; Tapponnier et al., 2001; Yin, 2000; Yin, 2010; ; Zhao et al., 2010; Zuza et al., 2018; Xiong et al., 2022). This collision has significantly impacted both the topography and climate of the surrounding regions (Zhang et al., 2004). However, the precise manner in which the Tibetan crust accommodates the convergence of the Indian and Eurasian plates, leading to the development of the Tibetan Plateau, remains a subject of intense debate (e.g., Cheng et al., 2019; Clark and Royden, 2000; Yin and Harrison, 2000; Tapponnier et al., 2001; van Hinsbergen et al., 2011; van Hinsbergen et al., 2012; Yakovlev and Clark, 2014; Ingalls and Stedman, 2016; Zuza et al., 2016; Zheng et al., 2017). Moreover, different segments of the Tibetan Plateau margins may exhibit varying deformation mechanisms, particularly across the northeastern Tibetan Plateau margin that stretches across strike for >1,000 km. (Lin et al., 2011; Duvall et al., 2013).

NQLS and the HXBS comprise the northernmost margin of the Tibetan Plateau. This contractional thrust-fold belt and foreland basin system represent an area of concentrated stress resulting from the convergence of the Indian and Eurasian plates (Xiong et al., 2019). Understanding the geometric and kinematic characteristics of the crustal structure in this region is essential for deciphering the deep interrelation between the Tibetan Plateau and the Eurasian plate. This knowledge holds the potential to significantly enhance our comprehension of continental tectonics (Yin et al., 2002; Molnar, 1988; Zhang et al., 2004; Zhou et al., 2006; Yin, 2010; Lease et al., 2012; Craddock et al., 2011; Craddock et al., 2014; Zuza et al., 2019).

Since the 1970s (Burke and Dewey, 1973; Molnar and Tapponnier, 1975), numerous end-member models have been proposed to explain the formation of the Tibetan Plateau, as a result of intensive study (Figure 1). In the northeastern Tibetan Plateau, these models include (1) Cenozoic underthrusting of the Asian lithosphere beneath the Tibetan Plateau (Willett and Beaumont, 1994; Kind et al., 2002; Zhao et al., 2011; Feng et al., 2014; Ye et al., 2015); (2) coupled distributed crustal shortening and underthrusting of the Alxa Block beneath Tibet, accounting for the observed crustal shortening and thickness distribution (England and Houseman, 1986; Dewey et al., 1997; Huang et al., 2018; Xiong et al., 2019; Zuza et al., 2019); (3) vertical inflation of the Tibetan crust due to lateral mid-lower crustal channel flow (Zhao and Morgan, 1987; Bird, 1991; Royden et al., 1997; Clark and Royden, 2000); and (4) intracontinental subduction facilitated by large-scale strike–slip faults (Zuza et al., 2016; Wu et al., 2023b).

These tectonic models offer specific predictions regarding the location, timing of deformation, and kinematics of crustal structures in the northeastern Tibetan Plateau (Lease et al., 2012; He, 2020; Zuza et al., 2016; 2019). Therefore, crustal-scale seismic imaging can provide insights into the dominant deformation mechanism along the northern margin of the Tibetan Plateau. Despite prolific geophysical research conducted over recent decades on the northeastern Tibetan Plateau, the mechanisms driving crustal deformation remain contentious. Given that the west-to-east-trending NQLS and HXBS stretch over a span of 1,000 km, the crustal deformation mechanism may indeed vary across different segments.

In the western segment, the Moho depth is notably deeper than that in the middle-eastern Qilian Shan fold-thrust belt, with a maximum depth of 73 km beneath Hala Lake, as substantiated by studies (Cui et al., 1995; Gao et al., 1999; Jolivet et al., 2001; Huang et al., 2021). This observation supports the hypothesis that the most concentrated convergent stress is accommodated through crustal deformation in the western segment. To address this question, we gathered seismic data from a linear array spanning 130 km, extending northward from NQLS, traversing the entire Jiuquan Basin, and reaching the Huahai Basin. The data were continuously recorded over approximately 30 days using the Smart Solo IGU-16 HR 3C geophones, which are three-component geophones with a dominant frequency of 5 Hz (Figure 2B). The dense seismic array comprises 261 short-period geophones spaced at 500-m intervals. Through receiver function processing of the seismic observations, we have unveiled the crustal structure beneath the North Qilian–Huahai Basin and gained insights into the crustal deformation mechanism within the western segment of NQLS–HXBS.

The northeastern Tibetan Plateau has an average elevation of ∼4.5 km, with the high topography decreasing rapidly to <1.5 km in HXBS to the north (Figure 2). Previous geophysical observations indicate that the Moho depth ranges from 55 to 73 km (Zhao et al., 2001; Yue et al., 2012; Gao et al., 2013; Yuan et al., 2013; Xu et al., 2018; Wu et al., 2023b). Additionally, the relatively stable crust beneath HXBS has a thickness of approximately 42–51 km (Cui et al., 1995; Xiong et al., 2019).

The NQLS formed as an accretionary complex due to the closure of the Qilian Ocean in the Early Paleozoic and was later reactivated following the Indian–Eurasian plate collision in the Cenozoic. The present-day tectonic belt, located between the northern Alxa Block–Beishan Block and the southern Middle Qilian Block, follows a narrow NW–SE-trending structure. This belt is truncated by the Altyn Tagh faults to the northwest (Figure 2). The Qilian Shan fold-thrust belt exhibits volcanic and magmatic rock distributions from the Paleoproterozoic to the early Paleozoic (Wu et al., 2023a). Since the Cenozoic era, NQLS has been involved in crustal thickening and deformation, marking its northernmost extension onto the Tibetan Plateau.

The Jiuquan Basin is situated at the western terminus of HXBS, and the Jiuquan Basin is bounded by the Jiayuguan fault, Altyn Tagh fault, Qilian Shan North Margin fault, and Heishan fault. Sedimentation in the Jiuquan Basin initiated in the Mesozoic, and this northwest-trending basin was subsequently influenced by the uplift of NQLS. Deep seismic reflection profiles reveal that the thickness of the Cenozoic sedimentary cover ranges between 1.5 and 2 km (Song et al., 2001; Huang et al., 2018).

The Huahai Basin lies to the north of the Jiuquan Basin bounded by the Kuantan Shan fault, Heishan fault, and Jintana Shan fault. Initially part of the Dunhuang Block, a relatively stable Precambrian block, the Huahai Basin is covered by Mesozoic–Cenozoic sediments with a sedimentary thickness of 5 km (Chen and Yang, 2010). Paleozoic granitoids have locally intruded into the sedimentary cover (refer to Figure 2).

In this study, we deployed a dense array of 261 Smart Solo 16 HR 3C nodal short-period seismometers for observations. These were positioned with an average spacing of 500 m in the field and operated for approximately 30 days. A receiver function is a time sequence derived from deconvolving the radial component of teleseismic P-waves. This process yields the vertical component that represents P-to-S wave conversion. Receiver function analysis is an effective approach for imaging crustal structures. By stacking multiple receiver functions with similar ray paths, the dense nodal array enhances the vertical resolution of receiver function images and boosts high-frequency signals (Langston, 1979; Xu et al., 2018; Ammon et al., 1990; Xu et al., 2020; Tian et al., 2021).

After obtaining the raw seismic data, several steps are typically followed to extract receiver functions. This process involves filtering, receiver function separation, and deconvolution. In practical research, receiver function studies usually focus on earthquakes with a magnitude of 5 or higher. For data preprocessing, earthquakes within 30°–90° epicentral distances (refer to Figure 3) were selected. We acquired a total of 47 seismic events with a magnitude (Ms) of ≥5 and epicentral distances between 30° and 90° from the IRIS website during the data acquisition period. Data for each station were extracted from 900 s after the origin time. Raw seismic data at stations contain not only information about the subsurface medium structure but also various types of noise from surface human activities, unrelated to underground structures and tectonics. Hence, noise filtering is essential, often using low-pass and band-pass filters. After removing the instrument response, the data were down-sampled from 500 Hz to 50 Hz (0.002 s–0.02 s), and band-pass filtering was conducted within the frequency range of 0.5 Hz–10 Hz. The raw three-component data coordinate system was rotated from NEZ to RTZ. Ray parameters and theoretical arrival times were computed. After rotation, P-wave energy concentrated mainly in the radial R component, while azimuthal energy in the tangential T component. Finally, time-domain deconvolution algorithm calculations were performed (Ligorria and Ammon, 1999), using a Gaussian filter factor 2. This iterative deconvolution subtracts the radial component of observed seismic records from calculated theoretical receiver functions, iteratively updating differences until either iterations exceed a threshold or differences fall within a given tolerance range (refer to Figure 4). Computed receiver functions were then selected, focusing on data with clear P-wave phases, distinct seismic phases, and good Ps–PpPs consistency.

The common midpoint (CMP) stacking concept, initially introduced by Dueker and Sheehan (1997) for reflection waves, was adapted for imaging converted waves in receiver functions. This method corrects and stacks slanting receiver functions dynamically to create a wavefield image of the interface below the array. Building on this, Zhu and Kanamori (2000) proposed the CCP stacking method for receiver functions. Based on a given velocity model, each receiver function is projected back into the spatial domain using ray paths, followed by stacking imaging in the spatial domain (Chen et al., 2022; Cheng et al., 2023). In our study, we applied the CCP stacking method to stack the processed receiver functions. The method assumes a horizontally layered medium beneath the station, mapping amplitude information of processed and deconvolved receiver functions along the ray path to the conversion point. The amplitudes of receiver functions at the conversion point are stacked to obtain the medium structure below the station (refer to Figure 4).

For time–depth conversion, we utilized the reference velocity model IASP91, with a rectangular shape and a maximum depth of 100 km (Chen et al., 2022). Bins were spaced at intervals of 5 km. The time-to-depth conversion used a depth interval of 1 km, with intervals of 2 km from 0 km to 100 km. After stacking receiver functions with these settings, the CCP profile image was generated (refer to Figure 5).

In the CCP image, red colors indicate positive amplitudes (increased velocity downward), while blue colors indicate negative amplitudes (decreased velocity downward). Our new seismic imaging reveals the following features.

A continuous positive amplitude stretches throughout the section at depths of 0–5 km; it deepens to 8 km in the contact zone of North Qilian and Jiuquan Basin and the center of the Huahai Basin; two distinct large negative amplitude bodies are visible at depths of 8–20 km below the NQLS and Jiuquan Basin; a relatively continuous positive amplitude extends at a depth of 12–20 km; the Moho, the most prominent positive velocity discontinuity, deepens from 45–50 km beneath the Huahai Basin to 55–60 km beneath NQLS; above the bottom Moho, a distinct positive amplitude at a depth of 45–50 km beneath NQLS, within 0–50 km distance, may represent the original Moho of the marginal NQLS (Figure 5).

Drawing on our receiver function image and the analysis presented previously, we have outlined a comprehensive tectonic representation on a crustal scale that encapsulates our observations and interpretations beneath the North Qilian–Huahai Basin (Figure 9). By amalgamating our seismic imaging with additional geological and geophysical data, we propose the existence of a potential decoupling plane situated at a depth of 12–20 km. Within this framework, south-dipping brittle thrust faults govern upper crustal deformation, while imbricate structures and duplexes operate in the lower crust. Notably, we conjecture that the underthrusting of the Beishan Block beneath NQLS might be accountable for the occurrence of a double Moho. In addition, we posit that a Moho ramp could be developing beneath the Jiuquan Basin to accommodate the pronounced compression. In the western stretch of NQLS and HXBS, we propose that coupled distributed shortening and the underthrusting of the Beishan Block beneath the Qilian Shan fold-thrust belt might collaboratively contribute to the intricate crustal deformation.