Quaternary Crustal Shortening of the Houyanshan Structure in the Eastern Chinese Tian Shan: Constrained from Geological and Geomorphological Analyses

: The Tian Shan is one of the most active intracontinental orogenic belts in the world. It has undergone complex deformation that has resulted in the formation of several fold-and-thrust belts (FTBs) in the piedmonts and intermontane basins. Investigating the deformation histories of these FTBs is important for understanding the tectonic propagation processes of the Tian Shan. Here, we gain insight into these crustal shortening processes by deciphering the Houyanshan structure, a typical fold-thrust belt in the eastern Chinese Tian Shan. We ﬁrst describe a curved thrust ramp and related fold pairs of the structure using high-resolution remote sensing photography, deformation of ﬂuvial terraces, and ﬁeld-based geological cross-section. Combined with deformed terrace records and optically stimulated luminescence (OSL) dating results, the kinematic style allows us to yield a geologic shortening rate of 1.6 ± 0.2 mm/a since ~52 ka. Second, to reduce uncertainty in the seismic interpretation and quantify the amount and time of crustal shortening, we interpret three seismic reﬂection proﬁles by using the theory of quantitative fault-related fold, area-depth-strain (ADS), and reverse modeling analyses. These proﬁles provide direct evidence that this structure connects by means of a listric thrust ramp to a shallow detachment level. ADS analysis reveals that the maximum shortening of the Huoyanshan structure is ~4.5 km, which is consistent with the result of quantitative inverse modeling. Each of the structural analysis methods gives similar parameters, and the high consistency of results greatly improves the soundness of a given geologic interpretation. Finally, the shortening rate and total shortening amount suggest that the structure may have formed at 1.8–3.7 Ma, which is nearly synchronous around the Tibetan Plateau. Together, these results indicate that this combined geological and geomorphological analysis provides greater insight into deformation information than can be achieved by any individual technique in studying fold-and-thrust belts worldwide.


The Tian Shan
The Tian Shan, a mountain belt more than 2500 km long, dominates the topography of central Asia (Figure 1a). From west to east, it extends through the Republics of Uzbekistan, Kazakhstan, Tajikistan, Kyrgyzstan, and China, generally rising to approximately 4000 m in altitude and with peaks exceeding 7000 m (Figure 1b). The Chinese Tien Shan Mountains are divided into two segments along 88 • E longitude: the eastern Chinese Tian Shan and the western Chinese Tian Shan. The western Chinese Tian Shan is further divided into the southern Tian Shan, central Tian Shan, and northern Tian Shan [63]. The Tian Shan has undergone complex deformation. The overall range was formed in the late Paleozoic [63][64][65][66][67][68], and was reactivated and intensely uplifted during the Oligocene and Miocene in response to the collision between India and Asia [42,[69][70][71][72]. As a consequence of regional S-N trending compression active since the late Cenozoic, deformation propagated basinward and several FTBs composed of subparallel structural belts formed in the range fronts. Different places have different stages of tectonic development. At the northern Tian Shan piedmont, at least three E-W trending structural belts have been formed; in Kuqa, there are four structural belts; Kashi has three structural belts; and in Turfan, there is mainly one structural belt. In general, the present-day Tian Shan topography, composed of a series of ranges with elevations >4000 m separated by large E-W-striking intermontane basins that are bounded by several thrust systems, reflects both its Paleozoic basement architecture and its Cenozoic reactivation history [73][74][75]. According to the results of GPS measurements, present-day geodetic shortening rates decrease from 20 mm/a at~75 • longitude [30,31], to less than 10 mm/a at 84 • longitude [76], and then to less than 5 mm/a at 87 • longitude [77].
The exact timing of the reactivation of the Tian Shan during the late Cenozoic remains controversial. On the basis of the unconformity at the base of the Oligocene conglomerates in the eastern Tarim and Turfan basins, Allen et al. [63] and Windley et al. [65] suggested that the late Cenozoic uplift of the Tian Shan was induced by the Indian-Asia collision. According to Avouac et al. [35], the uplift of the Tian Shan began in the early to middle Miocene (16 + 22/−9 Ma). Based on sediment flux acceleration detected from changes in the sedimentary logs of drill holes and cross sections from the Tarim and Jungar basins, Metivier and Gaudemer [78] found that sediment flux acceleration began at~17 Ma and enhanced again at~5 Ma. Abdrakhmatov et al. [30] estimated the shortening and shortening rate at 76 • E and concluded that the deformation initiated at~10 Ma. However, these pioneering studies contained a degree of uncertainty due to the absence of geologic cross sections across the range and precise chronological constraints. By utilizing low-temperature thermochronology, the uplift history of the Tian Shan has been better constrained. According to fission-track analyses, the unroofing commenced at 25-24 Ma at both the northern and southern margins of the Tian Shan [42,68]. Sobel and Dumitru [71] and Sobel et al. [72] studied sections of the northwestern margin of the Tarim Basin and indicated that exhumation started at~25-18 Ma, and pulsed events took place along the Remote Sens. 2023, 15, 1603 4 of 25 southern margin of the Tian Shan. Similarly, magnetostratigraphy also provides chronologic analyses of the Cenozoic deformation of the Tian Shan. Bullen et al. [44,79] developed magnetostratigraphy and detrital fission-track studies on the Kyrgyz Tian Shan and discovered pulsed deformation that occurred at~11 Ma and~3 Ma. Magnetostratigraphy of late Cenozoic sediments in the Kuqa Depressions, southern Tian Shan piedmont, revealed exhumation ages of 24-20 Ma, and deformation episodes around 16-15 Ma and 7 Ma were also found [41,80]. Based on the magnetostratigraphy of a section in the northern Tian Shan piedmont, Ji et al. [81] concluded that increases in the sedimentation rate at~26-22 Ma marked the onset of uplift [81]. However, based on increased sedimentation rates, changes in magnetic parameters and thermochronology along the northern and southern Tian Shan range fronts, Charreau et al. [82,83] implied that rock uplift did not begin until~11 Ma, which is a much younger age. Sun et al. [50,84] commented on a magnetostratigraphy study in the northern Tian Shan piedmont and concluded that the mountain building started at 7 Ma and lasted to~2. 6 Ma. The controversy in these conclusions has demonstrated the ongoing challenges associated with assessing the exact initial uplift. Currently, the deformation study of the Tian Shan is focused on several subparallel folds within the fold-and-thrust belt of the northern and southern margins, which documented the most recent records of the growth of the Tian Shan from the Miocene to the Holocene. Multiple research methods, including morphological and chronologic investigations of deformed geomorphological surfaces, structural analysis, and low-temperature thermochronology as well as magnetostratigraphy of terrestrial sediments of the fold in foreland basins, have been employed to constrain the deformation process [32,39,47,49,[85][86][87][88][89][90]. According to fission-track analyses, the unroofing commenced at 25-24 Ma at both the northern and southern margins of the Tian Shan [42,68]. Sobel and Dumitru [71] and Sobel et al. [72] studied sections of the northwestern margin of the Tarim Basin and indicated that exhumation started at ~25-18 Ma, and pulsed events took place along the southern margin of the Tian Shan. Similarly, magnetostratigraphy also provides chronologic analyses of the Cenozoic deformation of the Tian Shan. Bullen et al. [44,79] [41,80]. Based on the magnetostratigraphy of a section in the northern Tian Shan piedmont, Ji et al. [81] concluded that increases in the sedimentation rate at ~26-22 Ma marked the onset of uplift [81]. However, based on increased sedimentation rates, changes in magnetic parameters and thermochronology along the northern and southern Tian Shan range fronts, Charreau et al. [82,83] implied that rock uplift did not begin until ~11 Ma, which is a much younger age. Sun et al. [50,84] commented on a magnetostratigraphy study in the northern Tian Shan piedmont and concluded that the mountain building started at 7 Ma and lasted to ~2.6 Ma. The controversy in these conclusions has demonstrated the ongoing challenges associated with assessing the exact initial uplift. Currently, the deformation study of the Tian Shan is focused on several subparallel folds within the fold-and-thrust belt of the northern and southern margins, which documented the most recent records of the growth of the Tian Shan from the Miocene to the Holocene. Multiple research methods, including morphological and chronologic investigations of deformed geomorphological surfaces, structural analysis, and low-temperature thermochronology as well as magnetostratigraphy of terrestrial sediments of the fold in foreland basins, have been employed to constrain the deformation process [32,39,47,49,[85][86][87][88][89][90].

The Turpan Basin and the Huoyanshan Anticline
The Turpan Basin (Figure 1c), an intermontane basin located in the eastern Chinese Tian Shan, is considered one of the most hydrocarbon-prone areas in the world [91,92]. It extends approximately 500 km from east to west and 60-100 km from south to north and covers an area of approximately 53,000 km 2 . A minimum altitude of 154 m below sea level is located southwest of the Turfan Basin, known as Aiding Lake, which is the second lowest exposed land surface on Earth. Boundaries of the basin are defined by the Bogda Mountains to the north and the Juelotage Mountains to the south. To the north of the Turpan Basin, seismicity, active deformation slip estimates, and evidence of surface faulting and folding all demonstrate active deformation. This range has a history of infrequent but significant contemporary earthquakes (Figure 1c). The rate of shortening at the northern piedmont of the Bogda Mountains was calculated to be 1 mm/a [58]. At the southern margin of the Bogda Mountains, overthrusting of the Turpan Basin manifests as E-W alignments of active thrust-related folds, where a W-N-W tending FTB can be observed (Figure 1c). This foreland FTB includes (from west to east) the Hongshankou (HSK), Yanshan (YS), Kendeke (KDK), Huoyanshan (HYS), and Qiketai anticlines (QKT). These anticlines generally formed above the detachment in the early Jurassic sediments [37] and are thought to be a long-distance manifestation of the collision of India and Asia [63]. In contrast, Quaternary deformation in the southern margin of the basin is weak, and crustal shortening is not prominent in this range [57,60,93,94].
The Huoyanshan anticline, 90 km long west-east, and 10 km wide north-south, is the most active structure in the basin (Figures 1c and 2). The outcropped Jurassic, Cretaceous, Palaeogene, and Neogene deposits as well as the overlying late Quaternary alluvial terraces all have been warped by folding. The thrust fault lies under the Huoyanshan anticline outcrops at the south margin of the anticline, and the dip angle of the thrust is~30-50 • [59,95]. Based on deformed geomorphological features along the fold, Yang et al. [95] suggested that the Huoyanshan structure accommodates 2.0-3.2 mm/a of crustal shortening. Based on the surface deformation and bedding attitudes, Yang et al. [95] also established a kinematic model and indicated that the Huoyanshan structure is a fault-propagation fold over a listric thrust rooting into a flat detachment layer at depth. However, due to the lack of subsurface structure deformation, it is difficult to further retrieve the complete history of fold growth.

Stratigraphy
In the Turpan Basin, the sediment fill contains more than 7000 m of continental sediments, extending from the late Permian to the Quaternary, and the Huoyanshan anticline mainly exposes Mesozoic and Cenozoic deposits. From the Permian to the Jurassic, the Jueluotage Mountains were the primary provenance area for the basin, as compared to the locally uplifted Bogda Mountains, which provide little source rocks, forming sedimentary wedges with layers thinning southward. Then, with the gradual uplifting and folding of the Bogda Mountains since the latest Jurassic, the Bogda area built another provenance area for the basin. Until the late Cenozoic caused by the India-Asia collision, the intense uplift of the Bogda Mountains made it the dominant sediment source area [96,97], forming sedimentary wedges with layers thinning northward. Currently, the morphology pattern of stratigraphic units resembles a kind of seesaw ( Figure 3).

Geomorphology and Quaternary Chronology
Deformed fluvial terraces preserved along the Shengjingou Valley have recorded localized folding deformation across the Huoyanshan anticline [37,56,95]. In this study, the surficial features of these fluvial terraces have been extensively investigated through mapping and GPS surveying. High-resolution DEMs (~1 m) across the Shengjingou Valley were acquired by using UAV photography to correlate surface geomorphology across the entire transect of the Huoyanshan anticline ( Figure 5). A total of 5112 photos were captured, and 75 ground control points were placed by using a real-time kinematic GPS (RTK-GPS) with centimeter vertical and horizontal precision ( Figure 5) that covered all surveying areas to further improve the precision of the DEM. Details of this procedure can be found in Yang et al. [95]. Based on DEM interpretation and field surveying, we classified and mapped eight generations of terraces across the Shengjingou Valley (T1-T8, from young to old, Figure 6). Then, to obtain the spatial distribution of terraces and riverbeds, profiles along the tops of terrace treads and riverbeds coupled with the coordinates and elevations of survey points ( Figure 7a) were extracted from the DEM. To erase the riverbed gradient from the terrace profiles, elevations above the riverbed of the terraces were recalculated by subtracting the riverbed gradient using a polynomial that closely approximates the current riverbed profile (Figure 7b). To reflect the vertical uplift caused by folding, the undeformed terrace elevation was removed from the maximum terrace elevation ( Figure 8). Table 1 displays our evaluation of the uncertainties associated with these measures.

Geomorphology and Quaternary Chronology
Deformed fluvial terraces preserved along the Shengjingou Valley have recorded localized folding deformation across the Huoyanshan anticline [37,56,95]. In this study, the surficial features of these fluvial terraces have been extensively investigated through mapping and GPS surveying. High-resolution DEMs (~1 m) across the Shengjingou Valley were acquired by using UAV photography to correlate surface geomorphology across the entire transect of the Huoyanshan anticline ( Figure 5). A total of 5112 photos were captured, and 75 ground control points were placed by using a real-time kinematic GPS (RTK-GPS) with centimeter vertical and horizontal precision ( Figure 5) that covered all surveying areas to further improve the precision of the DEM. Details of this procedure can be found in Yang et al. [95]. Based on DEM interpretation and field surveying, we classified and mapped eight generations of terraces across the Shengjingou Valley (T1-T8, from young to old, Figure 6). Then, to obtain the spatial distribution of terraces and riverbeds, profiles along the tops of terrace treads and riverbeds coupled with the coordinates and elevations of survey points ( Figure 7a) were extracted from the DEM. To erase the riverbed gradient from the terrace profiles, elevations above the riverbed of the terraces were recalculated by subtracting the riverbed gradient using a polynomial that closely approximates the current riverbed profile (Figure 7b). To reflect the vertical uplift caused by folding, the undeformed terrace elevation was removed from the maximum terrace elevation ( Figure 8). Table 1 displays our evaluation of the uncertainties associated with these measures.       Topographic profiles of terraces T1-T8 that removed the modern river gradient. The distance is projected along the line C-C′ (location see Figure 5c), where C is the direction upstream. In both profiles, the vertical exaggeration is 10.  OSL dating of deformed surface sediments is proven suitable for constraining the chronology, yielding reliable absolute ages for terrace abandonment through the Holocene, under conditions of sufficient exposure of material being dated before deposition. We collected our sample SJK-T4 from a layer of interbedded silts and sands lying ~3 m beneath outwash gravel on a prominent terrace (T4) in the middle part of the Shengjingou Valley (Figure 5a). Two cores were collected by driving 50 cm aluminum tubes into the fluvial deposit, and light shielding materials were capped on the ends. Samples were processed at the Research Laboratory of Luminescence Dating at the Institute of Geology,  OSL dating of deformed surface sediments is proven suitable for constraining the chronology, yielding reliable absolute ages for terrace abandonment through the Holocene, under conditions of sufficient exposure of material being dated before deposition. We collected our sample SJK-T4 from a layer of interbedded silts and sands lying~3 m beneath outwash gravel on a prominent terrace (T4) in the middle part of the Shengjingou Valley (Figure 5a). Two cores were collected by driving 50 cm aluminum tubes into the fluvial deposit, and light shielding materials were capped on the ends. Samples were processed at the Research Laboratory of Luminescence Dating at the Institute of Geology, China Earthquake Administration. Sample preparation and analysis followed the standard procedures of Aitken [100,101]. The equivalent dose (De) was obtained from the Simplealiquot Regenerative Dose (SMAR) protocol of fine particle quartz samples (4~11 µm). The environmental dose rate was determined through activity concentration measurements of the U, Th, and K. The dose rate was computed from the measured activity concentration and attenuation values [100]. Both the cosmic ray contribution and the sample water content were used to assess the environmental dose [102]. The simple and weighted means and standard errors are at a 68% confidence level (Table 1).

Structural Interpretation from Outcrops and Seismic Profiles
The subsurface geometry of the Huoyanshan structure is constrained by the bedrock geology and seismic reflection profiles. Outcrops of the folded strata are exposed by the entrenchment along the Shengjingou Valley. Thus, 21 bedding attitudes were measured across the Shengjingou Valley (Figures 2b and 6). The results reveal a very uniform strike and a continuous and progressive variation in dip (Figure 8), indicating fold growth by limb rotation over listric faults. We acquired 2D pre-stack time migrated seismic reflection data along 3 lines across the eastern part of the Huoyanshan anticline (Figure 2a), which constrains the subsurface geometries to~7 km below sea level in detail (Figures 9-11). These seismic profiles were converted into depth by utilizing a multilayer velocity model [103] in the Move 2013 software (© 2011 Midland Valley Exploration Ltd., 144 West George Street, Glasgow, G2 2HG, UK). Field observations, regional geologic maps, and well data were used to constrain boundaries of horizons in seismic profiles. Several prominent reflections in the seismic sections, including the bottoms of lithological units J 2 x, J 2 q, J 3 k, K 1 s, Esh, N 1 t, N 2 p, and Q 1 x, can be readily traced and correlated, allowing a robust structural interpretation and analysis with little uncertainty. The kinematic model was inferred from the characteristics of surface and subsurface geometries. The crustal shortening was calculated using the trigonometric relationship in the kinematic model (Figure 8b), and when paired with the terrace abandonment age, the shortening rate was computed (Table 1).
China Earthquake Administration. Sample preparation and analysis followed the standard procedures of Aitken [100,101]. The equivalent dose (De) was obtained from the Simple-aliquot Regenerative Dose (SMAR) protocol of fine particle quartz samples (4~11 μm). The environmental dose rate was determined through activity concentration measurements of the U, Th, and K. The dose rate was computed from the measured activity concentration and attenuation values [100]. Both the cosmic ray contribution and the sample water content were used to assess the environmental dose [102]. The simple and weighted means and standard errors are at a 68% confidence level (Table 1).

Structural Interpretation from Outcrops and Seismic Profiles
The subsurface geometry of the Huoyanshan structure is constrained by the bedrock geology and seismic reflection profiles. Outcrops of the folded strata are exposed by the entrenchment along the Shengjingou Valley. Thus, 21 bedding attitudes were measured across the Shengjingou Valley (Figures 2b and 6). The results reveal a very uniform strike and a continuous and progressive variation in dip (Figure 8), indicating fold growth by limb rotation over listric faults. We acquired 2D pre-stack time migrated seismic reflection data along 3 lines across the eastern part of the Huoyanshan anticline (Figure 2a), which constrains the subsurface geometries to ~7 km below sea level in detail (Figures 9-11). These seismic profiles were converted into depth by utilizing a multilayer velocity model [103] in the Move 2013 software (© 2011 Midland Valley Exploration Ltd., 144 West George Street, Glasgow, G2 2HG, UK). Field observations, regional geologic maps, and well data were used to constrain boundaries of horizons in seismic profiles. Several prominent reflections in the seismic sections, including the bottoms of lithological units J2x, J2q, J3k, K1s, Esh, N1t, N2p, and Q1x, can be readily traced and correlated, allowing a robust structural interpretation and analysis with little uncertainty. The kinematic model was inferred from the characteristics of surface and subsurface geometries. The crustal shortening was calculated using the trigonometric relationship in the kinematic model (Figures 8b), and when paired with the terrace abandonment age, the shortening rate was computed (Table 1).

ADS Analysis
Area-depth-strain (ADS) analysis is a very useful technique for analyzing structural interpretations in contractional systems [104][105][106][107][108][109][110][111][112][113]. ADS analysis provides a simple method of quantitatively constraining the structural evolution history that would not be apparent from the interpretation alone, which in turn validates interpretations [114,115]. Based on the assumption of a constant change in the area during deformation, this method allows the determination of the magnitude of horizontal shortening for growth and pre-growth strata and detachment depth independent of kinematic assumptions. These parameters are determined by measuring the excess area (s) (against the regional level) and the depths (h) (against the reference level) of multiple marker horizons and plotting the corresponding area-depth graph. In this study, ADS analyses followed the method of Wang et al. [111], a generalized ADS method that can be used to handle contractional systems with strata thinning in one direction. The first step was loading three interpreted seismic sections into the program Move and then defining the boundaries where the marker horizons remain undeformed as the pin line and the moving line (regional limits) in each seismic section. Excess areas and depths against the reference level (zero-area line) were measured for each horizon within the program. In the area-depth diagram, x-axis represents the depth against the reference level, y axis represents the excess area, and each horizon has an area-depth value. Balanced interpretations within pre-growth strata yield lines with good linear trends in the area-depth diagram, and the slope of the line equals the overall boundary shortening. The x-intercept equals the distance of the predicted detachment against the reference level. Positive values indicate upward change, while negative values indicate downward change. Note that if the trend of the undeformed marker between two regional limits is not horizontal, the detachment depths at the two regional limits are distinct. Together, the ADS analysis method allows us to constrain the detachment depth and the overall shortening amount more precisely and meanwhile facilitate assessing the structural integrity in the seismic interpretation and geological modeling (Figures 9-11) [108,109].

Inverse Modeling
We adopted a listric fault-related fold kinematic model, including folding over a listric fault connected to a detachment layer at depth. Attention is only paid to the backlimb deformation in this model because this area is often incompletely imaged in seismic lines and it is difficult to model folding in front of propagating thrusts [116]. However, this problem can be solved by using trishear modeling [7,117,118], as it has the ability to explain many of the details of forelimb deformation and mimic those complexities [119]. Thus, in our study, a symmetric trishear was implemented in front of the listric fault to further optimize the kinematic model. The trajectory of listric fault propagation was determined by using the velocity description as described by Cardozo and Brandenburg [62]. The listric fault was modeled as a circular arc with a given center of curvature (CC), the radius of curvature (RC), and the maximum central angle (θmax). In the backlimb, parallel shear takes place and results in constant fault slip behind the trishear zone as well as backlimb geometry concentric to the listric fault. In the forelimb, linear and symmetric trishear operates in front of the fault tip [62,117], and two important parameters, the fault propagation to fault slip ratio (P/S) and trishear angle, mainly control the geometry of the forelimb (Figures 9-11). Several model parameters are relatively easy to determine if the final geometry of the structure is well interpreted, such as the CC, RC, θmax, and fault silp (S). However, the quantification of the trishear angle and P/S ratio has proven elusive [119], because trishear is an incremental rather than a graphical model and there is no mathematical method to quantify these two parameters from the final fold geometry [62]. The difficulty of defining these two parameters has been one of the greatest barriers to the widespread acceptance of trishear. To obtain the trishear angle and P/S ratio, evaluate the fitness of a given model parameter combination, and search for the optimal model parameter combination, the primary vehicle is to carry out inverse modeling. The inversion was carried out by using a script offered by Cardozo (https://thecdem.com/ (accessed on 1 March 2022)) that utilizes the simulated annealing MATLAB function (simulannealbnd in the MATLAB Global Optimization Toolbox) [120].
This inversion method is much faster in comparison to the standard grid search inversion proposed by Allmendinger [118]. Additionally, rather than providing a single solution, the algorithm gives a range of models that can fit the structure (Figures 12-14) [62]. mathematical method to quantify these two parameters from the final fold geometry [62]. The difficulty of defining these two parameters has been one of the greatest barriers to the widespread acceptance of trishear. To obtain the trishear angle and P/S ratio, evaluate the fitness of a given model parameter combination, and search for the optimal model parameter combination, the primary vehicle is to carry out inverse modeling. The inversion was carried out by using a script offered by Cardozo (https://thecdem.com/ (accessed on 1 March 2022)) that utilizes the simulated annealing MATLAB function (simulannealbnd in the MATLAB Global Optimization Toolbox) [120].
This inversion method is much faster in comparison to the standard grid search inversion proposed by Allmendinger [118]. Additionally, rather than providing a single solution, the algorithm gives a range of models that can fit the structure (Figures 12-14) [62].

Surface Geometry and OSL Age Constraint
Along the Shengjingou Valley, eight fluvial terraces were identified, T1 (youngest) to T8 (oldest) (Figures 5 and 6). Figure 7 shows the profiles of the terrace surfaces and riverbeds. The upper older suite (T5-T8) is not very continuous and is preserved primarily in both the terminal parts of the valley. T5-T8 are significantly elevated (~80-200 m) above the riverbed (Figure 7b). Among the lower set of terraces (T1-T4, <100 m above the modern riverbed), T4 is the most continuous terrace and forms an extensive surface along the valley; therefore, we focus on this prominent terrace (T4) in this study. The vertical uplift (subtracted undeformed terrace surface) of this terrace at the peak of the fold is~56 m ( Table 1). In general, terraces in the Shengjingou Valley are syntectonic sediments of folding. These geomorphic markers provide the most recent record of the folding deformation and fault slip activity of the Huoyangshan structure. Consistent and progressive rotation of the surfaces resulted in no parallelism between the younger and older terraces, which is consistent with backlimb deformation resulting from listric thrust-related folding. Further evidence from the outcrop of bedrock along the Shengjingou Valley is consistent with this observation. The attitudes of bedding along the Shengjingou Valley reveal a very uniform strike and a continuous and progressive shift in the dip. The gentler backlimb dips from 40 • to 9 • with northeast dip directions and no dip domain (Figure 8a). This characteristic also implies fold growth by limb rotation over listric faults.
Sample SJK-T4 was collected from fine silts within a 4.5 m section of fluvial gravel cap on the base rock of terrace T4 at the middle of the Shengjingou Valley (Figure 5a), which yields an age of 51.7 ± 5.1 ka. The position of this sample lies at the bottom of the outwash gravel section, indicating that this age is a maximum estimate of the time of terrace abandonment (Table 1).

Subsurface Geometry
Three seismic reflection profiles, Line 1, Line 2, and Line 3, located in the eastern segment of the Huoyanshan anticline (Figure 2a), were used to investigate the subsurface geological structure. These seismic lines have similar geological structures and sediment properties. In all three seismic reflection profiles, above~4.5 km depth, on the northern part of the profiles, we observed fairly continuous reflectors that are subhorizontal. In the central part, these reflectors gently dip northward and are subhorizontal again in the southern part (Figures 9-11). The detachment depth was estimated from the transition between undeformed and folded reflectors (a better constraint of the detachment depth will be presented in Section 4.2). The fault geometry was determined by the discontinuities of the reflectors, delineating a listric thrust boundary of the fold. These geometries outline that the Huoyanshan anticline is an asymmetric fold controlled by a listric thrust rooted in a subhorizontal detachment level. In great detail, as expected from the listric thrustrelated fold model characterized by limb rotation, smoothly continuous progressively rotated backlimb concentric to the listric thrust. The area near the surface is incompletely imaged; therefore, it is difficult to estimate the exact position of the fault tip. The thrust was interpreted to be modeled with a circular geometry, with a center (CC) and a radius (RC) of curvature (Figures 9-11). The interpreted maximum central angles are 55 • in Line 1, 54 • in Line 2, and 55 • in Line 3. Beyond the maximum central angles, the interpreted thrusts are planar and outcrop the ground with dip angle values the same as the maximum central angle values. In Section 4.3, the position of the fault tip in each profile will be expanded on through restorations of deformed reflectors. The locations of axial surfaces were identified by the transitions from the listric ramp to the upper planar ramp and detachment and were used to define the southern and northern fold limits, respectively. Collectively, the seismic interpretation results show direct and robust evidence that the Huoyanshan anticline is a listric fault-related fold, and that the thrust is mainly composed of a listric fault that connects the upper planar ramp with a detachment located~4.8 km under the surface (near the base of J 2 x). This interpretation is consistent with a previous study utilizing the combination of surface terrace deformation and bedding attitudes to infer the kinematics of the Huoyanshan structure [95].

Kinematic Model and Shortening Rate
The features of both the surface and subsurface geometries of the Huoyanshan structure are inconsistent with the traditional fault-related folding models [7,116,[121][122][123][124][125]. We adopted a listric thrust folding theory for the Huoyanshan anticline in this research [98,99,126]. In our listric thrust-related folding model, the thrust grows listrically by limb rotation and connects the upper planar ramp with a detachment located at depth (Figure 8b), and the dip angle of the planar segment equals the maximum central angle values. Shortening (S) was estimated from the dip angle of the thrust at the surface (θ) and the maximal vertical uplift (dz) of the terrace [95,99]. The trigonometric relationship between these parameters is as follows: S = dz sin θ This theoretical geometric relationship allows estimation of the horizontal shortenings associated with the different vertical uplifts and the surface fault dip angle (42 • ). The modelpredicted horizontal shortening is 81.7 ± 2.3 m for T4. Combining the horizontal shortening value with the abandoned age of T4, the horizontal shortening rate was estimated to be 1.6 ± 0.2 mm/a since~52 ka (Table 1, uncertainties were calculated by using Monte Carlo simulation) [127]. The terrace abandonment age is a maximum (described in Section 4.1.1), which results in the shortening rate becoming a maximum.

ADS Results
The upper horizons of the fold suffer from erosion, which in turn affects the integrity of these horizons. Therefore, in each seismic profile, four horizons, numbered 1-4, from the bottom up, which can be traced across the fold, were used to perform ADS analysis in this study (Figures 9-11). The regional level of each horizon was defined as the northward extrapolation of the undeformed horizon imaged on the southern part of the seismic lines, considering its dip as the regional level dip. Where the horizon intersects the move line and the pin line mark, the regional level depths for each horizon are set. The reference level was defined as the interpreted, predicted detachment level. Following this strategy, we measured the excess areas (S) for the four horizons and the distances (h) of the regional levels against the reference levels at the pin line as well as the moving line in each seismic profile (Figures 9-11). Figures 9c and 11c show the area-depth diagrams. Area-depth relationships on both move lines and pin lines all show prominent linear trends with R 2 values of >0.98 in all three seismic profiles, demonstrating that interpretations for all seismic reflection profiles are regarded as reasonably balanced and that the four horizons used to perform ADS analysis are pre-growth strata. Intercepts of the linear fits on the x-axis give the detachment adjustive values. For Line 1, the area-depth diagram gives the vertical distance between the intersections on the regional limits to be~0.4 km. The horizontal distance between the regional limits is~12 km. Thus, the angle between the reference level and the best-fit detachment level is calculated at 2 • , and the maximum central angle is adjusted to 52 • (Figure 9). For Line 2, the horizontal distance between the regional limits is 14 km, and the vertical distance between the intersections on the regional limits is~1.2 km. The angle between the reference level and the best-fit detachment level is calculated at 5 • . The maximum central angle is adjusted to 50 • (Figure 10). For Line 3, the vertical distance between the intersections on the regional limits is~1.3 km, and the horizontal distance between the regional limits is~15 km. The angle between the reference level and the best-fit detachment level is calculated at 5 • . The maximum central angle is adjusted to 49 • (Figure 11). According to Wang et al. [111], the crustal shortening determined from the area-depth diagram at the moving line has a very small error (~3%). Thus, in this study, we followed this strategy. The horizontal crustal shortenings for Line 1, Line 2, and Line 3 are 4.50, 4.44, and 3.10 km, respectively. The profile along the strike of the Huoyanshan fold shows that the displacement along the fold crest displays an asymmetric bow shape, and maximal displacement exists in locations around Line 1, with displacement decreasing eastward and westward (Figure 2c). In addition, from observation of the geological map of the Huoyanshan anticline (Figure 2b), it can be seen that the oldest strata are exposed at the locations of these three seismic lines, indicating that these locations have accumulated maximal shortening. Together, these two aspects demonstrate that the shortening of Line 1 can represent the overall shortening of the Huoayanshan anticline.

Best-Fit Model Parameters
Building on the information given by geometry interpretation and ADS analysis, the interpreted pre-growth horizons 1-4 and listric thrust geometries in Figures 9-11 were read back and digitized as a deformed section to invert for the combinations of best-fit model parameters. Six parameters, the CC, RC, θmax, P/S ratio, trishear angle, and fault slip, were searched in a constrained parameter space in each seismic line. The CC, RC, and θmax were well interpreted in the seismic line and the fault slip was constrained by the ADS analysis. Therefore, the main objective of the inversion is to search over reasonable ranges of P/S and trishear angles. The search space of these parameters has lower limits [CC x −10, The displacement identified by the inversion matches the shortening estimated from ADS analysis in each seismic line. These consistent results confirm the reliability of the adopted inverse procedure.

Estimation of Uncertainties in Structural Interpretation
Although high-quality seismic profiles provide an unprecedented level of detail in thrust settings, a good structural interpretation involves not only a well-imaged seismic profile but also a reliable depth conversion and an applicable theoretical model [128][129][130]. In our case, based on the features present in both surface and subsurface tectonic structures, we adopted a listric thrust-fold theory in our structural interpretation. The model we chose better accounts for features present in nature and avoids oversimplified kinematic interpretation. Then, we combined ADS analysis with reverse modeling of the interpreted structure to obtain a maximum amount of structural information. These structural analysis techniques quantified displacement and detachment depth, validated the original seismic interpretation, and furthermore gave the best-fit kinematic model parameters. However, there are still some uncertainties in each independent technique. For instance, despite the ADS analysis showing a good linear correlation in area-depth values for pre-growth horizons, which reveals a balanced interpretation, the detachment depth and overall shortening derived from ADS analysis may vary with the choice of seismic interpretation because seismic reflection profiles may not image well at the folded area where steeply dipping horizons exist. Similarly, with the inversion procedure, the question judging procedure works well or does not arise from the factor of the uniqueness of the optimal parameter solution. Evaluating the consistency of results obtained from ADS analysis and reverse modeling can help to answer this question and meanwhile estimate structural uncertainties present in interpretations. The inverse model-predicted displacement for each seismic line matches the ADS estimated displacement within a 3% deviation, which is likely within the uncertainty associated with seismic interpretation. In addition, when implementing the inverse modeling procedure, trishear is a continuous deformation. However, the computer models it as a series of discrete displacement steps, and this discrete displacement is very sensitive to inversion [119]. In this paper, we have used a sufficiently small slip increment (less than 0.8% of the total slip) to ensure that the inversion works well. In sum, this combined structure analysis approach provides greater insight into structural informa-tion than can be achieved by any independent technique and realizes mutual verification. This process has implications for guiding the development of fully quantified and robust constrained interpretations of complex structures.

Onset Time of the Huoyanshan Anticline and Its Reliability
For the Huoyanshan anticline, a total shortening of~4.5 km was constrained in Sections 4.2 and 4.3. When combined with the shortening rate of 1.6 ± 0.2 mm/a estimated by geomorphology and OSL dating in Section 4.1, it was inferred that the initiation of folding started at 1.8-3.7 Ma. This onset time of folding was determined from two aspects: the total shortening amount and the shortening rate. The total shortening value is quite robust because the combined structural analysis approach has quantified similar results, and the extent to which the results agree raises confidence in the reliability of this value. Whereas the use of the shortening rate is attached to a range of uncertainties, one must question the extent this rate can be linearly extrapolated in time. Although the assumption that geodetic slip rates are representative of Quaternary slip rates is challenged, there is still a good agreement between GPS and Quaternary slip rates [131]. The consistency between the long-and short-term horizontal shortening rates is broadly shown by the Kyrgyz Tian Shan [38], Pamir Frontal Thrust [53], and Longmen Shan piedmont [132,133]. Therefore, we begin with the comparison of GPS and geologic rates (here referred to as the N-S direction shortening rate) for faults in our study longitude. Results of present-day GPS measurements show that the N-S shortening rate across the eastern Chinese Tian Shan is 3-5 mm/a [76,77]. In the northern Bogda range, an~1 mm/a shortening rate since the late Quaternary was estimated on the northern Bogda fault [58]. To the south of our study range, Wang et al. [94] constrained the shortening rate of the Kumysh fault to~0.3 mm/a during the late Quaternary. Further south, in the Yanqi Basin, Huang et al. [93] estimated that the shortening rate of the north-edge fault is~0.4-0.5 mm/a, and Li et al. [134] suggested that the shortening rate in the Hejing thrust-fold belt has been~0.3 mm/a since the late Pleistocene. These available shortening rates of individual structures, taken together, indicate that the long-term rate of shortening across the eastern Chinese Tian Shan (except for the Houyanshan anticline) is~2 mm/a. If we use the total GPS rates minus this long-term rate, the left rate is~1-3 mm/a, which is a rough estimate. However, it is still consistent with the estimated shortening rate of the Huoyanshan anticline. Relying on this comparison, we conclude that the assumption of constant deformation rates since folding is acceptable. The onset time of the Huoyanshan anticline will be better constrained if we have the chronological constraint of the growth strata and/or long-term shortening rate. However, while these data are absent, our approach to the estimation of the onset time could still provide insight into the deformation process of the eastern Chinese Tian Shan range.

The 2 Ma Synchronous Growth Stage around the Tibetan Plateau
In addition, there seems to be a synchronous growth stage at~2 Ma at the margins around the Tibetan Plateau. For instance, the Dushanzi anticline, located in the third structural belt of the northern Tian Shan piedmont, and an angular unconformity exists between the early Pleistocene and the middle Pleistocene strata [50]. Previous studies all suggested that the folding started at 2.6-1.5 Ma [39,50,52]. Likewise, according to analyses of the growth strata in the Longmen Shan piedmont, middle to late Pleistocene activity at 2.5-1.8 Ma was revealed [133,[135][136][137]. At the deformed margin of the Himalayan foreland basin of northern India and Pakistan, the application of magnetic stratigraphy to constrain the chronology of sediment accumulation and unconformities shows that the thrust sheet began to form from 2.1 to 1.6 Ma [138,139]. At the Himalayan front of northwest India, the surface geology, interpretation of seismic reflection profiles, and magnetostratigraphic data suggest that a minimum shortening of 23 km has occurred since 1.9-1.5 Ma. At the Yarlung Tsangpo River, the eastern end of the Himalayan region, the chronology study of drill core sediments shows that the steepening of the Tsangpo Gorge started at approximately  [140]. At the northern Tibetan Plateau, detailed magnetostratigraphy in the western Kunlun Mountains reflects a stage of uplift that has occurred since 1.6 Ma [141]. In northeastern Tibet, the uplifts of the Heli Shan and Jintannan Shan were reported at 2 Ma [142][143][144]. As noted above, these studies indicate the~2 Ma synchronous growth stage around the Tibetan Plateau and further confirm the reliability of our estimation of the deformation onset time of the Huoyanshan anticline.

Conclusions
The key findings of this paper are as follows: (1) Based on the features present in both surface and subsurface deformation, we first describe a curved thrust ramp and related fold pairs of the structure. (2) By using the kinematic model, deformed record, and OSL dating results of terrace T4, we estimate a geologic shortening rate of 1.6 ± 0.2 mm/a since~52 ka. (3) We interpret three seismic reflection profiles using the theory of quantitative faultrelated fold, ADS, and reverse modeling analyses. These profiles provide direct evidence that this structure connects by a listric thrust ramp to a shallow detachment level. ADS analysis reveals that the maximum shortening of the structure is~4.5 km, which is consistent with a quantitative inverse model. (4) The above shortening rate and total shortening amount suggest that the structure may have formed at 1.8-3.7 Ma, a nearly synchronous growth stage around the Tibetan Plateau.

Data Availability Statement:
The seismic profiles presented in this study can be found in Supplementary Materials.