The physics-based deterministic scenarios for earthquake hazards and losses of the zhujiangkou fault in southern china

ABSTRACT Dynamic rupture directivity significantly impacts earthquake hazard and loss assessments; yet, it is often overlooked, leading to uncertainties. We developed high-resolution,

broadband, deterministic scenarios for earthquake hazards and losses of the Zhujiangkou Fault (ZF) in the Guangdong-Hong Kong-Macao Greater Bay Area of southern China. The calculations cover

the entire physical process from fault dynamic rupture through strong ground motion propagation to quantitative loss estimations. We fully incorporate the real 3D physical fields, including

fault structure, regional stress, friction conditions, inhomogeneous velocity, and topography. Results revealed higher earthquake hazards and losses triggered by the northwestern ZF

segment, affecting Hong Kong, Shenzhen, Dongguan, Macao, and Guangzhou. The rupture directivity effect, driven by seismic wave interferences in forward propagation, leads to stark

differences in hazards and losses across urban agglomerations, even under similar fault conditions. These findings help enhance the precision of earthquake risk assessments and provide

practical approaches for disaster prevention. SIMILAR CONTENT BEING VIEWED BY OTHERS SEISMIC HAZARD PREDICTION OF THE HUNHE FAULT IN THE SHEN-FU NEW DISTRICT Article Open access 25 June 2024

IRREGULAR RUPTURE PROPAGATION AND GEOMETRIC FAULT COMPLEXITIES DURING THE 2010 MW 7.2 EL MAYOR-CUCAPAH EARTHQUAKE Article Open access 17 March 2022 CONSECUTIVE RUPTURES ON A COMPLEX

CONJUGATE FAULT SYSTEM DURING THE 2018 GULF OF ALASKA EARTHQUAKE Article Open access 16 March 2021 INTRODUCTION The Guangdong-Hong Kong-Macao Greater Bay Area (GGBA), as one of the four

major bay areas globally, stands as one of the most modernized and economically vibrant regions in China1. In 2020, the economic output of the GGBA exceeded 1.6 trillion US dollars. This

output is comparable to the GDPs of Canada and South Korea, which ranked ninth and tenth globally that year, respectively2. Moreover, the remarkable economic strength of the GGBA has induced

a significant siphon effect on population dynamics. In the past decade, the permanent population surged by more than 21 million, reaching approximately 78.01 million in 20203. While the

GGBA experiences fewer earthquakes than China’s Sichuan-Yunnan region, it is still considered a seismically active area4,5. Should extreme earthquakes occur, they could profoundly impact the

GGBA. Other major bay areas globally have suffered extensive damage from earthquakes, offering painful lessons and crucial insights. In 1923, the Tokyo Bay Area endured the devastating Mw

7.9 Great Kanto Earthquake, which ignited massive fires, caused numerous building collapses, and triggered floods. This disaster led to over 100,000 fatalities6. The region faced further

calamity during the Mw 9.0 Tohoku Earthquake in 2011, significantly damaging its infrastructure, buildings, and economy and resulting in leakage at the Fukushima nuclear plant7. Similarly,

the Mw 7.9 San Francisco earthquake in 1906 struck the San Francisco Bay Area, which caused severe damage to the city, including building collapses, widespread fires, and the paralysis of

water and transportation systems, resulting in over 3000 fatalities and leaving a large number of residents homeless8. Moreover, the Mw 6.9 Loma Prieta Earthquake in 1989, though of lower

magnitude, still caused significant destruction of critical infrastructure such as bridges and highways and disrupted the social functions of urban agglomerations, resulting in 63 deaths,

thousands of injuries, and billions in property losses9. These seismic events have driven major bay areas globally to continually improve their earthquake preparedness and mitigation

strategies, enhancing their resilience and reducing the impacts of future earthquake disasters. The Zhujiangkou Fault (ZF) is located in the core region of the GGBA (Fig. 1), surrounded by

the most developed large urban agglomerations in China, including Guangzhou, Shenzhen, Hong Kong, and Macau. The fault extends over 90 km and is primarily distributed in the coastal area10.

This fault has exhibited a vertical annual slip rate ranging from 0.14 to 0.47 mm/a and presented a main feature of sinistral slip motion since the late Quaternary, indicating a certain

seismic potential4,5. Furthermore, the seismicity of the ZF has been further substantiated by combining experimental stress analysis techniques with rock mechanics, assessing the potential

maximum magnitude of approximately Mw 7.011. However, due to the historical rarity of earthquakes along the ZF and the scarcity of historical seismic records, there is a lack of scientific

seismic design guidelines for the surrounding regions. Moreover, infrastructures, socio-economic systems, and ecological environments are highly vulnerable to extreme seismic events. A

destructive earthquake could cause significant and unpredictable impacts on public life, socio-economics, and even national security. Therefore, it is essential to conduct comprehensive

research on the ZF and the associated potential earthquake hazards and risks in the GGBA. In recent decades, earthquake hazard analysis methods have become critical technologies for reducing

earthquake risk. Particularly, empirically-based Ground Motion Prediction Equations (GMPEs) are widely used to estimate ground motion or seismic intensity, often serving as reference

standards for seismic-resistant design in engineering projects and land planning12,13. However, the data supporting GMPEs frequently lack records from large earthquakes near faults. This

deficiency hampers the ability of GMPEs to accurately interpret the detailed impacts of ground motion near faults, which are often the primary areas of seismic damage14. Furthermore, GMPEs

involve an ‘ergodic assumption’, which leads to a lack of precision when predicting ground motion in specific regions based on multi-source data15,16. Failure to effectively consider the

‘spatial correlation’ between adjacent sites often results in underestimating earthquake hazard17,18,19,20. Additionally, existing GMPEs do not adequately address the fault rupture

directivity. This limitation results in the inability to fully represent the variability of specific seismic scenarios triggered by different nucleation zones along the same fault.

Physics-based deterministic simulation methods have been validated through recordings at seismic stations and empirical GMPEs21,22. Compared to the traditionally employed empirical GMPEs,

these models have shown great potential for accurately assessing earthquake hazards by comprehensively considering the complexities of the source and the multi-physical parameter fields of

the subsurface23,24,25,26. Nonetheless, these methods also exhibit certain limitations. Some can calculate only low-frequency ranges, while others employ stochastic methods for high

frequencies. Others overlook the impact of complex 3D physical parameters during simulation, such as fault geometry, friction conditions, stress field distribution, rupture directivity,

heterogeneous velocity structures, and undulating topography. Alternatively, they utilize kinematic models to characterize the seismic source. Therefore, the results may not be sufficient to

fully reflect the potential future earthquake hazards in the target area. In this study, we have developed physics-based, high-resolution, broadband, deterministic scenarios for earthquake

hazards and losses of the Zhujiangkou Fault in the Guangdong-Hong Kong-Macao Greater Bay Area of southern China. This work encompasses most of the crucial elements, from the dynamics of

fault rupture, seismic wavefield propagation, and surface ground motion damage to urban agglomeration losses quantitative assessment. The computation considers the actual 3D physical

parameter field, including fault geometry structure, focal mechanisms, regional stress fields, fault friction conditions, inhomogeneous velocity mediums, and undulating topography.

Specifically, the complex 3D structure of the fault is meshed, and the appropriate earthquake nucleation area is selected through stress analysis. Fault dynamics rupture is calculated using

the curved-grid finite difference method to obtain source models capable of characterizing complex source processes of earthquake scenarios. Through the refined strong ground motion

simulation, the broadband ground motion distribution has been obtained and validated by the empirical GMPEs. The earthquake loss estimation model based on historical data is used to

quantitatively estimate the human and economic losses under different earthquake scenarios. The uncertainty of the results of this study is further discussed through the calculation and

analysis of different nucleation depths and topography. The results indicate that the potential earthquake hazards and losses in the northwestern segment of the ZF are significantly higher

than those in the southeastern segment. Areas along the Bay in Hong Kong, Shenzhen, Dongguan, and most parts of Macao and Guangzhou face higher earthquake risks. Additionally, the

experimental results highlight the significant impact of dynamics' rupture directivity on earthquake hazard and loss assessment. For different earthquake scenarios initiated on the same

fault within the study area, despite identical input physical conditions and similar fault slip and magnitude, the directivity effect of fault rupture can give rise to vastly different

hazards and losses among urban clusters. This study provides important insights for improving the precision of earthquake risk assessments. Moreover, it presents a practical approach for

future earthquake disaster prevention efforts, promoting sustainable urban agglomeration development. RESULTS FAULT DYNAMIC RUPTURE In this study, ZF is identified as a vertical fault with a

main strike of 155°, extending approximately 93 km in length and 10 km in width along the dip (Fig. 2a). To comprehensively assess the potential earthquake hazards along the ZF, six

earthquake nucleation zones are set in areas of high shear stress on the fault plane (see Fig. 2e). The curved-grid finite difference method (CG-FDM)27,28 is employed to numerically simulate

the dynamic rupture process of the ZF. Figure 3 illustrates the results of the dynamic rupture simulation on the ZF, depicting the fault surface slip distribution and isochrones of rupture

timing (1 s intervals). Figure 3a–f correspond to Scenarios 1–6, respectively. Scenarios 1–3 are characterized by magnitudes close to Mw 7.0, with Scenario 2 exhibiting the most significant

magnitude at Mw 6.92 and a maximum slip displacement of 3.23 m. In contrast, Scenarios 4–6 show significantly lower magnitudes and lateral slip displacements. Scenario 4 features the lowest

magnitude at Mw 5.68, and its rupture arrests near the nucleation patch and does not reach the surface. The distribution of subsurface stresses can effectively explain these modes of fault

rupture. The ratio of shear stress to normal stress on the fault plane (_μ__0_) appears to be an intuitive measure of the potential rupture capability. Around the source location of Scenario

4, _μ__0_ is notably low (Fig. 2f), severely impedes its rupture propagation outward. The magnitude of stress projected onto the fault plane serves as a mapping expression of the fault

geometric structure. In this study, the principal stress direction is N115°, accompanied by a 20° change in strike direction in the central region of the fault. These geometric structures

cause variations in the projected values of normal and shear stresses. Ultimately, this leads to the formation of two areas with low _μ__0_ in the middle of the fault plane, making it

difficult for earthquakes triggered to rupture along the entire fault completely. Therefore, a bend in the ZF substantially reduces the hazards of triggering extreme earthquake scenarios.

Overall, the earthquake hazards posed by ruptures initiated on the northern side of the ZF (Scenarios 1–3) are significantly higher than those on the southern side. This disparity leads to a

notably increased earthquake risk for Guangzhou compared to other cities within the GGBA. The dynamic rupture processes described herein do not precisely predict actual surface damage, as

they are also closely linked to the effects of three-dimensional subsurface medium properties and topography on seismic wave propagation. Consequently, these factors necessitate further

calculations and analysis through strong ground motion simulations. STRONG GROUND MOTION We utilize results from dynamic rupture calculations as the input seismic source for strong ground

motion simulations. Figure 4a–f show snapshots of seismic wavefields and seismic intensity distributions for Scenarios 1-6, where seismic intensity is derived from the peak ground velocities

(PGV) according to the Chinese Seismic Intensity Scale. Overall, the southeastern segment of the ZF exhibits smaller scales of rupture and seismic intensity. Earthquake scenarios

originating from the northwestern segment of the fault demonstrate significantly greater potential earthquake hazards than those on the southeastern segment. The seismic intensity maps

clearly illustrate that seismic energy is primarily concentrated along the direction of fault rupture propagation. Scenarios 1–3 display markedly different intensity distribution patterns

while sharing similar magnitudes and slip distributions (Fig. 3a–c), display markedly different intensity distribution patterns. Similarly, Scenarios 5–6 exhibit the same phenomenon.

Scenario 4 features a smaller rupture scale and displays surface energy distribution patterns akin to a double-couple point source. The rupture directivity effect is mainly caused by seismic

waves emitted early in the rupture, which can interfere constructively with waves emitted later, particularly in the forward direction. The southeast-directed seismic wavefields well

demonstrate this in Fig. 4a, b and e, and the northwest-directed wavefields in Fig. 4c, f. The radiated energy from these wavefields aligns closely with the rupture direction, illustrating

how rupture directivity strongly correlates with surface energy distribution. In Scenario 1, energy radiation is concentrated in the southeast part of the GGBA, while in Scenario 3, it is

focused in the northwest. Scenario 2, which has its nucleation zone in the middle of the northern segment of the fault, shows a more uniform energy distribution along the fault. Similarly,

in Scenarios 5 and 6, where the nucleation zones are located at either end of the southeast segment of the fault, the directions of surface energy diffusion are opposite. These observations

emphasize the major influence of fault dynamic rupture directivity on the distribution of earthquake hazards and highlight the necessity of incorporating rupture directivity effects in

earthquake risk analysis. EARTHQUAKE LOSS AND RISK ANALYSIS Based on the Chinese earthquake population and economic loss assessment models29,30, and the spatial distribution of seismic

intensity, this study determines the fatality rate and direct economic loss rate for each grid point. By multiplying these loss rates by the population exposure and fixed capital stock

exposure at each grid point, the number of fatalities and direct economic losses for that grid point can be calculated. Spatial integration across the entire study area yields the total

expected fatalities and economic losses for the entire earthquake scenario. By incorporating the expected loss value and the standard deviation from the loss assessment model into the

cumulative normal distribution function, it is possible to evaluate the relative probability of occurrence that the total loss for this earthquake scenario falls within different loss

intervals (for detailed process, refer to the Methods section). Figure 5a–f present the earthquake loss assessment results for earthquake Scenarios 1-6. Each panel displays an assessment of

earthquake fatalities on the left and direct economic losses on the right. Pie charts illustrate the quantitative losses for different cities and their percentage of the total loss. The

center number in each chart indicates the total loss for that Scenario, and the colors correspond to city coding. Bar charts display the relative probability of the total loss for each

Scenario falling within specific loss intervals. The height of the bars indicates the relative probability, with colors transitioning from green (on the left) to red (on the right) to encode

the severity of losses. Scenario 1 exhibits the highest potential for earthquake losses, with estimated fatalities exceeding 12,000, and a 57.5% probability that the total fatalities will

exceed 10,000. In this Scenario, Shenzhen accounts for the highest proportion of deaths at 47.53%, followed by Guangzhou. Together, these two cities account for over 87% of the total

fatalities. The main reason for this phenomenon is that the population of Shenzhen is largely concentrated on the inner side of the Bay Area (many areas have population densities exceeding

20,000 people/km²), which is primarily affected by the seismic high-intensity areas of Scenario 1 (Fig. 4a). Although Guangzhou has a far lower population density in the high-intensity

impact areas than Shenzhen, it has a broader affected area and a larger total exposed population. The total economic loss for Scenario 1 is close to RMB 600 billion yuan (1 USD = RMB 7.13

yuan (at 21 August 2024)), with an 82.4% probability that the economic losses will exceed RMB 100 billion yuan. The most probable range for economic losses lies at RMB 100–1000-billion-yuan,

accounting for 42.9%. Hong Kong suffers the greatest economic losses, accounting for nearly 40% of the total, contrasting sharply with its share of fatalities (about 6%). This discrepancy

is attributed to the extremely high exposure of social assets in highly developed Hong Kong, resulting in a higher percentage of economic loss than fatality in all scenarios. Scenario 4

exhibits the smallest scale of earthquake damage losses (Fig. 5d), with a 96% probability that fatalities will be less than 100% and a 94% probability that economic losses will be under RMB

100 billion yuan. This can be considered a highly probable scenario for minor losses compared to other scenarios. The probability that fatalities will fall within the 1–10 range is the

highest, at 47.1%, with an expected value of 9 fatalities. Similarly, economic losses are most likely to fall within the RMB 1–10 billion-yuan range, with a probability of 45.5% and an

expected value of RMB 7.2 billion yuan. Overall, within the same study area, the number of fatalities and economic losses across different earthquake scenarios notably correlate with the

magnitude of the earthquakes. For Scenarios 1-3, which have magnitudes ranging from Mw 6.90 to 6.92 (Fig. 3), the probability of fatalities exceeding 1000 people is greater than 92%, and the

probability of economic losses exceeding RMB 100 billion yuan is over 78%. However, due to the presence of rupture directivity effects, significant differences in losses still exist between

different earthquake scenarios under conditions of the same physical parameters and similar magnitudes. Scenario 2 exhibits the highest magnitude, yet its most likely fatality range is

103–104, comparatively lower than the 104–105 range for Scenarios 1 and 3. However, its anticipated economic loss falls within a moderate range compared to the other two scenarios. These

findings indicate that the correlation between earthquake losses and magnitude is not linear, even within the same study area. This highlights the importance of conducting a comprehensive

analysis of earthquake hazards and losses considering ground motion and the quantity and value of exposure (population and fixed assets). Furthermore, the variability in exposure across

different earthquake scenarios triggered by the same fault is due to the effects of rupture directivity. This difference is even more pronounced in Scenarios 5 and 6. Despite sharing

identical physical parameters and nucleation depths, as well as similar magnitudes, slip distributions, and proximate nucleation locations, the directivity effects of dynamic rupture result

in vastly different earthquake losses. The expected fatalities in Scenario 6 are more than three times that in Scenario 5, and the economic losses exceed twice those in Scenario 5. In

addition, the interval of fatalities and economic losses corresponding to the potential maximum probability of Scenario 6 is one order of magnitude larger than that of Scenario 5. The

spatial variability inherent in these scenarios in the same fault cannot be adequately expressed when calculated for ground motion using traditional empirical GMPEs, with most parameters

being identical. Therefore, this study can conduct more accurate and reasonable comprehensive assessments of earthquake risks for specific faults in target study areas by considering the

directivity effects of fault rupture. With its smaller overall population and economic volume, Macao experiences relatively lower loss values than other cities. However, its extremely high

population density and exposure to social assets make it highly vulnerable to seismic damage. Therefore, spatial risk analysis of earthquake scenarios along the ZF is necessitated. The

analysis in Fig. 6 presents the spatial distribution of earthquake losses and risks for the seven cities bordering the Bay in the GGBA. The six panels correspond to Scenarios 1-6, with the

left side of each panel displaying the distribution of risk for earthquake fatalities and the right side depicting the economic losses risk. Based on the loss severity per unit grid area,

the earthquake risks are categorized into five levels: very low, low, middle, high, and very high, each represented by color coding. Areas of ‘very high’ risk for earthquake fatalities are

primarily concentrated on the eastern to northern sides of the GGBA, particularly near the Bay in Hong Kong, Shenzhen, Dongguan, and most parts of Guangzhou. The Macao Peninsula also falls

within the ‘high’ to ‘very high’ risk levels in Scenarios 1-2. Areas of ‘very high’ economic loss risk often correspond to highly urbanized and modernized city areas suffering from

high-intensity impacts, specifically near the Bay in Shenzhen, Hong Kong, and Macao. If an earthquake occurs on the northwestern segment of the ZF, most areas of the seven cities near the

Bay would be at ‘high’ to ‘very high’ risk levels. Conversely, if an earthquake is nucleated in the southeastern segment, only a few areas would be affected by ‘high’ risks. Overall, future

efforts should focus on strengthening earthquake-resistant construction and emergency management levels in the areas from the eastern to northern shores of the Bay. Attention should also be

paid to most parts of Guangzhou from the central to southern areas, as regardless of a future earthquake initiated on any segments of the ZF, the proportion of areas in GGBA affected by

‘middle’ or higher risk levels is relatively significant. DISCUSSION In comparing Scenarios 1-3, it is observable that earthquakes initiated in different nucleation zones along the fault

exhibit markedly different patterns of ground motion distribution despite having similar magnitudes. This variation is primarily due to the directivity effects of fault rupture. However,

empirically-based GMPEs do not adequately represent these phenomena. To further discuss the advantages of physics-based simulation methods in capturing the directivity of rupture and

enhancing the rationality of the present computations, we compared the numerical simulation results of ground motion from Scenarios 1-3 with those predicted by two widely recognized

empirical ground motion prediction equations, as depicted in Fig. 7. Grey dots represent the PGV at the surface computational grid points. Blue and red lines represent the predicted PGV and

their standard deviations for the BSSA1431 and CB0832 models, respectively. These models employ RJB (the closest distance to the surface projection of the rupture plane) and RRUP (the

closest distance to the rupture plane itself) to denote distance. The distances are equal because the fault in this study is a strike-slip fault, and the ruptures occurred at the surface.

Figure 7 shows that the discrepancy between the simulations and the GMPEs increases with the distance from the fault. This trend may be attributed to the physics-based simulations

effectively capturing the directivity effects of fault dynamics during rupture. As the propagation distance increases, the impact of rupture directivity becomes more pronounced. Furthermore,

the disparity in Scenario 2 is notably less than in Scenarios 1 and 3. This discrepancy may stem from the fact that the epicenter in Scenario 2 is located precisely at the central area of

the rupture zone (Fig. 3b), which significantly mitigates the effects of rupture directivity. In Scenario 2, the rupture process occurs symmetrically from both sides of the nucleation zone.

The rupture distance at each end is shorter than the unilateral rupture distance observed in Scenarios 1 and 3 (Fig. 3a, c), which means that the impact of rupture directionality is less

pronounced. In Scenarios 1 and 3, the nucleation areas are positioned at the edges of the rupture zone, amplifying the rupture directivity effects and increasing the disparity between the

numerical simulation results and the GMPEs. Despite these differences, the numerical simulation results and the empirical GMPEs demonstrate good consistency, underscoring the advantages of

physics-based simulations in capturing rupture directivity and cross-validating the rationality of the numerical computations. Using Scenario 2 as a reference, we calculated the earthquake

scenario based on the flat topography model, as shown in Fig. 8a. Compared with the flat topography model, the peaks in the western region of Hong Kong show obvious energy concentration in

Scenario 2, considering the undulating terrain model. Since the terrain in this study area is relatively flat (Fig. 1), with the highest altitude of about 1 km, the overall difference in the

calculated results is very small. The maximum seismic intensity of both scenarios is X, and the calculated PGV means differ by only 3.52%. This study established the nucleation depth at 7.5

 km based on the hypocenter depths of small to moderate earthquakes in the surrounding region5. These earthquakes are primarily located within a 5 km to 10 km depth range. Using Scenario 2

as a reference, we calculated additional earthquake scenarios with varying nucleation depths of 5.5 km, 6.5 km, and 8.5 km. The results of the dynamic rupture process are illustrated in Fig.

S1. Our findings indicate that the nucleation depth has little impact on the dynamic rupture process in this study area. Consequently, the final slip distribution and the moment magnitude

of the fault remain nearly identical across the different nucleation depths. Figure 8b illustrates the seismic intensity distribution for the earthquake scenario with the shallowest

nucleation depth of 5.5 km (see Fig. S2 for other scenarios with different nucleation depths). The maximum seismic intensity for all earthquake scenarios is X. The mean PGV values calculated

for the two scenarios with nucleation depths of 5.5 km and 7.5 km differ by only 6.7%. It is important to note that the ZF has no recorded history of destructive earthquakes, so the results

of this study can be viewed as representing potential extreme scenarios of destructive earthquakes along the fault. Some research suggests that the ZF may not be continuous, which could

mean that its maximum magnitude might not reach Mw 7.05. Furthermore, our simulations do not account for the viscoelastic attenuation of seismic waves, which could lead to overestimating the

earthquake hazards and losses. Other factors that could lead to inaccuracies in assessing earthquake hazards and risks include neglecting the plastic deformation of material responses,

failing to incorporate near-surface low-velocity media, and idealizing model setups such as fault geometry, stress configurations, and friction criteria. The onsite geology conditions, such

as lithological contacts and geotechnical properties, greatly impact the ground motions and seismic response, especially the shear wave speeds near the ground surface (Vs30). However, we

don’t have detailed geology information, such as the detailed velocity model near the ground surface, the soil type or stratification, etc.; it is beyond the scope of this work and our

capability. We will discuss the effects of local detailed geological information by cooperating with the relevant department in the future. Additionally, this study provides only a

preliminary earthquake risk assessment based on physical simulations, which include ground motion, fatalities, and economic losses. It does not consider more detailed issues such as time

exposure, building type, social conditions, etc., which is also the direction of our future work. This study employs physical-based methods to construct high-resolution, broadband,

deterministic scenarios for earthquake hazards and losses of the ZF. It fully considers various factors such as 3D fault geometry structures, stress fields, friction criteria, heterogeneous

media, and irregular topography. A comprehensive analysis of earthquake losses and risks in the Guangdong-Hong Kong-Macao Greater Bay Area was conducted. The results indicate that the

potential earthquake hazards of the southwestern segment of the ZF are significantly higher than those of the southeastern segment, attributable to the constraints imposed by the fault

geometry and stress distribution on its dynamic rupture. Additionally, the study finds that the fault rupture directivity effect, compared to traditionally emphasized factors like site

effects and propagation paths, plays a crucial role in determining the distribution patterns of seismic energy on the surface, thereby critically influencing the earthquake hazard analysis.

This study comprehensively assesses potential future earthquake hazards, losses, and risks. The research thoroughly considers the entire physical seismic generation process, from dynamic

rupture to ground motion propagation. Compared to traditional methods, this approach enables more realistic and efficient computations. Theoretically, given sufficient computational hardware

resources, simulating earthquake wavefields and seismic motions at any bandwidth and spatial resolution is possible. This study offers practical methodologies for regional earthquake

prevention and mitigation planning and may provide oppotunities for substantially enhancing urban resilience. METHODS FAULT GEOMETRY AND NUCLEATION CONDITION In this study, we construct a 3D

nonplanar ZF model based on a study on crustal structure in the Pearl River Delta by the Guangdong Earthquake Agency10. In this model, the fault extends approximately 93 km with a main 155°

strike (Fig. 2a). This fault model is set to be vertical, and the fault width along the dip is 10 km according to geological surveys, gravity, and magnetic prospecting results

interpretation5,33. Moreover, the fault plane is embedded in an elastic and heterogeneous media model. Figure 2b shows the interpolated shear wave speed distribution on the ZF plane. To

initiate the dynamic rupture, we place a circular nucleation patch on the ZF plane, and the radius of the nucleation patch is 2 km. Within the nucleation patch, the initial shear stress is

slightly higher (0.5%) than the fault strength. To comprehensively evaluate the earthquake hazards and losses of the ZF, we anchored six nucleation zones (Fig. 2a) according to the regions

with large shear stress of the fault, all of which are 7.5 km deep. We used MATLAB R2020b to complete the modeling (fault geometry, stress, velocity medium, etc.). INITIAL STRESSES We adopt

the triaxial stress scheme to configure the initial stress on the fault plane, which consists of three principal compressional stresses, namely, the maximum horizontal principal stress

(_σ__H_), the minimum horizontal principal stress (_σ__h_), and the vertical principal stress (_σ__V_). We set _σ__V_ to be depth-dependent: it increases with depth from 0 to 5 km (_σ__V_ = 

0.8·_ρgh_, where _ρ_ is the density, _g_ is the gravitational acceleration, and _h_ is the depth) and then remains constant (_σ__V_ = 0.8·_ρg_·5000). _σ__H_ and _σ__h_ are determined by the

relative magnitudes of the three principal stresses. The ratio between these stresses ((_σ__V_ -_σ__H_)/(_σ__h_ -_σ__H_)) in the research area is 0.434. Moreover, the lateral pressure

coefficient (_σ__h/__σ__V_) is assigned as 0.64 according to the stress states in the upper crust35. Therefore, we have _σ__H_ = 1.24·_σ__V_ and _σ__h_ = 0.64·_σ__V_. These depth-dependent

principal stresses are schematic in Fig. 2c (left column). The initial normal stress and shear stress are revealed by projecting the triaxial stress onto the fault plane, and the projecting

results depend on the angle between the fault strike and the azimuth of the maximum principal stress. The focal mechanisms inverted show that the azimuth of the maximum horizontal principal

stress in the research area is approximately N115°5,34. According to this stress scheme, the ZF is characterized by left-lateral strike-slip movement (red arrows in Fig. 2a). By applying the

above stress parameters, the resolved normal stress (_T__n_) and shear stress (_T__s_) on the ZF are presented in Fig. 2d, e. FRICTION PARAMETERS We use the slip-weakening law36 to control

the frictional behavior of the ZF during the dynamic rupture simulation, which consists of three key parameters: the static frictional coefficient (_μ__s_), the dynamic frictional

coefficient (_μ__d_), and the critical distance (_D__c_). The static and dynamic frictional coefficients are 0.4 and 0.24, respectively. These frictional coefficients are obtained by tuning

for a maximum moment magnitude of approximately 7. We adopt a depth-dependent critical distance, in which the _D__c_ decreases from 1.2 m at the surface to 0.4 m at a depth greater than 5 km

(right column in Fig. 2c). Figure 2f shows the ratio between the shear and the normal stress (_μ__0_) on the ZF, which reveals the fault potential rupture ability. The closer _μ__0_ is to

the _μ__s_, the easier it is for the rupture to initiate. Moreover, with the friction parameters, we can calculate the potential dynamic stress drop distribution on the fault plane using the

initial shear stress minus the residual fault strength (_T__s_ - _μ__d_·_T__n_). The calculated potential dynamic stress drop is presented in Fig. 2g. CG-FDM The study utilized the

curved-grid finite difference method (CG-FDM)27 to conduct dynamic rupture and strong ground motion simulations. The CG-FDM effectively addresses traction-free boundary conditions in

simulations of seismic wave propagation. It employs curved grids to discretize the computational volume and, therefore, can incorporate irregular topography of the simulated research area.

The split-node method was introduced to implement the fault boundary conditions and then perform the spontaneous dynamic rupture simulation of the nonplanar faults28. The CG-FDM is

distinguished by its computational efficiency and versatility in simulating complex fault geometries, including faults with stepovers and surface roughness. It has been validated against the

benchmark problems of the Southern California Earthquake Center-U.S. Geological Survey (SCEC-USGS) dynamic rupture code validation project37. Furthermore, the CG-FDM has been widely

validated by simulating the dynamic rupture process of natural earthquakes38,39,40,41 and scenario earthquakes42,43,44. In fault dynamic rupture simulations, the computational domain is

configured in 3D, aligned with the fault strike, normal, and dip directions, with dimensions of approximately 100 km, 20 km, and 14 km, respectively. The grid spacing is set at 50 m, with a

time step of 0.004 s over 8000 steps, thus totaling a simulation duration of 32 s. In strong ground motion simulations, the 3D computational domain spans 184 km, 223 km, and 40 km in

latitude, longitude, and vertical depth, respectively. The spatial resolution is set at 100 m, resulting in a computational grid scale of 1840 × 2230 × 400. The time step is established at

0.0065 seconds, with a total of 20,000 steps, culminating in a simulation duration of 130 seconds. The calculations fully incorporate the impact of the 3D heterogeneous velocity structure

and undulating topography. Topography data are sourced from the Shuttle Radar Topography Mission (SRTM90), with an approximate resolution of 90 m45. Bathymetry data are sourced from the

General Bathymetric Chart of the Oceans (GEBCO). We employed a 3D high-resolution crustal shear wave velocity (Vs) model developed by the Guangdong Earthquake Agency, China, based on ambient

noise tomography in the Pearl River Delta46. The datasets comprise approximately 30 days of seismic ambient noise recorded by 124 land seismic stations and 28 Ocean Bottom Seismometers

(OBS). The minimum shear wave velocity at the surface of the model domain is 2445 m/s, which increases with depth to 4550 m/s. According to the requirements of the CG-FDM for defining the

effective wavenumber with six points per wavelength, the upper limit of the effective frequency for ground motion simulation in this study is approximately 4.1 Hz. Figure 1 shows a

cross-section of the velocity structure. Figure 2b illustrates the interpolated shear wave speed distribution on the ZF plane. EARTHQUAKE LOSSES ESTIMATION This study utilizes gridded

population datasets produced by WorldPop, with a spatial resolution of 100 m47. Provincial fixed capital stock data are calculated using the perpetual inventory method29,30 and are gridded

based on nighttime light data48. The research employs the earthquake fatality estimation model49 and the earthquake economic loss estimation model50 based on the development of regression

analysis of earthquake disaster loss datasets in China from 1949 to 201951. The following equations are used for estimating fatalities or economic losses:

$$\begin{array}{l}{E}_{F\,or\,E}=\mathop{\sum}\limits_{grid}{r}_{F\,or\,E}(Int)\times {\left(\frac{HD{I}_{regression-year}}{HD{I}_{event-year}}\right)}^{\!\!h}\\\qquad\qquad\times\,

Exposure(gird,event-year,Pop\,or\,Fcs)\end{array}$$ (1) $$P(a \,<\, losses\le b)=\Phi \left[\frac{{\mathrm{ln}}\,b-\,{\mathrm{ln}}({E}_{F\,or\,E})}{\zeta }\right]-\Phi

\left[\frac{{\mathrm{ln}}\,a-\,{\mathrm{ln}}({E}_{F\,or\,E})}{\zeta}\right]$$ (2) where \({E}_{F\,or\,E}\) represents fatalities estimated (\({E}_{F}\)) or economic losses estimated

(\({E}_{E}\)) for an earthquake. \({r}_{F\,or\,E}(Int)\) represents a fatality or economic loss ratio function based on seismic intensity (_Int_). The Human Development Index (HDI) is the

indicator the United Nations Development Programme proposed for comprehensively evaluating the social development level. \(Exposure(gird,event-year,PoporFcs)\) represents the current year’s

total population (_Pop_) or fixed capital stock (_Fcs_) at the grid point within the calculation domain. Φ is the normal cumulative distribution function. \({\mathrm{ln}}({E}_{F\,or\,E})\)

is the expectation of the distribution, calculated by the loss estimation model with Eq. (1). _ζ_ is the distribution’s standard deviation, representing the uncertainty of the loss

estimation model. The reasonableness of the model has been confirmed through normality tests and has been subject to uncertainty analysis49,50. For the temporal dimension, which encompasses

potential impacts on earthquake disaster loss rates due to societal development factors such as building fortification, emergency preparedness, and public awareness, the model adjusts using

the HDI. The model employs a zonal strategy to address the differences in seismic structures and economic development that might affect earthquake disaster loss rates in the spatial

dimension. This strategy integrates seismological and sociological factors to construct a set of loss rate functions with distinct mathematical expressions for various sub-regions in China.

In this study, based on the location of the ZF and the applicable intensity range, Scenarios 1-6 are categorized under the South China fatality assessment model. Scenarios 1-3 fall under the

East China economic loss assessment model, while Scenarios 4-6 are grouped under the South China economic loss assessment model. The corresponding standard deviation (_ζ_) and HDI term

coefficients (_h_) are 1.3539, 1.0; 1.9267, 1.7; and 1.6914, 1.6, respectively36,37. The fatality and economic loss ratio functions are as follows: $$\log ({r}_{F})=-12.625+1.145\cdot Int$$

(3) $${r}_{E}=\varPhi \left[\frac{1}{0.251}\cdot \,\mathrm{ln}\left(\frac{Int}{10.853}\right)\right]$$ (4) $$\log ({r}_{E})=-57.904\cdot \exp (-0.540\cdot Int)$$ (5) Where \({r}_{F}\) and

\({r}_{E}\) represent the fatality ratio and economic loss ratio, respectively. Equation (4) applies to Scenarios 1-3, and Eq. (5) applies to Scenarios 4-6. This approach has been

successfully applied to studies of natural earthquakes and setting earthquakes in China in recent years41,44,52,53. DATA AVAILABILITY The topographic data (SRTM90) are available from

https://srtm.csi.cgiar.org. The bathymetry data are available from https://www.gebco.net. The gridded population datasets produced by WorldPop are available from https://www.worldpop.org/.

The gridded nighttime light data are available from https://doi.org/10.7910/DVN/YGIVCD. Data on the faults, velocity medium, and fixed capital stock are available from the authors upon

request. CODE AVAILABILITY The codes of CG-FDM are available from https://zenodo.org/records/13358124. REFERENCES * Huang, Q. F., Wang, L. H., Jia, B. H., Lai, X. & Peng, Q. Impact of

Climate Change on the Spatio-Temporal Variation in Groundwater Storage in the Guangdong-Hong Kong-Macao Greater Bay Area. _Sustainability-Basel_ 15, 10776 (2023). * World Bank Group. World

Bank Open Data. https://data.worldbank.org/ (2024). * Zhu, H. A Comparative Study of Population Development in the Nine Cities of the Guangdong-Hong Kong-Macao Greater Bay Area over Ten

Years under New Urbanization — An Analysis Based on the Seventh Census Data. _Future City Design and Management_, 46-48 (2022). * Song, F. et al. Late Quaternary Vertical Dislocation Rate on

Several Faults in the Zhujiang Delta Area. _Seismology And Geology_, 203–210 (2003). * Liang, G. _Active Fault Mapping and Seismic Hazard Evaluation in Guangzhou_. (Science Press, 2013). *

Scawthorn, C., Nishino, T., Schencking, J. C. & Borland, J. Kantō Daikasai: The Great Kantō Fire Following the 1923 Earthquake. _Bull. Seismological Soc. Am._ 113, 1902–1923 (2023).

Article  Google Scholar  * Maly, E. & Suppasri, A. The Sendai Framework for Disaster Risk Reduction at Five: Lessons from the 2011 Great East Japan Earthquake and Tsunami. _Int J. Disast

Risk Sc._ 11, 167–178 (2020). Article  Google Scholar  * Scawthorn, C., O’Rourke, T. D. & Blackburn, E. T. The 1906 San Francisco earthquake and fire-enduring lessons for fire

protection and water supply. _Earthq. Spectra_ 22, S135 (2006). Article  Google Scholar  * Haynes, B. E. et al. Medical Response to Catastrophic Events - California Planning and the

Loma-Prieta Earthquake. _Ann. Emerg. Med_ 21, 368–374 (1992). Article  CAS  Google Scholar  * Ye, X. et al. Crustal structure and tectonic attribute revealed by a deep seismic sounding

profile of Dinghu-Gaoming-Jinwan in the pearl river delta. _Chin. J. geophysics_ 63, 1959–1969 (2020). Google Scholar  * Yan, Y., Li, Y., Liu, T., Ouyang, L. & Guan, Y. Experimental

Research of Active Fault in Guangdong Province and Its Adjacent Region. _SOUTH China J. Seismol._, 14–25 (2006). * Mcguire, R. K. Probabilistic Seismic Hazard Analysis and Design Earthquakes

- Closing the Loop. _Bull. Seismol. Soc. Am._ 85, 1275–1284 (1995). Article  Google Scholar  * Bommer, J. J. Deterministic vs. probabilistic seismic hazard assessment: An exaggerated and

obstructive dichotomy. _J. Earthq. Eng._ 6, 43–73 (2002). Article  Google Scholar  * Douglas, J. & Aochi, H. A Survey of Techniques for Predicting Earthquake Ground Motions for

Engineering Purposes. _Surv. Geophysics_ 29, 187–220 (2008). Article  Google Scholar  * Infantino, M., Mazzieri, I., Özcebe, A. G., Paolucci, R. & Stupazzini, M. 3D Physics‐Based

Numerical Simulations of Ground Motion in Istanbul from Earthquakes along the Marmara Segment of the North Anatolian Fault. _Bull. Seismol. Soc. Am._ 110, 2559–2576 (2020). Article  Google

Scholar  * Stupazzini, M., Infantino, M., Allmann, A. & Paolucci, R. Physics-based probabilistic seismic hazard and loss assessment in large urban areas: A simplified application to

Istanbul. _Earthq. Eng. Struct. Dyn._ 50, 99–115 (2021). Article  Google Scholar  * Erdik, M. Earthquake risk assessment. _Bull. Earthq. Eng._ 15, 5055–5092 (2017). Article  Google Scholar 

* Chen, Y. L. & Baker, J. W. Spatial Correlations in CyberShake Physics-Based Ground-Motion Simulations. _Bull. Seismol. Soc. Am._ 109, 2447–2458 (2019). Article  Google Scholar  * Du,

W. Q. & Ning, C. L. Modeling spatial cross-correlation of multiple ground motion intensity measures (SAs, PGA, PGV, Ia, CAV, and significant durations) based on principal component and

geostatistical analyses. _Earthq. Spectra_ 37, 486–504 (2021). Article  Google Scholar  * Infantino, M., Smerzini, C. & Lin, J. Y. Spatial correlation of broadband ground motions from

physics-based numerical simulations. _Earthq. Eng. Struct. Dyn._ 50, 2575–2594 (2021). Article  Google Scholar  * Rodgers, A. J., Pitarka, A., Petersson, N. A., Sjögreen, B. & McCallen,

D. B. Broadband (0-4 Hz) Ground Motions for a Magnitude 7.0 Hayward Fault Earthquake With Three-Dimensional Structure and Topography. _Geophys. Res. Lett._ 45, 739–747 (2018). Article 

Google Scholar  * Withers, K. B., Olsen, K. B., Day, S. M. & Shi, Z. Q. Ground Motion and Intraevent Variability from 3D Deterministic Broad-band (0-7.5 Hz) Simulations along a Nonplanar

Strike-Slip Fault. _Bull. Seismol. Soc. Am._ 109, 229–250 (2019). Article  Google Scholar  * Denolle, M. A., Dunham, E. M., Prieto, G. A. & Beroza, G. C. Strong Ground Motion Prediction

Using Virtual Earthquakes. _Science_ 343, 399–403 (2014). Article  CAS  Google Scholar  * Yu, H., Zhang, W., Zhang, Z., Li, Z. & Chen, X. Investigation on the Dynamic Rupture of the

1970 Ms 7.7 Tonghai, Yunnan, China, Earthquake on the Qujiang Fault. _Bull. Seismol. Soc. Am._ 110, 898–919 (2020). Article  Google Scholar  * Wang, X. C., Wang, J. T. & Zhang, C. H.

Deterministic full-scenario analysis for maximum credible earthquake hazards. _Nat. Commun._ 14, (2023). * Yu, H., Zhang, Z., Hu, F., Xu, D. & Chen, X. Estimation of the Nucleation

Location and Rupture Extent of the 1850 Xichang, Sichuan, China, Earthquake by Dynamic Rupture Simulations on a Multi-Segment Stepover Structure. _Earth Space Sci._ 10, e2022EA002775 (2023).

Article  Google Scholar  * Zhang, W. & Chen, X. Traction image method for irregular free surface boundaries in finite difference seismic wave simulation. _Geophys. J. Int._ 167, 337–353

(2006). Article  Google Scholar  * Zhang, Z., Zhang, W. & Chen, X. Three-dimensional curved grid finite-difference modelling for non-planar rupture dynamics. _Geophys. J. Int._ 199,

860–879 (2014). Article  Google Scholar  * Shan, H. J. Reestimating the capital stock of China: 1952-2006. _J. Quant. Tech. Econ._ 10, 17–31 (2008). Google Scholar  * Zhang, J. Estimation of

China’s provincial capital stock (1952–2004) with applications. _J. Chin. Economic Bus. Stud._ 6, 177–196 (2008). Article  Google Scholar  * Boore, D. M., Stewart, J. P., Seyhan, E. &

Atkinson, G. M. NGA-West2 equations for predicting PGA, PGV, and 5% damped PSA for shallow crustal earthquakes. _Earthq. Spectra_ 30, 1057–1085 (2014). Article  Google Scholar  * Campbell,

K. W. & Bozorgnia, Y. NGA ground motion model for the geometric mean horizontal component of PGA, PGV, PGD and 5% damped linear elastic response spectra for periods ranging from 0.01 to

10 s. _Earthq. Spectra_ 24, 139–171 (2008). Article  Google Scholar  * Ren, Z. _Deep Fault System and Seismotectonics in the Coastal Region of South China_ (Guangdong Science and Technology

Press, 2016). * Kang, Y., Yang, X., Chen, X., Chen, G. & Zheng S. Inversion of stress field in Guangdong and its adjacent area. _Acta Seismol. Sin._, 59–66+105 (2008). * Wang, Y., Cui,

X., Hu, X. & Xie, F. Study on the stress state in upper crust of China mainland basedon in situ stress measurements. _Chin. J. Geophys._ 55, 3016–3027 (2012). Google Scholar  * Ida, Y.

Cohesive force across the tip of a longitudinal-shear crack and Griffith’s specific surface energy. _J. Geophys. Res. (1896–1977)_ 77, 3796–3805 (1972). Article  Google Scholar  * Harris, R.

A. et al. A Suite of Exercises for Verifying Dynamic Earthquake Rupture Codes. _Seismol. Res. Lett._ 89, 1146–1162 (2018). Article  Google Scholar  * Zhang, Z., Zhang, W., Xin, D., Chen, K.

& Chen, X. A Dynamic-Rupture Model of the 2019 Mw 7.1 Ridgecrest Earthquake Being Compatible with the Observations. _Seismol. Res. Lett._, (2020). * Chen, K., Zhang, Z., Liang, C., Xue,

C. & Liu, P. Kinematics and Dynamics of the 24 January 2020 Mw 6.7 Elazig, Turkey Earthquake. _Earth Space Sci._ 7, e2020EA001452 (2020). Article  Google Scholar  * Xu, D. et al. The

2016 Menyuan Earthquake: The Largest Self-Arrested Crustal Earthquake Ever Observed. _Geophys. Res. Lett._ 50, e2023GL103556 (2023). Article  Google Scholar  * Wang, W. et al. Rapid

estimation of disaster losses for the M6.8 Luding earthquake on September 5, 2022. _Sci. China Earth Sci._ 66, 1334–1344 (2023). Article  Google Scholar  * Zhang, Z., Zhang, W., Chen, X.,

Li, P. & Fu, C. Rupture Dynamics and Ground Motion from Potential Earthquakes around Taiyuan, China. _Bull. Seismol. Soc. Am._ 107, 1201–1212 (2017). Article  Google Scholar  * Gu, Y.,

Zhang, Z., Wang, W. & Wang, Z. Dynamic rupture simulations based on interseismic locking models—taking the Suoerkuli section of the Altyn Tagh fault as an example. _Geophys. J. Int._

234, 1737–1751 (2023). Article  Google Scholar  * Li, Y., Wang, Z., Zhang, Z., Gu, Y. & Yu, H. A Physics-Based Seismic Risk Assessment of the Qujiang Fault: From Dynamic Rupture to

Disaster Estimation. _Int J. Disast Risk Sc._ 15, 165–177 (2024). Article  CAS  Google Scholar  * Reuter, H., Nelson, A. & Jarvis, A. An Evaluation of Void-Filling Interpolation Methods

for SRTM Data. _Int. J. Geogr. Inf. Sci._ 21, 983–1008 (2007). Article  Google Scholar  * Lv, Z. Y., Ye, X. W., Zhang, Y. X., Qian, Y. P. & Wen, G. G. 3-D crustal shear wave velocity and

azimuthally anisotropic structure of the Pearl River Delta onshore-offshore area and its tectonic implications. _Tectonophysics_ 868, (2023). * Bondarenko, M., Kerr, D., Sorichetta, A.

& Tatem, A. Census/projection-disaggregated gridded population datasets, adjusted to match the corresponding UNPD 2020 estimates, for 51 countries across sub-Saharan Africa using

building footprints. https://www.worldpop.org/doi/1...ON/WP00683 (2020). * Chen, Z. et al. An extended time-series (2000-2023) of global NPP-VIIRS-like nighttime light data. V5 edn. Harvard

Dataverse (2020). * Li, Y., Xin, D. & Zhang, Z. A rapid-response earthquake fatality estimation model for mainland China. _Int J. Disast Risk Re_ 66, 102618 (2021). Google Scholar  * Li,

Y., Xin, D. & Zhang, Z. Estimating the economic loss caused by earthquake in Mainland China. _Int J. Disast Risk Re_ 95, 103708 (2023). Google Scholar  * Li, Y., Zhang, Z. & Xin, D.

A Composite Catalog of Damaging Earthquakes for Mainland China. _Seismol. Res. Lett._ 92, 3767–3777 (2021). Article  Google Scholar  * Li, Y., Zhang, Z., Wang, W. & Feng, X. Rapid

Estimation of Earthquake Fatalities in Mainland China Based on Physical Simulation and Empirical Statistics—A Case Study of the 2021 Yangbi Earthquake. _Int. J. Environ. Res. Public Health_

19, 6820 (2022). Article  Google Scholar  * Wang, Z., Li, Y., Wang, W., Zhang, W. & Zhang, Z. Revisiting Paleoearthquakes with Numerical Modeling: A Case Study of the 1679 Sanhe–Pinggu

Earthquake. _Seismol. Res. Lett._ 94, 720–730 (2022). Article  Google Scholar  Download references ACKNOWLEDGEMENTS This work is supported by the Guangdong Basic and Applied Basic Research

Foundation (Grant No. 2023A1515110278), the National Natural Science Foundation of China (Grant No. U1901602), National Key R&D Program of China (2023YFC3007404), Guangdong Provincial

Pearl River Talents Program (2019QN01G801), and Guangdong Provincial Key Laboratory of Geophysical High-resolution Imaging Technology (2022B1212010002). The simulation is supported by the

Center for Computational Science and Engineering at the Southern University of Science and Technology. We thank the GEBCO Compilation Group (2024) GEBCO 2024 Grid

(doi:10.5285/1c44ce99-0a0d-5f4f-e063-7086abc0ea0f). We thank Dr. Zuoyong Lv (Guangdong Earthquake Agency) for providing the high-resolution shear velocity model that was used in our

simulation. AUTHOR INFORMATION AUTHORS AND AFFILIATIONS * Joint International Research Laboratory of Catastrophe Simulation and Systemic Risk Governance, Beijing Normal University, Zhuhai,

China Yilong Li, Ming Wang & Kai Liu * School of National Safety and Emergency Management, Beijing Normal University, Zhuhai, China Yilong Li * Department of Earth and Space Sciences,

Southern University of Science and Technology, Shenzhen, China Yilong Li, Houyun Yu, Zhenguo Zhang, Wei Zhang & Xiaofei Chen * Guangdong Provincial Key Laboratory of Geophysical

High-resolution Imaging Technology, Southern University of Science and Technology, Shenzhen, China Houyun Yu & Xiaofei Chen * School of National Safety and Emergency Management, Beijing

Normal University, Beijing, China Ming Wang & Kai Liu * China Earthquake Administration Key Laboratory of Earthquake Monitoring and Disaster Mitigation Technology, Guangdong Earthquake

Agency, Guangzhou, China Xiuwei Ye * Guangdong Science and Technology Collaborative Innovation Center for Earthquake Prevention and Disaster Mitigation, Guangzhou, China Xiuwei Ye * Key

Laboratory of Disaster Forecast and Control in Engineering of the Ministry of Education, Urban Earthquake Safety Laboratory, Guangzhou, China Xiuwei Ye Authors * Yilong Li View author

publications You can also search for this author inPubMed Google Scholar * Houyun Yu View author publications You can also search for this author inPubMed Google Scholar * Ming Wang View

author publications You can also search for this author inPubMed Google Scholar * Xiuwei Ye View author publications You can also search for this author inPubMed Google Scholar * Kai Liu

View author publications You can also search for this author inPubMed Google Scholar * Zhenguo Zhang View author publications You can also search for this author inPubMed Google Scholar *

Wei Zhang View author publications You can also search for this author inPubMed Google Scholar * Xiaofei Chen View author publications You can also search for this author inPubMed Google

Scholar CONTRIBUTIONS Y.L. and H.Y. conceptualized and designed the research, conducted the analyses, prepared the figures, and wrote the initial manuscript. M.W., X.Y., K.L., Z.Z., W.Z.,

and X.C. contributed to the revision process and commented on different manuscript versions. CORRESPONDING AUTHOR Correspondence to Houyun Yu. ETHICS DECLARATIONS COMPETING INTERESTS The

authors declare no competing interests. ADDITIONAL INFORMATION PUBLISHER’S NOTE Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional

affiliations. SUPPLEMENTARY INFORMATION SUPPLEMENTARY INFORMATION RIGHTS AND PERMISSIONS OPEN ACCESS This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives

4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original

author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted

material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise

in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the

permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/. Reprints and

permissions ABOUT THIS ARTICLE CITE THIS ARTICLE Li, Y., Yu, H., Wang, M. _et al._ The physics-based deterministic scenarios for earthquake hazards and losses of the Zhujiangkou fault in

southern China. _npj Nat. Hazards_ 2, 6 (2025). https://doi.org/10.1038/s44304-024-00054-3 Download citation * Received: 18 September 2024 * Accepted: 19 December 2024 * Published: 17

January 2025 * DOI: https://doi.org/10.1038/s44304-024-00054-3 SHARE THIS ARTICLE Anyone you share the following link with will be able to read this content: Get shareable link Sorry, a

shareable link is not currently available for this article. Copy to clipboard Provided by the Springer Nature SharedIt content-sharing initiative