Contemporary crustal kinematics in the Guangdong-Hong Kong-Macao Greater Bay Area, SE China: Implications for the geothermal resource exploration

Fault kinematics plays an important role in estimating the stress state and permeability of faults which are controlling factors in the formation of a geothermal system. However, there are very few studies on the kinematic characteristics of the major faults in the Guangdong-Hong Kong-Macao Greater Bay Area (GBA), SE China. To obtain a better understanding of the fault kinematics, we establish a comprehensive 3D geomechanical model for the GBA. Our results show that the NE-trending faults in the west of the Pearl River Estuary (PRE) usually slip faster than the faults striking in the ENE-WSW direction to the east of the PRE. The NW-trending faults have the lowest modeled fault slip rate. Slip rates of the faults are generally low with a maximum value of 0.12 mm/a occurring on the northeastern segment of the Wuchuan-Sihui fault. The NE-trending faults display sinistral motion, while the ENE-trending faults are dextral. The opposite slip senses on these two fault groups are inferred to be caused by the lateral variations in the crustal stress. Based on the analysis of contemporary kinematics and the heat flow in the GBA, we suggest that the fault segments with relatively high slip rates, such as the north-eastern segment of the Wuchuan-Sihui fault, the Kaiping fault, the Enping fault, and the middle segments of the Wuhua-Shenzhen, Zijin-Boluo, and Heyuan faults, have a high prospect for geothermal resources. The intersections of the NW-trending extensional faults and the NE-/ENE-trending faults could also be potential areas of interest for future geothermal exploration.


Introduction
The South China Block, located at the southeast margin of the Eurasian Plate (Fig. 1), has attracted widespread interest in the field of geodynamics due to its extensive magmatism and super-giant ore deposits (Dong et al., 2020;Li et al., 2020;Li, 2000;Zhou et al., 2006).Geochronological studies suggest that this block was initially formed by the Neoproterozoic collision of the Yangtze and Cathaysia blocks and then experienced three tectonic-magmatic events, which occurred sequentially in Silurian, Early to Middle Triassic and Late Mesozoic periods (Shu et al., 2021).Among these events, the Late Mesozoic tectonism-magmatism is considered to be caused by the northwestward subduction of the Paleo-Pacific Plate with a slab roll-back (Dong et al., 2020;Zhou et al., 2006;Zhou and Li, 2000).The back-arc extensional setting and the asthenosphere upwelling in response to the roll-back of the subducting Paleo-Pacific plate jointly promoted the magma emplacement and formed the famous Yanshanian (Middle Jurassic-Middle Cretaceous, 180-80 Ma) granitoids in the South China Block (Dong et al., 2020;Zhou et al., 2006).The widely distributed granitoids provide stable heat sources to make the South China Block a region rich in geothermal resources, especially in its southeastern part (i.e.Cathaysia Block).The latest terrestrial heat flow map of China shows that the Cathaysia Block has a high average heat flow of 69.4 mW/m 2 with local anomalies over 100 mW/m 2 (Jiang et al., 2019), which is generally consistent with the distribution of the Yanshanian granitoids in SE China (Zhang et al., 2018a).
As an important economic center in SE China, Guangdong Province is characterized by high terrestrial heat flow values and abundant geothermal resources, and it has the first Chinese geothermal power station that was built in 1970 (Huang, 2012).Previous geothermal surveys have identified over 300 hot springs within the province, with low to moderate temperatures mainly ranging from 25 to 118 • C (Wang, 2018;Wei et al., 2021).In 2019, an administrative region called the Guangdong-Hong Kong-Macao Greater Bay Area (GBA) (Fig. 1), including Hong Kong, Macao, and nine cities in Guangdong, was established by the Chinese government to accelerate the economic development of this region.This has led to a strong demand for energy in the area.However, the abundant geothermal energy resources in the GBA have not been fully exploited yet, and one of the barriers to the exploration of this resource is that the kinematic characteristics of faults in this region are still not clear.
Fault kinematics plays an important role in estimating the fault stress state, which further governs the fault permeability and eventually affects the migration of fluids through the subsurface (Jolie et al., 2021;Siler et al., 2018).For example, for a right-stepping fault system, a dextral slip will produce extensional structures such as normal faults and extension fractures in the step-over (e.g.Fig. S1a), while a sinistral slip produces contractional structures (Fossen, 2010).Fluids in an extensional regime can migrate more easily from deeper to shallower depths to form a geothermal field than in a contractional regime.In such a setting, compression would reduce the permeability of faults and tend to close flow channels to impede fluid migration (Tanikawa et al., 2010;Evans et al., 1997).In addition, fault permeability is also affected by the fault activity (Tanikawa et al., 2010).High fault slip rates would be favorable for the generation of new fractures or make pre-existing ones remain open.Since the geothermal resources in the GBA depend mainly on deep hydrothermal circulation systems controlled by faults (e.g.Lu and Liu, 2015;Wang et al., 2018;Wei et al., 2023;Wei et al., 2021), a good understanding of its fault kinematics is essential to provide a solid foundation for better geothermal resource assessment in the GBA.

Fault kinematics in the GBA
Active faults are widely developed in the GBA and can be roughly divided into three groups based on their orientations: the NE-SW group located in the western part of the study area, the ENE-WSW group in the eastern and southern parts, and the NW-SE group in the central part, specifically the Pearl River Estuary area (Fig. 1).Several previous qualitative studies suggest that the activity of these faults is relatively low.For example, Chen et al. (2002) measured the soil gas emission data along the Guangzhou-Conghua fault (F25), Luofushan fault (F12), Xijiang fault (F24), and Baini-Shawan fault (F23) and pointed out that Fig. 1.Map of active faults, earthquakes and heat flow measurements of the Guangdong-Hong Kong-Macao Greater Bay Area and its adjacent area.Fault traces are from the Seismic Active Fault Survey Data Center (activefault-datacenter.cn), thick and thin lines denote lithospheric scale and crustal scale faults, respectively; earthquake data are from the National Earthquake Data Center (data.earthquake.cn);heat flow measurements are from the International Heat Flow Commission (IHFC, ihfc-iugg.org)and Tang et al. (2014).The inset shows the tectonic setting of the study area, in which the black box represents the location of Fig. 1 and the spatial extent of the model region.
these faults have slight activities.This understanding is further strengthened by comprehensive geological surveys and soil gas measurements carried out specifically along the Baini-Shawan fault (F23) (Dong et al., 2016).The low fault activity presents a challenge in acquiring precise fault displacements in the field, thereby impeding the quantitative determination of fault slip rates.Different from these aforementioned studies, Song et al. (2003) used differences in elevation within the same stratum on both sides of the fault in the GBA, which are revealed through borehole analysis, and preliminarily calculated vertical fault slip rates ranging from 0.14 to 0.47 mm/a.Yu et al. (2016) also compiled stratigraphic profiles from boreholes and seismic reflections and suggested that the ENE-and NW-trending faults are the primary active faults.Although a few studies on fault activities have been conducted quantitatively in the region (e.g.Lei et al., 2018;Sun et al., 2007;Wang et al., 2011), they have focused only on two faults, namely, the Wuhua-Shenzhen fault (F10) and the Guangzhou-Conghua fault (F25) (see Fig. 1).The detailed kinematic characteristics of the various other faults remain unclear.
Numerical modeling provides a powerful method to study the largescale kinematics of complex fault systems.It can integrate the observational data from different geoscience fields such as geology, geophysics, and geodesy to derive physics-based fault slip rates and crustal deformation considering the acting forces and mechanical properties of faults and rocks.Thus, the modeling approach can provide detailed and consistent estimates for the kinematics even of a complex fault system.Several modeling studies on fault slip rate have been conducted in the Marmara Sea region, northwestern Turkey (Hergert and Heidbach, 2010;Hergert et al., 2011) and the eastern Tibetan Plateau (Li et al., 2022;Li et al., 2021).The modeled rates in these studies have shown to be in good agreement with observations.However, in the GBA, to our knowledge, no similar work has been conducted on the kinematics of the faults except a couple of two-dimensional (2D) (Chen et al., 2014;Wen et al., 2001) or three-dimensional (3D) models with simple fault geometries (Dong et al., 2021) focusing on crustal stress.
In this paper, we present a 3D geomechanical model for the GBA incorporating a complex 3D fault system and inhomogeneous rock properties as well as topography and gravity.Detailed contemporary kinematics of the faults and the crust in the study area are derived from the numerical simulation and validated by comparison with modelindependent data.Finally, based on the model, a comprehensive analysis of the implications for geothermal exploration in the study area is given.

Model concept and input
The modeling process in this study is composed of five steps, which are model design (including a 3D fault system and model geometry), model discretization, computation, results calibration, and analysis of results (Fig. 2).In the beginning, we construct a 3D geometry of the study area with a complex fault system by using the CAD software Rhinoceros®v7.The detailed geometrical parameters used to constrain the model in this step are extensively compiled from published geological and geophysical literature and reports.The 3D model volume is discretized into ~ 4 million linear tetrahedral elements with a resolution of 1-2 km on the faults and ~ 5 km near the model boundaries using the meshing software Hypermesh®v2019.After assigning inhomogeneous material properties to the finite element model and applying values of initial stress and gravity to the model volume, the numerical simulation was carried out with the finite element analysis software Abaqus®v2019 using displacement boundary conditions.The modeled displacements at the surface and on the faults are then compared with model-independent data, such as the GPS observations and geological or geodetic fault slip rates.Finally, we analyze the modeled kinematic results and discuss their implications on the geothermal resource exploration in the GBA.

Model geometry
The model geometry has a cuboid shape with a spatial extent of 572 km in E-W (110.5 • -116 • E) and 490 km in S-N (20.5 • -25 • N) direction.The model thickness is ~ 38 km (Fig. 3).From top to bottom, the model contains four units, namely the upper, middle, and lower crust, and the top layer of the upper mantle.These units are bounded by five surfaces: the topography, three interfaces of the upper, middle, and lower crust, Fig. 2. Model setup and workflow of the geomechanical modeling in this study (modified from Hergert et al., 2015).The partial differential equations of the equilibrium of forces in the 3D model are solved using the finite element method (σ ij stress tensor, x j Cartesian coordinates, ρ density, and ρx i body forces).
X. Li et al. and the bottom of the model at 38 km depth.The relief of the topography and crustal interfaces are taken from the GTOPO30 global digital elevation model (USGS) and the CRUST1.0model (Laske et al., 2013), respectively.
For the 3D fault system embedded in the model, the upper terminations of the faults are constrained by the fault traces from the Seismic Active Fault Survey Data Center (https://www.activefault-datacenter. cn), which incorporates the new findings of active fault investigations in recent years and forms the latest active fault database in China.The fault geometries at depth are determined from deep seismic sounding profiles, focal mechanism solutions, seismic tomography, earthquake relocations, etc.For an individual fault, if none of the references  mentioned above is available to constrain its subsurface geometry, we use the shallow fault dip generally obtained by surface fault surveys to extend the fault plane to a greater depth.Based on geological and geophysical data in the study area, the faults implemented in the model can be divided into three types according to their depth extent.One is the lithospheric scale faults that penetrate the entire model down to a depth of ~ 38 km.The other two are both crustal scale faults but terminate differently at the bottom of the upper (~10 km) and middle crust (~20 km), respectively.The detailed geometrical parameters of the faults and the corresponding references are summarized in Table 1.

Rock properties and friction coefficient
In this study, the mechanical behavior of the model units is set to be linear and of static elasticity, an approach frequently used in geomechanical modeling (Ahlers et al., 2021;Hergert et al., 2011;Li et al., 2021;Reiter and Heidbach, 2014).The corresponding rock properties (i. e. Young's modulus E, Poisson's ratio ν and density ρ) are referenced from the CRUST1.0model (Laske et al., 2013), in which the Young's modulus is further converted to the static ones by empirical relations (Brotons et al., 2016).According to the seismic velocity structures obtained by Cao et al. (2014), the model domain can be divided into the normal continental crust of the Cathaysia Block and the thinned continental crust of the South China Sea separated by the Littoral fault (F1, Fig. 1).Table 2 lists all the material parameters used in the model.
The faults implemented in the model are defined by pairs of contact surfaces that can slide relative to each other and are constrained by the Coulomb friction law.In each contact pair, two surfaces are in direct contact without any additional thickness between them, and the elements on both sides of the contact pair are discontinuous.After testing a series of effective friction coefficients ranging from 0 to 1, both 0.01 and 0.02 are optimum values that can minimize the misfit between the measured velocities at the surface and the modeled velocity field (Fig. 4), resulting in nearly identical simulated fault slip rates.The value of 0.02 is used to continue the analysis in the following model.Although the friction coefficient of crustal rocks is usually within the range of 0.6-0.85(Byerlee, 1978), observations at large faults (e.g.Li et al., 2015;Zoback et al., 1987), as well as lots of numerical simulations (He and Lu, 2007;Hergert and Heidbach, 2010;Li et al., 2021), suggest that the effective friction coefficient for large active faults in nature is generally much lower.The low friction coefficient is likely caused by specific materials distributed along faults.For instance, Carpenter et al. (2015) reported a static friction coefficient as low as 0.1 on the San Andreas fault and suggested the abundant magnesium-rich clay localized within the faults leads to the low frictional strength.Additionally, the effect of pore pressure can also contribute to the reduction of the effective friction coefficient.

Gravity and initial stress state
Gravity is applied to the whole model with an acceleration of 9.81 m/s 2 .An appropriate initial stress state is also accounted for in the modeling approach, which is important to simultaneously model both the kinematics and crustal absolute stresses.The initial stress state is calibrated by the semi-empirical horizontal-to-vertical stress ratio k proposed by Sheorey (1994), where E is Young's modulus (GPa) and z is the burial depth (m).Vertical k-z profiles at three test sites (see Fig. 1) are extracted for the calibration.These selected sites are far away from the faults and model boundaries to reduce potential stress perturbations and boundary effects.The good consistency between the modeled initial k-value and the theoretical initial k-value predicted by Eq. ( 1) is shown in Fig. 5.

Kinematic boundary conditions
For the kinematic boundary conditions, we use the hitherto most comprehensive Global Position System (GPS) data compiled by Wang and Shen (2020) to drive the horizontal movement of the model and simulate the tectonic loadings (Fig. 6).Since all the GPS stations in this dataset are in continental China, they cannot constrain the motions of the southern and southeastern model boundaries which are located in the South China Sea.
Although a previous geodetic study has reported a motion rate (10.7 mm/a, relative to stable Eurasia Plate) of the Dongsha Islands, located near the southeastern corner of the model (Yu et al., 1999), the rate is  abnormally high and would drag the model southeastward dramatically and produce a NE-SW-trending maximum horizontal stress (S H ) in the central part of the study area (Fig. S2a), which would contradict the NW-SE-trending S H revealed by focal mechanism solutions (Kang et al., 2008) and in situ stress measurements (Heidbach et al., 2016) (Fig. S2a).
In this study, we assign a relatively low displacement rate (6.9 mm/a) to the southeastern corner, which is inferred from the variation trend of the continental GPS to ensure a suitable stress orientation.This inferred rate, together with the other GPS measurements near the model boundary, can generate a stress field with predominant NW-SE-trending S H , which is consistent with the model-independent stress observations (Fig. S2b).
To check the effect of changing this inferred displacement rate on the modeled kinematic results, two tests with ± 10 % perturbation to the 6.9 mm/a have been conducted.The results indicate that the slip rates and sense of slip on the onshore faults in both tests are almost the same (Fig. S2c and d), which are also similar to the modeled results when the value of 6.9 mm/a is used (see Fig. 7).Therefore, the modeled kinematic results on the onshore faults are marginally affected by the perturbation of the displacement rate at the southeastern corner and the displacement rate of 6.9 mm/a adopted in this study can be considered reasonable.
Then the inferred motion rate and the GPS measurements are jointly used to interpolate the velocities at nodes located on the boundary of the model.A time span is selected to convert the velocities at the boundary nodes to displacement boundary conditions to drive the lateral motion of the model (see the black arrows in Fig. 6).The time span is set to be 600 ka to allow the accumulated displacements at the boundaries to propagate into the model thus generating a proper contemporary state of stress and deformation.In addition, good correlations between the polarization directions of fast shear waves and the GPS measurements in the South China Block (Wu et al., 2007) suggest that the crustal strain in the study area is consistent over the whole depth range of the crust, which provides a reasonable basis for using the GPS measurements at the surface to constrain the model at all depths.The model bottom is fixed only with respect to movements in the vertical direction while allowed to move horizontally (see Fig. 3).The model surface is unconstrained.

Model results
In this study, the GPS data used to validate the model mostly represent the interseismic strain accumulation (Wang and Shen, 2020).Therefore, in order to compare the modeled velocities with GPS measurements, we first lock the faults in the model by assigning to them an infinitely large effective coefficient of friction.The modeled velocities at the exact positions of the GPS stations are extracted from the model and plotted together with the GPS data in Fig. 6.The comparison shows that the modeled crustal velocities in most places are in good agreement with the GPS data.Slight deviations are mostly within the uncertainty ranges of the GPS observations.After a good fit between the modeled velocities and GPS measurements is achieved (Fig. 6), the faults are changed to be unlocked with a low effective coefficient of friction (0.02), and the same boundary conditions as in the previous locked-fault model are imposed to obtain the long-term kinematics of the study area.Fig. 7 shows the modeled horizontal surface velocities and fault slip rates.Overall, the crust in the GBA moves ESE-ward at an average surface velocity of ~ 7.2 mm/a (Figs. 6 and 7).The northern part of the study area exhibits a higher velocity of ~ 7.5 mm/a, gradually decreasing to ~ 7.0 mm/a towards the south, indicating a low velocity gradient, which implies relatively weak crustal deformation within this region.
The faults in the study area can be divided into three groups based on their orientations: NE-SW, ENE-WSW, and NW-SE, located in the western part, the eastern and southern part, and the Pearl River Estuary, respectively.These faults exhibit generally low slip rates (Fig. 7), which is consistent with the previously mentioned weak crustal deformation.Among these three fault groups, the NE-trending faults have the highest slip rates.For example, a prominent high slip rate of approximately 0.12 mm/a is observed on the northeastern segment of the Wuchuan-Sihui fault (F32), gradually decreasing towards the southwest.To the east, the Kaiping and Enping faults (F26, F28) have slip rates ranging from 0.05 to 0.06 mm/a.To the west, the Dashiding fault (F33) slips at a rate of approximately 0.07 mm/a.Furthermore, the discontinuous velocity contour lines intersected by the NE-oriented faults (Fig. 7) indicate that these faults have a significant influence on adjusting crustal deformation due to their relatively high slip rates.In the east of the Pearl River Estuary, the faults are distributed nearly parallel in the ENE-WSW direction with slip rates of < 0.04 mm/a, which are lower than the slip rates inferred for the NE-trending faults.The Littoral fault (F1) and its subfaults located in the South China Sea also exhibit very low slip rates (<0.03 mm/a).The lowest modeled fault slip rate (<0.02 mm/a) occurs on the NW-trending faults that are mainly located in the Pearl River Delta region.Compared with the northwestern segments of these faults, the slip rates on the southeastern segments are extremely low to almost no slip (Fig. 7).
For the fault slip sense, our modeled results show that faults with different orientations exhibit different slip senses.As shown in Fig. 7, the NE-trending faults located in the western part of the study area display sinistral motion, while the ENE-trending faults in the eastern part are dextral.The opposite slip senses on these two fault groups are speculated to be caused by the changes in the crustal stress field, which will be discussed in detail in section 5.2.Unlike the previous two fault groups, the slip sense of the NW-trending faults in the central part of the study area is variable.For example, the Baini-Shawan fault (F23) is dextral, while the Xijiang fault (F24) to its west changes to sinistral.The extremely low fault slip rates are most likely responsible for the different fault slip senses.

Discussion
Faults and fractures play an important role in the circulation of geothermal fluids in the crust (Egger et al., 2014).They can not only channel fluids from hot deep depths of the crust to shallow depths to form extractable geothermal resources, but also facilitate the downward flux of meteoric fluids to deeper crust and the sustaining of broad convective geothermal systems (Jolie et al., 2021).However, such fluid movements are strongly affected by fault permeability, which varies according to the current stress field, structural setting, fault slip rate, etc. (Cox, 2010;Egger et al., 2014;Tanikawa et al., 2010).As mentioned in the Introduction, the stress state of a given structural setting is closely related to fault kinematics.Estimation of fault kinematics therefore is essential to understand the stress regime and, ultimately, the fault permeability, which helps to locate hidden geothermal systems and provides insight into their genesis.In the following sections, we will first compare the results with model-independent observations, such as reported fault slip rates and fault slip senses, to validate the model and then discuss in detail several research findings revealed from the kinematic characteristics of the study area.

Comparison with reported fault slip rates, earthquakes and slip senses
In order to validate our modeled results, published modelindependent observational data in the study area have been extensively collected (such as reported fault slip rates) and compared to our modeling results.In recent years, there have been a few investigations of fault activities conducted in the GBA.For example, Lei et al. (2018) used a dynamic deformation instrument to monitor the ENE-trending Wuhua-Shenzhen fault (F10) for four years and found that the fault slip rate is 0.01-0.06mm/a.Sun et al. (2007) analyzed the geological characteristics of the Wuhua-Shenzhen fault (F10) and suggested that the fault's activity is very low.Based on shallow seismic explorations combined with trench profiles, Wang et al. (2011) argued that the NE-trending Guangzhou-Conghua fault (F25) is almost inactive since Late Quaternary times.In summary, it is widely accepted that the slip rates on the faults in the GBA are generally very low, which is consistent with our modeled results.
To explore the relation between fault slip rate and the earthquake locations in the study area, we also plot earthquakes with M ≥ 3 together with fault slip rates, as shown in Fig. S3.As can be seen from the figure, earthquakes in the study area are generally sparsely distributed and there are no earthquake clusters in fault segments with relatively high slip rates even when the magnitude threshold is decreased to a much lower value (e.g., Xia et al., 2022).In fact, most of the earthquakes in the study area are located in the intersection zone of faults with NW and NE orientations (Sun et al., 2012).For example, the 1969 M6.4 Yangjiang earthquake occurred in the intersection area of the NE-trending Pinggang fault (F29) and the NW-trending Yangbianhai fault (F30) (Fig. S3b); the 1962 M6.1 Xinfengjiang earthquake also occurred in the intersection zone of the NNW-trending fault (namely the Shijiao-Xingang fault) and the ENE-trending Heyuan fault (F13) (Fig. S3c).Detailed studies on the seismogenic structures of small earthquakes located in the Pearl River Estuary also indicate that these earthquakes are in the intersection area of NW-trending and NE-trending faults (e.g.Chen et al., 2021;Xia et al., 2018).We consider this is because our study area is in an intraplate setting of weak deformation rates, within which the fault intersection is a typical stress concentrator place, where tectonic movement can cause a localized build-up of stresses and, ultimately, earthquakes (Gangopadhyay and Talwani, 2003).
For the sense of fault slip, Fan et al. (2022) used remote sensing interpretation as well as field investigation to estimate the kinematic characteristics of the branch of the Zijin-Boluo fault (F11) and found that this ENE-trending fault has an obvious dextral slip component (Fig. S1a); Huang and Zheng (2001) analyzed a series of microstructures in the NE-trending Wuchuan-Sihui fault zone (F32) and suggested that this fault is sinistral.These senses of fault slip revealed by field investigations match well with our model results (Fig. 7).
Although not all the faults implemented in the model can be validated against model-independent data due to the limited number of published fault activity data in the study area, the good consistency of fault slip rates and fault slip senses between the model-independent measurements and the modeled results on the typical faults indicate that our modeled kinematic results are generally reliable.

Mechanism of the opposite slip sense on ENE-and NE-trending faults
Two groups of faults oriented in NE-SW and ENE-WSW directions are developed in the west and east of the Pearl River Estuary, respectively (Fig. 1).Our modeling results show that the sense of slip on these two groups of faults is opposite to each other: the NE-trending faults are sinistral, while the ENE-trending faults show dextral movements (Fig. 7).Yu et al. (2016) first pointed out the difference in activity between these two groups of faults but gave no further explanation.
In order to gain a detailed understanding of the relationship between fault slip sense and fault strike, we design a series of sub-models with generic faults oriented in different directions.A submodeling technique is applied in the computation, during which the boundary nodes of the sub-models are driven by the modeled displacement of the GBA model.The generic faults in the sub-models are simplified to be vertical planes, consistent with the large dip angles of the faults in the study area (Table 1).The fault orientations range from 0 • to 170 • in steps of 10 • .Since the strike of a vertical fault is bidirectional, the chosen strike range can represent the full range of 0-360 • .The sub-models are cylindrical in shape with a diameter of 100 km and a height of 10 km, representing the thickness of the upper crust.Since the focus of this test is on the fault slip sense, the friction coefficient of the scenario faults in the sub-models is set to 0 to make the faults slide more easily.We then use our kinematic model to drive this series of sub-models and obtain the slip senses for the faults in different orientations.In view of the different dominant orientations of the faults on both sides of the Pearl River Estuary, this test is performed on the east and west sides of the Pearl River Estuary separately (as outlined by the two dashed circles in Fig. 7).The results are summarized in Fig. 8.
In the west of the Pearl River Estuary, faults with strikes between 110 • and 170 • exhibit dextral slip, while faults with strikes between 10 • and 80 • exhibit sinistral slip (Fig. 8a).Compared with the results in Fig. 8a, the pattern of fault slip senses in the east of the Pearl River Estuary (Fig. 8b) has an anticlockwise rotation, since the strike of dextral and sinistral faults is concentrating on the range of 50 • -120 • and 140 • -210 • , respectively (Fig. 8b).This anticlockwise rotation of the fault slip sense pattern reflects the lateral variation of the crustal stress field.Wen et al. (2001) calculated the modern crustal stress field in the southern part of the South China Block by the finite element method and found that the compressional stress axis in the west of the Pearl River Estuary trends nearly N-S and rotates to NW-SE in the east, which is consistent with the results of our sub-model simulations.Therefore, we suggest that the opposite slip senses on ENE-and NE-trending faults are caused by the lateral variation of the crustal stress field.The changing compressional stress axis is oblique to the ENE-and NE-trending faults and promotes the ENE-trending faults to slip in the dextral direction while the NE-trending faults show sinistral movements, as illustrated in Fig. 8c.
It should be noted that the angle between the compressional stress axis and the ENE-and NE-trending faults is usually larger than 45 • (Fig. 8c).This means that these faults are not oriented in the preferred direction, which is typically ± 20 • -30 • (Fossen, 2010) relative to the maximum principal stress axis σ 1 .The NE-and ENE-trending faults are pre-existing faults inferred to have been initiated in Triassic times which have experienced a complex history of deformation, including alternating episodes of shortening and extension, during the northwestward subduction of the Paleo-Pacific Plate (Li et al., 2020;Qiu et al., 1991).Since the Late Cenozoic, the stress regime of the study area has changed into the NW-SE-trending compression and converted the NE-and ENEtrending faults into transpressional structures (Xie et al., 2016;Zhang et al., 1999).
The shortening component of the transpressional movements will reduce the fault permeability and make it difficult to circulate geothermal fluids.On the contrary, the NW-trending faults, especially those located in the east of the Pearl River Estuary, are parallel to σ 1 (see Fig. S2b), i.e. they are extensional and can be most conducive to channel X. Li et al. geothermal fluids (e.g.Jolie et al., 2021).Therefore, we suggest that the NW-trending faults as well as their intersection zones with other faults are favorable sites for the exploitation of geothermal systems.In addition, more attention should also be paid to local structures, such as stepovers or fault bends along the NE-and ENE-trending faults, where pullapart basins can be formed due to the strike-slip component of the deformation (e.g.Fan et al., 2022).Normal faulting stress regimes in these local areas would facilitate the flow of fluids and host a geothermal system.

Spatial correlation of hot springs, heat flow, and fault slip rate
The extensive magmatism and widespread granitoids (Zhou et al., 2006), the gradually thinning crust towards the sea (Huang et al., 2015), and the uplift of the Moho interface (Zhang et al., 2018b) make the GBA and its vicinity an area of high heat flow values (Jiang et al., 2019) and abundant hot springs (Wang, 2018).Fig. 9 shows the spatial distribution of heat flow, hot springs, and fault slip rates in the GBA region.There is a significant high heat flow anomaly covering most of the western part of the study area with heat flow values reaching up to 95 mW/m 2 .Another high heat flow anomaly zone is located in the eastern part of the study area and extends in the NE-SW direction, nearly parallel to the dominant strike of the faults in this region.In addition, there is also a local high thermal anomaly (with heat flow as large as 80 mW/m 2 ) situated at the southwestern end of this thermal anomaly zone, around the Shenzhen-Hong Kong area.
Hot springs are often taken as indicators of geothermal anomalies.When underground fluids are heated by deep heat sources (such as recently intruded dikes and plutons), the fluid density decreases, and the buoyant forces promote the upward movement of fluids along faults to the surface, thus forming hot springs.Therefore, the presence of hot springs indicates the presence of heat sources below as well as the existence of fluid pathways and high permeability of the faults.There are at least 106 hot springs in the study area (Wang, 2018) primarily distributed in the southwestern coastal region and the eastern part of the study area, where there are also high heat flow anomaly zones (Fig. 9).However, the formation of hot springs is also influenced by specific hydrogeological conditions.For example, the presence of a lower water table in the geothermal system may prevent the formation of hot springs (Jolie et al., 2015).Therefore, areas without hot springs can also be strong geothermal anomaly zones (Zhong and Zhan, 1991).
In general, faults can increase permeability (e.g.Evans et al., 1997) and promote the generation of new fractures or make pre-existing ones to remain open, favorable for the formation of hydrothermal fields (Curewitz and Karson, 1997).As aforementioned, the granitoids are distributed widely in the study area.Laboratory experiments suggest that the permeability of granite increases significantly at high slip rates, due to the microcracks and mesoscale fractures forming near the slip surface (Tanikawa et al., 2010).Jolie et al. (2019) measured soil gas emission data in a geothermal field in the East African Rift System and found that the increased values of soil gas efflux are mainly located in the area with the maximum displacement of faults, indicating that at the regional scale in the field, higher fault slip rates can enhance the permeability of the fault and facilitate the migration of geothermal fluids.Furthermore, hydrothermal alteration processes in the fault often make the fault tight by minerals precipitated from hot fluids in the crust, such as the giant quartz reef 40 km long and > 75 m wide formed in the core of the Heyuan fault (F13) (Tannock et al., 2020).The activity of the fault would break the solidified fault cores and allow the fluid to flow.
Geothermal projects installed worldwide are located in various geological settings and commonly require highly permeable faults (Moeck, 2014), not only because these faults can facilitate the development of a circulation system (e.g. by providing pathways for hot fluids to ascend from greater depth), but also can provide high flow rates for commercial usage.In Germany, for example, fault zones are the primary targets of all geothermal projects, no matter whether the geothermal reservoir is of aquifer type in the foreland basin (Munich area in the Molasse) (Cacace et al., 2013) or of crystalline basement type in the Upper Rhine graben (Frey et al., 2022).
Our modeled results show that the highest slip rate of approximately 0.12 mm/a is observed along the northeastern segment of the Wuchuan-Sihui fault (F32).The Kaiping and Enping faults (F26, F28) as well as the middle segments of the Wuhua-Shenzhen fault (F10), Zijin-Boluo fault (F11), and Heyuan fault (F13) in the eastern study area also have a relatively high slip rate ranging from 0.03 to 0.06 mm/a.These active fault segments are believed to have higher fracture densities due to their higher slip rates, which can also form a pronounced normal faulting stress regime in releasing fault bends or step-overs (e.g.Fan et al., 2022, Fig. S1).Therefore, the permeability of these segments would be high.Considering the high background heat flow values, we infer that these segments could be interesting prospects for future geothermal exploration.Additionally, in the Shenzhen-Hong Kong area, the NW-SE compression of the background stress regime will increase the dilation tendency of the NW-trending fault (F16) (see Fig. 8c), which extensionally intersects the ENE-trending fault (F10) and further increase the permeability there.The expected elevated permeability and high background heat flow values in the intersection area suggest that geothermal exploration in that area could also be prospective.

Comparison with existing geothermal systems
Detailed geothermal investigations found that many hot springs or hot water wells are located along the middle segment of the Heyuan fault (F13) (Tannock et al., 2020), of which the fault slip rate is relatively high.The Huangshadong geothermal field, where a borehole of 3000 m produces geothermal water of 118 • C at the rate of 137 m 3 /h (Fan et al., 2022) is also located near the middle segment of the Zijin-Boluo fault (F11) with a relatively high slip rate.In the Tangkeng geothermal field, where the first Chinese geothermal power station was established, hot springs and geothermal wells are mainly distributed along the NW-trending faults that cross-cut the NE-trending faults (Luo et al., 2022), which means that the intersection of the faults is favorable to host promising geothermal resources.
Fault stress can also provide valuable clue for the potential geothermal resources.We extracted the stress from the 3D geomechanical model and calculated the slip and dilation tendency on the faults (Fig. S4), which is similar to the distribution pattern of the slip rate, indicating that the delineated potential geothermal resources based on stress analysis are consistent with those based on slip rate analysis.
Some limitations on the model primarily stem from uncertainties in the model-input parameters due to sparse subsurface or surface information.For example, the fault geometries employed in this study are highly simplified and some fault branches are omitted in the model.These simplifications are appropriate to obtain the first-order characteristic of slip rates on these faults, but more realistic and complex fault geometries and fault distributions would certainly improve the applicability and reliability of the model.Moreover, the absence of data on GPS measurements at the southern and southeastern model boundaries could potentially affect the reliability of the fault kinematics within the South China Sea.More high-quality records of the movement rate of the South China Sea are necessary to increase the reliability of the modeled kinematics of the faults in the sea area.

Conclusions
In this study, we establish a 3D geomechanical model of the Guangdong-Hong Kong-Macao Greater Bay Area.The model considers a complex 3D fault system, inhomogeneous rock properties, topography, initial stress as well as gravity and provides the spatially continuous contemporary kinematics of the faults and the crust for the study area.
The average velocity of the ESE-directed movement of the crust in the GBA is ~ 7.2 mm/a.The northern part of the study area has a higher velocity of ~ 7.5 mm/a, gradually decreasing to ~ 7.0 mm/a towards the south, indicating a low velocity gradient.Slip rates on the faults are generally low.The NE-trending faults slip at the highest rates of 0.05-0.12mm/a in the study area, which is generally higher than the rates on the ENE-trending faults (<0.04 mm/a).The NW-trending faults have the lowest modeled fault slip rates (<0.02 mm/a).
For the fault slip sense our modeling results show that the NEtrending faults located in the western part of the study area display sinistral motion, while the ENE-trending faults in the eastern part always undergo dextral slip.It is inferred that the opposite slip senses on these two fault groups are caused by changes in the orientation of the regional compressional stress axis.
From the perspective of geothermal resource exploration, our results indicate that the fault segments that have relatively high slip rates, e.g., the northeastern segment of the Wuchuan-Sihui fault (F32), the Kaiping and Enping faults (F26, F28) as well as the middle segments of the Wuhua-Shenzhen fault (F10), Zijin-Boluo fault (F11) and Heyuan fault (F13), are in areas with relatively high background heat flow.The apparent strike-slip component of the deformation will promote the local structures, such as releasing step-overs or fault bends along the NEand ENE-trending faults to form pull-apart basins.Normal faulting stress regimes and high fracture densities inferred to develop in these local areas would facilitate the flow of fluids and host a geothermal system.Additionally, the intersections of the NW-trending extensional faults and the NE-/ENE-trending faults could also be potential areas of interest for future geothermal exploration.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 3 .
Fig. 3. Geometry of the 3D geomechanical numerical model and the implemented 3D fault system.The white arrows and circles indicate the displacement boundary condition applied to the model.

Fig. 4 .
Fig.4.Variability of the misfit between the modeled crustal velocities and GPS measurements with the effective friction coefficient.For details of the method to calculate the misfit, it is referred toHergert and Heidbach (2010).The value of 0.02 is selected as the best friction coefficient.Note that GPS measurements within ~ 5 km of faults are excluded in the misfit calculation due to the strong perturbations in the proximity to faults.

Fig. 6 .
Fig. 6.Displacement boundary conditions applied to the model and comparison between the GPS observations and modeled horizontal velocities.The black arrows represent the horizontal boundary movements.The red and blue arrows represent the modeled velocities at the surface and the GPS observations (relative to stable Eurasia), respectively.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. Modeled horizontal fault slip rates (color-coded) and fault slip senses (black arrows).The background contours and the corresponding labels denote the magnitude of the crustal surface velocities (in mm/a) relative to stable Eurasia.The dashed circles represent the location of Fig. 8.

Fig. 8 .
Fig. 8.The preferred slip senses for different orientations of scenario faults in (a) the western and (b) the eastern part of the study area (detailed locations are outlined in Fig. 7).(c) Scheme interpreting the opposite slip senses on the NEand ENE-trending faults, which are inferred to be caused by the lateral variations of the crustal stress.

Fig. 9 .
Fig. 9. Spatial distribution of heat flow, hot springs and the modeled fault slip rates in the study area.Background color contours denote the heat flow, which is interpolated from the discrete data of IHFC and Tang et al. (2014) (see Fig. 1).White dashed ellipses represent the sites of prospective geothermal systems inferred by this study.

Table 1
Geometry of the faults implemented in the model.

Table 2
Elastic parameters and densities in the 3D geomechanical model.