Dynamics of the 1057 M 6 earthquake rupture and seismic hazard implications for the Beijing region, China

Xie, Zhangdi

doi:10.3389/feart.2025.1669495

ORIGINAL RESEARCH article

Front. Earth Sci., 17 October 2025

Sec. Solid Earth Geophysics

Volume 13 - 2025 | https://doi.org/10.3389/feart.2025.1669495

Dynamics of the 1057 M 6 earthquake rupture and seismic hazard implications for the Beijing region, China

Zhangdi Xie^1,2*

¹Beijing Gangzhen Technology Co., Ltd., Beijing, China
²Key Laboratory of Optoelectronic Technology and Systems (Ministry of Education of China), College of Optoelectronic Engineering, Chongqing University, Chongqing, China

Seismic activity frequently occurs in the North China Plain seismic belt. Historical earthquake records for the fault southwest of Beijing are sparse, with only one significant event: the 1057 M 6 earthquake. This study focuses on the southern extension of the Nanyuan-Tongxian Fault. Using dynamic numerical methods, we simulate the rupture process of this earthquake to investigate its impact on strong ground motions in Beijing and surrounding areas. First, on the basis of the fault geometry and hydraulic fracturing stress data, we set the dynamic parameters of the fault plane and construct a curved fault model. Numerical simulation of the dynamic rupture process yielded the following results: a maximum slip of 0.75 m, an average slip of 0.49 m, a rupture duration of approximately 10 s, and a moment magnitude (M_w) of 6.5. The rupture propagated bilaterally; influenced by fault geometric bending, the peak slip rate reached 0.7 m/s on the SW side, whereas the slip rate on the NE side was smoother. The strong ground motion distribution showed a maximum intensity of VIII, which was concentrated along the fault trace with a significant hanging wall effect. Notably, owing to the combined effects of rupture directivity and fault bending, anomalously high-intensity zones appeared in areas distant from the epicenter, such as Shunyi and Tongzhou. This pattern differs significantly from the historical seismic impact pattern of the Sanhe-Pinggu Fault on the northeastern side. Further simulations of the epicenter within the fault indicate that if a future earthquake’s epicenter is located west of the fault center, Tongzhou would still face high-intensity ground motion risks. The results from both planar-fault and 3-D velocity-model simulations further demonstrate that strong ground-motion characteristics are influenced by multiple factors. This study reveals the unique seismogenic mechanism of the Nanyuan-Tongxian Fault, providing a key basis for seismic hazard assessment southwest of Beijing.

Plain language summary

This study investigates the 1057M 6 earthquake in the Beijing region via dynamic numerical simulations to understand its rupture process and potential seismic hazards. This research focused on the southern extension of the Nanyuan–Tongxian Fault and constructed a curved fault model on the basis of fault geometry and stress data. The simulation results reveal a maximum slip of 0.75 m, an average slip of 0.49 m, and a rupture duration of approximately 10 s, with a moment magnitude (M_w) of 6.5. The rupture propagated bilaterally, with significant slip rate peaks on the southwest side due to geometric bending. The strong ground motion distribution indicates a maximum intensity of VIII along the fault trace, influenced by rupture directivity and fault geometry. It is noteworthy that areas far from the epicenter, including Shunyi and Tongzhou, recorded unusually high ground-motion intensities as a result of the combined effects of rupture directivity and fault bending, and fault geometry and velocity model significantly controls rupture propagation and ground motion distribution. This study highlights the importance of fault geometry in controlling rupture propagation and ground motion distribution, providing critical insights for seismic hazard assessment in the Beijing area.

Keypoints

The 1057M 6 earthquake rupture process was simulated via dynamic numerical methods, revealing significant slip and rupture characteristics.

Strong ground motion results in high intensity in areas far from the epicenter due to rupture directivity and fault geometry.

Fault geometry and velocity model significantly controls rupture propagation and ground motion distribution, highlighting the need for detailed fault studies in seismic hazard assessments.

1 Introduction

The North China Plain seismic belt has experienced multiple earthquakes above magnitude 6, particularly in the Beijing area, where active faults are densely distributed and where seismic activities are frequent. Historical seismic records exist for most active faults. For example, the Nankou–Sunhe Fault, with a dip angle of 70°–80° (Che, 1994; Liang, 2019), has paleoearthquake records dating back 4,000 years (Jiang et al., 2001). Since the 294 AD Yanqing M 6.0 earthquake, M ≥ 6.0 earthquakes have successively occurred in Daxing, Yanqing, Tongxian, and western Beijing (Gu, 1984), with the most recent being the 1730 western Beijing M6.5 earthquake (Gu, 1984; Xie and Cai, 1983). Subsequently, no M ≥ 6.0 earthquakes have been recorded in the Beijing area (Jiang and Yin, 2021). Related seismic and geophysical research (He et al., 2020; Jiao et al., 2005; Li et al., 2021; Qin et al., 2014; Wu et al., 2013; Zhang et al., 2017a; Zhang et al., 2017b; Zhang and Zhu, 2019) has revealed the seismic geological conditions of the Beijing area, providing a solid foundation for analyzing its historical earthquakes and potential future large earthquakes.

To investigate the seismogenic mechanisms and rupture processes of historical earthquakes in Beijing and its surrounding areas, some studies have employed numerical methods to simulate the rupture process and the resulting strong ground motions, such as the 1,679 Sanhe–Pinggu M8.0 earthquake (Ba et al., 2022; Fu et al., 2015; Gao et al., 2002; Liu et al., 2007; Pan et al., 2006; Wang et al., 2022). They have also studied the effects of ground motions generated by different strata, topographies, and surrounding earthquakes in the Beijing area (Ding et al., 2004; Liao, 2017; Yang et al., 2008). The results of these studies indicate that areas in and around Beijing, such as the Babaoshan–Huangzhuang–Gaoliying Fault and the Yanhuai Basin, still possess the potential for significant seismic activity (Jiao et al., 2005; Qin et al., 2014; Wu et al., 2013; Yang et al., 2008; Zhang and Zhu, 2019). The above research focused primarily on faults in the northeastern part of Beijing, which have frequent historical earthquake records. In contrast, the southwestern part of Beijing has fewer historical earthquake records and is characterized mainly by the development of the eastern margin fault of the Daxing Uplift (Daxing Fault), the Nanyuan-Tongxian Fault, and the Yongdinghe Fault. The eastern margin fault of the Daxing Uplift (Daxing Fault), located in the intermediate zone between Beijing’s Daxing District and Langfang city, Hebei Province, is also known as the northern margin fault of the Gu’an Basin (Gao and Ma, 1993). It formed under the Paleogene SE‒NW extensional stress regime of the North China Plain (Zhao and Liu, 2002) and has been notably active since the Quaternary (Li et al., 1985). This fault is approximately 90 km long and 20–40 km wide, trending NE and dipping SE. Its northern section connects to the Xiadian Fault, which generated the 1,679 Pinggu M 8.0 earthquake, and the junction area experienced the 1536 M 6.0 earthquake (Li et al., 2021). Simultaneously, this fault separates the Daxing Uplift from the Langgu Depression, controlling and influencing the tectonic evolution of the Langgu Depression (Gui et al., 2012).

North of the eastern margin fault of the Daxing Uplift lies the Nanyuan–Tongxian Fault, which is approximately 110 km long. This fault serves as the boundary between the Daxing Uplift and the Beijing Depression (Zhao and Zhu, 2003). It is S shaped, dips NW at angles of 50°–70° (Jiao and Qiu, 2006; Xu et al., 2015), and is divided from north to south into three segments: the northeastern segment, central segment, and southwestern segment (Lei et al., 2021; Zhao et al., 2021). Its southern segment extends to the Nankou–Sunhe Fault. Because the historical earthquakes on the series of parallel faults northeast of Beijing (Huangzhuang--Gaoliying Fault, Shunyi-Liangxiang Fault, Nanyuan-Tongxian Fault, and Xiadian Fault) occurred less than 400 years ago, whereas the faults southwest of Beijing (Yongdinghe Fault, southern extension of the Nanyuan-Tongxian Fault, and eastern margin fault of the Daxing Uplift) had only one recorded M 6.0 earthquake in 1,057, the accumulated regional stress in this area may have been released less than that in the northeast. Moreover, the specific location and seismogenic structure of this earthquake remain undetermined (Li et al., 2021). To investigate the rupture process of the 1057 M 6.0 earthquake on the Nanyuan–Tongxian Fault and its potential impact on strong ground motions in the Beijing area, a curved fault model is constructed here, referencing the southern extension of the Nanyuan–Tongxian Fault (Figure 1). Spectral-element methods are employed to simulate the rupture process under different scenarios and calculate the resulting ground motions and their effects on Beijing and surrounding areas.

Figure 1

Map highlighting seismic activity near Beijing. Red lines indicate fault lines, with a marked epicenter labeled

Figure 1. Active faults and focal mechanisms of historical earthquakes in Beijing and surrounding areas. The blue line indicates the southern extension of the Nan Yuan–Tongxian fault, which has a NE strike, SE dip, and dip angle of 70°. The red circle represents the reference epicenter location of the 1057M6 earthquake. The active fault data are from (Xu et al., 2016); The focal mechanism solutions are taken from the 2009–2021 dataset for mainland China and adjacent regions (Guo et al., 2022).

2 Data and methods

2.1 Fault geometry

During the Eocene–Oligocene extensional rifting in the Cenozoic Neogene, a series of NE-trending normal faults developed in the Beijing region (Huang et al., 1991). Among them, the Nanyuan-Tongxian fault strikes NE, dips NW, and has a dip angle of 50°–75°. It extends from the vicinity of Nan’gao village in the south, passes through Nanyuan, Nanmofang, east of Dingfuzhuang, and reaches Ping’an Tong in the north, with a total length of 110 km. In the Tertiary, this fault served as the boundary between the Beijing downfault basin and the Daxing uplift (Jiao and Qiu, 2006). Since the Quaternary, its roughly EW-trending normal-fault pattern has been consistent with the regional tectonic stress field (Bai et al., 2018; Feng, 2014; Zhang et al., 2014). The central segment of the Nanyuan-Tongxian fault was active in the Holocene, potentially offsetting Quaternary strata to the surface and/or providing stress transfer and accumulation for the formation of ground fissures (Fang et al., 2016; Zhao et al., 2018). Niu et al. (2023) integrated geophysical profiling, multi-borehole cross-sections, and geochronology to determine a fault dip of 56°–75° for the Nanyuan-Tongxian fault. According to the relationship between magnitude and fault dimensions, and statistical data on the source-rupture lengths of earthquakes in North China since 1965 (Cheng et al., 2019; Long et al., 2006; Ran, 2011), the rupture length for M 6 events is generally less than 70 km. Because the SW segment of the Nanyuan-Tongxian fault (between the Nankou-Sunhe fault and Zhuozhou) is a bedrock fault that shows no Quaternary activity (Huang et al., 2012; Jiang et al., 2001), we exclude this inactive SW segment. Referring to previous surface traces of central and NE segments and earlier fault-attribute studies (Jiao and Qiu, 2006; Niu et al., 2023; Xu et al., 2016), the fault model is configured as a 70°-dipping, 20-km-deep normal fault with a curved surface length of approximately 68 km (Figure 2).

Figure 2

3D models labeled (a) and (b) represent geological structures. Panel (a) shows terrain with fault traces marked in yellow and red. Panel (b) depicts a cross-section indicating a dip of seventy degrees and depth of twenty kilometers. Coordinates are 40.7 degrees north, 38.6 degrees north, 115.2 degrees east, and 117.3 degrees east.

Figure 2. Numerical simulation region and fault-plane geometry. (a) Simulation domain incorporating real topography; (b) spatial position of the fault geometry with a 70° dip, a curved surface length of ∼68 km, and a width of ∼20 km.

2.2 Dynamic parameters

Investigating the constitutive relationships that pervasively govern rupture processes is the physical basis for understanding complex rupture behavior; these relationships define how stress, slip, slip rate, and other relevant physical parameters are linked on the fault plane. Friction laws are crucial because they control how earthquakes are triggered, how rupture propagates after initiation, and how it ultimately arrests. To date, most friction laws have been derived from laboratory rock-fracture experiments (Dieterich, 1972; Dieterich, 1978; Dieterich, 1979; Ruina, 1983). Numerous low-velocity rock-friction measurements (Dieterich, 1978; Weeks and Tullis, 1985) can be explained by rate- and state-dependent friction laws (Dieterich, 1979; Dieterich, 1981; Ruina, 1983), whereas high-velocity slip experiments (Okubo and Dieterich, 1986) reveal that material strength weakens with slip (Dieterich, 1981; Ohnaka et al., 1987; Okubo and Dieterich, 1984). Many researchers have proposed various slip-weakening friction laws on the basis of theoretical considerations or laboratory tests (Campillo et al., 2001; Ida, 1972; Matsu’ura et al., 1992; Ohnaka and Yamashita, 1989). Ida (1972) assumed several types of stress evolution and first introduced the slip-weakening friction law into dynamic rupture models; as an approximation of rate- and state-dependent friction, the simple slip-weakening law can explain the mechanism of dynamic strike-slip rupture propagation (Okubo and Dieterich, 1984; Okubo and Dieterich, 1986). The results from studies of actual earthquakes (Ide and Takeo, 1997; Xie et al., 2024; Zhang et al., 2003) indicate that stress decreases with increasing slip, broadly satisfying the slip-weakening friction law (Ida, 1972); hence, the slip-weakening friction law is adopted in this study (Figure 3).

Figure 3

Graph depicting stress versus slip with a triangular area labeled $G_c$ shaded in blue. The y-axis indicates stress from $\sigma_y$ to $\sigma_s$, while the x-axis represents slip ending at $D_c$. The slope decreases from $\sigma_y$ to $\sigma_s$, where a constant stress level is maintained.

Figure 3. Slip-weakening friction law (Ida, 1972) $σ_{y}$ : yield stress; $σ_{0}$ : initial stress; $σ_{s}$ : residual stress; $D_{c}$ : slip-weakening distance; $G_{c}$ : fracture energy.

In the slip-weakening friction law, the slip-weakening distance Dc is also a critical parameter. First proposed theoretically by Ida (1972) as a key parameter of the slip-weakening friction law (Palmer and Rice, 1973), Dc is defined as the slip required for the local stress to drop to the residual static value. It reflects the physical state of the fault and primarily governs the fracture energy (Abercrombie and Rice, 2005; Cocco et al., 2006; Tinti et al., 2004). Research by Peyrat et al. (2001) demonstrated that similar dynamic rupture processes and numerical results that are consistent with observations (Guatteri and Spudich, 2000) can be obtained via different yield stresses and slip-weakening distances. While studying the dynamic rupture process of the 2018 Hokkaido–Japan earthquake, Xie et al. (2024) reported that the slip-weakening distance controls both the rupture propagation speed and the slip rate on the fault plane. Furthermore, the rupture duration for M6.0 earthquakes typically ranges between 10 and 15 s (Xie et al., 2024; 2025; Xu et al., 2023). To constrain the slip rate and rupture speed on the fault plane, a slip-weakening distance is constructed, as shown in Figure 4.

Figure 4

Three-dimensional graph depicting depth in kilometers versus latitude and longitude. Color gradient from blue to red indicates varying Dc values in meters, with a red star marking a significant point.

Figure 4. Spatial distribution of slip-weakening distances on the fault plane The red star marks the epicenter.

In dynamic simulation processes, stress on the fault plane is difficult to measure directly and is typically determined via methods such as borehole measurements or hydraulic fracturing (Chen et al., 2012; Hast, 1969; Li and Wang, 1979; Zoback et al., 1980). These methods can establish the background stress state on the fault plane and are therefore widely applied in studies of paleoearthquakes (Sunilkumar et al., 2024; Wang et al., 2022). Alternatively, another approach to determine stress on the fault plane involves calculations based on slip amounts on the fault surface (Ide and Takeo, 1997; Xie et al., 2024; Zhang et al., 2003). This method more readily yields dynamic rupture processes that are consistent with observations but requires first performing a kinematic rupture process inversion. The shear stress time history is subsequently determined by solving the equations of elastodynamics (Zhang et al., 2003) or using frequency-domain methods (Andrews, 1980; Xie et al., 2024), from which the dynamic parameters on the fault plane are derived. Given the lack of seismic observation data for paleoearthquakes, a suitable approach is to set the fault plane stress on the basis of existing hydraulic fracturing or borehole data. Wang et al. (2012) investigated the stress state in the upper crust via 1,780 hydrofracturing and stress-measurement datasets from mainland China:

S_{H} = 22.66 z + 6.59

S_{h} = 16.35 z + 3.59

where $S_{H}$ and $S_{h}$ are the maximum and minimum horizontal principal stresses (MPa) in mainland China, respectively, and z is depth (km); Feng (2014) using hydraulic-fracture measurements, obtained the present-day crustal stress variation with depth in the Beijing area. The results show that within the 0–40 km crustal depth range, the three principal compressive stresses are related as follows: in the 0–15 km depth range $S_{H} > S_{v} \geq S_{h}$ , a strike-slip stress state; in the 15–35 km depth range $S_{H} > S_{v} > S_{h}$ , a strike-slip stress state; and in the 35–40 km depth range $S_{v} > S_{H} > S_{h}$ , a normal-fault stress state. Because the simulation results in this paper focus mainly on 0–40 km, especially within the 0–20 km range of the fault plane, and referring to the above results and Wan (2010), R = 0.3 at 40°E, 117°N in the study area, the vertical principal compressive stress $S_{v}$ is set (Figure 5). Since a linear relationship exists between the stress on the fault plane and the resulting slip amount (Ide and Takeo, 1997; Xie et al., 2024; Zhang et al., 2003), the dynamic and static friction coefficients were adjusted multiple times and finally set to 0.4 and 0.15, respectively, to ensure that the moment magnitude (M_w) obtained from the dynamic simulation approximates that of the 1,057 earthquake.

Figure 5

Graph depicting principal stress versus depth. The y-axis represents depth in kilometers, ranging from zero to thirty. The x-axis shows stress in megapascals, from zero to one hundred. Three lines indicate different stresses: dashed for $ S_H $, solid for $ S_h $, and dotted for $ S_v $.

Figure 5. Variation in the principal stresses with depth. The maximum horizontal principal stress S_H, the minimum horizontal principal stress S_h, and the vertical principal stress S_v increase linearly with depth and then remain constant; modified from Wang et al. (2012).

For the numerical simulations, the layered velocity model (Figure 6) at the central longitude and latitude was extracted from the USTClitho 1.0 (Unified Seismic Tomography models for continental China lithosphere) lithospheric velocity model for mainland China (Xin et al., 2018) and used as the input velocity structure.

Figure 6

Graph of a velocity model showing velocity versus depth. The blue line represents Vp (primary wave velocity), and the red line represents Vs (secondary wave velocity). Depth is in kilometers, and velocity is in kilometers per second. The chart indicates depth levels from 0 to -60 kilometers and velocity from 2 to 8 kilometers per second.

Figure 6. One-dimensional layered velocity model. Derived from the USTClitho 1.0 velocity model at the epicenter’s central longitude and latitude.

2.3 Numerical method

The simulation of seismic source rupture processes can employ various methods. For earthquakes recorded at densely distributed stations, kinematic methods can be used to invert the slip distribution on the fault plane from observational data (Guo et al., 2024; Liu et al., 2013). Such kinematic source methods cannot reveal the physical mechanisms of rupture propagation. Consequently, dynamic numerical methods based on rock physics experiments are increasingly being used to simulate seismic source rupture processes (Wang et al., 2022; Weng and Ampuero, 2019; Xie et al., 2024). Unlike kinematic source methods, this approach requires setting stress parameters on the fault plane and constructing a physically consistent dynamic rupture process of the seismic source on the basis of the equations of elastodynamics.

Because the dynamic rupture process is influenced by many factors, analytical solutions are generally unattainable; numerical methods are therefore routinely employed. A variety of numerical techniques capable of simulating seismic wave propagation have been successfully applied to dynamic rupture modeling, such as the finite-difference (FD) method introduced by Madariaga (1976) and Andrews (1976) for earthquake rupture studies (Day, 1982a; 1982b; Miyatake, 1992); the boundary‒integral equation method, first applied to dynamic rupture by (Das and Aki, 1977; Andrews, 1985; Das and Aki, 1977); and the spectral‒element method (a high-order finite‒element‒spectral method). The spectral-element method combines the accuracy of spectral techniques with the flexibility of finite-element methods (Patera, 1984; Seriani and Priolo, 1991) and is therefore a high-order numerical scheme (Komatitsch et al., 2005). Initially, proposed by Patera (1984) for solving the Navier–Stokes equations, this method incorporates domain decomposition and Gauss–Lobatto–Legendre (GLL) nodes on spectral elements, enhancing numerical precision while allowing unstructured meshes outside the rupture zone to reduce computational cost (Festa and Vilotte, 2005; Galvez et al., 2020).

In this study, the spectral-element code Specfem3D (Komatitsch and Tromp, 1999; Komatitsch and Tromp, 2002a; Komatitsch and Tromp, 2002b) is employed to simulate the dynamic rupture process. Following previous work (Huang et al., 2012; Jiang et al., 2001; Jiao and Qiu, 2006; Niu et al., 2023), the upper edge of the fault plane extended to the surface. The computational domain measures 2.1 ° × 2.1 ° × 40 km; the fault plane is discretized into 38,077 (377 columns × 101 rows) subfaults, each 200 m × 200 m. Absorbing boundary conditions (Clayton and Engquist, 1977) are applied on all outer boundaries except the free surface to suppress reflections that would otherwise contaminate the rupture process. Each element uses 5 GLL nodes; the time step is 0.0001 s and the total simulation duration is 20.0 s.

3 Results

On the basis of the above dynamic parameters, the dynamic rupture process of the 1057M6 earthquake in the Beijing area can be obtained. The dynamic simulation yields a maximum slip of approximately 0.75 m on the fault plane (Figure 7), an average slip of approximately 0.49 m, a rupture duration of ∼10 s, and a moment magnitude of M_w 6.5 (Figure 8), which aligns closely with M6¾ reported in earlier studies (Li et al., 2021). The maximum slip on the fault plane is located northeast of the epicenter. Owing to fault geometry, slip gradually increases to ∼0.75 m in this area. At the epicenter location, slip is not distributed uniformly in regions with significant geometric deformation (Figure 7). Simultaneously, influenced by fault geometry, the rupture propagation process on the fault plane undergoes changes (Figure 9): rupture propagates bilaterally from the epicenter; at ∼2 s, it bends on the SW side of the fault plane; at 3 s, it propagates to the upper and lower boundaries of the fault plane and reaches the maximum slip rate (∼0.7 m/s); at 4 s, rupture crosses the bend on the SW side of the fault plane, with its slip rate exceeding that of the smoother propagation process on the NE side; especially at 5.5 s, an arcuate rupture front propagating at the maximum slip rate is distinctly observable on the SW side. Moreover, affected by the greater slip-weakening distance at the upper boundary (Figure 4), the slip rate at the upper boundary during rupture is significantly lower than that at the lower section of the fault plane, and the rupture timing is delayed; at 6.5 s, the rupture propagates to the SW boundary of the fault plane, whereas the NE side achieves the maximum slip rate; at ∼8 s, it also propagates to the NE boundary of the fault plane; at ∼9 s, the rupture fully traverses the entire fault plane, with energy essentially released, with only minor slip velocities exhibited in partial areas owing to surface reflections.

Figure 7

Three-dimensional graph showing slip distribution in meters along a fault plane, with depth in kilometers. Color gradient ranges from blue at lower slip values to red at higher slip values, with a red star indicating a significant point.

Figure 7. Spatial distribution of fault slip obtained from dynamic simulation. The maximum slip of 0.75 m is located on the NE side of the epicenter, with an average slip of approximately 0.49 m, and it has a nonuniform distribution due to the geometric curvature of the fault; The red star marks the epicenter.

Figure 8

Graph illustrating the source time function with moment rate measured in units of ten to the power of seventeen newton-meters per second on the vertical axis, and time in seconds on the horizontal axis. The graph peaks between three and four at approximately six seconds. A box notes the seismic moment as 8.17 times ten to the power of eighteen newton-meters and moment magnitude as 6.5.

Figure 8. Moment magnitude function and rupture time. The curve of the seismic moment release with time shows a moment magnitude of M_w6.5, with a rupture duration of approximately 10 s. The main energy release is completed in the first 8 s, and only residual slip remains after 9 s.

Figure 9

Twelve panels show seismic slip rate evolution over time from zero to fifteen seconds, in increments of 0.5 to 3 seconds. Each panel displays a 3D graph with time labels, depth, and geographical coordinates, with a color gradient indicating slip rates from zero to 0.6 meters per second. A red star marks the origin point in each panel.

Figure 9. Dynamic rupture propagation process. The spatiotemporal evolution of the fault slip rate; the rupture propagates bilaterally: owing to geometric bending, the SW side has a higher slip velocity peak and propagates faster than the NE side does; the upper boundary rupture is inhibited by high D_c; The red star marks the epicenter.

As a historically earthquake-prone region, Beijing has been extensively studied with respect to the impacts of historical or scenario earthquakes on strong ground motions in Beijing and surrounding areas. The peak ground acceleration in the North China Plain indicates greater seismic hazards in Tangshan, Taiyuan, and Beijing than in other regions (Yang et al., 2008). For example, the 1,679 Sanhe-Pinggu M8 earthquake occurred on the nearby Sanhe-Pinggu Fault. Various strong ground motion simulations (Ba et al., 2022; Fu et al., 2015; Gao et al., 2002; Liu et al., 2007; Pan et al., 2006) have shown that strong ground motions from M8 earthquakes on the Sanhe-Pinggu Fault are influenced by basin amplification effects and surface sedimentary layers in the Beijing area; simulations of the 1998 Zhangbei earthquake revealed greater amplification effects in Quaternary sedimentary layers than in Tertiary caprocks (Zhou and Chen, 2008). Additionally, strong ground motions are simultaneously affected by rupture processes (Xie et al., 2025) and topography (Wang et al., 2018). Owing to significant geometric bending and kinks in this earthquake fault, its generated strong ground motions (Figure 10) have a maximum intensity of VIII, with the VIII-degree zone primarily distributed along the fault trace and exhibiting a strong hanging wall effect (Xie et al., 2025; Zhang et al., 2003). However, influenced by the rupture process (Figure 9), higher slip velocities occur in the NE-trending rupture, combined with SE-ward fault bending, creating rupture front effects. Consequently, areas distant from the epicenter (Shunyi, Tongzhou, Sanhe) along the rupture front experience substantial strong ground motions; Wang et al. (2022) simulated another fault in this area (Sanhe‒Pinggu Fault Zone), which did not generate VIII-degree strong ground motions in Shunyi and Tongzhou.

Figure 10

Seismic intensity map of a region near Beijing, China, with colored gradients representing intensity levels from low (blue) to high (red). Major locations like ChangPing, HaiDian, and TongZhou are marked, with black-and-white circles indicating seismic activity. Red lines depict fault lines.

Figure 10. Intensity distribution of strong ground motion. The intensity distribution map converted from the peak ground acceleration (PGA) of the simulated seismic motion. The Ⅷ intensity zone is concentrated along the fault trace, with a significant hanging wall effect. Owing to the influence of rupture directivity and fault bending, abnormally high intensities appear in far-field areas such as Shunyi and Tongzhou; The focal mechanism solutions are taken from the 2009–2021 dataset for mainland China and adjacent regions (Guo et al., 2022).

4 Discussion

4.1 The impact of the epicenter location on strong ground motion

This earthquake (the 1057 M6 earthquake) occurred nearly a thousand years ago, and the accumulated stress was relatively high. Determining the location of the next earthquake is difficult. To improve the seismic fortification requirements, two earthquakes are set on this fault to simulate the dynamic rupture process. Figure 11 shows the strong ground motions simulated at different epicentral locations. The results indicate that there is a significant directivity effect in the strong ground motions generated whether the initial rupture point is on the left side of the fault plane (Figure 11a) or the right side of the fault plane (Figure 11b) (Wang et al., 2018; Yang and Zhao, 2010). When the epicenter is on the left side of the fault plane (Figure 11a), the slip rate reaches the maximum value when the rupture propagates to the right side of the fault plane, and strong ground motions are generated in the Tongzhou area following the bending direction of the fault plane. When the epicenter is on the right side of the fault plane, the slip rate reaches the maximum value when the rupture propagates to the left side, and larger intensity areas are generated in the Gaobeidian and Zhuozhou areas. Combined with Figure 10, owing to the influence of the fault geometry and the rupture propagation process, strong ground motions are generated in the Tongzhou area if the epicenter is on the left side of the center of the entire fault.

Figure 11

Comparison of two earthquake intensity maps. Both panels (a) and (b) display intensity levels from zero (blue) to eight (red) across a region with marked locations. Red lines indicate fault lines, while black and white circles represent seismic stations. A red star and circle denote epicenter and maximum intensity areas.

Figure 11. Comparison of seismic motion scenarios at different epicentral locations. The red star indicates the set epicenter location, and the red circle represents the original epicenter location. (a) When the epicenter is on the left side of the southern segment of the fault, Tongzhou experiences high-intensity shaking. (b) When the epicenter is on the right side, the intensity increases in Gaobeidian and Zhuozhou; The focal mechanism solutions are taken from the 2009–2021 dataset for mainland China and adjacent regions (Guo et al., 2022).

4.2 Influence of the fault geometry

Dynamic rupture is also governed by fault geometry (Dong and Zhang, 2019; Ma et al., 2021; Ma et al., 1996; Zhang et al., 2016); for large earthquakes, irregular fault geometry can even lead to supershear rupture (Bhat et al., 2007). To investigate how fault geometry influences the strong ground motion of the 1057 M 6 earthquake, a planar fault with the same dip (70°) (yellow line in Figure 12) was constructed across the epicenter for dynamic rupture simulation, and the resulting ground motions were compared. To minimize the influence of other factors, the previously determined dynamic parameters were retained. Figure 12 shows the ground motions generated by the planar fault: the maximum intensity is IX, and the intensities on both sides of the fault are symmetrically distributed. Because the epicenter is closer to the southwestern end of the fault, the rupture propagates more slowly than on the northeastern side, so no pronounced directivity effect is produced in the southwest. The planar rupture (Figure 12) produces intensities of VII or higher in XiangHe, whereas the nonplanar rupture does not produce comparable intensities (Figure 10), which is consistent with the findings of Zhang et al. (2019) when studying rupture intensities for the Wenchuan earthquake using different fault geometries: planar and nonplanar ruptures yield markedly different intensities in different regions.

Figure 12

Seismic intensity map showing regions with varying colors representing intensity levels from zero to eight, indicated by a vertical color scale. Red lines show fault lines. Marked cities include HuaiLai, ChangPing, and MiYun. A red circle indicates a significant seismic event near HaiDian. Latitude and longitude coordinates are marked along the edges.

Figure 12. Trace of the planar fault and the spatial distribution of the resulting strong ground motion. The yellow line indicates the constructed planar fault; the blue line represents the original curved fault. The red circle denotes the epicenter—i.e., the initial rupture point. The color scale represents the intensity converted from the peak ground velocity; The focal mechanism solutions are taken from the 2009–2021 dataset for mainland China and adjacent regions (Guo et al., 2022).

4.3 Influence of the velocity model

In addition to the location of the initial rupture point and the fault geometry, the medium parameters also control the strong ground motion generated during dynamic rupture. To verify the ground motions of the 1,057 earthquake under different velocity models, we constructed a 3-D velocity model for the study area on the bases of USTClitho 1.0 (Deng et al., 2003; Li et al., 2013; Xin et al., 2018). USTClitho1.0 is a 3-D lithospheric velocity structure for mainland China obtained from seismic travel-time tomography Xin et al. (2018); its horizontal resolution is ∼50–100 km and is therefore suitable for building a 3-D data model for our study region. Using this 3-D velocity model and the dynamic parameters determined with the curved fault, we simulate the resulting strong ground motion (Figure 13a).

Figure 13

Map and cross-section charts displaying geological data from a specified region. Panel (a) is a topographic map with faults and labeled locations like HuaiLai, MiYun, and ZhuoZhou. Panels (b) and (c) show longitude-based depth and velocity data, while panels (d) and (e) present latitude-based depth and velocity data. Color scales indicate intensity and velocity variations.

Figure 13. 3-D velocity model and velocity cross-sections near the epicenter (a) The blue line marks the surface trace of the curved fault; the yellow lines indicate the locations of the velocity cross-sections (b-e) Vp and Vs. velocity cross-sections along the yellow lines shown in (a).

The 3-D velocity model simulation revels that the hanging-wall side (southeast side) of the fault exhibits significantly greater ground motions than those produced by the 1-D layered model (Figure 10), particularly within the intensity-VII zone. Gao et al. (2005) used finite-element modeling to examine the relationship between the ground motions of the 1999 Chi-Chi earthquake and the stiffness of the fault-zone medium, and reported that the peak ground acceleration is approximately proportional to the fault-zone stiffness, whereas the peak ground velocity scales roughly with the square root of the stiffness. To investigate velocity contrasts across the fault, we extracted depth profiles along lines AB and CD (yellow lines in Figure 13) passing through the epicenter; the resulting Vp and Vs. cross-sections (Figure 13b–e) revealed pronounced differences across the Nanyuan-Tongxian fault. At depths down to 20 km and in the shallow crust, the hanging-wall (southeast) presents consistently lower Vp and Vs. than the hanging wall (northwest), a disparity caused mainly by different deformation styles across the fault (Le Pichon et al., 2005; Thurber et al., 2006).

Our study area lies within the Quaternary sedimentary basin. Ding et al. (2004) demonstrated through synthetic seismograms that Quaternary sediments strongly amplify ground motions in the Beijing region. Likewise, Zhou and Chen (2008), using the discrete wavenumber method to model two-dimensional ground motions in Beijing induced by the 1998 Zhangbei earthquake, reached conclusions similar to those of Ding et al. (2004). Moreover, larger contrasts in material properties lead to more pronounced spatial asymmetry in ground motion (Yuan and Zhu, 2016; Zhu and Jie, 2016). Therefore, the larger ground motions on the southeast side of the fault in our 3-D velocity simulations most likely resulted from the combined effects of fault-zone heterogeneity and amplification by Quaternary sedimentary cover.

5 Conclusion

This study reconstructed the rupture process of the 1057M6 earthquake through dynamic numerical simulation, revealing the unique faulting mechanism and disaster characteristics of the fault southwest of Beijing. The simulation results based on the geometric structure of the Nan Yuan–Tongxian fault and the parameters of the regional stress field reveal that the earthquake rupture lasted approximately 10 s, with a maximum slip of 0.75 m and a moment magnitude of Mw6.5, which is consistent with the historically recorded magnitude of M6³/⁴. Rupture propagation was significantly controlled by the fault geometry. The SW side experienced a sharp increase in slip velocity due to bending (peak value of 0.7 m/s), whereas the NE side experienced more stable propagation. The difference in the attenuation distance of slip at the fault boundaries further delayed the rupture process at the upper boundary.

The results of strong ground motion show that the maximum intensity of Ⅷ is distributed along the fault and that there is a significant hanging wall effect. Unlike the seismic effects of faults northeast of Beijing (such as the Sanhe–Pinggu fault), this simulation revealed abnormally high intensities in areas far from the epicenter, such as Shunyi and Tongzhou. This is the result of the combined action of the NE-dipping bend of the fault and the directivity of the rupture, forming a “rupture front effect”. A comparison of the simulation results at different epicentral locations reveals that if the fault point is located to the left of the center of the fault, the Tongzhou area will experience more intense shaking, whereas the right epicenter will experience Gaobeidian and Zhuozhou high-risk areas.

The planar-fault simulation shows that, when the rupture nucleates away from the fault’s geometric center, a pronounced directivity effect develops toward the northeastern side of the fault, producing intensities of VII or greater in XiangHe; in contrast, the curved-fault rupture does not generate a comparable intensity pattern, so the intensities produced by planar and nonplanar ruptures differ markedly in different regions. The 3-D velocity-model results further demonstrate that strong ground motion depends not only on the fault’s dynamic parameters and geometry but also on the medium encountered during rupture propagation: the greater the contrast in material properties is, the stronger the spatial asymmetry of ground motion. The greater ground motion on the southeastern side of the fault relative to the northwestern side is most likely due to the combined effects of heterogeneous media on either side of the fault and amplification by Quaternary sediments.

This earthquake occurred nearly a thousand years ago. The accumulated stress released by the southwest fault is low, and there is still a risk of recurrence in the future. The study results show that fault geometry is the key factor controlling rupture propagation and the spatial distribution of ground motion. The amplification effect of the sedimentary layer of the Beijing Basin may further exacerbate the disaster. This result provides a dynamic basis for the seismic fortification of Beijing and its surrounding areas (especially potential high-risk areas such as Tongzhou and Shunyi) and highlights the need to focus on the stress accumulation state of the Nanyuan–Tongxian fault.

Data availability statement

The focal mechanism solution data are from Guo et al. (2022). The active fault data are from Xu et al. (2016). The numerical simulation program Specfem3D is from Komatitsch et al. (2024). Figures were created via Matplotlib and GMT Wessel et al. (2013). The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

ZX: Conceptualization, Data curation, Formal Analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review and editing.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. We are grateful to Beijing Gangzhen Technology Co., Ltd. for providing the research topic “The Impact of Strong Ground Motion on Buildings” (Fund No. YF-202520). The funder was not involved in the study design, collection, analysis, interpretation of data, the writing of this article, or the decision to submit it for publication.

Acknowledgments

We are deeply grateful to the editors and the two reviewers for their thorough evaluation of our manuscript and for providing numerous detailed and constructive suggestions, which have greatly improved its clarity, accuracy, and overall quality. We also thank the Computational Infrastructure for Geodynamics (CIG, http://geodynamics.org).

Conflict of interest

Author ZX was employed by Beijing Gangzhen Technology Co., Ltd.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Abercrombie, R. E., and Rice, J. R. (2005). Can observations of earthquake scaling constrain slip weakening? Geophys. J. Int. 162 (2), 406–424. doi:10.1111/j.1365-246X.2005.02579.x