The data used in this study can be mainly divided into two parts: 1) surface traces of active faults, and 2) deep structures of active faults. The surface traces are derived from the accurate and reliable geological mapping of active faults (Deng et al., 2002; Zheng et al., 2016, Zheng et al., 2019). The deep structures of active faults come from published research of seismic reflection profiles, geological profiles, magnetotelluric profiles, tomographic imaging profiles, receiver function profiles, accurate location of small earthquakes, and focal mechanism solutions in the arc-shaped tectonic belt in the northeastern margin of the Tibetan Plateau. The geological mapping of active faults provides accurate fault surface traces for this study, which is not only the highest priority data for limiting fault surface locations but also the prerequisite for constraining the deep structural profile data of faults. With current technology, it is impossible to directly observe the deep structure of the earth. The indirect observation of the deep geometric structures of the faults through numerous 2D profiles published by researchers is the main way to constrain the deep structures of active faults in this study.

Before the establishment of the 3D geological models, this study fully collected a series of data revealing the deep geometric structures of active faults in the study area, mainly including digital image data of the profile locations and the 2D profiles. Digital image data involves multiple disciplines, technical equipment updates, and different processing institutions, making the collected data have the properties of multiple sources, multiple resolutions, and multiple standards. Therefore, through the evaluation and weighting of multi-source data, this research takes the data closer to the true deep geometric structure of the fault as the main body, extracts control points at a density of about 2 km horizontally, and extracts control points of low-priority data at a horizontal interval of about 4 km, to obtain the 3D space coordinate sets corresponding to the faults and detachment layers.

In the past decades, to analyze the potential seismic hazards, a series of formulas of the moment magnitude-rupture parameters have been proposed by researchers, among which the most widely used is the empirical formulas proposed by Wells and Coppersmith (1994) that the earthquake magnitude corresponds to the surface rupture length (SRL), subsurface rupture length (RLD), rupture area (RA), and average displacement (AD) respectively. Cheng et al. (2020) collected seismic data from 91 earthquakes in mainland China and adjacent areas and used general orthogonal regression to derive the magnitude-rupture scaling relations. Their results were close to the statistical results given by global data without considering the five major earthquakes in eastern China before 1900 (Blaser et al., 2010; Hanks and Bakun, 2002; Wells and Coppersmith, 1994). Therefore, in this study, Eqs 1–5, which are more representative of the empirical seismic magnitude-rupture scale relationship in China, were used to estimate the seismic moment magnitudes (Cheng et al., 2020).

For strike-slip earthquakes M w = 3.77 + 0.98 L o g ( R A ) ( M w < 6.7 ) (1) M w = 1.86 + 1.64 L o g ( R A ) ( M w ≥ 6.7 ) (2) L o g ( R L D ) = ( − 2.45 ± 0.17 ) + ( 0.61 ± 0.03 ) M w (3)

For thrust earthquakes M w = 4.41 + 0.83 L o g ( R A ) (4) L o g ( R L D ) = ( − 3.27 ± 0.66 ) + ( 0.72 ± 0.10 ) M w (5)

At the early stage of this study, the initial 3D fault model V1 of active faults in the arc-shaped tectonic belt in the northeastern margin of the Tibetan Plateau was established (Figure 4, hereinafter referred to as the “V1 model”). The V1 model involves the deep structures of four large-scale active faults, namely, the HYF, the XS-TJSF, the LPSF, and the WQLF, and a set of detachment layers with continuously changing depths is also set up. The surface part of the initial V1 model is controlled by the main 3D control points extracted from the fault surface trace. The deep structure setting of the model is as follows. For the HYF, the WQLF, and the northwestern segment of the XS-TJSF, which are left-slip in nature, the dip angle is uniformly set to 90°. For the LPSF and the southeastern segment of the XS-TJSF with a thrust nature, the dip angle variation along the strike of the faults is determined according to the actual geological surveys. In this model, it is considered that the detachment layer with a depth of about 24 km controls the development of all faults in this region.

In the initial V1 model, firstly, the density of control points extracted from the surface traces of each fault is insufficient, which cannot accurately reflect the surface strike change of each fault, and secondly, the deep structures of each fault are only constrained by the geological dip angles. Therefore, the constraints at all levels of the V1 model are not accurate enough, and only the structure of active faults in the region can be briefly displayed. The model is still far from reaching our modeling goal, but it verifies for us the feasibility of establishing the 3D fault model of major active faults in the arc-shaped tectonic belt in the northeastern margin of the Tibetan Plateau.

Given the shortcomings of the V1 model mentioned in Section 3.1 above, rational optimization of this model is the main direction of this research. The most effective way to optimize the fault surface trace is to encrypt the control points, to improve the accuracy of the fault surface trace. Therefore, we encrypt and extract the latitude and longitude coordinates of the control points at the turning points and curvature changes of the trace, then project them onto DEM to obtain the complete 3D spatial coordinates with depth information. For the deep structure information of each fault, this study combines published geological, geophysical, and seismic data to interpret the deep structures of main active faults and detachment layers in 7 magnetotelluric profiles and 7 auxiliary profiles (Figure 5) in the study area. Then extract 3D spatial coordinates of the deep structures, and obtain the accurate location of small earthquakes in the study area simultaneously (Zhan et al., 2020). Finally, with the aid of the SKUA-GOCAD software modeling function, the 3D fault model V2 of the active faults in the arc-shaped tectonic belt in the northeastern margin of the Tibetan Plateau was obtained (Figure 6, hereinafter referred to as the “V2 model”).

The modeling scope of the V2 model has been extended to more active faults in the study area, such as the MXSF, the YTSF, the XGSF, and the QS-MZF, which are relatively small-scale faults. According to the published research, two sets of detachment layers with depths of 12 and 24 km were installed in the V2 model to replace the detachment layer with a depth of 24 km in the V1 model (Du et al., 2018; Tian et al., 2021; Xia, 2019). In addition, the accurate locations of small earthquakes are added to the model to visualize the relationship between the deep geometry of active faults and the distribution of small earthquakes, which provides model support for studying the relationship between the fault structure and earthquakes and evaluating the regional seismic risk.

The optimization direction of the third model is mainly about the deep structure of the relevant active faults. To achieve this purpose, we searched a large number of published literature and obtained a series of 2D profiles passing through the arc-shaped tectonic belt in the northeastern margin of the Tibetan Plateau. The surface traces, accurate location of small earthquakes, and detachment layers distribution of the third model follow the extracted information of the V2 model. For the deep structures of the faults, this study adopts 7 seismic reflection profiles, 3 geological profiles, 22 magnetotelluric profiles, 9 accurate location profiles of small earthquakes, 4 tomographic imaging profiles, and 1 receiver function profile from many collected profiles (Figure 7, Table 1), aiming to use multidisciplinary data to finely characterize the 3D subsurface structures of active faults. Based on the above data, we built the 3D fault model V3 of the active faults in the arc-shaped tectonic belt in the northeastern margin of the Tibetan Plateau (Figure 8, hereinafter referred to as the “V3 model”).

Through three versions of the models established successively, it is not difficult to see that our results have been gradually transformed from simple modeling to a 3D fine model under multiple constraints. First, from the view of the whole model, the major active faults in the arc-shaped tectonic belt in the northeastern margin of the Tibetan Plateau have been all involved in the modeling process from the V1 model to the V3 model, whether it is a large-scale active block boundary fault across the study area or a small-scale fault that is only under the control of local neotectonic movement. Secondly, starting from the V2 model, referring to the accurate geological mapping of active faults, the surface traces of main active faults in the model are constrained by more 3D control points along the fault strike. Meanwhile, the accurate locations of small earthquakes are considered, and the relative positions are visualized in the model. Finally, the most important deep geometric structure, their constraints are gradually strengthened. The V1 model only depends on the geological profiles, the dip angles of left-slip faults are set to 90° uniformly, and thrust faults refer to the surface dip angles to infer the deep extension of faults. The V2 model mainly refers to a series of magnetotelluric profiles and auxiliary profiles to reveal the deep geometry of active faults. The deep structure of each fault in the V3 model is constrained by multi-source 2D profiles. In general, our current V3 model has meticulously depicted the distribution along the strike and extension along the longitudinal direction of the main active faults in the arc-shaped tectonic belt in the northeastern margin of the Tibetan Plateau from surface to subsurface.

During our model optimization process, the interpolated formed fault planes have also changed dramatically with more detailed changes as the constraint data increased. Given the fact that in the V1 model only the dip angles from the geological field survey can hardly outline the deep geometric structures of each fault, the following fault plane comparisons only involve the V2 model and the V3 model.

The WQLF with left-slip nature is a deep boundary fault that divides the Tibetan Plateau block from the West Qinling block in the southern part of the study area. The fault cuts through the Moho and extends much deeper than our modeling scale. The top of the WQLF plane constructed by fault surface traces, geological profiles, and DEM data shows the elevation variation of the fault surface trace along the strike. The WQLF planes of the V2 model and V3 model are established from different deep geometric information and the plane of the V3 model has more abundant changes in dip angles by contrast (Figures 9A,B). By analyzing and comparing the normal distribution of dip angles on the two fault planes, the dip angles are mainly distributed in the high-angle area (Figure 9C). For the dip angles of the WQLF plane in the V2 model, the average value is 88.8°, the standard deviation is 1.15°, and 99.7% of the dip angles are distributed in the range of 85.4°–90.0°. It shows that the distribution of dip angles is concentrated the nearly vertical. For the dip angles of the WQLF plane in the V3 model, the average value is 77.4°, the standard deviation is 6.99°, and 99.7% of the dip angles are distributed in the range of 60.0°–87.3°. Although the dip angles of the WQLF plane in the V2 model have a more concentrated distribution, the dip angles of the WQLF plane in the V3 model constrained by multi-source data are more in line with the geological surveys (Shao et al., 2011).

The LPSF with a thrust nature is the boundary fault, which divides the Tibetan Plateau block from the Ordos stable block in the east of the study area. Similar to the analysis of the WQLF plane above, the LPSF surface traces of the V2 model and V3 model show consistent elevation relief, and the shape of fault planes is significantly different (Figures 10A,B). The normal distribution of dip angles shows that the distribution center of dip angles in the V3 model is about 16° higher than that of the V2 model, but the distribution widths of dip angles are extremely similar (Figure 10C). For the dip angles of the LPSF plane in the V2 model, the average value is 37.3°, the standard deviation is 12.17°, and 99.7% of the fault dip angles are distributed in the range of 27.9°–70.0°. For the dip angles of the LPSF plane in the V3 model, the average value is 53.6°, the standard deviation is 12.05°, and 99.7% of the fault dip angles are distributed in the range of 19.1°–72.3°.

In the V2 model, the near-surface dip angles of northern LPSF are about 40°, and those of southern LPSF are about 60°. The dip angles of the LPSF gradually increase southward along the strike. Meanwhile, as the LPSF extends deeper, the dip angles decrease, and the LPSF plane gradually merges into the detachment layer with a depth of about 12 km (Figure 10A). However, there is a completely different situation in the V3 model. The dip angles of LPSF gradually decrease southward along the strike. The near-surface dip angles of northern LPSF are about 60° and the northern LPSF plane gradually merges into the detachment layer with a depth of about 22 km. The near-surface dip angles of southern LPSF are about 40°, and the southern LPSF plane gradually merges into the detachment layer at a depth of about 12 km. According to geological investigations (Shi, 2011b; Wang, 2018), the shallow dip angles in the V3 model are more in line with geological facts.

Through the establishment from the V1 model to the V3 model, we have the following insights into the main active faults in the arc-shaped tectonic belt in the northeastern margin of the Tibetan Plateau.

(1) The WQLF is a large left-slip fault cutting the Moho, which has a depth of about 50 km, a total length of about 443.94 km, an overall strike of 278°–287°, a tendency of SSW, and the overall dip angles of 60°–90°. Along the strike, the near-surface dip angles show the characteristics of low in the west and high in the east. Vertically, the dip angles do not change much and the fault plane extends deeper maintaining the original trend.

We suggest that the major active faults in the arc-shaped tectonic belt in the northeastern margin of the Tibetan Plateau are largely controlled by two sets of detachment layers at depths of about 12 and 22 km. Except for the WQLF in the scale of super crust and the HYF and the MXSF in the scale of at least lower crust, thrust faults terminate between 12 and 22 km in depth, indicating that they are mainly controlled by the above two sets of detachment layers. And the QS-MZF and the northwestern segment of the XS-TJSF with left-slip nature extend to about 12 and 22 km, respectively. In this paper, our research is limited to the HYF (narrow sense), through which most of the 2D profiles show that it has not yet cut through the Moho. However, some researchers have suggested that the western segment of the HYF (broad sense) may be a deep fault cutting through the Moho (Shi et al., 2017; Zhan et al., 2008; Zhao et al., 2019, Zhao et al.,2020). The LPSF continues southeastward along the tail of the HYF, of which the northwestern segment has a similar extended depth to the tail of the HYF and surface dip angles of about 60°, while the southeastern segment has an extended depth of only about 12 km and gentler dip angles of about 40°. The change of the deep structure of the LPSF from northwest to southeast probably reflects the fact that the strike-slip component of the near E-W trending HYF is absorbed by the near NW-SE trending LPSF, and the nature changes from predominantly left-slip to predominantly thrust (Zhang et al., 1991; Zheng et al., 2013, Zheng et al., 2016).

To a certain extent, the V3 model has visualized the geometric structure of main active faults in the study area from the overall model. However, there are still a series of uncertainties in the V3 model, which affect the accuracy of the V3 model and determine whether it has practical application value. First of all, the multi-source 2D profiles, which are the basis for establishing the V3 model, have different geological interpretations for the same fault, and the resolution of profiles obtained by different methods is different. Through unified evaluation and screening of these data in the early stage, we can slightly reduce the impact of these uncertainties. Secondly, as the ordinate of some seismic profiles is the two-way travel time, we uniformly adopt the empirical seismic wave velocity of 5 km/s to make the time-depth conversion, instead of using the actual seismic wave velocity. This approach will inevitably distort the converted fault plane in the model and cannot make it closest to the real structure of the fault plane. Finally, during the modeling process, we found that the density of control profiles in different parts of the study area is different, such as 2D profiles through the convergence area of the LPSF, HYF, and southeastern segment of the XS-TJSF are extremely dense. In contrast, the control profiles through the northwestern segment of the XS-TJSF and the western HYF are relatively sparse. As a result, the control points in different areas of the V3 model are unevenly distributed, and different parts of the fault plane constructed by interpolation also have different accuracy.

Through the visual V3 model, we can easily obtain the cutting profiles of the arc-shaped tectonic belt in the northeastern margin of the Tibetan Plateau in any direction, which reflects the deep geometry of regional active faults (Figures 11, 12). Therefore, the V3 model can provide pre-study results for subsequent researchers who plan to carry out relevant geological research in this area, such as collecting long-distance 2D profiles in a certain direction, helping them to arrange and adjust subsequent fieldwork.

In addition, the arc-shaped tectonic belt in the northeastern margin of the Tibetan Plateau is the front of the northeastward expansion of the plateau with active tectonic movements, developing a series of active faults. The evolutionary history of this region since the late Cenozoic is of great significance for studying the expansion of the Tibetan Plateau. At present, the geological research of this area is only limited to the regional surface geological evolution and the deep structures of several active faults (Lei et al., 2016; Shi et al., 2013; Wang et al., 2013; Wang et al., 2014; Zheng et al., 2013), and have not combined the deep structures of the overall active fault and surface geology to jointly explain the regional geological evolution history. The 3D fault model may serve as a bridge to combine the two, effectively linking regional surface geology, deep structures, most of the active faults, and so on. However, if we want to study the small-scale geological structures and the finer structure of a single fault in this area, our model is far from reaching the corresponding resolution and accuracy. The actual distribution of active faults in this region is far more complex than the V3 model. To reveal the deep segmental and geometric features of active faults more precisely in the 3D fault model depends on the improvement of the regional 2D profiles coverage density.

In the modeling process of this study, considering that the resolution and density of the modeling data cannot identify the segmental properties of active faults, the discontinuous surface traces of the faults are fitted into continuous traces. However, in the earthquake occurrence, not the entire fault is active at the same time, but only the energy release of the partial fault. Given this, it is necessary to clarify the segmentation of faults before quantifying the potential seismic hazard (Ding, 1993; Ding, 1995). For strike-slip faults, this study segments the fault planes according to the segmentation of fault surface traces (Howarth et al., 2021; Wesnouky, 1988; Wesnouky, 2006). For thrust faults, the dip angle contour of the fault planes was generated based on SKUA-GOCAD software, and the changing interfaces of dip angles are used to segment the fault planes (Qiu et al., 2016) (Figure 13). The accurate locations of small earthquakes indicate that there are still a considerable number of small-earthquake distributions within the region at least 30 km deep (Du et al., 2018; Sheng, 2015; Xin et al., 2010; Zhan et al., 2020; Zhang B., 2020), whereas our model sets a uniform cut-off of the fault plane extension depth to 30 km, for which the statistical fault plane area can be used to quantify the potential seismic moment magnitude.

Assuming that the planes and surface traces of the seismogenic faults all rupture during an earthquake, we can calculate the potential seismic moment magnitude of each fault segment according to the empirical formulas proposed by Cheng et al. (2020) (Figure 14). If the major active faults in the arc-shaped tectonic belt in the northeastern margin of the Tibetan Plateau produce earthquakes, the estimated moment magnitudes are generally higher than magnitude 6.0.

(1) The WQLF can be divided into five segments. From W to E, the moment magnitudes of each segment are 6.4, 6.6, 6.6, 6.8, and 6.0 deduced from the RAF (RAF = rupture area-moment magnitude empirical formula), and 6.7, 7.1, 7.0, 7.6, and 7.0 according to the RLF (RLF = rupture length-moment magnitude empirical formula).

The moment magnitudes respectively estimated from the rupture area and rupture length of thrust faults are in good agreement, with the difference basically within ±0.3, while for left-slip faults the moment magnitudes calculated from the rupture length are mostly larger than those calculated from the rupture area, with significant deviation. In combination with the deep structures of the main active faults shown in the V3 model, the thrust faults are relatively shallow in extended depth. The potential seismic hazards of thrust faults in this region are greater than those of left-slip faults with high dip angles. Therefore, we pay more attention to the possible seismic hazards caused by thrust faults in this region.

For thrust faults, the estimated moment magnitudes are all less than or equal to 7.3. Except for the southernmost part of the southeastern segment of the XS-TJSF, all moment magnitudes estimated from the RAF are greater than or equal to those estimated from the RLF. These thrust faults, which have experienced many historical strong earthquakes (Figure 1), have the potential for large earthquakes of magnitude 7 or greater. Therefore, the estimations based on the RAF are closer to reality and more reliable. Among thrust faults, the LPSF and the southeastern segment of the XS-TJSF, where the maximum seismic moment magnitude can be respectively up to 7.3 and 7.2. Moreover, the historical records of strong earthquakes in the middle and south segments of the LPSF are few and the passing time is long, so there is a risk of strong earthquakes (Du et al., 2018; Lei et al., 2007; Mei et al., 2012; Wang et al., 2018).

At the same time, the seismic hazards caused by left-slip faults, which also have the potential for a strong earthquake of magnitude 7 or more, should not be ignored. For the WQLF, the QS-MZF, the MXSF, and the northwestern segment of the XS-TJSF, all moment magnitudes estimated from the RLF are greater than those estimated from the RAF. The eastern segment of the HYF is consistent with this phenomenon, but the western segment is opposite. We suggest that this may be because the western segment of the HYF, which is constrained by seismic reflection profile, has smaller dip angles of about 60° compared to other left-slip faults (Guo et al., 2016). Therefore, the estimations based on the RLF are more reliable for the left-slip faults with steep planes, and conversely, the estimations based on the RAF are more reliable. Among the left-slip faults, the westernmost segment of the HYF and the WQLF are at relatively higher seismic risks. The potential moment magnitudes of the westernmost segment of the HYF can reach up to 7.3. Considering that the WQLF is beyond the modeling range, it extends continuously to the west of the model, and the surface trace is fitted into a continuous trace. Therefore, the V3 model underestimates the length of the westernmost segment of the WQLF, and the moment magnitude calculated by the RLF is small, which is most likely larger than 7.0. All segments of the WQLF have potential seismic risks of moment magnitude 7.0 or greater, and one segment can reach a maximum of magnitude 7.6.

Compared with the historical earthquake records, we find that the moment magnitudes estimated by the V3 model are small. We believe that this error mainly comes from the following two reasons. Firstly, the fault traces are simplified in the model, and the rupture length and rupture area of the seismogenic fault are underestimated. Secondly, we calculate the moment magnitudes under the assumption that only a single fault segment ruptures, without considering the cascade rupture of active faults, and the seismogenic structure generally involves more than a single fault segment. At present, the accuracy of the V3 model is not yet sufficient to solve the problem. However, our result reveals the relative seismic hazard and the minimum seismogenic potential of each rupture in the region, which is of the reference value.