Tectonic belts along active tectonic block boundaries comprise one or more active faults; along which, large earthquakes recur. Therefore, it is important to establish the recurrence behavior of large earthquakes along such boundary zones for studying their characteristics and developments. Many paleoearthquake studies make it possible to investigate the recurrence behavior of large earthquakes along the northern boundary of the Ordos block (NBOB). Based on the previous studies, data from 52 trenches were collected to reconstruct prehistoric earthquakes using an improved multiple trench constraining method. This method is based on paleoearthquake indicators and trench location distribution to constrain the rupture time and length, thereby reducing the selection bias of fixed rupture length to construct additional rupture scenarios. The results suggest that the NBOB comprises four normal faults (from west to east): the Langshan Piedmont Fault (LPF), Sertengshan Piedmont Fault (SPF), Wulashan Piedmont Fault (WPF), and Daqingshan Piedmont Fault (DPF); along which, six, seven, eight, and six paleoearthquakes have occurred within approximately 15,000 yr, respectively. In addition, recurrence behaviors of the individual faults exhibit remarkable periodicity. The regional fault network along the NBOB reveals clustered characteristics with six clusters propagating either westward or eastward and a recurrence time of approximately 1,300 yr. Large earthquake events have occurred along the LPF, WPF, and DPF according to the most recent cluster; however, earthquakes were absent along the SPF, and no evidence of large earthquakes was observed along the NBOB after the 849 CE earthquake. Thus, we discuss the possibility of occurrence of large earthquakes along the SPF after the 849 CE earthquake based on earthquake recurrence and cluster migration behavior. Additional research is required to assess the potential risk of the occurrence of a large earthquake along the SPF in the future.

The Chinese mainland is divided into several active tectonic blocks by boundary fault zones. Large magnitude earthquakes occur frequently along these boundaries, implying that there is a link between the two [18]. Therefore, it is critical to understand the formation mechanism and occurrence of such earthquakes along active tectonic block boundaries, based on their long recurrence cycles and exceedingly destructive characteristics [911]. Paleoearthquake studies can correct flaws in instrumental seismic records (temporary and infrequent) [12, 13] and historical earthquakes records (inaccurate and omissive) [14, 15] by characterizing and dating major prehistoric earthquakes based on geologic evidence [1619]. The methodology allows paleoearthquake records to extend back thousands of years and offers a foundation for an improved understanding of earthquake behavior [18, 2027]. Evidence of several prehistoric earthquakes has been excavated [2833] from trenches dug at particular sites, sections, or faults along the northern boundary of the Ordos block (NBOB) in North China, which has indicated that different fault segments of these regions were active. However, whether large earthquakes occur in a regular spatiotemporal pattern along the NBOB or more randomly in time and space remains unknown. Moreover, the activities of different earthquakes along single faults and regional fault networks remain unclear. To resolve these knowledge gaps, paleoearthquake records, which are fundamentally important for seismic risk assessment, must be completed. Fortunately, many paleoearthquake studies with age constraints on the NBOB have provided a solid basis for addressing such issues.

Numerous paleoearthquake indicators have been observed in trench profiles, including morphological and sedimentological evidence of ground deformation across boundary fault zones, and have been used to optimize data on the time and length characteristics of paleoearthquake rupture [16, 18]. However, trenches may not record all the evidence of an earthquake, as observed in studies along strike-slip and normal faults such as the San Andreas Fault [16, 21, 25, 34, 35] and Wasatch Fault [21, 36, 37], respectively. Trench exposures may be shallow or may record for a shorter period or fewer events [3841]. In contrast, multiple trenches located on different fault sections can reconstruct a complete sequence of recent paleoearthquakes [42, 43], as the events observed in a single trench are either reoccurring or supplement each other [44]. Earthquakes are generally restricted within a fault segment [45] and are subject to the characteristic earthquake model [21] based on the elastic rebound theory [46], which states that the rupture length and magnitude (limited by the rupture length) are fixed [21, 4749]. However, actual earthquake records show that characteristic earthquakes may be an artifact of selection bias [48, 5053] that have occurred sparingly [5456]. Some of these events have occurred erratically and appear to support a wide range of earthquake models with various features [53, 57, 58], indicating that previous rupture scenarios were inadequately constrained based on the fixed fault segment rather than on the location distribution of the trench. Herein, we provide a solution for paleoearthquake record limitations, which is based on paleoearthquake indicators to constrain the occurrence time and rupture length of the event.

Based on data accumulated in recent years, it is necessary to investigate the occurrence and implications of long-lasting seismic quiescence in the vicinity of the NBOB. In this study, we collected data from 52 trenches in the NBOB region to understand its seismic behavior. We used an improved multiple trench constraining method to reconstruct large earthquakes in the region and explored more reasonable rupture forms, built a relatively complete paleoearthquake series, and established the spatiotemporal rupture history pattern. We then compared the most recent and historical paleoearthquakes and summarized their characteristics according to the large earthquake recurrence mode parameters of single faults and clusters along the boundary fault zone. Furthermore, we integrated relevant parameters and earthquake recurrence modes for the NBOB, discussed the potential of large earthquakes since the historical 849 CE earthquake, and estimated the maximum potential magnitudes of future earthquakes.

The Ordos block is a unique uplifted region and an active tectonic block located in the northern part of North China that undergoes little interior deformation [5, 59]; it is surrounded by a series of fault basins and faults that experience frequent seismic activity [60]. The Ordos peripheral fault system comprises of the Shanxi and Yinchuan Grabens to the east and west, respectively, which are characterized by extension and right-lateral shear, and the Weihe and Hetao Grabens to the south and north, respectively, which are characterized by left-lateral shear. In contrast with other grabens that have records of historical M 8 earthquakes (Figure 1(a)), including the 1303 M 8 Hongdong earthquake in the Shanxi Graben [6062], the 1556 M 8.5 Huaxian earthquake in the Weihe Graben [60, 6366], and the 1739 M 8 Yinchuan-Pingluo earthquake in the Yinchuan Graben [60, 67, 68], the Hetao Graben does not contain known records of M 8 earthquakes.

The NBOB comprises of the Yinshan uplift, the Hetao Graben, and the Ordos uplift. The Hetao Graben is 440 km long and 40–80 km wide, and it is orientated E–W between two uplifts. The northern margin of the Hetao Graben is bounded by the Sertengshan Piedmont Fault (SPF), the Wulashan Piedmont Fault (WPF), and the Daqingshan Piedmont Fault (DPF), while the southern margin of the graben is bounded by the Northern Ordos Fault [30, 60, 69]. The Langshan Piedmont Fault (LPF) is located to the west [29, 31, 70], and the Helingeer Fault is located to the east (Figure 1(b)). The normal faults along the western and the northern margins of the graben constitute the boundary fault zone of NBOB and have been highly active during the Holocene and control the main tectonic activity and geomorphological structures in the Hetao Graben [59, 60]. Moreover, these boundary faults can be divided into smaller but relatively independent segments across intersegment zones (known as stepovers or bends) due to their differing geometries, geomorphologies, and kinematics [21, 71, 72]. Previous studies have established the segmentation of the abovementioned four active faults (Table S1). In the Hetao Graben, the Linhe, Baiyanhua, and Huhe Sag are separated by Xishanzui and Baotou uplift (from west to east), with an alternating arrangement of depressions and uplifts [60, 73].

Previous studies have indicated that active deformation still occurs along the NBOB. According to historical records, two large earthquakes may have occurred in 7 BC and 849 CE in the Hetao Graben [74, 75]. In addition, moderate earthquakes, including the 1976 M 6.2 Bayan earthquake, the 1979 M 6.0 Wuyuan earthquake, and the 1996 M 6.0 Baotou earthquake, and other small M 4–5 earthquakes have frequently occurred along the Hetao fault depression zone since the twentieth century, indicating that intense tectonic activity has occurred recently in this region (Figure 1(b)).

3.1. Trench Data Collection

Data were collected from 52 effective trenches excavated from 1988 to 2020 CE to reconstruct paleoearthquakes, investigate long-term seismic processes, and assess seismic risk [2832, 73, 7683]. We collected as much trench data as possible from all available documents, and this collection is the most complete so far. The LPF, SPF, WPF, and DPF contained 14, 12, 6, and 20 trenches, respectively. Based on the completeness of collected data, the trenches were divided into three categories. The trench that is the most complete consists of the profile, detailed profile interpretation, and useful dating data, and it is highly reliable. The moderately reliable trench data consists of detailed profile interpretation and useful dating data without the profile. The trench data with low reliability only has dating data. It should be noted that this classification of data represents only its completeness for usage and does not represent its quality. Among these (Table S1), 42 trenches were highly reliable, while an additional six trenches were moderately reliable. The remaining four trenches had low reliability. More than 90% of the trench data, including the required event indicators and corresponding dating data, is available.

In normal faults, precise ages are determined by either directly dating the soil in the colluvial wedge, dating the soil profile that developed after stabilization, dating the uppermost deposit that underlies the colluvial wedge, or directly dating the deformation feature that displaces the surface [18, 2426, 8487]. The ages collated in this study were obtained using the first three methods stated above. A total of 46 radiocarbon (14C), 23 optically stimulated luminescence (OSL), and 56 thermally stimulated luminescence (TL) dated samples were synthesized in this study (Table S2). The dating methods used for the collected data were all within the applicable range, and all stratigraphic age reversal data were eliminated. The 14C dates of Holocene samples were reliable and consistent with those obtained using OSL dating [88, 89] and are preferred as they have fewer analytical errors. As large inaccuracies might result in 14C dating data that are acquired from typically conventional ages corrected by isotopic fractionation rather than by tree rings and different radiocarbon calibration methods used in previous studies might produce different stratigraphic units ages, we performed consistent recalibration (Table S2) for radiocarbon measurements of all the dated samples using the OxCal 4.4 software with the IntCal 20 curve at a confidence interval of 95.4% (2σ) [90]. To sum up, the collected data has good integrity, strong validity, and high accuracy.

3.2. Improved Multiple Trench Constraining Method

A paleoearthquake can be recognized if the individual wedge where soil had developed on a colluvial unit was subsequently buried by sediment (Figure 2(a)). The number of colluvial wedges can be used to determine the number of events that have occurred [21]; however, events may be missed if a wedge is located in an unfavorable location for paleoearthquake preservation. In general, it is unlikely that a single trench will contain a complete sequence of events [42]. The estimation of the ages of paleoearthquakes in which the upper and lower bounds of events that are interpreted from multiple trenches are arranged by time and the boundaries of some events can be refined stepwise [42, 43] is known as the multiple trench constraining method.

Primary paleoearthquake evidence observed in trenches on normal faults mainly comprises of colluvial wedges and faulted strata produced by tectonic deformation, followed by soil horizons that underlie and overlie each wedge [18]. The age of a fault is lower than that of the highest faulted strata and higher than that of the lowest unfaulted strata (Figure 2(a)). Therefore, a colluvial wedge partially represents the occurrence of a paleoearthquake, whereas faulted and unfaulted strata represent the upper and lower limits of a paleoearthquake, respectively.

The multiple trench constraining method effectively addresses the problems with incomplete and unreliable paleoearthquake sequences; however, the rupture is considered an independent segment, and the length of multiple events on the same fault segment is fixed [21, 4749]. Multiple earthquakes, including the 1887 Mexico Earthquake (Mw 7.4) [91] and the 2008 Wenchuan Earthquake (M 8.0) [92], have indicated that seismic rupture can breakthrough steps and gaps between fault segments [93]. Therefore, we improved the aforementioned method by weakening the segmentation while constraining paleoearthquakes. For each trench location, we determined whether the rupture broke through the bedrock, stepover, or barrier along the fault; that is, we considered that each paleoearthquake was not entirely constrained by the segment but by the evidence of the rupture at the corresponding time. In addition to the original segment location/region, we interpreted the subsegment location/region (smaller scale stepover), which showed that the rupture may have terminated either at the segment or subsegment region of the trench. This improved method reduces the selection bias caused by the characteristic earthquake model, enriches the rupture behavior between different faults, and is conducive to the limitation of time and length of large earthquakes with multisegment rupture and small earthquakes with partial segment rupture.

For this method, we developed three qualifying criteria. First, three trench qualities were used to determine the order of precedence for constraining. Moderately reliable trenches lagged behind highly reliable trenches and took precedence over trenches with low reliability (Figure 2(b)). Second, colluvial wedge data were preferentially considered when constraining the age of a paleoearthquake (Figure 2(c)). Third, among the three main dating methods, 14C-derived ages showed precedence over OSL-derived ages, which in turn had precedence over TL-derived ages (Figure 2(d)). Last but not the least, change in strike, step zone, and bedrock deformation was used to determine the segment or subsegment, which indicated the end of the rupture. A rupture was considered to have occurred within a segment as long as there is no trench along the extension of the fault segment (Figure 2(e)).

3.3. Event Indicator Analysis

The improved multiple trench constraining method was used to comprehensively analyze and constrain the ages of paleoearthquakes. For example, 20 trenches that constrained the events on the DPF were previously reported [33, 60, 75, 78, 8082, 94, 95]. Event D1 was constrained by evidence from three trenches, all of which were dated using the TL dating (Figure 3(a)). In trench D-TC19, the event occurred between 15,310±1,060 yr BP and 13,770±1,060 yr BP. In trench D-TC18, the same earthquake occurred between 15,320±1,180 yr BP and 13,890±1,070 yr BP. Both the two trenches were highly reliable. In the low-reliability trench D-TC2, the earthquake occurred before 13,460±1,040 yr BP. These constraints define the age of event D1 to be 14,600±710 yr BP, which was a full-length rupture (Figures 3(b) and 3(c)). The age of event D5 was determined according to the evidence from nine trenches. The event was constrained well by trenches D-TC4 and D-TC7, wherein the colluvial wedge and 14C dating were used for identification and dating, respectively, and it was indicated that the event occurred at 3,630±260yr BP (Figure 3(a)). Moreover, moderately reliable trenches D-TC3, D-TC9, and D-TC11, in which the colluvial wedge and 14C dating were used, revealed a lower priority than a highly reliable trench. Apparently, the event D5 was a multisegment rupture that occurred from the Tuyouqixi segment to the Huhehaote segment. Furthermore, if a piece of sufficient trench evidence was found in the Baotou segment, event D5 would be a full-length rupture (Figure 3(c)). Because of the absence of early paleoearthquake indicators, evidence is inadequate (Figure 4(a)) and the rupture was that of a partial segment that ended at the subsegment of stepover for event L1 (Figures 4(b) and 4(c)), which was limited by trench L-TC11. These examples show that the improved method constrains the rupture time of the events along a single fault using paleoearthquake indicators and then obtains the rupture length using the locations of trenches. The improved method can obtain the multimode rupture behavior and more accurate rupture time and length with a specific limit when paleoearthquake indicators are more complete.

The identification indicators were initially classified and then used to evaluate the reliability of each paleoearthquake. In previous studies, all evidence within a suspected seismic event layer were evaluated, the highest and average grades were identified, and the total number of identification indicators was calculated for each layer, from which a semiquantitative description of the probability of an event based on the aforementioned analysis was obtained [25, 26]. However, paleoearthquake indicators along the NBOB are mainly colluvial wedges or fault strata, both of which are highly reliable. Therefore, we redefined levels based on the indicator type and the number of trenches (Table 1).

3.4. Earthquake Magnitude and Energy Estimation

Although various types of primary evidence have been used to infer paleoearthquakes magnitudes, the surface rupture length (SRL) on continental fault traces is by far the most commonly used [96, 97]. The inferred SRLs for paleoearthquakes can be compared to worldwide rupture length data for historic earthquakes with known magnitudes to estimate potential paleoearthquake magnitudes. Extensive research has been conducted regarding the relationship between earthquake magnitude and SRL [9799]. Equation (1) can be used with historic data to obtain the paleoearthquake magnitude [97], and it can also be used to estimate the released energy of an earthquake based on the Gutenberg–Richter magnitude-energy relation (Equation (2)) [12].
(1)M=4.86+1.32logSRL,(2)logE=4.8+1.5M,
where M is the earthquake magnitude, SRL is the surface rupture length (km), and E is the energy (J) released by the earthquake. Due to the limitations of paleoseismology, SRL does not correspond to the rupture observed on the surface, but it is the length between the estimated endpoints of a single earthquake rupture indicated by trench evidence.

3.5. Earthquake Recurrence Mode

For paleoearthquakes with three elements (time, magnitude, and rupture range), the earthquake recurrence mode can be obtained using paleoearthquake parameters. The temporal difference between two adjacent events is the time interval Xi, and sorting n+1 events by time yields n intervals (X1,X2,Xn). The intervals between each event can then be used to determine the mean recurrence interval (Ta) and its standard deviation (σ) by Equations (3) and (4), respectively. Furthermore, the ratio of the coefficient of variation (COV) between the standard deviation and the mean recurrence interval (Equations (5))) is used to measure the level of regularity in a dataset [100, 101]. By considering the time error, in the absence of more evidence, paleoearthquakes can be limited to a normally distributed time range, with a median of μ and a time error of 2σ uncertainties. This indicates that the probability distribution function (PDF) of each event is normally distributed. The COV is calculated by dividing the PDF of each event into 10,000 pieces of equal probability and then using a random sampling approach to draw possible event histories for the site, which requires that the events are in chronological order until 1,000 successful event histories are achieved. In this study, we use two models to calculate the mean recurrence interval and COV directly (regardless of error) for model 1 and indirectly (considering error) for model 2. The COV is a simple metric for estimating the variability among intervals and can be used to quantitatively determine seismic recurrence patterns. A COV of 0 indicates perfect periodic behavior, while a COV of 1 is random. COV values between 0 and 1 are called quasiperiodic [52, 100]. A COV>1 reflects clustered behavior, with multimodal variability in the interval length.
(3)Ta=X1+X2++Xnn=i=1nXin,(4)σ=1n1i=1nXiTa2,(5)COV=σTa.

The results of our study indicate that six, seven, eight, and six events were observed on the LPF, SPF, WPF, and DPF (Figures 36), respectively (Table 2). Among these events, L6, S3, S4, S7, W4, D1, D3, and D4 were full-length ruptures, whereas the others were single or mixed segment ruptures. Notably, the number of trenches for events L1 and L2 were too small to constrain the events well, and they may have terminated in a subsegment (overstep that divides the east segment (LE) of the LPF into LE1 and LE2) or may have had shorter rupture lengths (Figure 4(c)). Similarly, event S6 may also have terminated in a subsegment (stepover that divides the Dashetai segment (SD) of the SPF into SD1 and SD2) (Figure 5(c)). Based on the existing evidence, the reliabilities of the paleoearthquakes on all faults, except the WPF, were very high. Possible ruptures were derived from the locations of the trenches, the ensemble magnitude was inferred, and the released energy was then calculated by empirical Equations (1) and (2) (Table 2).

5.1. The Most Recent Earthquake along the LPF and DPF

In this study, we collected as much data as possible to refine the results of constrained paleoearthquakes and improve the completeness of existing records. Based on historical records, two large earthquakes (in 7 BC and 849 CE) may have occurred in the Hetao Graben [74, 75]. Evidence has indicated that the epicenter of the 849 CE earthquake was located near Baotou in Inner Mongolia, but the magnitude is disputed with a potential magnitude in the range of M 6.7–8.26 [60, 8082, 94]. There is little discussion about the 7 BC earthquake [75], and there are doubts about its rupture location and magnitude.

Based on data from seven trenches and samples measured using 14C dating, surface ruptures and dislocated strata, since 2,000 yr, were found in the Baotou, Tuyouqixi, and Tuzuoqixi segments of the DPF (Figure 3(a)). Available paleoearthquake data and evidence from trench D-TC4 indicate that the latest disturbance occurred 1,280 yr BP ago. Evidence of activity has not been identified on other faults in the Hetao Graben. In addition, the earthquake intensity map of the 849 CE earthquake was compiled according to the earthquake deformation trails, building damages using historical records, and the control of fault zones, strata, and landforms on the attenuation of seismic intensity [102]. Both the range of damage and the overall trend indicate that the earthquake occurred near Baotou, and rupture range and time are consistent with the most recent earthquake event D6. Collectively, dislocations along the DPF most likely occurred due to the seismogenic fault of the 849 CE earthquake, which ruptured from the Baotou segment to the Tuzuoqixi segment and was the same as that of event D6, and the empirical Equation (1) gives a magnitude of 7.7 (Figure 3(a)).

According to historical records, the 7 BC earthquake affected a wide area, including Gansu, Ningxia, Shaanxi, Inner Mongolia, Shanxi, Hebei provinces, and Beijing, and it was one of the most destructive earthquakes in history. Due to the uncertainty of historical records, there are differences in understanding the earthquake, and the epicenter is believed to be in the Hetao Graben by some scholars [102]. Based on documented evidence, the strike of the damage generally runs parallel to the trend of LPF. Moreover, existing paleoearthquake evidence (faulted strata, colluvial wedges, and stable strata) indicates that the most recent earthquake event L6 on the LPF was identified in seven trenches with an age of approximately 2,000 yr (Figure 4(a)), which is similar to the 7 BC earthquake in terms of time. Adequate evidence is available to show that event L6 ruptured the entire LPF with a full-length rupture of 160 km. Using the empirical Equation (1), a magnitude of M 7.8 was arrived at for event L6, which is consistent with the 7 BC earthquake based on data obtained from the trenches. Thus, we can tentatively associate events D6 and L6 with the 849 CE and 7 BC historical earthquakes, respectively.

5.2. Recurrence of Large Earthquakes along Single Faults

Earthquake recurrence behavior can be qualitatively described by comparing the recurrence intervals between adjacent events and the mean recurrence interval. The results of paleoearthquake constraining show that six or more surface ruptures have occurred on the four normal faults on the NBOB, which are suitable to calculate the paleoearthquake parameters using Equations (3))–(5)). The mean recurrence intervals calculated by the two models in the present study were the same. Model 1 yielded mean recurrence intervals and standard deviations of 2,290±580, 1,790±820, 1,680±910, and 2,700±340yr, for LPF, SPF, WPF, and DPF, respectively (Figure 7). These values were determined by all the observed surface ruptures, including those that were unlikely expressed. Notably, the intervals on the DPF were more concentrated, followed by those on the LPF, whereas the WPF exhibited relatively random intervals. The interval (1–2) of SPF is 3,440 yr long (13,920–10,480 yr BP), almost two to three times longer than those of the others. This particular value may represent an early missed event due to the absence of event indicators when all other intervals are continuous and similar.

At the same time, the recurrence behavior can be described more quantitatively using COV. COVs of LPF and DPF were 0.25 and 0.12 in model 1 and 0.28 and 0.16 in model 2, respectively. However, a reversal may have occurred on SPF and WPF, for which the COVs in model 1 were 0.46 and 0.54 and those in model 2 were 0.45 and 0.52, respectively. In general, due to the compensation for temporal errors in model 2, its COVs were larger than those of model 1. The COVs of the four faults for the ground-rupturing earthquakes since 15,000 yr were <1, which is consistent with the quasiperiodic behavior even though the values differed marginally (Figures 7 and 8). The interval and frequency histogram of events calculated by random sampling in model 2 showed that the distribution curve of the DPF with good periodicity had one peak, which is close to the average recurrence interval, and it is similar to that of the normal distribution, and in contrast, the WPF with poor periodicity has complex peaks in the distribution curves (Figure 8), and the other two faults have better periodicity with two peaks. Nevertheless, a dominant peak was observed for each fault, which is similar to the recurrence time that corresponds to the median of the normal distribution (Figure 8). Furthermore, the quasiperiodic earthquake recurrence conforms to the elastic rebound theory, wherein faults rupture in response to a specific amount of accumulation strain [27, 34, 46, 103]; however, the strain accumulation process can be disturbed by a complex tectonic setting, resulting in a recurrence interval that is shorter or longer than the average.

Characteristics that may be good indicators of quasiperiodic recurrence on faults can be identified based on the periodicities exhibited by the four faults. The linear regression to fit the COVs and the mean magnitudes of the paleoearthquakes on the four faults defined in this study show a negative correlation. A weak link was observed for approximately 15,000 yr, with R2<0.6. Nevertheless, R2>0.9 was observed since ~10,000 yr (Figure 9). When the magnitude characteristics of the earthquakes were considered, the magnitudes of the DPF earthquakes were larger than M 7.5, with an average of M 7.9, and a significant periodicity. The magnitudes of the earthquakes on the WPF were < M 7.5, with an average of M 7.1 and poor periodicity. SPF and the LPF had uneven magnitude distributions; however, their average earthquake magnitudes were negatively correlated to the COVs (Figure 9). Available evidence presents further proof that larger earthquakes on the DPF and SPF were relatively more periodic.

Regional geology shows that the active tectonic block boundary of the Hetao Graben along the NBOB is located in several depressions, including the Linhe, Baiyanhua, and Huhe depressions. SPF and DPF correspond to the largest (Linhe and Huhe) depressions, respectively, indicating that these two faults may dominate the region. Previous research has found that larger earthquakes occur where multiple segments break; however, their frequencies are lower when more strain is relieved [10]. The frequency of earthquake occurrence is also related to the regional structures. For example, the recurrence intervals of the San Andreas Fault [34] and Altyn Tagh Fault [84] are both on a centennial scale, while the Central Longmen Shan Thrust Belt [104] and the Haiyuan Fault [24, 105] have recurrence intervals on a millennial scale with relatively large magnitudes; however, their recurrence frequencies are relatively different. Geophysical evidence indicates that the NBOB is affected by (from east to west) the compression of the northeastern margin of the Qinghai-Tibet Plateau, the upwelling of hot material, and the deep subduction of the Pacific plate [60], which provided the stress conditions necessary for the initiation and occurrence of strong earthquakes in this region. The earthquake recurrence cycle occurs on a millennial scale with obvious periodicity. The results indicate that larger earthquakes produced by primary faults exhibited better periodicities with lower COVs, whereas secondary faults were affected by primary faults in the same tectonic zone and produced smaller earthquakes, thereby exhibiting poor periodicities and higher COVs. This conclusion is different from previous studies in that the larger the magnitude, the longer the recurrence period. The recurrence period, that is, the frequency of the earthquake, depends on the complexity of the geological background.

5.3. Large Earthquake Cluster along the NBOB

Although we determined that the recurrence behaviors of these four faults are quasiperiodic, the absence or errors in the timings of earthquakes highly affect their behavior. If the outlying interval of 1–2 on the SPF is removed, the recurrence interval is reduced to 1,460±140 yr, with a COV of 0.09 (as estimated in model 1) (Figure 7). Thus, other characteristics that may be effective for assessing earthquake risk should be identified. Six earthquake clusters in the NBOB with a cluster interval of approximately 1,330±300 yr (since approximately 15,000 yr) were divided by integrating a time window of 500 yr and the conditional probability distribution curve of the event modeled times (Figure 10). Clusters 1 (C1), 2 (C2), and 3 (C3) comprise multiple temporal groups. Moreover, seismogenic markers have presumably degraded due to denudation, or disturbances in other strata may have influenced the dates of earlier strata, thereby degrading the older seismogenic event clusters. By analyzing the paleoearthquakes in different clusters, we found that each cluster contained surface ruptures of the LPF, SPF, WPF, and DPF except for C6, which indicates that this cluster may not have terminated yet. Earthquakes that occurred during the clusters were dissimilar from sequence to sequence in terms of their locations, magnitudes, and rupture sequences. However, if only SPF and DPF were considered, cluster migration occurred in two types, the westward-propagating (DPF–SPF type and DS type: C1, C2, and C4) and the eastward-propagating (SPF–DPF–SPF and SDS type: C3 and C5).

In addition to the characteristics mentioned above, it is possible to compare the integrity of different clusters through the accumulated energies released by regional earthquakes. In Figure 10, the x-axis indicates time and the y-axis indicates accumulated energy; the total cumulative energy released by the earthquake is 7.7×1016 J, 5.1×1016 J, 15.0×1016 J, 8.8×1016 J, 10.4×1016 J, and 4.9×1016 J for C1, C2, C3, C4, C5, and C6, respectively. The average energy released by the different earthquake clusters was 8.6×1016 J (Figure 10). C2 and C6 are significantly below average. The events in C2 are relatively complete, but the rupture length of events may be shorter due to the lack of adequate evidence, resulting in an assessment of a smaller magnitude. Comprehensive analyses for C6 indicate incomplete surface rupture and low energy release, suggesting that this cluster may not have terminated yet. Among these, the main controlling faults, SPF and DPF, released more energy, whereas the secondary faults, LPF and WPF, released moderate amounts of energy.

Clusters are widespread in the study area, but the migration and rupture patterns in the clusters are different. The material and segmentation of the fault zone may play a pivotal role in controlling rupture propagation and termination [106], while the sequence of individual earthquakes in a specific cluster may result from a combination of stress evolution along the fault system [107]. It can be interpreted that the occurrence of adjacent earthquakes can trigger or inhibit subsequent earthquakes. The 1992 M 7.3 Landers earthquake led to the early occurrence of the 1999 M 7.1 Hector Mine earthquake due to the viscoelastic loosening of the lower crust and upper mantle after the Landers earthquake [108, 109]. The 2008 Mw 7.9 Wenchuan earthquake triggered the 2013 Mw 6.6 Lushan earthquake, but the rupture gap between the two earthquakes was inhibited [110]. If the viscoelastic medium model and fault dislocation model of the study area are known, we can use GeoTaos based on PSGRN/PSCMP [111] to calculate the Coulomb stress changes of coseismic and postseismic events [110]. However, it is difficult to obtain the slip angle and coseismic horizontal dislocation of events, which makes it difficult to establish a fault dislocation model. Moreover, the complicated geological background made the changes of the cluster more difficult to understand. More evidence and data are needed to support the interaction of static stress changes between paleoearthquakes.

5.4. Seismic Risk Assessment

We considered a richer set of seismic rupture modes, in which the rupture length was not fixed; however, the temporal characteristics were similar to those of characteristic earthquakes. Notably, the simple renewal model of an elastic rebound-driven earthquake cycle requires expansion to accommodate variations that span multiple earthquake cycles. Existing data can be explained by highly periodic temporal ruptures (Figure 11, explanation 1: characteristic earthquake model), wherein the dominant peak replaces the mean rupture recurrence interval (explanation 2: considers outliers based on the characteristic earthquakes) and random ruptures are generated by meeting parameter conditions (explanation 3: conforms to the actual scenario). Each model implies a different risk level in the NBOB and supports a different physical faulting model. As mentioned above, C6 may still be continued. Explanation 1 indicates that the SPF was ruptured, whereas explanations 2 and 3 suggest that the clustering has finished (Figure 11). Thus, an earthquake may have already occurred; however, no evidence of the event is available. Moreover, based on the earthquake catalog, an earthquake with a magnitude of approximately 63/4 occurred in the NBOB in 1934 CE, and analyzing its relationship with predicted events may help evaluate the risk of an earthquake along the SPF. In any case, C7 did not show any potential indication for an earthquake to occur along the DPF, which suggests that the activity of the DPF may be decreasing, and it may have a low seismic risk. LPF and WPF are at a high seismic risk to compensate for their accumulated stress. Collectively, the analysis indicates that the most urgent degree of seismic risk is along the SPF, wherein possibly a full-length rupture with a magnitude of 7.8 might occur. This might be followed by a multisegment rupture with a magnitude of 7.4 and a full-length rupture with a magnitude of 7.8 in WPF and LPF, respectively. The seismic risk along the DPF has been relatively low since 849 CE.

This study provided novel insights on the recurrence behavior of the NBOB by integrating evidence obtained from paleoearthquake trenches along its boundary faults since approximately 15,000 yr. Our conclusions are as follows:

  • (1)

    In general, six, seven, eight, and six events have been constrained on the LPF, SPF, WPF, and DPF, respectively, in which surface ruptures may have occurred as full-length, single-segment, and multisegment rupture events. The most recent paleoearthquakes along the DPF and LPF correspond to the 849 CE and 7 BC earthquakes, respectively

  • (2)

    For a single fault, the recurrence behaviors exhibited significant periodicity with recurrence intervals of 2,290±580, 1,790±820, 1,680±910, and 2,700±340 yr for the LPF, SPF, WPF, and DPF, respectively. The behavior of the regional fault network along the NBOB was clustered, with six clusters that had a recurrence time of approximately 1,300 yr. Two cluster migration types were observed: the westward-propagating DS type and the eastward-propagating SDS type

  • (3)

    According to the earthquake recurrence behavior and cluster migration described in the present study, earthquakes on the SPF were noticeably absent and no large earthquakes appeared to have occurred on the NBOB since the 849 CE earthquake. Therefore, the risk of a future large earthquake along the SPF should be considered

The data used to support the findings of this study are included within the article.

The authors declare no conflict of interest.

This work is supported by the National Key Research and Development Program of China (2017YFC1500104), National Natural Science Foundation of China (41774049, 42174062, and 41972228), and Guangdong Province Introduced Innovative R&D Team of Geological Processes and Natural Disasters around the South China Sea (2016ZT06N331). We gratefully acknowledge Yanxiu Shao, Xingwang Liu, and Gan Chen for their assistance and support.

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).