Crustal Apparent Density Variations in the Middle Segment of the North Tianshan Mountains and Their Tectonic Significance

Abstract

In this study, we collected mobile gravity observations in the middle segment of the North Tianshan Mountains from August 2016 to July 2022 and carried out classical adjustment calculations under the constraint of the absolute gravity datum to obtain the spatiotemporal variation pattern of the local gravity field. We used equivalent source inversion to obtain the spatiotemporal variation characteristics of crustal apparent density. We also extracted the coseismic deformation field from SAR data, using the 2016 Hutubi earthquake as an example, and constructed a model of the seismogenic fault. The gravity monitoring network in the study area performed well in resolving the earthquake source parameters. Both the time-varying gravity field and equivalent apparent density variation pattern show prominent zoning characteristics with a smoothly evolving spatial distribution over time. The variation trends of the gravity field and equivalent apparent density are in line with the orientation of tectonic structures, and their anomalous signals can be detected before and after an earthquake. The constructed seismogenic structure of the 2016 Hutubi earthquake indicates a typical thrust earthquake, probably occurring on a north-dipping blind fault beneath a region with intense crustal deformation. The subsurface tectonic system reflected by this earthquake can be informatively extended to the entire middle segment of the North Tianshan Mountains by subsurface configuration. Our findings can serve as a reference for analyzing the source characteristics of the time-varying gravity field and interpreting anomalous pre-seismic signals, and aid in understanding earthquake preparation zones and the mode of crustal tectonic movements.

1. Introduction

The earth’s gravity field is a basic physical field reflecting subsurface material migration that directly reflects the process and nature of deep tectonic movements and surface mass redistribution [1]. Earthquakes are closely related to the gravity field through the tectonic movements and variations in the density of subsurface media [2]. There is a long history of research on the relationship between seismic activities and gravity changes, with a landmark achievement obtained by Barnes in 1966. He proposed that before and after an earthquake, the density of subsurface media along the fault varies due to the changes in the regional crustal stress field triggered by fault movement. At the same time, medium deformation leads to the formation of new cracks or expands pre-existing fractures in the media. As a result, fluid media, such as groundwater or volcanic magma, may directly flow in or out, resulting in fluid medium density alterations around the observation site and directly causing temporal variations in the regional gravity field [1,3,4].

Time-varying gravity anomalies are reliable earthquake precursors [5]. In recent years, the rapid development of high-precision absolute gravimeters has made available time-varying gravity data obtained by the terrestrial hybrid repeated gravity observation (THRGO) system, which have been applied to investigate crustal deformation, deep crustal fluid migration, and crustal apparent density variations in tectonically active regions [6,7]. Recently, relatively accurate medium-term predictions of several large earthquakes have been made by interpreting the data of anomalous gravity changes with reference to historical local earthquake records [8,9,10]. Many scholars therefore believe that it is possible to extract the gravity signals that are probably associated with deep mass transfer processes from THRGO records [11], which is the basis for applying time-varying gravity observations to in-depth explorations of earthquake mechanisms and seismic risk prediction [12]. With the advancement of technology, analyzing earthquake-associated surface deformations using interferometric synthetic aperture radar (InSAR) has become a research hotspot, and this field is among those where InSAR technology has been most successfully applied [13]. With the gradual maturity of InSAR technology, increasingly detailed coseismic deformation fields can be extracted to inverse coseismic slip distribution patterns [14]. Due to their large spatial coverage and high-accuracy deformation monitoring capacity, short-period InSAR data have facilitated the rapid acquisition of earthquake-induced surface deformation, post-earthquake field investigations, surface rupture zone identification, and seismic regime prediction, thus providing essential scientific support for earthquake disaster mitigation [15].

In the present study, we used mobile gravity observation data collected in the middle segment of the North Tianshan Mountains from August 2016 to July 2022 to obtain time-varying gravity values on an annual scale based on the repeated ground-based observations from the monitoring network in the study area. We used the classical constrained adjustment method, with an idea of integrating the field and seismic sources. We further obtained the local gravity field variation features and then performed equivalent source inversion based on the time-varying gravity signals to reveal the crustal apparent density changes in the study area. These provide a geophysical basis for a thorough investigation of the process of deep crustal material migration and the generation mechanism of moderate to strong earthquakes. We used the differential InSAR (DInSAR) technique [16,17] to obtain the LOS coseismic deformation field for the 2016 Hutubi earthquake (magnitude 6.2) that occurred during the observation period, and inverted the seismogenic fault based on the elastic dislocation model [18] to study the regional earthquake generation mechanism and fault activity. The software used for gravity field inversion was GEOIST1.0, the software used for InSAR data processing was GMTSAR6.2, and the calculation program used for seismic fault inversion was the Steepest Descent Method (SDM). The outcomes may aid in understanding the tectonic background of earthquake preparation and the tectonic activity pattern around the faults in the study area. Therefore, the purpose of this work was to interpret the seismogenic structure and focal mechanism of the 2016 Hutubi earthquake, with the aim of exploring the geological significance of the deep tectonic features in an extended spatial range covering the entire middle segment of the North Tianshan Mountains.

The research objective was to discuss the seismicity characteristics and earthquake risk areas of the study area by analyzing the characteristics of the time-varying gravity field and the apparent density of the crust in the study area, and combining it with a comprehensive analysis of the meaning of the underground structure obtained from the inversion.

2. Materials and Methods

2.1. Gravity Monitoring Network in the Middle Segment of the North Tianshan Mountains

Mobile gravity monitoring has been conducted by the Seismological Bureau of the Xinjiang Uygur Autonomous Region since the 1980s, with two periods of repeated observation at each station every year. The gravity monitoring network in Northern Xinjiang had been modified twice by 2015, with the number of monitoring stations increased to 97 and network coverage extended to the middle segment of the Tianshan Mountains and the Yili Basin. This enhanced the effectiveness of surface gravity measurements and allowed for better capture of the gravity anomalies in the Xinjiang region. The study area (Figure 1) is located at the southern margin of the Junggar Basin and the middle segment of the North Tianshan Mountains, and features complex geological structures and frequent seismic activities [19]. The monitoring network covers densely populated cities, including Urumqi, Changji, Hutubi, Shihezi, and Dushanzi, which has enhanced the seismic monitoring capability of the region and laid a solid foundation for obtaining data to reveal the spatiotemporal variation characteristics of the local gravity field. Mobile gravity observation in the study area was carried out using CG5 and CG6 high-precision gravimeters, and the personal digital assistant (PDA) system provided by the Institute of Seismology, China Earthquake Administration, was used for data recording and processing. The observations featured a self-deviation within 25 × 10⁻⁸ m·s⁻², a mutual deviation within 30 × 10⁻⁸ m·s⁻², and a unit weight variance of 10.0 × 10⁻⁸ m·s⁻², and therefore fulfilled the specification requirements. In the classical constrained adjustment calculations, segment differences with large errors were rejected or down-weighted. After determining the instruments’ precision, the measurement results were re-weighted accordingly to ensure data accuracy. The mobile gravity records were then gridded to obtain the dynamic spatiotemporal variation patterns in the study area.

Seismic activities are frequent in the study area, with more than one hundred felt earthquakes occurring during our research period. The largest event, with magnitude 6.2, occurred in Hutubi on 8 December 2016. Here, we used mobile gravity observations from August 2016 to July 2022 to explore gravity field changes at the annual scale during these six years. We performed equivalent source inversion based on the time-varying gravity field to obtain the characteristics of the crustal apparent density distribution in the study area and further explored the physical variation features of seismic sources within the crust. Figure 1 was drawn by the Generic Mapping Tools (GMT), in which the topographic relief data are the data of GMT, and the seismic data are from the China Earthquake Networks Center (CENC). Fault data and ground gravity point data are from the Earthquake Agency of the Xinjiang Uygur Autonomous Region.

2.2. Tectonic Background

The Tianshan Mountains are characterized by active tectonic movements and frequent seismic activities. Under the influence of the collision between the far-field Indian plate and the Eurasian plate and the direct squeezing of the Tarim and Junggar basins on both sides, the Cenozoic tectonic changes in the Tianshan area are mainly characterized by the bidirectional thrusting of the mountain and the shortening of the Tianshan crust. In the Late Paleozoic, the Junggar, Tuha, Bogda, and Harlik blocks accreted during the collision between the Siberian and Tarim blocks to form the east–west trending Tianshan orogenic belt. This is a typical reactivated orogenic belt with a basic structure built in the Late Paleozoic and reactivated during the Cenozoic [20]. In the Cenozoic, the Tarim block subducted northward, compressing the Tianshan Mountains under the remote impacts of the India–Eurasia collision and Indian plate subduction, which resulted in the reactivation and uplift of the Tianshan orogenic belt and the formation of the present-day topographic pattern with alternating large mountain belts and basins in the Tianshan region [21]. Plate collision and subduction activities have altered the geological structure of the Tianshan Mountains and the adjacent basins. According to GPS-derived horizontal and vertical displacement data, the Tianshan orogenic belt still has an overall uplift trend [22]. Our study area lies in the middle segment of the North Tianshan Mountains (Figure 2), principally including the structural units of the Urumqi piedmont depression and Tianshan intensive uplift belt. The study area can be roughly divided at 44° N into a subsidence zone in the north, with a maximum subsidence rate of about −3 mm·a⁻¹, and the North Tianshan uplift zone in the south, where the Haxilegen Danban in the Yilianhabierga Mountain area west of Urumqi is the region with the strongest uplift throughout Xinjiang and displays an uplift rate of up to 17.8 mm·a⁻¹ [23].

A series of active faults developed along the foothills of the North Tianshan Mountains. From the south to the north, these include the Yamate fault, Junggar Southern Margin fault, Bogda arcuate fault zone, and Horgos–Manas–Tugulu fault. Multiple types of tectonic structures are co-developed in this area, forming complex and superposed fault systems [24]. The crust in the North Tianshan region gradually shortens from the west to the east and from the south to the north [25,26]. However, this shortening is not uniformly absorbed, indicating that the thrust force gradually weakens from the west to the east along the Tianshan Mountains, in a direction away from the plate boundary; the intensity of folding in the Tianshan Mountains decreases accordingly [27]. Due to the thrusting of the Indian plate [28], the Tianshan Mountains were intensely uplifted in the Late Cenozoic, forming a regenerated orogenic belt [29]. Heavily deformed compressive structures were formed in the piedmont depressions on both sides of the mountain range. These are piedmont thrust sheets inherited from the structures formed during the Late Mesozoic tectonic activities, which continuously developed from the Mesozoic to the Cenozoic, manifesting as a series of thrust fault–anticline belts. The several sets of thrust fault–anticline belts in the Urumqi piedmont depression north of the Tianshan Mountains are the most remarkable features of this range [30].

The Urumqi piedmont depression is located at the northern foot of the North Tianshan Mountains. It is a Mesozoic–Cenozoic depression stretching from the east of Urumqi City in the east to Wusu in the west, with a southern margin separated from the North Tianshan Mountains by the Junggar Southern Margin fault [30]. The maximum thickness of the Mesozoic–Cenozoic sediments in the Manas region within the depression reaches 12,000 m [31]. Subsidence centers of different ages have migrated from the east to the west. The Early Tertiary subsidence center was located around Shawan and Anjihai, and the deposition thickness was the greatest around Wusu in the Late Tertiary. The subsidence center in the Quaternary had migrated to Wusu and west of this location [30]. Three sets of nearly east–west trending thrust fault–anticline belts developed within the Urumqi piedmont depression. From south to north, these are the Qigu piedmont, Horgos–Manas–Tugulu, and Dushanzi–Halaande–Anjihai thrust fault–anticline belts, separated by synclines [30]. In the following, these are referred to as the Q, M, and D thrust fault–anticline belts, respectively. The Q anticlines are parts of an arcuate structure, consisting of six anticlines in an echelon arrangement (i.e., the Qigu, Changji, Kalazha, Xishan, Qidaowan, and Nanfukang anticlines), which were formed by the northward thrusting of the Bogda block [32]. The north limb of the Xishan anticline lies in the Wangjiagou–Toutunhe area, where the Upper Pleistocene gravel bed is tilted with the folds. Its southern limb is faulted by the Xishan fault, with the Badaowan anticline folded into the Upper Pleistocene strata. The Nanfukang folds in the North Bogda Mountain in the east were formed under intense folding activities involving the Upper Pleistocene strata [33]. The M anticlines consist of the Horgos anticline, Manas anticline, Turpan anticline, and Beisantai uplift, which are associated with more intense folding activities than the Q anticlines, and the Upper Pleistocene sequences are generally deformed and folded. The hinge sections or two limbs of these anticlines are typically faulted, and these folds are all active in the Holocene. The intense folding activities have lasted from the end of the Early Pleistocene to the present, controlling the occurrence of large earthquakes, e.g., the Manas earthquake of magnitude 7.7 in December 1906 in the North Tianshan region [34,35,36]. The D anticlines are younger folds involving the Upper Pleistocene strata, and the Holocene active faults are developed along the anticline structures. Although the D thrust fault–anticline belt was formed much later, roughly from the end of the Early Pleistocene to the beginning of the Middle Pleistocene [37], the associated tectonic activities were very intense, with an average shortening rate of 1.7–1.9 mm·a⁻¹ [30,38].

The thrust fault–anticline belts in the Urumqi piedmont depression are mainly asymmetrical anticlines with steep north limbs and gentle southern limbs (Figure 2), with active vertical faults parallel to the fold axes on their north limbs and hinge zones and north-dipping backthrust faults on their south limbs. The activity of the multiple sets of piedmont fold belts gradually increases from the south to the north (Figure 2). These folds mostly have a thin-skinned structure with a shallow depth, forming a pattern of closely packed stripes trending east–west [23].

The study area is one of the most seismically active regions, with frequent moderate-to-strong earthquakes. It has been shown that the frequent earthquakes in this area may be caused by differential block movements [39]. Since the twentieth century, the seismicity in the study region first increased, followed by a calm period, and then increased again [23]. Two earthquakes of magnitude 7 or above have occurred in this area, with the largest one being the Manas earthquake of magnitude 7.7 in December 1906, followed by a magnitude 7.2 event at Xinyuan in 1944 and a magnitude 6.6 earthquake northeast of Urumqi in 1965. No earthquakes of magnitude six or above occurred during the last few decades, and the faults remained locked until the calm period was broken by the 2012 Xinyuan–Hejing earthquake of magnitude 6.6. The most recent strong earthquake was the 2016 Hutubi earthquake of magnitude 6.2, occurring in the piedmont of the North Tianshan Mountains. The seismogenic fault is a thrust fault, and no surface rupture was detected in the vicinity of the epicenter. Due to the lack of understanding of the seismogenic mechanism of intracontinental earthquakes and inadequate research on deep tectonic structures in this region, it is difficult to clarify the local earthquake generation mechanism, and no unified interpretation of the seismogenic fault of the Hutubi earthquake has not yet been obtained [40,41,42].

3. Results

3.1. Equivalent Source Inversion

3.1.1. Basic Principles of the Inversion

Changes in the gravity field are not only the results of crustal movements, but they are also affected by density variations in the materials inside the earth [43]. The time-varying gravity field obtained by repeated surface gravity observations has the advantages of a high spatial resolution, observation accuracy, and data reliability, which provides important support for extracting anomalous signals during the earthquake preparation process. Therefore, this field can be inverted using the “equivalent source” model to obtain the characteristics of apparent density changes and infer the dynamic density variations in the subsurface media [44].

The experimental method can be summarized as the tuning and inversion of time-varying gravity field signals. First, the resolution of the source model is assessed and adjusted through a checkerboard test. The source term is discretized into hexahedral grids in a spherical domain [45], and spatiotemporal constraints are introduced. Then, the “equivalent source” inversion model is optimized and solved following Akaike Information Criterion and Bayesian Information Criterion (ABIC) [46].

The equation of equivalent source inversion [47] can be expressed as

$$\Phi =W_0{\left|Gm-f\left(x,y,t\right)\right|}^2+{\Phi }_s\left(m\right)+{\Phi }_T\left(m\right)$$

(1)

The three terms on the right-hand side of the equation represent the observations, spatial constraints, and temporal constraints, respectively, the details of which depend on the a priori assumptions of the model; $W_0$ is a weight parameter; $G$ is the kernel function matrix; $m$ is the source model; and $f\left(x,y,t\right)$ denotes the three-dimensional space–time model, where $x$, $y$, and $t$ represent the latitude, longitude, and time, respectively.

The spatial constraint term ${\Phi }_s\left(m\right)$ can be expressed as

$${\Phi }_s\left(m\right)=\iint \left[W_1{\left(\frac{\partial m}{\partial x}\right)}^2+W_2{\left(\frac{\partial m}{\partial x\partial y}\right)}^2+W_3{\left(\frac{\partial m}{\partial y}\right)}^2\right]dxdy$$

(2)

$W_1$, $W_2$, and $W_3$ are the weight coefficients (hyperparameters) related to the spatial smoothness of the model in each direction, whose physical implication is that the source model to be solved is expected to have second-order smoothness in horizontal directions of $X$, $XY$, and $Y$.

The temporal constraint term ${\Phi }_T\left(m\right)$ can be expressed as

$${\Phi }_T\left(m\right)=\iint W_4{\left(\frac{\partial m}{\partial t}\right)}^{2}dt$$

(3)

$W_4$ is the weight coefficient (hyperparameter) related to the temporal smoothness of the model, whose physical implication is that the source model to be solved is expected to have second-order smoothness on the temporal scale.

During the inversion, it is necessary to determine all the hyperparameters in the equation. Therefore, the ABIC criterion needs to be introduced to optimize the hyperparameters [48]. It is defined as follows:

$$ABIC=-2\log L+2P$$

(4)

where $L$ is the likelihood function of the model, and $P$ is the number of hyperparameters. The criterion can be minimized using multi-parameter and nonlinear optimization methods, such as Newton’s method or the simplex method.

3.1.2. Checkerboard Test

Before the inversion of the surface time-varying gravity field, it is necessary to examine the source resolution capacity of the gravity monitoring network in the study area. To this end we applied the checkerboard test, which is commonly used in seismic tomography studies [49]. As we aimed to conduct equivalent source inversion for crustal apparent density, we only tested the resolution in the horizontal directions. We assigned alternatively distributed sources of positive and negative anomalies of the same size in the study area, measured the theoretical gravity anomalies at the surface monitoring stations by forward gravity modeling, and then used the inversion method to solve the source parameters and compare the results with the designed values to evaluate the differences in spatial resolution caused by the uneven distribution of the existing monitoring stations [50]. We designed discretized source bodies with a size of 0.5° × 0.5°, a burial depth of 10 km, an apparent density of ±1.0 kg/m³, and a thickness of 1 km to examine the resolution capacity of the gravity monitoring network in the middle segment of the North Tianshan Mountains with the hypothetical equivalent source model [51]. The range of surface gravity anomalies obtained by forward modeling was about ±30 μGal.

Figure 3 shows the theoretical gravity anomalies at the ground surface associated with a checkerboard model with a resolution of 0.5° × 0.5°, as well as the spatial distribution of the monitoring stations in the study area. The locations of the stations are marked by grey rectangles, and the color of the rectangular patches indicates the measured gravity anomalies at the ground surface. The monitoring stations are unevenly distributed in the study area, and the employed resolution was chosen after considering the source resolution capacity of the stations, the optimal design and layout of the monitoring network, and the uncertainty of the source signals [52]. Then, the theoretical gravity anomalies with a checkerboard pattern were inverted to obtain the crustal apparent density using the equivalent source inversion method. The results are presented in Figure 4, where the grey rectangles mark the locations of the surface monitoring stations. The obtained source parameters can provide a reliable reference for the gravity variations in different regions. The figure shows that the predetermined source parameters cannot be accurately recovered in regions far away from the monitoring stations, and that there are certain areas that are not covered by the monitoring network. At the edges of the covered area, the inverted source body densities therefore show transition features, mainly manifested as the constraints of spatial smoothing regularization. The test shows good resolution in the area covered by the monitoring network, and the results were representative of the surface gravity changes induced by the regional crustal material migration. This indicates that the network’s resolution would not have significant impacts on the subsequent inversion of the actual mobile gravity data.

3.1.3. Time-Varying Gravity Field

In this study, we selected the mobile gravity observations recorded by the monitoring network in the study area from August 2016 to July 2022 and further adjusted and inverted the gravity data to analyze the dynamic changes at the annual scale in the gravity field and equivalent crustal apparent density. Figure 5 demonstrates the dynamic changes in the gravity field in the study area, with the rectangles of different colors representing the gravity differences at the monitoring stations after adjustment. The contours of the gravity field changes are also plotted. The focal mechanisms of the 2016 Hutubi earthquake, 2018 Urumqi earthquake, 2018 Hutubi earthquake, and 2020 Toksun earthquake are indicated by black beach balls. Solid black lines delineate faults, and yellow circles of different sizes denote felt earthquakes during the research period.

Figure 5a shows the gravity change pattern in the study area from August 2016 to August 2017. The gravity field in the entire region exhibited apparent zoning characteristics. The Urumqi piedmont degression region was dominated by negative gravity variations, and the foothills of the North Tianshan Mountains all showed positive gravity changes. On 8 December 2016, a magnitude 6.2 earthquake occurred in Hutubi. As shown in the figure, this earthquake occurred at the junction of positive and negative gravity changes, confirming the migration of subsurface materials in this region. The foothills of the North Tianshan Mountains are subjected to intense uplift, and the gravity field change should thus be negative. The Urumqi piedmont degression is undergoing mild subsidence, and its gravity field change should be a positive value of relatively small magnitude. However, we found exactly the opposite (Figure 5a), with the strongly uplifting foothills of the North Tianshan Mountains showing dramatic positive variations, indicating that subsurface medium migration did occur in this region. The regional tectonic stress around the epicenter had changed, and the subsurface media were migrating from the north to the south, eventually triggering the Hutubi earthquake. We believe that these factors probably constituted the main causes of the Hutubi earthquake. Figure 5b shows the gravity change pattern in the study area from August 2017 to August 2018. Compared with the pattern from August 2016 to August 2017 (Figure 5a), the regions with positive gravity changes increased throughout the study area, and the gravity field variation around the epicenter turned from negative to positive. Negative gravity anomalies appeared in the eastern part of the Horgos–Manas–Turpan fault, north of Urumqi City, and in other local areas, suggesting that the tectonic stress field and subsurface medium distribution in the Urumqi piedmont depression region had recovered after the earthquake. On 20 January 2018, a magnitude 4.9 earthquake occurred at the junction of positive and negative gravity changes in the Urumqi area, presenting apparent correlations with the positive gravity variation anomalies at the foothills of the North Tianshan Mountains. Figure 5c shows the gravity change pattern in the study area from August 2018 to July 2019, which presents as completely changed, with the foothills of the North Tianshan Mountains and the Urumqi piedmont depression region showing negative gravity variations and the southern Junggar Basin turning into a positive gravity change zone. On 18 August 2018, another magnitude 4.9 earthquake occurred in Hutubi; at this time the gravity field anomaly was changing from positive to negative, which further indicates that subsurface medium migration or exchange might occur in the study area. It has been frequently reported that moderate to strong earthquakes tend to occur in the four quadrants of gravity changes, anomalous regions showing alternating positive and negative gravity changes, or areas with a strong gravity gradient [1,53,54].

Figure 5d demonstrates the gravity change pattern in the study area from July 2019 to August 2020. Compared with Figure 5c, most of the study area was dominated by negative gravity changes, and the region with positive change was migrating southward in the horizontal direction. The degree of negative gravity variation increased in the western part of the study area, including the Cholma region, indicating that migration of subsurface materials from the west to the east was again occurring in the study area. There was an isolated point of positive gravity change in the Shihezi area, suggesting that tectonic stress might be accumulating in this region. Figure 5e shows the gravity change pattern in the study area from August 2020 to July 2021. Regions with positive gravity changes were distributed in Dushanzi, Shihezi, and Cholma and extended along the Horgos–Manas–Turpan fault to the north of Urumqi. The foothills of the North Tianshan Mountains still showed negative gravity changes. On 8 August 2020, a magnitude 4.9 earthquake occurred in Toksun County, a transition zone between positive and negative gravity anomalies. This seismic signal was once again correlated with the gravity field anomaly. Figure 5f shows the gravity change pattern in the study area from July 2021 to July 2022. The gravity field in the study was dominated by positive changes, with negative gravity variations only locally distributed in Urumqi. The enhanced degree of gravity field changes in the study area indicates that the energy level in this region increased and that tectonic stress was accumulating.

The time-varying gravity field in the study area exhibited clear zoning features with smooth changes. Although the magnitude of the gravity variations was not drastic, the site values showed alternating positive and negative changes. Comparison with the seismic events and fault activities in the study area suggests an apparent correlation between the earthquake clusters and gravity anomalies. These earthquakes mostly occurred in the regions with a high gravity gradient or the anomalous zones with alternating positive and negative gravity changes [55]. Considering the regional tectonic background, it can be concluded that the study area is still at risk of moderate to strong earthquakes.

3.1.4. Apparent Density Changes of the Sources

Based on the adjustment results of the mobile gravity observations in the study area, we used the equivalent source inversion method to obtain crustal apparent density values (Figure 6). The geographical elements in this diagram are the same as those in Figure 5. Strictly speaking, the gravity variations are closely related to environmental factors local to each monitoring station, such as the migration of subsurface fluids, soil moisture, vehicle vibration, and elevation changes [12]. These factors, which are difficult to determine, can be uniformly defined as external factors of the sources. They are highly contingent on small-scale drivers and difficult to quantitatively eliminate. Reducing any of them to a fixed value will affect the reliability of the results. In this study, we therefore adopted the equivalent source inversion method and introduced the a priori condition of space–time smoothing into the inversion model to suppress anomalous signals, such as local disturbances and source errors. This allowed us to more accurately extract the anomalous time-varying gravity field signals and determine the parameters of the source model.

For the inversion, we selected equivalent crustal density anomalies of 1 km thickness at a depth of 10 km in the study area that could be effectively detected by repeated surface gravity measurements. The characteristics of apparent density variations obtained by inversion are similar to the evolution trends of the gravity values, and the major anomalous regions are consistent. Data (crust1.0) show that the average crustal density in the study area is approximately 2700 kg·m⁻³, and the apparent density variation of the sources in the study area is roughly ±3.0 kg·m⁻³. The inverted density variation is about 1‰ of the crustal density. The time scale of each panel in Figure 6 is the same as that in Figure 5, and the color bar for the time-varying gravity field variation is consistent with that for crustal apparent density variations to allow for a better illustration of the anomalous features revealed by the inversion.

Figure 6a shows the apparent density variation pattern in the study area from August 2016 to August 2017. Apparent density anomalies were generally distributed along the trend of the tectonic structures, with distinct zoning features. The Urumqi piedmont depression was a region with negative changes, while the foothills of the North Tianshan Mountains showed positive variations. These density changes verified our inference based on the gravity changes described above, i.e., that the crustal media were migrating from the north to the south. There was an incongruous negative anomaly at the epicenter location, where the density variation contour bends southward, out of line with the tectonic trends. We believe that this feature was not simply caused by noise but very likely constitutes a manifestation of deep crustal material migration before and after the earthquake. According to the theory of isostasy, the regional uplift of the crust will lead to a decrease in its density, thus achieving a state of mass balance among different regions [56]. Therefore, the anomaly shown in Figure 6a may reflect the influence of post-seismic surface deformation and crustal uplift on subsurface media; that is, the seismic activity led to crustal uplift, which induced a decrease in crustal density. The possible reasons for this anomaly are further discussed and verified in combination with the InSAR data in the next section. Figure 6b shows the apparent density change pattern in the study area from August 2017 to August 2018. The region with positive apparent density variation, which was previously distributed along the trend of tectonic structures, significantly expanded toward the north, and the negative anomalies were only distributed in a small area northeast of Urumqi. This phenomenon suggests that the crustal density deficit at the Urumqi piedmont depression was restored after the earthquake. On 20 January 2018, a magnitude 4.9 earthquake occurred at the junction of positive and negative gravity field changes in the Urumqi area, with an epicenter also located at the junction of positive and negative changes in apparent density. We suggest that the migration of subsurface materials altered the tectonic stress field, thereby triggering the fault activity. Figure 6c shows the apparent density change pattern in the study area from August 2018 to July 2019. The regions with positive and negative apparent density changes in the study area were completely interchanged, with the foothills of the North Tianshan Mountains and the Urumqi piedmont depression showing a decrease in apparent density and the north of Urumqi City and the Hutubi area showing an increase. During the crustal density change, a magnitude 4.9 earthquake occurred again in Hutubi on 18 August 2018. Figure 6d shows the apparent density change pattern in the study area from July 2019 to August 2020. The zoning features in the study area became more prominent, with an enhanced negative apparent density change in the Dushanzi area, and the region with an increase in apparent density shifted to the south, indicating that crustal material migration occurred again from the west to the east in the study area. Figure 6e shows the apparent density change pattern in the study area from August 2020 to July 2021, displaying obvious zoning features. The regions with an increase in apparent density expanded around Dushanzi, Cholma, and Shihezi but shrunk around Urumqi. The degree of apparent density decrease at the foothills of the North Tianshan Mountains was reduced. On 8 August 2020, a magnitude 4.9 earthquake occurred in Toksun County, where the apparent density change was turning from positive to negative, once again corresponding to the crustal material migration feature in this region. Figure 6f shows the apparent density change pattern in the study area from July 2021 to July 2022, displaying a predominantly increasing trend, while an apparent density decrease was only locally distributed in Urumqi. This increase in apparent density in the study area suggests that the tectonic stress had been accumulating and that the energy level in the study area may be increasing.

In summary, the positive and negative equivalent apparent density anomalies in the study area show zoning characteristics, and their spatial distribution pattern has been evolving and repeating. The source variation features correspond well to the seismicity pattern, and their measurement accuracy is sufficient for further studies on the magnitude of equivalent apparent density changes obtained by inversion. For our study area, the inversion results of the sources’ apparent density can serve as the basis for analyzing local seismicity and determining the sources of future seismic risk.

3.2. Tectonic Significance of the Middle Segment of the North Tianshan Mountains

3.2.1. The Hutubi Earthquake

Regional tectonic movements are closely related to seismicity. Analyzing the seismogenic structures and focal mechanisms of moderate to strong earthquakes in an area can help in the interpretation of the area’s tectonic significance. According to the China Earthquake Networks Center, a magnitude 6.2 earthquake occurred in the Hutubi area on 8 December 2016, Beijing time, with the mainshock epicenter located at 43.83° N, 86.35° E, and with a focal depth of 6 km. The earthquake has been the subject of much subsequent research. Kong et al. [57] reported out that the aftershocks were concentrated to the west of the mainshock, and they inferred, based on the focal mechanism and the spatial distribution of the aftershock sequence, that the Junggar South Margin fault was the seismogenic structure of the Hutubi earthquake. Liu et al. [58] suggested that the seismogenic structure of this earthquake was a north-dipping thrust fault, probably formed by fault inversion. Lu et al. [41] concluded that the focus of this earthquake was at a depth of 16.5 km, located near the Horgos–Manas–Tugulu fault; they consequently suggested that both the Hutubi earthquake and the 1906 Manas earthquake occurred on the Horgos–Manas–Tugulu fault. Yang et al. [42] performed cut-and-paste inversion to obtain the focal mechanism of the earthquake and speculated that the seismogenic fault might be the south-dipping Qigu fault. As shown by the inverted equivalent apparent density changes described above, there existed an anomaly around the epicenter of the Hutubi earthquake between August 2016 and August 2017 (Figure 6a). We believe that this feature could be a direct manifestation of the crustal density change induced by the movement of the subsurface tectonic system before and after the earthquake; this is further explored below. Due to the scarcity of GPS stations in the study area, no coseismic signals of this earthquake were recorded, and the brief observation period prevented obtaining a significant transient response of the time-varying gravity field to the earthquake. Therefore, we chose to use the SAR data to comprehensively study the earthquake with the D-InSAR technique in combination with the time-varying gravity field and crustal apparent density with the aim of interpreting the tectonic significance of the entire middle segment of the North Tianshan Mountains.

3.2.2. InSAR Coseismic Deformation Field

We employed C-band Sentinel-1A TOPS-mode SAR images from both the ascending and descending tracks provided by the European Space Agency (ESA) to obtain the line-of-sight (LOS) coseismic ground deformation field using the D-InSAR technique. The parameters of the SAR images are listed in

Table 1. Parameters of the Sentinel-1A and Sentinel-1B SAR images used in the analysis of the 2016 Hutubi earthquake.

Satellite	Track (No.)	Subdivision Number	Observation Time	Vertical Baseline of Interference Pairs (m)	Azimuth (°)	Incidence Angle (°)
Sentinel-1A	Ascending track (114)	12	Master image: 19 November 2016 Sub-image: 13 December 2016	−0.6	346	39.4
Sentinel-1B	Descending track (92/3)	3	Master image: 24 November 2016 Sub-image: 18 December 2016	125.4	194	33.6

.

We used the GMTSAR software (UC San Diego: Scripps Institution of Oceanography) for data processing and made use of ESA high-precision orbit data, SRTM DEM data published by NASA (with a spatial resolution of 30 m), and GACOS atmospheric data for orbital correction, topographic phase removal, and atmospheric phase weakening, respectively. The interferograms were processed by Gaussian filtering to raise their signal-to-noise ratio. Then, we applied the minimum-cost flow method for phase unwrapping [59] and established bilinear trend surfaces for the orbits based on the far-field data in the coseismic deformation field to remove the trend of residual phases generated by the errors in the high-precision orbital data [60]. After geocoding and phase transformation, the coseismic deformation fields from the ascending and descending tracks covering the entire study area were obtained (Figure 7).

Figure 7 shows that the coseismic deformation field of the Hutubi earthquake had an irregular elliptical shape. It was located on the north side of the Junggar South Margin fault, crossed through by several fold belts. The InSAR images from both the ascending and descending tracks show positive deformation signals in the earthquake region. After removing the noise interference, the maximum displacement in the LOS direction was 30 mm in the ascending track image (Figure 7a) and 25 mm in the descending track image (Figure 7b). The surface deformation caused by the earthquake was relatively small, and the fault rupture did not extend to the ground surface, indicating that the earthquake-induced ground deformation was mainly controlled by the uplift movement. Thus, the surface displacement associated with the Hutubi earthquake might have been caused by fault thrusting, which agrees with the focal mechanisms provided by the USGS and GCMT.

3.2.3. Inversion of the Coseismic Fault Slip Model

Using the acquired InSAR ascending and descending track data as constraints, we obtained the parameters of the seismogenic fault using the simulated annealing inversion method based on the theory of dislocations in an elastic half-space. We then conducted linear inversion to obtain the fault slip characteristics, finally deriving the geometric position of the slipping fault and the fault movement model. The spatial continuity and wide coverage of the ascending and descending track data selected for this study resulted in a large data volume, which had to be downsampled to improve computational efficiency. To ensure the reliability of the results, we first masked the deformation field based on the coefficient of coherence to remove values with low accuracy. Then, the unmasked area was downsampled using the quadtree sampling method [61] to retain the completeness of deformation characteristics in the study area. Because the quadtree sampling method can effectively compress the data while keeping the deformation information well and can achieve the purpose of noise elimination to a certain extent, it has good universality. And, compared with other methods, the quadtree sampling method has the advantages of being a simple algorithm, and having fast calculation efficiency and wide application. The results are shown in Figure 8.

The relocating of the Hutubi earthquake and its aftershocks [57,62] revealed that the aftershocks were mainly concentrated to the west of the mainshock. Based on the characteristics of the coseismic deformation field and the spatial distribution of the aftershock sequence and in view of the wedge-shaped structural characteristics of the fault system in the fold belts north of the North Tianshan Mountains, we preliminarily constructed a north-dipping wedge-shaped fault model with a length of 40 km to describe the seismogenic structure of the Hutubi earthquake. The search intervals for the optimal strike, dip angle, and focal depth were set as [260°, 280°], [5°, 80°], and [−20 km, −5 km]. The optimal geometric parameters determined using simulated annealing consisted of a strike of 276°, a dip angle of 70°, a depth of −17.85 km, a slip angle of 98.77°, and an earthquake magnitude of Mw5.98. A comparison to focal mechanism parameters sourced from different institutions and researchers (

Table 2. Parameters of the 2016 Hutubi earthquake derived from various sources.

Data Source	Strike (°)	Dip (°)	Slip Angle (°)	Focal Depth (°)	Magnitude (MW)
USGS	269	71	93	17.6	6.0
GCMT	272	68	92	20.3	6.0
SCANDEC	276	70	95	20	6.0
CENC	277	69	88	19	6.0
IGCEA	273	70	108	17.6	5.95
Present study	276	70	98.77	17.85	5.98

) suggests that these strike, dip angle, slip angle, and focal depth values are within a reasonable range.

After obtaining the geometric parameters of the fault plane, we discretized the fault plane into small rectangular blocks of 2 km × 2 km and inverted the motion vectors of the fault plane under the constraints of the downsampled deformation fields from the ascending and descending tracks, in order to obtain a detailed fault slip distribution. The resulting reverse fault slip distribution of the Hutubi earthquake is shown in Figure 9. The length, width, and depth of the fault plane were approximately 40 km, 40 km, and 35 km, respectively; the arrows represent the motion vectors of the discretized fault slips. The coseismic fault slip was mainly distributed at a depth of 10–30 km, in a range of 40 km long along the strike. The maximum slip was approximately 10 cm at a depth of approximately 15–20 km. The model suggests that the seismogenic structure of the Hutubi earthquake was a north-dipping thrust fault with a small slip magnitude and that the fault rupture did not extend to the ground surface. Research based on teleseismic data has shown that the spatial coherence of the north-dipping fault with the aftershocks was significantly higher than that of the south-dipping fault, indicating that the 2016 Hutubi earthquake likely resulted from the movement of a north-dipping reverse fault [58]. Therefore, here, we only discuss the north-dipping fault, and conclude that the parameters obtained from the experiments are in line with the actual situation. Combining the data and experimental results, we believe that the north-dipping fault obtained by the inversion carried out in this study was the seismogenic fault of the 2016 Hutubi earthquake.

3.3. Tectonic Significance of the Study Area

Figure 10 shows the geological profile inferred from the seismic reflection data of the 2016 Hutubi earthquake. Figure 10 is modified from Lu et al. [41]. The location of this profile is marked by the solid cyan line in Figure 1. Prominent fold-and-thrust belts developed in the North Tianshan foreland basin, with the solid lines indicating the interpreted faults and the dashed lines representing the axial planes. F1 denotes the Yamate fault, F2 the Junggar South Margin fault, and F3 the Horgos–Manas–Tugulu fault. The solid white circles indicate the aftershocks of the Hutubi earthquake, and the solid black lines passing through them are possible seismogenic faults. The yellow star marks the epicenter of the 1906 Manas earthquake, and the solid white circles indicate the aftershocks, whose spatial distribution implies that the seismogenic fault should be a north-dipping fault. However, as no north-dipping faults have been detected in this region, we suggest that the mainshock of the Hutubi earthquake might have occurred on a north-dipping blind fault at a depth of 5–20 km. This structure is extremely similar to the fault model obtained by the inversion described above; this also verifies that the inverted fault model fulfils the research requirements and expectations. F2 and F3 extend downward and finally intersect the detachment plane at a depth of about 24 km. Between them, strata of different ages overlap to form a wedge structure. The geological structures depicted in this figure are extremely complex, and the deformation characteristics of these faults and their related folds are similar to those of large-scale wedge structures [63,64,65], being composed of small wedges formed under forward and backward thrusts [66]. Blind faults usually occur as individual faults, sets of faults, or thrust systems. The overlapping strata may respond to a blind fault system by back-thrusting, coupling, forward-thrusting, or sedimentation [67]. The thrust structure of a large-scale tectonic wedge contains many blind faults, but F2 and F3, rupturing at the ground surface, are not blind faults. The Hutubi earthquake did not cause any surface rupture as its epicenter was deep, and its magnitude was not particularly great.

The 1906 Manas earthquake occurred at a much greater depth, but historical records show that it caused considerable damage. According to its deep structural features and the regional tectonic background, this may have been because energy was released at a deep blind fault under the reactivated uplift of the Tianshan orogenic belt, which triggered tectonic activities at a shallow depth that resulted in violent crustal deformation. The young fault scarps along the Horgos–Manas–Tugulu thrust fault–anticline belt may be surface ruptures associated with the M 7.7 Manas earthquake [33,68]. Therefore, we believe that the triggering mechanism of the 2016 Hutubi earthquake may be similar to that of the 1906 Manas earthquake, both of which were initiated by deep fault activities. The 2016 Hutubi earthquake might have been induced by the slip of a blind fault at a depth of 5–20 km under tectonic stress, which in turn caused the movement of deep tectonic systems that manifested as the uplift of the ground surface at the epicenter.

Our study area is located in the middle segment of the North Tianshan Mountains, between the northern margin of the Tianshan Mountains and the Junggar Basin. The Tianshan orogenic belt is one of the largest orogenic belts within the Asian continent and was formed by multiple phases of plate subduction and collision [69,70,71]. The 2016 Hutubi earthquake occurred right on a fold-and-thrust belt within the Urumqi piedmont depression, a region that absorbs about 50% of the shortening from the collision between the Indian and Eurasian plates [72,73,74]. This earthquake was the largest to occur in the North Tianshan fold-and-thrust belt in the past century (1916–2023), and it released some of the accumulated strain in the northern part of the Tianshan fold-and-thrust belt [75]. The formation of fault-related folds can be explained by the upper crustal folding associated with brittle processes [64,76,77]; these fold structures have played an important role in coseismic deformation [78,79,80,81,82]. The occurrence of the 2016 Hutubi earthquake revealed a macroscopic tectonic background that suggests that the back-thrusting of the wedge structure in the Tianshan area offset the tectonic displacement from the Tianshan Mountains. Conclusive evidence can be obtained from the surface deformation pattern and the inverted fault model, both of which indicate that the Hutubi earthquake was caused by blind thrusting and the associated folding of the overlying strata in the North Tianshan fold-and-thrust belt.

The geological structures in the study area are mainly manifested as east–west trending faults, and the risk of fold-type earthquakes is high in the northern part of the Tianshan fold-and-thrust belt [36,83]. Our analysis of the earthquakes associated with deep blind thrust faults and fault-related folds suggests that attention should be paid to the shallow folding activities within active fold-and-thrust belts, since the faults developed in these shallow fold belts may be ruptured during moderate to strong earthquakes occurring on deep faults and in some cases generate more severe consequences than the mainshocks. At the same time, the activities of deep blind faults and the folding within thrust belts should also be monitored, as the macroscopic tectonic background they reflect constitutes an important basis for studying the tectonic significance of the study area and the entire Tianshan orogenic belt.

4. Discussion

The time-varying gravity field in the study area over the past 6 years was characterized by apparent zoning of positive and negative anomalies, and the transition between anomalous zones was found to be smooth. Positive and negative gravity changes appeared alternately but at a moderate magnitude. A comparison with the seismic activities in the study area shows an apparent correlation between earthquake clusters and gravity anomalies; earthquakes mostly occurred around areas in the high strain zones, in the four quadrants of gravity anomalies, or in the anomalous regions showing alternating positive and negative gravity changes. Before the 2016 Hutubi earthquake, the time-varying gravity field around its epicenter appeared as a transition zone between positive and negative gravity anomalies.

The apparent density variations obtained by equivalent source inversion are matched well with the time-varying gravity field changes. The positive and negative equivalent apparent density anomalies show a prominent zoning characteristic, with an evolving and repeating spatial distribution pattern (Figure 6). The source variation features exhibit a good correlation with the seismicity pattern. Seismic activities mostly occurred in regions with dramatic apparent density changes, such as the transition between positive and negative changes and rapid apparent density increases or decreases. The crustal apparent density inversion in this study fulfils the research objective. The examined evolutions of the time-varying gravity field and spatial and temporal crustal apparent density distribution can serve as a basis for analyzing the seismic activity and determining future seismic risk in the study area.

We used the 2016 Hutubi earthquake as a case study to thoroughly explore the tectonic significance of the deep part of the study area by examining the geometric parameters of the earthquake’s seismogenic fault. The fault model of the earthquake was obtained by inversion based on the SAR-derived coseismic deformation field, and the focal mechanism of this earthquake was investigated to explore the transient response of the subsurface tectonic system to the earthquake. Finally, the tectonic structures were extended to the entire study area. This was done to discuss the regional tectonic background and interpret the tectonic significance of the area’s deep layers. The aim was to provide a basis for elucidating local seismogenic mechanisms and enabling early earthquake warnings. The regional time-varying gravity field was obtained by THRGO to extract and analyze the long-term features involved in the earthquake preparation process. Further, the coseismic deformation field was obtained from the SAR data to invert the fault model and synthetically reflect the transient response of the tectonic system to the occurrence of an earthquake. This approach can not only make up for the shortcomings of the THRGO system, which is not able to capture the transient features of earthquake events, but also provides detailed background information on regional seismic risks and tectonic evolutions for coseismic deformation studies. For example, Zhang conducted inversion algorithms and simulation experiments on a spherical model of the time-varying gravity field, and Wang used the gravitational field changes before and after the Lushan earthquake in 2013 to reveal the transfer of deep matter.

Previous studies have shown that the Junggar block is subducted beneath the Tianshan block under near north–south compressive stress, forming a multi-layered wedge-shaped thrust system in the shallow crust of the Tianshan piedmont in the study area [84]. Although earthquakes in the study area have mainly occurred on shallow blind thrust faults or fault-related folds, moderate to strong earthquakes on deep faults can also trigger ruptures at shallow depths, resulting in even more serious consequences than the mainshock. The deep blind faults and fault-related folds present in the study area should therefore also receive attention. The spatiotemporal variations in the gravity field can reflect the movements of major tectonic faults in the monitored area [2], and deep tectonic activities are also capable of altering the surface gravity field [85]. This enables the extraction and analysis of anomalous gravity field signals to evaluate the seismic risk of a region. Because the inversion of crustal apparent density can reflect deep material migration and medium flow, the comprehensive analysis of the apparent density variation trend and the tectonic activity background can help in the identification of possible anomaly signals before an earthquake, which is essential for seismic risk assessment and early earthquake warnings.