Changes in Fault Slip Potential Due to Water Injection in the Rongcheng Deep Geothermal Reservoir, Xiong’an New Area, North China

: The Xiong’an New Area is abundant in geothermal resources due to its unique geological structure. To address whether large-scale deep geothermal exploitation will induce a fault slip, we ﬁrst determined the initial in situ stress ﬁeld using shallow (~4000 m) in situ stress measurements from the North China plain. After characterizing the in situ stress ﬁeld, we analyzed the initial stability of the main active faults in the sedimentary strata of the Rongcheng deep geothermal reservoir based on the Mohr–Coulomb failure criteria. Assuming that this area will be subjected to forty years of continuous ﬂuid injection, we calculated the excess pore pressure in the deep geothermal reservoir and, subsequently, estimated the fault slip potential of the main active faults in this region from 2021 to 2060. Our results indicated that both the in situ stress ﬁeld in the shallow crust of the Xiong’an New Area and the Middle-Late Pleistocene active faults will initially maintain a stable state. With constant ﬂuid injection for forty years at six geothermal wells in the Rongcheng deep geothermal reservoir, the maximum superposed excess pore pressure at a single well is 18 MPa; this excess pore pressure value impacts the stress state of faults within 8 km of the well location. These pore pressure perturbations heavily impact the F5-10, F5-11, and F9-2 segments of the Rongcheng uplift boundary fault, with FSP values of 92%, 23%, and 47% in 2060, respectively. Porosity exacts little impact on the fault slip potential on the boundary fault segments of F5-10 and F9-2 in the Rongcheng deep geothermal reservoir, while an enhanced permeability can weaken the FSP values for these faults. of forty years can be up to M w 5.0 with a 5% ﬂuid loss in the Rongcheng deep geothermal reservoir. Long-term water injection may increase the ambient thermoelastic stress to the point where faults in a critical (or subcritical) stress state become unstable. The results can provide a reference for geothermal development in terms of injection rate and locations of geothermal wells.


Introduction
Previous work has established that wastewater proposals [1,2], oil and gas operations [3,4], Carbon Capture and Storage (CCS) [5,6], and geothermal energy development [6][7][8] are capable of inducing earthquakes [9][10][11]. Earthquakes occur when the crustal The location of Xiong'an New Area in China and regional division of the recent tectonic stress field in China (modified from [26]). (b) The geothermal fields distributed in the Jizhong Depression (modified from [22]). (c) The Middle-Late Pleistocene active faults, historical earthquakes (since 1970), and existing geothermal exploration boreholes in the Xiong'an New Area (modified from [21][22][23][24]   The location of Xiong'an New Area in China and regional division of the recent tectonic stress field in China (modified from [26]). (b) The geothermal fields distributed in the Jizhong Depression (modified from [22]). (c) The Middle-Late Pleistocene active faults, historical earthquakes (since 1970), and existing geothermal exploration boreholes in the Xiong'an New Area (modified from [21][22][23][24]

Geological Setting
Our study area, the Xiong'an New Area, is located in the Jizhong Depression of North China plain (Figure 1a) [13,14]. The destruction of the North China Craton began in the late Mesozoic, when it was destabilized by the subduction of the Pacific slab; the tectonic regime transitioned from a NS-oriented transpression to a nearly EW-oriented transtension [27]. As a result, extensional structures such as metamorphic core complexes, detachment faults, and fault-basins developed in the Early Cretaceous [28]. In the early Cenozoic, the Jizhong Depression entered the fault-depression stage and the Xiong'an New Area began experiencing normal faulting depression in the Early Paleogene [29]. Normal faults, such as the Xushui-Dacheng (F3), Rongcheng (F5 and F9), and Niudong faults (F2), and secondary structural units, such as the Xushui Depression (XS-D), the Rongcheng uplift (RC-U), and the Baxian Depression (BX-D), developed in the Paleocene-Eocene. While the Taihang Mountain fault became quiescent in the Middle Eocene, there is evidence of ongoing seismic activity on the Niudong fault (F2). In the Late Eocene-Oligocene, normal faulting depression continued in the Xushui Depression (XS-D) and the Baxian Depression (BX-D). The normal faulting depression tapered off between the Neogene to the Late Quaternary; since the Middle-Late Pleistocene, nearly all of the major faults in the area (with the exception of the Niudong fault (F2)) have become inactive. In general, the Middle-Late Pleistocene active faults have been inactive since the Holocene, and the geological structure of the Xiong'an New Area is relatively stable [22][23][24].
Historical seismic data indicate that there have been no earthquakes with magnitudes M > 3.0 in the Xiong'an New Area since 1970. While a moderate number of earthquakes with magnitudes M ≥ 3.0 have occurred within 50 km of Baiyangdian Lake (BYD Lake; Figure 1c), many of these events were located in the Renqiu and Wen'an counties; of these seismic events, the largest earthquake (M 5.1) occurred in Wen'an county on 4 July 2006 (Figure 1c) [30]. Most of these earthquakes occurred near the southeastern segments of Xushui-Dacheng Fault (F3), and some small earthquakes that occurred within Xiongxian (XX) may be related to the Niudong Fault (F2).
To analyze the initial faulting stability in the Rongcheng deep geothermal reservoir, we gathered data from the main active faults surrounding the Rongcheng deep geothermal reservoir. Due to the difference in fault strike, the main active faults were simplified to facilitate our analysis of the fault slip potential ( Figure 2). In our study, there were 5 active faults consisting of 25 fault segments. The fault strikes and dip angles had uncertainties of ±5 • and ±10 • , respectively. The parameters of faults are listed in Table 1 [19,21,23,24,[31][32][33][34].  [24,31]). (b) Geological profile of the Rongcheng deep geothermal reservoir (modified from [31]).   [24,31]). (b) Geological profile of the Rongcheng deep geothermal reservoir (modified from [31]).

Mohr-Coulomb Failure Criteria
The Mohr-Coulomb failure criteria is a useful framework for understanding how increasing the pore pressure via fluid injection can trigger a slip along faults. Due to the critically stressed nature of the crust, a fault will remain in a locked state as long as the applied shear stress is lower than the strength of the contact between the rocks on either side of the fault. The critical shear stress on an earthquake fault under static friction is expressed as [35,36]: where τ c is the shear stress (MPa), µ is the coefficient of friction, σ n is the normal stress (MPa), P f is the total pore pressure (MPa), P 0 is the initial pore pressure or hydrostatic pressure (MPa), and ∆P is the excess pore pressure caused by fluid injection (MPa). Under ambient conditions, the effective normal stress, which is oriented normal to the plane of the fault, effectively clamps the fault closed and reduces the likelihood of a slip occurring on the fault. During fluid injection, as the pore fluid pressure increases, the effective normal stress decreases proportionally; this reduction in the normal stress unclamps the fault and may result in a slip along pre-existing subcritical ruptures [25].

Hsieh and Bredehoeft Hydrological Model
Injection of fluids into a porous medium causes an increase in the pore pressure that decays exponentially with radial distance from the injection source. This pressure change radiates away from the well radially as injection continues; as such, our model includes a radially symmetric pressure profile for each injection well at a given point in time. This profile, known as the Hsieh and Bredehoeft hydrological model, is expressed as [37]: where ∆h is the vertically averaged buildup of hydraulic head above the initial head (m), ∆P is the excess pore pressure (MPa), T is the principal value of the transmissivity (m 2 /s), S is the storage coefficient of the reservoir, x and y represent the spatial location of the water injection wells (m), Q(t) is the injection rate (L/s), t is the injection time (s), and r is the specific weight of the fluid (N/m 3 ). These groundwater flow equations describe the two-dimensional (2D) radial flow in a vertically confined aquifer containing a well with a variable injection rate. In order to compute the injection-induced pressure buildup and subsequent falloff, our model relies on several simplifying assumptions [25,[37][38][39]: (1) the porous medium is fully saturated and has a uniform pressure distribution, (2) the hydraulic head is constant in all wells prior to the onset of the fluid injection, (3) injection wells are treated as point sources in the 2D grid, (4) the permeability and porosity are constant and isotropic, and (5) the total effect of pressure plumes interacting with one another is the linear sum of the individual pressure plume effects in that area.
As there are few fractures in the upper (~2580-2600 m) and bottom boundaries (~2820-2850 m) of the water injection layer (~2600-2820 m), the Hsieh and Bredehoeft hydrological model is suitable for calculating the pore pressure perturbations caused by water injection in the Rongcheng deep geothermal reservoir.

Initial In Situ Stress Field in the Shallow Crust of the Xiong'an New Area
Currently, there are no public in situ stress measurements of the Xiong'an New Area. However, Huang et al. (2013) characterized the averaged stress field in North China plain using 1017 points of in situ stress data from the "Fundamental database of crustal stress environment in continental China"; their results showed that the linearly increasing gradients of the two horizontal principal stresses (σ H and σ h ) are 0.0233 MPa/m and 0.0162 MPa/m at shallow depths (~4000 m), respectively [40,41]. To gather more specific data on our study area, we conducted 16 additional in situ stress measurements using hydraulic fracturing in Shunping county in 2018, which is located~70 km away from the Xiong'an New Area ( Figure S1), and the results of hydraulic fracturing in situ stress measurement are shown in Table S1 [42]. The linearly increasing gradients of σ H and σ h at Shunping county are 0.0252 MPa/m and 0.0164 MPa/m, respectively, which are slightly higher than the average in North China. To remedy the lack of in situ stress information in the Xiong'an New Area, we obtained the integrated increasing gradients of σ H and σ h in the Xiong'an New Area by refitting the 16 in situ stress measurements in Shunping county and the 1017 points of in situ stress data from North China. The integrated fitting results revealed that the maximum (σ H ) and minimum (σ h ) horizontal principal stresses have linearly increasing gradients of 0.0233 MPa/m and 0.0162 MPa/m, respectively ( Figure 3). We can learn that the linearly increasing gradients of the σ H and σ h after refitting do not change significantly, indicating a relatively stable stress level in the shallow crust of the North China region. The integrated stress field approximately represents the ambient stress field of Xiong'an New Area. The magnitude of the vertical principal stress σ V is approximately equal to the weight of the overburden. The density of the dolomite unit varies between 2.51 and 2.69 g/cm 3 with an average value of 2.60 g/cm 3 at the depth of 4000 m, and then the linearly increasing gradient of σ V is set to be 0.0260 ± 0.001 MPa/m [43]. Furthermore, Figure 3 suggests that below a depth of 1286.48 m, the stress regime is characterized by normal faulting; this stress state is consistent with the extensional tectonic dynamic history of the Jizhong Depression [27].  Observations of borehole breakouts in the existing boreholes of the central Jizhong Depression (at distances of 28-35 km from Rongcheng thermal storage area, with depths of 1000-4086 m) indicate that the σH orientation is N 77°-86° E [44]; moreover, the σH orientation obtained from the focal mechanism inversion of the M 5.1 Wen' an earthquake (at distances of ~42 km from Rongcheng thermal storage area) is N 68° E [45] (Figure 1c). These stress data show that the orientation of the maximum principal stress in Xiong'an New Area is ENE (N 68°-86° E, with an average value of N 77° E ± 9°), an orientation that is consistent with the regional stress field (N 82° E) [26].
The logging data from representative deep geothermal wells show that the relatively stable static water level is 27-65 m (e.g., Borehole D18, located in the Rongcheng deepgeothermal reservoir) [32]. The initial pore pressure (P0) is approximately equal to the static water pressure. With these data, we estimated that the initial pore pressure gradient with depth in the Xiong'an New Area is ~0.0094-0.0097 MPa/m, with an average value of 0.00955 ± 0.00015 MPa/m.

Friction Coefficients of the Main Active Faults
As shown in Equation (1), it is necessary to determine the friction coefficient (μ) of Observations of borehole breakouts in the existing boreholes of the central Jizhong Depression (at distances of 28-35 km from Rongcheng thermal storage area, with depths of 1000-4086 m) indicate that the σ H orientation is N 77 • -86 • E [44]; moreover, the σ H orientation obtained from the focal mechanism inversion of the M 5.1 Wen' an earthquake (at distances of~42 km from Rongcheng thermal storage area) is N 68 • E [45] (Figure 1c). These stress data show that the orientation of the maximum principal stress in Xiong'an New Area is ENE (N 68 • -86 • E, with an average value of N 77 • E ± 9 • ), an orientation that is consistent with the regional stress field (N 82 • E) [26].
The logging data from representative deep geothermal wells show that the relatively stable static water level is 27-65 m (e.g., Borehole D18, located in the Rongcheng deepgeothermal reservoir) [32]. The initial pore pressure (P 0 ) is approximately equal to the static water pressure. With these data, we estimated that the initial pore pressure gradient with depth in the Xiong'an New Area is~0.0094-0.0097 MPa/m, with an average value of 0.00955 ± 0.00015 MPa/m.

Friction Coefficients of the Main Active Faults
As shown in Equation (1), it is necessary to determine the friction coefficient (µ) of the faults. The friction coefficient is influenced by many factors such as stress, temperature, and the fault material [46]. In our study, only the static friction coefficient was considered. The results from numerous rock mechanical tests suggest that the friction coefficient of brittle rock in frictional equilibrium lies between 0.6 and 1.0 [47]. When it comes to assessing the fault slip potential, an empirical friction coefficient of 0.6 is typically invoked as the critical value [25,48,49]. For example, Walsh and Zoback (2016) evaluated the injection-induced faulting instability and seismicity in northern and central Oklahoma using a frictional coefficient of 0.6 [25]. Zhai and Shirzaei (2018) also used a frictional coefficient of 0.6 to explore the relationship between the high-volume deep fluid injection and the increasing seismicity in the Barnett Shale in Texas [49]. As such, we employed a critical friction coefficient of 0.6 in our evaluation of the initial fault slip potential in our study area.

Initial Stability of the Main Active Faults
We utilized the FSP v.1.0 software package to estimate the fault slip potential of the main active faults in the sedimentary strata of the Rongcheng deep geothermal reservoir. This software allows users to generate either a deterministic or a probabilistic geomechanical model of the fault slip potential. Both deterministic and probabilistic geomechanical models rely on three simplifying assumptions [25,[37][38][39]: (1) both the initial pore pressure and the in situ stress tensor are uniform across the study area and linearly increase in magnitude with depth; (2) one of the principal stress vectors is oriented vertically; and (3) the stress state is determined by the relative magnitude of the vertical stress vector (maximum, intermediate, or minimum). With the FSP software, we estimated the likelihood that the planar fault segments in question will be critically stressed within a given stress field.

Deterministic Geomechanical Assessment of Fault Slip Potential
The deterministic geomechanical assessment of the slip potential of the main active faults in the Rongcheng deep geothermal reservoir in the absence of fluid injection was first calculated using the parameters below.
The increasing gradients of σ H , σ h , σ V , and P 0 at a reference depth of 2600 m (i.e., the intermediate depth of Gaoyuzhuang Formation seen at D16 borehole) are 0.0233 MPa/m, 0.0162 MPa/m, 0.0260 ± 0.001 MPa/m, and 0.00955 ± 0.00015 MPa/m, respectively. As the increasing gradients of σ H and σ h were obtained using the data of the averaged stress level in North China, which can characterize the horizontal principal stress distributions of the shallow crust in the study area, we only considered the uncertainties of σ V controlled by gravitational loading. As mentioned earlier, we used a critical friction coefficient (µ) value of 0.6. The orientation of σ H is N 77 • E.
Based on initial stress boundary conditions confined by the integrated stress field of Xiong'an New Area, we calculated the effective normal stress and shear stress along the main active faults, and assessed the initial stable state based on Mohr-Coulomb failure criteria (Equation (1)). Figure 4 shows the results of our deterministic geomechanical assessment of the pore pressure required to generate a fault slip at a given fault segment. We found that the main active faults in the Rongcheng deep geothermal reservoir are relatively stable and that it is very unlikely that they will slip in the present stress field (Figure 4a). The pore pressure required to produce a fault slip varies for each fault segment (Figure 4b). For example, for segments F5-1-F5-7 and F5-10-F5-14, the pore pressure required to cause a fault slip ranges from 5.58 MPa to 32.32 MPa, while fault segments F5-5, F5-6, and F5-10 are relatively close to the critical stress state. For segments F9-1-F9-5, the pore pressure required to cause a slip ranges from 6.34 MPa to 20.24 MPa, and 15.31 MPa, where the segments F9-4 and F9-5 are close to the critical stress state.   Figure 4a also shows that the minimum critical friction coefficient of the reservoir rock mass that would result in a slip is 0.48. A lower critical friction coefficient of fault can reduce its frictional resistance to sliding and accelerate the process of faulting instability (Equation (1), Figure 4a). To explore the worst-case scenario, namely, to assess the maximal influence of fluid injection on the stable state of faults, we used 0.48, instead of 0.6, as the critical friction coefficient of the main active faults in subsequent calculations.

Probabilistic Analysis of Fault Slip Potential
As the deterministic model ignores some uncertainties that are often present in the strike, dip, in situ stress field, and the coefficient of friction, the deterministic geomechanical results are not entirely reliable. To minimize these uncertainties, we evaluated the possibility of a slip along these faults in response to an increase in pore pressure using a Monte Carlo method to randomly sample the specified uniform uncertainty distributions of the input parameters [25,50]. The Monte Carlo method is an extremely useful method of assessing the model error when there are uncertainties in the model parameters and little to no information on historical fault slips in our study area [25]. For example, the resampling of the probabilistic distribution of the model parameters for segment F5-10 is shown in Figure 5. Figure 4a also shows that the minimum critical friction coefficient of the reservoir rock mass that would result in a slip is 0.48. A lower critical friction coefficient of fault can reduce its frictional resistance to sliding and accelerate the process of faulting instability (Equation (1), Figure 4a). To explore the worst-case scenario, namely, to assess the maximal influence of fluid injection on the stable state of faults, we used 0.48, instead of 0.6, as the critical friction coefficient of the main active faults in subsequent calculations.

Probabilistic Analysis of Fault Slip Potential
As the deterministic model ignores some uncertainties that are often present in the strike, dip, in situ stress field, and the coefficient of friction, the deterministic geomechanical results are not entirely reliable. To minimize these uncertainties, we evaluated the possibility of a slip along these faults in response to an increase in pore pressure using a Monte Carlo method to randomly sample the specified uniform uncertainty distributions of the input parameters [25,50]. The Monte Carlo method is an extremely useful method of assessing the model error when there are uncertainties in the model parameters and little to no information on historical fault slips in our study area [25]. For example, the resampling of the probabilistic distribution of the model parameters for segment F5-10 is shown in Figure 5. Using the data discussed in Section 4.3.1, we produced a probabilistic fault slip analysis in the absence of fluid injection for 2021 ( Figure 6). Many results have a lower initial fault slip potential that are less than 5%, while the fault slip potentials of segments F9-4, F9-5, and F5-6 are slightly higher, being 6%, 11%, and 9%, respectively. Overall, we found that, in the initial stress field, the main active faults in the Rongcheng deep geothermal reservoir are in a stable stress state. Using the data discussed in Section 4.3.1, we produced a probabilistic fault slip analysis in the absence of fluid injection for 2021 ( Figure 6). Many results have a lower initial fault slip potential that are less than 5%, while the fault slip potentials of segments F9-4, F9-5, and F5-6 are slightly higher, being 6%, 11%, and 9%, respectively. Overall, we found that, in the initial stress field, the main active faults in the Rongcheng deep geothermal reservoir are in a stable stress state.

Injection-Induced Changes in the Pore Pressure in the Rongcheng Deep Geothermal Reservoir
In general, the risk of faulting instability increases as the distance decreases between a fluid injection well and a fault [25,51]. To calculate the excess pore pressure that occurs over a forty

Injection-Induced Changes in the Pore Pressure in the Rongcheng Deep Geothermal Reservoir
In general, the risk of faulting instability increases as the distance decreases b a fluid injection well and a fault [25,51]. To calculate the excess pore pressure tha over a forty-year time scale (2021-2060), we selected six representative deep geo wells (D11, D12, D15, D18, D20, and D21) that are within 1 km of the Rongchen boundary fault.
According to the Hsieh and Bredehoeft hydrological model (Equations (3) a the physical and mechanical properties (e.g., reservoir thickness, rock density, pe ity, porosity) and the geothermal well parameters (e.g., the location of wells, the i depth, the injection rate, the injection time) directly affect the distribution of exc pressure caused in the presence of fluid injection. Previous studies on geotherm D16 indicated that the density of the reservoir rock is 2675 kg/m 3 [43], the aquifer th of the karst fissure geothermal reservoir is ~220 m [34], the porosity of the deep mal reservoir is 1.34-4.08% [21], the reservoir permeability is 1.33-2.92 mD [21], According to the Hsieh and Bredehoeft hydrological model (Equations (3) and (4)), the physical and mechanical properties (e.g., reservoir thickness, rock density, permeability, porosity) and the geothermal well parameters (e.g., the location of wells, the injection depth, the injection rate, the injection time) directly affect the distribution of excess pore pressure caused in the presence of fluid injection. Previous studies on geothermal well D16 indicated that the density of the reservoir rock is 2675 kg/m 3 [43], the aquifer thickness of the karst fissure geothermal reservoir is~220 m [34], the porosity of the deep geothermal reservoir is 1.34-4.08% [21], the reservoir permeability is 1.33-2.92 mD [21], the density of water is 1000 kg/m 3 , the coefficient of water viscosity is 0.0008 Pa·s (the default value in the FSP software), and the injection rate (i.e., the maximum pumping rate) is 171 m 3 /h [33]. Herein, we selected the maximum values of porosity and permeability to calculate the possible excess pore pressure caused by fluid injection. Figure 7 shows the projected pore pressure perturbations generated by fluid injection into the six injection wells from 2021 to 2060.

Injection-Induced Changes in the Fault Slip Potential in the Rongcheng Deep Geothermal Reservoir
Based on the initial stability of the main active faults (Figure 4), as well as the projected pore pressure changes caused by forty years of fluid injection (Figure 7), we used the FSP v.1.0 software package to calculate the changes in the fault slip potential of the main active faults between 2021 and 2060 ( Figure 9).

Injection-Induced Changes in the Fault Slip Potential in the Rongcheng Deep Geothermal Reservoir
Based on the initial stability of the main active faults (Figure 4), as well as the projected pore pressure changes caused by forty years of fluid injection (Figure 7), we used the FSP v.1.0 software package to calculate the changes in the fault slip potential of the main active faults between 2021 and 2060 ( Figure 9).

Injection-Induced Changes in the Fault Slip Potential in the Rongcheng Deep Geothermal Reservoir
Based on the initial stability of the main active faults (Figure 4), as well as the projected pore pressure changes caused by forty years of fluid injection (Figure 7), we used the FSP v.1.0 software package to calculate the changes in the fault slip potential of the main active faults between 2021 and 2060 ( Figure 9). Figure 9 indicates that forty years of fluid injection at wells D11, D12, D15, D18, D20, and D21 increases the fault slip potential at segments F5-10, F5-11, and F9-2, while the fault slip potential at other fault segments is largely unchanged. As shown in Figure 10, the fault slip potential increases exponentially with injection time. For example, due to injection at wells D12, D18, D20, and D21, the fault slip potential of segment F5-10 increases from 1% in 2021 to 92% in 2060. The fault slip potentials of segments F5-11 and F9-2 increase from little risk in 2021 to higher risk with the FSP values of 23% and 41% in 2060, respectively.  Figure 9 indicates that forty years of fluid injection at wells D11, D12, D15, D18, D20, and D21 increases the fault slip potential at segments F5-10, F5-11, and F9-2, while the fault slip potential at other fault segments is largely unchanged. As shown in Figure 10, the fault slip potential increases exponentially with injection time. For example, due to injection at wells D12, D18, D20, and D21, the fault slip potential of segment F5-10 increases from 1% in 2021 to 92% in 2060. The fault slip potentials of segments F5-11 and F9-2 increase from little risk in 2021 to higher risk with the FSP values of 23% and 41% in 2060, respectively.
We compared the fault parameters of F5-10 with F9-2 (Table 1), and we found that both fault segments have the same friction coefficient and a similar strike and well distance (Figures 1c and 2a, Table 1). However, the dip of F5-10 (70°) is greater than that of F9-2 (53°). The higher FSP value (92%) for F5-10 in 2060 may indicate that some steeply dipping faults are more prone to be reactivated with the same in situ stress and pore pressure perturbations caused by fluid injection. For example, recent widespread seismicity in Oklahoma (USA) is attributed to the reactivation of pre-existing basement structures facilitated by steeply dipping basement-rooted faults [52,53].

Effects of Porosity and Permeability on the FSP Values of Segments F5-10 and F9-2
The main water yielding stratum logging results of the Gaoyuzhuang Formation i We compared the fault parameters of F5-10 with F9-2 (Table 1), and we found that both fault segments have the same friction coefficient and a similar strike and well distance (Figures 1c and 2a, Table 1). However, the dip of F5-10 (70 • ) is greater than that of F9-2 (53 • ). The higher FSP value (92%) for F5-10 in 2060 may indicate that some steeply dipping faults are more prone to be reactivated with the same in situ stress and pore pressure perturbations caused by fluid injection. For example, recent widespread seismicity in Oklahoma (USA) is attributed to the reactivation of pre-existing basement structures facilitated by steeply dipping basement-rooted faults [52,53].

Effects of Porosity and Permeability on the FSP Values of Segments F5-10 and F9-2
The main water yielding stratum logging results of the Gaoyuzhuang Formation in D16 (at depths of 2600-2820 m) showed that the porosity and permeability range between 1.34-4.08% and 1.33-2.92 mD, respectively [21]. The influences of the porosity and permeability were not considered in the previous analysis. Thus, we recalculated the probabilistic FSP on segments F5-10 and F9-2 in 2030, 2040, 2050, and 2060, respectively, with porosity and permeability values ranging from 1.34 to 4.08% in 0.865% increments, and from 1.33 to 2.92 mD in 0.53 mD increments, respectively.
With a constant permeability (2.92 mD), we first calculated the effect of porosity on the FSP values. The results showed that the FSP values of segments F5-10 and F9-2 in the same year do not have obvious changes (∆FSP < 8%), with increasing porosity (Figure 11).

Effects of Porosity and Permeability on the FSP Values of Segments F5-10 and F9-2
The main water yielding stratum logging results of the Gaoyuzhuang Formation in D16 (at depths of 2600-2820 m) showed that the porosity and permeability range between 1.34-4.08% and 1.33-2.92 mD, respectively [21]. The influences of the porosity and permeability were not considered in the previous analysis. Thus, we recalculated the probabilistic FSP on segments F5-10 and F9-2 in 2030, 2040, 2050, and 2060, respectively, with porosity and permeability values ranging from 1.34 to 4.08% in 0.865% increments, and from 1.33 to 2.92 mD in 0.53 mD increments, respectively.
With a constant permeability (2.92 mD), we first calculated the effect of porosity on the FSP values. The results showed that the FSP values of segments F5-10 and F9-2 in the same year do not have obvious changes (∆FSP < 8%), with increasing porosity (Figure 11).  With a constant porosity (4.08%), secondly, we calculated the effect of permeability on the FSP values. The results showed that an enhanced permeability could weaken the FSP values for segments F5-10 and F9-2 in 2030, 2040, 2050, and 2060, respectively. The changes in permeability have largely influenced the FSP values for F9-2 after 2040. For example, in 2060, the FSP values will decrease by 39% with the permeability increased from 1.33 mD to 2.92 mD (Figure 12b). Cappa et al. (2018) suggested that permeability enhancement has an important effect on the pressure diffusion and seismic slip growth during fluid injection [54]. Their results revealed that a more pronounced permeability enhancement results in a larger seismic slip zone. As such, the permeability should be considered when conducting the seismic hazard assessment of a given region due to fluid injection.

The Predicted Maximum Moment Magnitude of Injection-Induced Seismicity
Previous studies have suggested that faults with a higher fault slip potential are more prone to induce earthquakes [1,4,25]. Some earthquakes induced by fluid injection in the Rongcheng uplifts may occur near the segments F5-10 and F9-2 with higher FSP values (92% and 47%, in 2060, respectively). McGarr (2014) predicted the possible maximum magnitude of injection-induced earthquakes by simulating a fully saturated reservoir with critically stressed and ideally oriented faults in the vicinity of an injection well [55]. The model generates a linear relationship between the maximum magnitude and the net injected volume (∆V): where M 0 (max) is the maximum seismic moment, G is the modulus of rigidity, and ∆V is the net injected volume.
Water 2022, 14, x FOR PEER REVIEW 16 of 21 enhancement has an important effect on the pressure diffusion and seismic slip growth during fluid injection [54]. Their results revealed that a more pronounced permeability enhancement results in a larger seismic slip zone. As such, the permeability should be considered when conducting the seismic hazard assessment of a given region due to fluid injection.

The Predicted Maximum Moment Magnitude of Injection-Induced Seismicity
Previous studies have suggested that faults with a higher fault slip potential are more prone to induce earthquakes [1,4,25]. Some earthquakes induced by fluid injection in the Rongcheng uplifts may occur near the segments F5-10 and F9-2 with higher FSP values (92% and 47%, in 2060, respectively). McGarr (2014) predicted the possible maximum magnitude of injection-induced earthquakes by simulating a fully saturated reservoir with critically stressed and ideally oriented faults in the vicinity of an injection well [55]. The model generates a linear relationship between the maximum magnitude and the net injected volume (ΔV): where M0(max) is the maximum seismic moment, G is the modulus of rigidity, and ∆V is the net injected volume. For the parameters used in Equation (5), according to the results of experiments examining the triaxial compressive strength of carbonate rock from the Gaoyuzhuang Formation (at the depths of 4195-4422 m) in Xiong'an New Area, the maximum Young's modulus E of the dolomite is equal to 65.03 GPa [43]. Then, the modulus of rigidity G (G = E/2 (1 + v)) is equal to 54.19 GPa.
In the calculation, we explored the possible earthquake magnitudes induced by fluid injection using the maximum injection rate of 171 m 3 /h for forty years of fluid injection (from 2021 to 2060), which corresponds to a monthly injection volume of 1.23 × 10 5 m 3 . We also considered the net injected volume (∆V) as 5% V [17], due to fluid loss injected into the reservoir. In addition, we examined the injection-production time in Xiong'an New Area, and the production logs suggested that the productions and injections usually start from November of that year to March of the second year. Thus, the injection time can be considered as 5 months per two years.
The predicted maximum moment magnitudes are shown in Figure 13. For the parameters used in Equation (5), according to the results of experiments examining the triaxial compressive strength of carbonate rock from the Gaoyuzhuang Formation (at the depths of 4195-4422 m) in Xiong'an New Area, the maximum Young's modulus E of the dolomite is equal to 65.03 GPa [43]. Then, the modulus of rigidity G (G = E/2 (1 + v)) is equal to 54.19 GPa.
In the calculation, we explored the possible earthquake magnitudes induced by fluid injection using the maximum injection rate of 171 m 3 /h for forty years of fluid injection (from 2021 to 2060), which corresponds to a monthly injection volume of 1.23 × 10 5 m 3 . We also considered the net injected volume (∆V) as 5% V [17], due to fluid loss injected into the reservoir. In addition, we examined the injection-production time in Xiong'an New Area, and the production logs suggested that the productions and injections usually start from November of that year to March of the second year. Thus, the injection time can be considered as 5 months per two years.
The predicted maximum moment magnitudes are shown in Figure 13.  Figure 13 shows that the estimated maximum moment magnitude increases with the injection time, namely, injection volume. The maximum moment magnitude of earthquakes induced by continuous water injection for 10, 20, 30, and 40 years is up to Mw 4.6, Mw 4.8, Mw 4.9, and Mw 5.0, respectively, which are smaller than that of the largest natural earthquake with Mw 6.3 that may occur in the Baoding seismic area [56,57].  Figure 13 shows that the estimated maximum moment magnitude increases with the injection time, namely, injection volume. The maximum moment magnitude of earthquakes induced by continuous water injection for 10, 20, 30, and 40 years is up to M w 4.6, M w 4.8, M w 4.9, and M w 5.0, respectively, which are smaller than that of the largest natural earthquake with M w 6.3 that may occur in the Baoding seismic area [56,57].

Injection-Induced Changes in the Thermoelastic Stress and Faulting Instability
The cumulative stress changes in the horizontal principal stress (∆σ horizontal ) caused by the thermal contraction of a reservoir can result in large-magnitude seismic events [53]. Assuming an isotropic, porous, elastic, and laterally extensive reservoir, that is thin relative to its lateral extensiveness; then, the change in the horizontal principal stress as a result of changes in the thermoelasticity can be expressed as [58]: where ∆σ horizontal is the change in the minimum horizontal principal stress due to fluid injection (MPa), υ is Poisson's ratio, α d is the drained thermoelastic effective stress coefficient (MPa ( • C) −1 ), ∆T is the change in temperature throughout the reservoir ( • C), K is the intrinsic material constants, and β d is the volumetric expansion coefficient (( • C) −1 ). For the parameters used in Equations (5) and (6), we set Poisson's ratio to be 0.233, which was carried out from conventional triaxial compression and tensile experiments of carbonate rock in D34, being 25 km away from Southeastern Rongcheng (at depths of 4422 m) [43]; α d was set to 0.26 MPa ( • C) −1 [59]; and, based on a twenty-year study of fluid injection in Geysers, California, ∆T was set to −6 • C [60].
Once these parameters were defined, we used Equations (5) and (6) to calculate ∆σ horizontal within the reservoir. We determined that the overall change in the minimum horizontal principal stress is −1.08 MPa in the normal faulting-stress regime ( Figure 14). As shown in Figure 14, this change in ∆σ horizontal may result in Mohr-Coulomb failure and, subsequently, faulting instability at a number of critical (or subcritical) fault segments (e.g., segments F2-2, F5-5, and F9-4 in Figure 14b). The faulting stability of fault segments such as F2-4, F5-3, and F5-14 remain largely unaffected by this change in the thermoelastic stress.

Conclusions
This paper quantitatively describes how the stable state of the main active faults in Rongcheng deep geothermal reservoir would change in response to a long-term fluid injection in the future. Here are the main conclusions.
The linearly increasing gradients of σH, σh, and σV in the shallow crust of the Xiong'an New Area are 0.0233 MPa/m, 0.0162 MPa/m, and 0.0260 MPa/m, respectively; these gradients are representative of a normal-faulting stress regime at depths greater than ~1200

Conclusions
This paper quantitatively describes how the stable state of the main active faults in Rongcheng deep geothermal reservoir would change in response to a long-term fluid injection in the future. Here are the main conclusions.
The linearly increasing gradients of σ H , σ h , and σ V in the shallow crust of the Xiong'an New Area are 0.0233 MPa/m, 0.0162 MPa/m, and 0.0260 MPa/m, respectively; these gradients are representative of a normal-faulting stress regime at depths greater thañ 1200 m. The σ H orientation is ENE with an average of N 77 • E ± 9 • . In the present in situ stress field, the main active faults in the Rongcheng deep geothermal reservoir will not slip instantaneously.
With an injection rate of 171 m 3 /h, the maximum superposed excess pore pressure caused by continuous injection at a single deep geothermal well for forty years is 18 MPa; the magnitudes of these pore pressure perturbations decay exponentially with distance from the well and exert little influence on the surrounding area beyond a distance of 8 km from the well.
With an injection rate of 171 m 3 /h, forty years of fluid injection heavily impacts the stability of the southwestern Rongcheng uplift boundary fault (i.e., fault segments F5-10, F5-11, and F9-2), while the impact on the other Quaternary fault segments is minimal. For forty years of water injection (from 2021 to 2060), the FSP values of fault segments F5-10, F5-11, and F9-2 exponentially increase with injection time. By 2060, the fault slip potential of these three segments is 92%, 23%, and 47%, respectively. Under the same conditions of in situ stress, fault strike, and injection activity, the southwestern segment of Rongcheng Boundary fault (F5-10) with a greater dip (70 • ) has a higher fault slip potential value than that for F9-2 with a smaller dip (53 • ), indicating that some steeply dipping faults are more prone to be reactivated by fluid injection.
In the Rongcheng geothermal reservoir, the fault slip potential on main boundary faults (such as F9-2 and F5-10) induced by long-term fluid injection changes slightly in response to variations in porosity. However, the fault slip potential values on these faults decrease with the unit permeability during the fluid injection.
With an injection rate of 171 m 3 /h, the magnitudes of injection-induced earthquakes increase with the injection volume, and the maximum moment magnitude can be up to M w 5.0 for continuous fluid injection from 2021 to 2060, considering a 5% fluid loss. The predicted maximum magnitude of injection-induced seismicity would be smaller than that of the largest natural earthquake with a magnitude of 6.3 in the Baoding seismic region.
The thermoelastic stress change produced by long-term water injection may reduce the minimum horizontal principal stress, which in turn creates a high risk of faulting at fault segments at critical (or subcritical) stress states. Therefore, precise assessment of the initial stability of the main active faults is a necessity when assessing the fault slip potential in a given study area. The results of the present analysis shed light on the strong correlation between fluid injection into deep geothermal reservoirs and faulting instability. However, the mechanism of faulting instability due to fluid injection still remains concealed. To remedy this lack of assessment of fault slip potential produced by fluid injection in Xiong'an New Area, in this paper, we predicted the changes in FSP values that may arise from large-scale geothermal exploitation, and this approach can provide important implications for deep geothermal reservoir exploitation in China, especially for other carbonate geothermal reservoirs that are close to the Quaternary active faults. In the future, we will focus on building an exact model of the three-dimensional stress field of the Xiong'an New Area, and assessing the effects of poroelasticity on the stability of nearby faults due to pore pressure diffusion.  Table  S1: Results hydraulic fracturing in situ stress measurements in Shunping county of Hebei province.