Linear geomorphic markers (e.g., rivers, mountain ridges, and terrace risers) crossing a fault are effective records of the cumulative offset of a fault, from which the offset history can be reasonably inferred (Klinger et al., 2011; Korjenkov et al., 2012; Ren et al., 2015; Tibaldi et al., 2015; Jiang et al., 2017). We analyzed Google Earth images to determine the surficial fault traces in the study area and used a DJI Phantom 4 RTK UAV to map areas with offset gullies in the field. The UAV is connected to a continuous operational reference system (CORS) station by a 4G communication system. The UAV is equipped with a 20-million-pixel camera with a complementary metal oxide semiconductor (CMOS) sensor and a GPS that captures photographs with a 3-cm/pixel accuracy at 120 m above the ground. Each photograph provides a high-precision geographic position (vertical 1.5 cm + 1 ppm and horizontal 1 cm + 1 ppm, where 1 ppm indicates an error of 1 mm over 1 km of movement) via real-time kinematic (RTK) technology. Orthoimages and a digital elevation model (DEM) of these areas were derived using Pix4Dmapper photogrammetry software based on SfM survey photographs with a vertical and horizontal accuracy of ∼ 6 cm. Subsequently, gully offsets were measured using LaDiCao_v2, a professional analysis software program developed by Zielke (Zielke et al., 2010; Zielke and Arrowsmith, 2012; Zielke et al., 2015) on the MATLAB platform. Considering the offset measured in the field, the measurement error range of LaDiCao_v2 was within the 95% confidence interval.

Seismic waves can be produced by an artificial seismic source and propagated in media, such as strata and rocks. Differences in the elasticity of media cause changes in the characteristics (e.g., the speed, path, frequency, and strength) of the seismic waves as they propagate. An analysis of these differences can be used to estimate the parameters (e.g., properties, structure, and geometric locations) of the subsurface strata or rocks. The seismic reflection method was adopted in this study to ascertain traces of buried faults. The main equipment used consisted of an Aries 2.66 digital seismograph and a 20-ton controllable vehicle-mounted vibrator. The observation system used 66 times coverage and was equipped with 400 receiver channels, with a 2-m channel spacing and a 6-m shot spacing. The seismic source was triggered in the middle of the study area.

Ground penetrating radar imagery (GPR) was also a kind of high-resolution geophysical method (Daniels, 1996) and mainly used to investigate the locations and activities of the buried fault in recently studies (McClymont et al., 2009; Roberts and Raithatha, 2010; Beauprêtre et al., 2012). Electromagnetic waves were transmitted by GPR to penetrate subsurface media. They will penetrate to different depths in different media because of variations in the dielectric constants and electrical conductivities of media. Therefore, the characteristics of the subsurface rock-soil masses can be indirectly inferred from the reflection times of electromagnetic waves recorded by a receiver (Daniels, 1996). A third-generation MALA Professional Exploration digital radar system with MALA RTA25 antennas (25 MHz) was used in this study. This system can image a subsurface area within 50 m of the surface in detail. The point measurement method was adopted with a 50-m spacing between survey points.

The dip slip between two walls of a fault could result in the vertical offset of existing stratums along the fault and then an unequal thickness layer on each side of the fault due to the rapid deposit in the negative terrain (Schwartz and Coppersmith, 1984; McCalpin, 1996). Thus, seven boreholes were dug across the S₃ fault to obtain the joint-drilling geologic section and analyse the vertical displacement of the fault from the difference in the depths of layers. The uppermost point was determined based on the age of the latest strata, which is dislocated. An organic carbon sample was extracted from borehole D2 and was dated using ¹⁴C at the Beta Analytic test laboratory (USA). The radiocarbon age was calibrated using BetaCal 3.21 with the IntCal 13 atmospheric model.

Fault scarps and offset gullies along the fault are visible in Google Earth (Figures 2A,B). Field geological survey results show that the NAJF trends 20–30° N and begins at Juxian in the south, extending approximately 135 km to Anqiu in the north. There are several distinct outcrops along the fault (e.g., Figure 2C). Our survey showed that the NAJF consists of three discontinuously exposed fault segments that form right-stepping faults (Figure 1B). Between the exposed fault segments, the NAJF is buried in river deposit plains with relatively thick sediments.

High-resolution images of the offset streams in representative areas of the NAJF were captured by the UAV. Then high-accuracy orthoimages and DEM (Figure 3) data were obtained from these images. All the gully offsets were interpreted in a high-resolution DEM (Figure 3), verified through orthoimages, and measured by LaDiCao_v2 software. Along the fault traces, almost every offset gully was identified on both sides of the fault scarps with dextral strike-slip. Finally, 87 dextral strike-slip offsets of gullies were collected, including five data points from a previous study (Wang et al., 2015) (see Appendix Supplementary Table SA), ranging from 4.0 (−0.9/+0.1) to 50 (+2/−0.5) m (Figure 4).

The southern stepover is located between S₂ and S₃ and north of Juxian (Figure 1B). Few paleoearthquake relics have been preserved in this area because of human activities. Thus, two seismic sections, I-I’ (1700-m long) and II-II’ (3326-m long), were placed in the survey area to locate the buried fault (Figure 5A).

These two seismic sections reveal unequivocal signs of the fault (Figures 6A,B): in both sections, there is a high-energy reflection layer at depths of 20–150 m (T₁), which is vertically offset at a site at approximately 800 m in section I-I′ and at sites at 425 and 2,325 m in section II-II’. These phenomena were inferred to result from fault activities. The fault in section I-I′ and the left-branch fault in section II-II′ are S₃, dipping to the west at ∼ 65°; the right-branch fault in section II-II’ was inferred to be S₂ and dips to the east at approximately 70°. Previous geological data (Wang et al., 2015) show that T₁ is the boundary between Quaternary sediments and bedrock. In addition, a fault outcrop was found at approximately 14 km from the survey line (Figure 5B), where there is clear fault contact between the Cretaceous amaranth sandstone and Cretaceous yellow sandstone with a high dip angle (65°), further demonstrating the reliability of the sections.

Seven boreholes were drilled to reveal the uppermost point of S₃ (Figure 4E). Based on the deposit thicknesses above the bed rock, the borehole drilling depths varied from 20 to 39 m. First, the fault plane was found in the core of borehole D5 at a depth of ∼ 32 m (Figure 6F). The joint-drilling geologic section showed that all strata at the borehole D3 site are higher than those at the borehole D5 site by 1–7.8 m and that the offset of the strata increases gradually from top to bottom, suggesting that the main fault is located between boreholes D3 and D5 and that older strata may have recorded more paleoearthquake events. The light-black clay stratum U2 is the latest offset stratum at the borehole sites and is a Holocene stratum formed at approximately 9,740 ± 30 cal BP according to ¹⁴C dating (see Appendix Supplementary Figure SA). Combined with the previous paleoearthquake research (He et al., 2005), S₃ was thus inferred to be active during the Holocene.

Five GPR survey lines (L1-L5) were placed along the Qu River plain (Figure 5C). Survey lines L3 and L5 were placed along the line extending from the northern end of S₂, and survey lines L1, L2 and L4 were placed along the line extending from the southern end of S₁ (Figure 5C).

The two GPR sections in S₂ both yielded relatively good stratification information (Figure 7), including two clear stratigraphic boundaries (T₀ and T₁). Based on a previous study (Wang et al., 2015), T₀, at a depth of 2–4 m, represents the bottom of the Holocene sediment and corresponds to U2 in the drilling section (Figure 6E), and T₁, at a depth of 10–20 m, represents the bottom of the Quaternary deposits above the bedrock. The signals of T₁ reflected by the upper media are relatively disordered compared with those of T₀, and there are marked differences across the depth of T₁—the west side of T₀ is ∼ 10 m higher than the east side, from which the sudden changes in T₁ (at the 100-m site on L3 and the 200–300-m site on L5 in horizon) were inferred to be fault locations. In addition, the top surface of bedrock gradually lowers from west to east, consistent with the topography for a transition from hills to plains. Moreover, T₀ breaks at the fault locations on L3 and L5, which suggests that the uppermost point of the fault may be located ∼ 2 m underneath the surface and that S₂ of the NAJF dislocated the Holocene strata.

The GPR sections in S₁ (Figure 8) show two similar reflection layers, T₀ and T₁, representing the bottom of the Holocene sediment and the bottom of the Quaternary deposits above the bedrock, respectively. The fault is also located where there are clear changes in the layer depth (the 240–260-m site on L1, 290–360-m site on L2 and 200–250-m site on L4 in horizon). In addition, the burial depths of the Quaternary sediment on both sides of the fault of 17 and 21 m, respectively, the height of 4 m is significantly larger than that of the Holocene sediment, indicating that the fault has been continuously active since the Quaternary. These results were corroborated by the discovery of a fault outcrop (the fault plane dips to 292° at ∼ 62°) in S₁ within the basin (Figure 5D).

Characteristic earthquakes refer to multiple earthquakes that are similar in rupture length, offset distribution and magnitude and occur periodically on a fault during its long-term active period (Schwartz and Coppersmith, 1984). The seismic activity of the AJF has been characterized by primarily strong earthquakes with similar magnitudes and a lack of moderate and small earthquakes, where a linear relation with a low-b-value between high-magnitude earthquakes has been found (Chao et al., 1994). In studies on the SAJF, a characteristic displacement of approximately 9 m has been obtained (Jiang et al., 2017), with a recurrence interval of 3,000–4,000 a and a magnitude of ∼ 8.5 (Lin and Gao, 1987; Huang, 1993; Wang, 1996; Chao et al., 1997). Therefore, the earthquake events of the SAJF may obey the characteristic earthquake model.

In this study, we obtained 87 offset values, all of which were projected onto the NAJF based on distance and statistically analyzed to determine the cumulative offset probability distribution (COPD) (Figure 4). Large offsets tend to correspond to a long activity history, and offsets >25 m do not have statistical significance because of the large time error involved. The coseismic displacement of the NAJF has similar features to that of the SAJF. The horizontal offsets are mostly concentrated in four intervals, for which the COPD peaks are 5.1, 10.3, 15.8, and 20.3 m along the NAJF (Figure 4), i.e., multiples of ∼ 5 m. Despite the relatively few data points, the COPD in S₃ exhibits two similar offset peaks at ∼ 5 and ∼ 10 m. These data are consistent with previous research showing that the gully offsets on the AJF are concentrated at 4–9, 14–20, and 28–34 m (Huang, 1988). Among these offsets, the largest density distribution corresponds to the 5 m offset, and the number of gullies gradually decreases as the gully offset increases, suggesting that the smallest offset of 5 m represents the latest earthquake of the NAJF and that 5 m may be the coseismic displacement of this event. Moreover, the larger offsets represent cumulative displacements associated with the relatively early earthquake events of the NAJF. Thus, similar to the SAJF, the NAJF may also have a characteristic offset of ∼ 5 m.

The rupture segmentation of a fault can be used to evaluate its future seismic risk to serve as an important reference (Ding, 1992; Ding, 1993; Ding, 1995). The NAJF has previously been separated into S1–S3 segments based on the spatial distribution of the exposed faults (Wang et al., 2015). The geophysical detection results show that S₁, S₂, and S₃ spatially form right-stepping faults, and relatively small widths (approximately 1 km) were found for the stepovers among these faults in this study. It is generally believed that stepovers wider than 3–5 km (Wesnousky, 2006; Wesnousky, 2008; Wesnousky and Biasi, 2016), and even those over 8 km (Wang et al., 2018), obstruct surface rupture propagation. Thus, the stepovers of the NAJF theoretically will not affect the overall rupture segmentation. This conclusion is also corroborated by our results. The coseismic offset and cumulative offset values for the three secondary faults are similar and are all multiples of approximately 5 m, suggesting similar rupture behavior for these segments.

The following empirical equations for the magnitude, coseismic displacement (D_co), and surface rupture length (SRL) have been established from the previous study (Wells and Coppersmith, 1994) on strike-slip faults: logSRL = a + blogDco (1) M = a + blogDco (2) where M is the magnitude, SRL is the surface rupture length, D_co is the coseismic displacement caused by a single seismic event, and a and b are parameters.

The possible SRLs calculated using a D_co of 5 m for a single seismic event (Table 1) range from 93.5 to 115 km, averaging 103.3 km. The data in Figure 4 show that the SRLs in S₁, S₂, and S₃ are 45, 73, and 42 km, respectively, which are all far shorter than 103.3 km. However, the combined length of these three segments of 130 km is relatively consistent with the calculated value. Thus, we infer that the earthquake that led to a D_co of 5 m was caused by a cascade rupture of the three faults. The rupture segmentation method (Ding, 1992; Ding, 1993; Ding, 1995) shows that the three secondary faults are combined into one segment.

Some debate remains regarding the seismogenic fault and magnitude of the Anqiu earthquake in 70 BC (Li, 2014; Wang et al., 2015). The epicenter of the Anqiu earthquake was located near the NAJF based on the earthquake intensity and the disaster distribution estimated from the “History of the Han Dynasty” records (Li, 1981; Zhu and Sun, 1991) (Figure 1B). Therefore, although no geological evidence has been found to verify the seismogenic fault of this earthquake in the NAJF, most researchers consider that this earthquake probably ruptured the NAJF. In this study, the D_co on the NAJF was found to be approximately 5 m and used to calculate the magnitude of the Anqiu earthquake (Table 1). Most results show that the magnitude of the Anqiu earthquake was above 7.5, which is far higher than the historical record of M7 (Department of Earthquake Damage and Defense, CEA., 1999). In addition, historical records (Li, 2014) show that the Anqiu earthquake was felt sufficiently strongly in Xi’an City, the capital of the Han Dynasty, located ∼960 km to the west of Anqiu, that the emperor had to leave the royal palace. By consulting the earthquake intensity scale developed by the China National Standardization Management Committee (The China National Standardization Management Committee & General Administration of Quality Supervision, Inspection and Quarantine of the People’s Republic of China, 2008), the intensity of the Anqiu earthquake in Xi’an City is inferred to have been IV–V. To verify the calculated results, the intensity distributions of several recent earthquakes with magnitudes ≥7.0 in North China, including the 1976 Tangshan M7.8 earthquake, the 1966 Xingtai M7.2 earthquake, the 1975 Haicheng M7.1 earthquake, and the 1937 Heze M7 earthquake, were compared with that of the Anqiu earthquake (Table 2). The results show that the intensity distribution of the Anqiu earthquake was far higher than that of the Heze M7 earthquake and close to that of the Tangshan M7.8 earthquake. Thus, we propose that the magnitude of the 70 BC Anqiu earthquake has been underestimated and was above M7.5. As the seismogenic structure of the Anqiu earthquake remains in dispute (Li, 2014; Wang et al., 2015), we have only presented one possible deduction based on empirical formulas for the magnitude, D_co and SRL, and more data are needed to verify this inference.

Based on previous research (Lin and Gao, 1987; Huang, 1993; Wang, 1996; Chao et al., 1997; Jiang et al., 2017), the SAJF has a recurrence interval of 3,000–4,000 a and a characteristic displacement of ∼ 9 m, suggesting that the SAJF follows the characteristic earthquake model. The results of this study show that the NAJF likewise has a D_co of ∼ 5 m. Although the recurrence interval remains unclear because of insufficient research, the latest earthquake on the NAJF may have occurred in 70 BC with a magnitude >7.5. Obviously, there are large characteristic differences between these two AJF segments. Moreover, both the fault structure and numerical simulation results indicate that the Juxian Basin between the NAJF and the SAJF is the permanent termination point for rupture propagation. The stepover between the NAJF and SAJF is approximately 8 km wide, which exceeds the general width (5 km) that allows rupture propagation (Wesnousky, 2006, Wesnousky, 2008; Wesnousky and Biasi, 2016), and the surface rupture of the Tancheng M8.5 earthquake in 1,668 also did not break through this stepover (Jiang et al., 2017). The simulation results of GPS data (Li et al., 2020) show a low probability of simultaneous rupture of both faults because the NAJF and the SAJF have different fault slip deficit rates and are in different locking states. This result indicates that the Anqiu–Juxian fault can be reliably divided into two permanent segments for earthquake rupture: the NAJF and the SAJF.