Modelling of ﬂuid pressure migration in a pressure sensitive fault zone subject to cyclic injection and implications for injection-induced seismicity

injection veriﬁes the accuracy of the numerical model. The parametric study reveals that the injection pressure attenuation, quantiﬁed by the amplitude ratio and phase shift, is enhanced by a lower initial fault permeability, a smaller stress sensitivity coefﬁcient and a shorter period of pressure cycle (i.e. a higher frequency). Besides, the amplitude of the pressure cycle has a negligible effect on the injection pressure attenuation. We also discuss the implications of our results for the less amenable far-ﬁeld seismic hazard and post shut-in seismicity.


I N T RO D U C T I O N
Induced seismicity associated with hydraulic stimulation treatments in geo-energy reservoirs can be used to map the stimulated volume. However, large magnitude induced earthquakes impede the development of geo-energy systems, such as deep geothermal systems (e.g. Grigoli et al. 2017;Kim et al. 2018;Zang et al. 2014;Hofmann et al. 2018;Schultz et al. 2020). To achieve permeability enhancement with minimized seismic hazards, cyclic injection has been proposed to replace conventional monotonic injection schemes (Zang et al. 2013;Hofmann et al. 2018;Zimmermann et al. 2019). The effectiveness of cyclic injection in mitigating seismic risks has been proven at various scales. In the laboratory scale, cyclic water injection can reduce the peak rate of injection-induced fault slip and the resulting seismic moment release (Ji et al. 2021a,b, Ji et al. 2022b. Cyclic hydraulic fracturing of core samples results in lower maximum acoustic emission (AE) amplitudes compared to monotonic hydraulic fracturing (Zhuang et al. 2019). At theÄspö Hard Rock Laboratory (Sweden), representing the mine scale, cyclic injection is characterized by significantly larger b-values of AEs, suggesting a smaller portion of large magnitude seismic events and thus a reduced seismic hazard (Niemz et al. 2020). In the field case of the Pohang Enhanced Geothermal System (EGS) in South Korea, a tailor-made cyclic injection into well PX-1 combined with a traffic light system has restricted the maximum magnitude below M w 2 during the stimulation between 2017 August 7 and 14 and during the subsequent flowback period (Hofmann et al. 2019). It must be noted that two months after further stimulation operations in well PX-2, an M w 5.5 earthquake occurred on 2017 November 15 near the geothermal site, which was linked to pore pressure change caused by the fluid injection into well PX-2 (Yeo et al. 2020). However, the optimization of cyclic injection parameters, including cyclic frequency and amplitude, considering reservoir hydraulic properties for minimizing seismic risks remains challenging.
In typical EGSs in crystalline reservoirs, numerous fractures may intersect with the reservoir and form several dominant fault zones (Baria et al. 1999;Genter et al. 2000;Jung 2013). Fluid is commonly injected near the fault zones, which normally act as major fluid pathways due to their high permeability compared to the surrounding rock matrix (Evans et al. 2005;Anyim & Gan 2020). The injected fluid migrates along the permeable fault zone and controls the occurrence of induced seismicity. The causal link between fluid diffusion and the expanding seismic cloud has been established (Shapiro et al. 1997), developed (Shapiro & Dinske 2009) and acknowledged (e.g. Katayama et al. 2018;De Barros et al. 2022). The fluid diffusion process is nonlinear because the local fault permeability is a function of local fluid pressure (Pedrosa 1986;Zhang & Mehrabian 2020. Particularly, pore pressure increase reduces the effective stress and thus increases the permeability along the fault, which in turn affects the fluid diffusion and the associated seismic hazard. The cyclic change of injection pressure could further complicate this nonlinear process. Therefore, it is crucial to understand the fluid pressure diffusion process during injection pressure cycling into a fault zone with pressure-dependent permeability. In this study, we numerically simulate the fluid pressure migration in a fault zone subject to cyclic injection considering pressuredependent permeability. Meanwhile, we also present an analytical solution for cyclic fluid injection into a pressure-sensitive fault for the first time. The numerical model is first described in detail and then validated by laboratory results and the newly proposed analytical solution. We then examine the influence of the intrinsic hydraulic fault properties (i.e. initial fault permeability and stress sensitivity coefficient) and the extrinsic injection parameters (i.e. period and amplitude of the pressure cycle) on injection pressure attenuation by performing a comprehensive parametric study. We further explore the temporal and spatial evolution of fluid pressure along the fault zone during a typical phase shift and finally discuss the implications for injection-induced seismicity.

Model configuration
The model construction is based on the geological conditions of typical EGS sites, such as Basel and Soultz-sous-Forêts, where there is a subvertical fault zone near the injection well, and the permeable fractures allow the injected fluid to migrate into the fault zone (Evans et al. 2005;Häring et al. 2008;Baisch et al. 2010). This can be simplified to obtain a conceptual 3-D reservoir model displayed in Fig. 1(a) (adapted from De Simone et al. 2017), including a subvertical fault zone directly intersected by the vertical injection well drilled from the centre of the reservoir with a side length of 1000 m. Due to the permeability contrast between the fault zone and the rock matrix, fluid is assumed to flow only along the fault zone. Thus, a 1-D fault model is sufficient to represent a horizontal cross-section at 4000 m depth of the 3-D reservoir model (Fig. 1b). Considering the symmetry of the 1-D fault with respect to the injection well, we simulate fluid pressure migration along the half-length of the fault (500 m). The 1-D linear flow of fluid along the fault zone is justified by the large ratio of fault length to fault width, which is normally on the order of 1/50 (Zang & Stephansson 2010). Note that the width of the fault will not influence the modelling results because of the linear flow assumption; here the fault width is set to a field-measured value of 15 m (Evans et al. 2005). In the numerical model, we assign a variable fluid pressure along the fault width and the other three sides were no-flow boundaries bounded by granite matrix (Fig. 1b). The finite element model is built in COMSOL Multiphysics (COMSOL Multiphysics 1998;Li et al. 2009;Ji et al. , 2022a with the minimum and maximum triangular element sizes of 0.01 and 5 m, respectively.

Governing equations
Along the fault length direction, the fluid flow equation in the Darcy flow regime derived from the mass balance equation can be written as where x is the distance along the fault length; p is the pore pressure; μ is the fluid viscosity; C is the total compressibility; is the porosity; t is the time; and k(p) is the pressure-sensitive fault permeability expressed as (Pedrosa 1986) where k i is the initial fault permeability; p i is the initial fluid pressure; σ ' is the change of effective normal stress; σ is the normal stress; σ i is the initial normal stress; and γ is the stress sensitivity coefficient.
Since the surrounding granite matrix will not be perturbed by the fluid injection in the fault zone due to the ultralow permeability, the normal stress on the fault is assumed to be constant at the depth of 4000 m (Fig. 1a), and thus the local fault permeability changes only with the local fluid pressure according to (3)

Validation against laboratory results
We performed a cyclic fluid injection test using the MTS 815 rock mechanics test system on a sawcut rectangular fault to validate the numerical model. The fault was created by first halving a 100mm-long and 50-mm-diameter granite core parallel to the sample axis and then grinding the sawcut surfaces using silicate-powder of 9 μm particle size (Fig. 2a). The granite core was drilled from a Luhui granite block from Shandong Province, China. The granite is mainly composed of 43 per cent feldspar, 28 per cent quartz, and 25 per cent illite (Zhao et al. 2017) with a low porosity ranging from 0.3 to 0.55 per cent and a low permeability in the order of 10 −19 m 2 ). After the sample was placed into the triaxial cell, the sample was fully saturated under a hydrostatic confining (silicon oil) pressure of 2 MPa and fluid (distilled water) pressure of 0.2 MPa. After full saturation, the sample was pre-conditioned by cycling the hydrostatic confining pressure from 20 to 40 MPa to remove any misalignment and plastic deformation (e.g. Kohli & Zoback 2013;Ji et al. 2021c). Afterwards, fluid was injected from the bottom inlet and the pressure was monitored in the top outlet (connected to a dead volume) when the fault was under a constant hydrostatic pressure of 40 MPa (equivalent to normal stress on the  fault). The injected fluid pressure was cycled between ∼7.5 and ∼12.5 MPa with a cycle period of ∼100 s for 10 cycles (Fig. 2b). The low permeability of the granite matrix ensures that the fluid flows dominantly through the fault. Using COMSOL Multiphysics, a finite element numerical model was built based on the laboratory conditions to simulate fluid migration in a planar fault with a length of 100 mm and a width of 50 mm. The minimum and maximum triangular element sizes are 1e-6 m and 1e-3 m, respectively (inset in Fig. 2b). Because the fluid pressure difference between the inlet and outlet was so small in the cyclic fluid injection test that we could not consider the stress sensitivity coefficient and thus it is set to zero (e.g. Nicolas et al. 2020). For the laboratory sawcut fault, the fault porosity approaches unity. Here, we set it to 0.9 (Cornelio et al. 2020). The initial fault permeability was measured as 3.27e-12 m 2 under a hydrostatic confining pressure of 40 MPa before the pore pressure cycling. The viscosity of the injected water at room temperature is 1e-3 Pa s. We directly used the experimentally applied injection pressure as the boundary condition at the injection side and adjusted the total compressibility to match the outlet pressures obtained from the numerical model and cyclic fluid injection test. Fig. 2(b) shows that the numerically derived and experimentally measured outlet pressures agree well with each other. The final total compressibility is 1.5e-5 Pa −1 , which is larger than the sum of typical laboratory fault compressibility (1e-9 to 1e-6 Pa −1 ; Tan et al. 2019) and water compressibility (4.5e-10 Pa −1 ; Farkas et al. 2021). This is primarily because of the additional compliance contributed by the fluid system including fluid pumps, dead volume fluid and tubing (e.g. Brantut & Aben 2021;Ji et al. 2021c). The input parameters in the numerical model for laboratory validation are listed in Table 1.

Validation against analytical results
Due to the near-uniform fluid pressure distribution over the laboratory fault in the cyclic injection experiments (Fig. 2b), the stress 1e-3 1e-3 a The compressibility in the laboratory scale is the total compressibility of the fault, fluid, fluid pumps and tubing. b The compressibility in the field scale for analytical validation is the total compressibility of the fault and the injected water. sensitive coefficient is not included in the laboratory validation (e.g. Nicolas et al. 2020). Therefore, we solved a mathematical model for fluid pressure evolution in a fault zone with pressure-dependent permeability subject to cyclic injection. The detailed steps for the described analytical solution are shown in the Appendix. Again, in the analytical model we neglect the normal stress change caused by fluid injection due to the low permeability of the surrounding host rock (eq. 3). For validation, we considered typical values for input parameters as listed in Table 1. Specifically, the fault porosity is set as 0.3 (Gan et al. 2021). The compressibility of the fault and fluid are both 4.5e-10 Pa −1 (Farkas et al. 2020), leading to a sum of 9e-10 Pa −1 , and thus we set the total compressibility here to 1e-9 Pa −1 . We set the initial fault permeability to 1e-13 m 2 (Gan et al. 2021), and the stress sensitivity coefficient to 2e-8 Pa −1 (Zhang & Emami-Meybodi 2020). The hydrostatic initial pore pressure at 4000 m depth is 40 MPa. In this study, we simulate cyclic injection pressures using sine waves. The base injection pressure is 60 MPa with a cycle period of 4 days and a pressure amplitude of 2 MPa. The viscosity of the injected water is constant at 1e-3 Pa · s. The numerically derived fluid pressures at 100 and 500 m from the injection well match well with the analytical solutions to the same problem (Fig. 3a). The spatial fluid pressure distribution after 1-day, 6-day and 12-day cyclic injection obtained from the numerical model also agree with the analytical solution (Fig. 3b).

R E S U LT S A N D A N A LY S I S
The advantage of the analytical solution is that it can offer more mechanistic explanations of the problem through mathematical terms without entailing numerical errors, while the numerical solution may be slightly inaccurate due to numerical errors. However, the developed analytical solution (see the Appendix) in our study includes two nested numerical integrations in eqs (A24) and (A26), and thus its computational efficiency is much lower compared with the numerical model. Because Sections 4 and 5 call for plotting numerous data points, we used the numerical model here for the parametric study so that a significant amount of time can be saved without sacrificing the accuracy since the numerical and analytical solutions are mutually validated (Fig. 3). We conducted a comprehensive parametric study to examine the influence of various variable parameters on fluid pressure migration along a fault zone, including hydraulic fault properties (i.e. initial fault permeability and stress sensitivity coefficient), and cyclic injection parameters (i.e. period and amplitude of pressure cycle). The other constant parameters, including fault porosity, total compressibility, initial pore pressure, base injection pressure as well as fluid viscosity, are set following the case for the analytical validation as listed in Table 1. Because the parameters used for the parametric study here are typical realistic values from the field and injection design, we present the results below using the absolute values of these parameters for  better understanding and to provide more direct guidance to the field operation.

Effect of hydraulic fault properties
We investigate the effect of initial fault permeability k i varying from 1e-15 m 2 to 1e-12 m 2 . The period and amplitude of the cyclic injection pressure are 2 days and 4 MPa, respectively. The remote pressure at 500 m from the injection well, denoted as p re , changes cyclically when k i is relatively high (Fig. 4a). However, the cyclical characteristic disappears within the 12 days of investigation when k i is lower than 1e-14 m 2 . We may refer to the critical initial fault permeability below which the cyclic characteristic of the remote pressure disappears as the cyclic breakdown permeability. Moreover, an extremely low initial fault permeability of 1e-15 m 2 results in an unchanged remote pressure. Note that theoretically there is no cyclic breakdown pressure because the remote pressure can finally exhibit the cyclic pattern given sufficient time of injection and monitoring, but with a much weaker amplitude. Nevertheless, we presume it is useful to introduce this term for scenarios of realistic fault hydraulic parameters, injection designs, operating time and time after active injection. We define two parameters, i.e. amplitude ratio α and phase shift ϕ, to characterize the attenuation of the cyclic injection pressure. The amplitude ratio is defined as the ratio between amplitudes of remote pressure and injection pressure. The phase shift is obtained by normalizing the phase shift between injection pressure and remote pressure with respect to the cycle period of injection pressure. Thus, a larger amplitude ratio and a smaller phase shift signify a smaller attenuation of the injection pressure after migrating through the fault. The remote pressure is not exactly a sine wave as the injection pressure due to the nonlinear diffusion process modelled in this study (Fig. 4a). Here we use a sine function to fit the remote pressure once it reaches a quasi-periodic pattern to derive the amplitude ratio and phase shift. For example, as shown in Fig. 4(b), in the case with an initial fault permeability of 5e-13 m 2 , the amplitude ratio and phase shift are estimated as 3.39/4 = 0.85 and 0.18/2 = 0.09, respectively. The amplitude ratio and phase shift are plotted in Fig. 4(c) as a function of the initial fault permeability ranging from 5e-14 m 2 to 1e-12 m 2 . It is found that a higher initial fault permeability favours a higher amplitude ratio and a smaller phase shift, reducing the attenuation of the injection pressure. The cases with the initial fault permeability varying from 1e-15 m 2 to 1e-14 m 2 are not plotted in Fig. 4(c) because the remote pressures in these cases show no cyclical features, which further confirm that a lower initial fault permeability could enhance the injection pressure attenuation.
Considering that the initial fault permeability plays a dominant role in controlling the remote pressure, the effect of the stress sensitivity coefficient is investigated by simulating two scenarios, one with a higher initial permeability of 1e-12 m 2 and the other with a lower fault permeability of 1e-13 m 2 . Again, the period and amplitude of the cyclic injection pressure is 2 days and 4 MPa, respectively. The stress sensitivity coefficient ranges from 2e-8 Pa −1 to 8e-8 Pa −1 and it exerts a smaller effect on the remote pressure when the initial fault permeability is higher (Figs 5a and b). This is quantitatively shown in Fig. 5(c), which demonstrates that the amplitude ratio is higher, and the phase shift is smaller in the case with a higher initial fault permeability. The amplitude ratio increases, and the phase shift decreases with higher stress sensitivity coefficient, suggesting that the injection pressure attenuation is suppressed when the pressure-induced fault permeability change is more significant.

Effect of cyclic injection parameters
To explore the effects of injection parameters (i.e. cycle period and pressure amplitude), we also consider two cases, one with a lower (1e-13 m 2 ) and the other with a higher (1e-12 m 2 ) initial fault permeability. The stress sensitivity coefficient is set to 2e-8 Pa −1 . To assess the effect of the cycle period, we set the amplitude of the injection pressure to 4 MPa and vary the cycle period from 1 day to 4 days. Here by setting the cycle period to several days, we intend to examine the effect of long-term cycles on fluid pressure migration along a pressure-sensitive fault zone (Hofmann et al. 2018), which is particularly important for understanding large-magnitude induced earthquakes. Consistent with the trend observed in Section 4.1, a smaller initial fault permeability promotes the injection pressure attenuation, signified by a smaller amplitude ratio and a higher phase shift (Fig. 6). The increase of amplitude ratio and decrease of phase shift with longer cycle period (Fig. 6c) indicate that the cyclic injection pressure with a longer cycle period is less vulnerable to attenuation after passing through a fault zone. The effect of pressure amplitude is examined by varying the pressure amplitude from 2 MPa to 8 MPa in the two cases with different initial fault permeability (Figs 7a and b). The amplitude ratio and phase shift remain relatively unchanged with increasing amplitude of injection pressure (Fig. 7c), suggesting that the injection pressure amplitude has a negligible effect on pressure attenuation.

Implications for injection-induced seismicity
For cyclic reservoir stimulation treatments, our study suggests that the cyclic injection parameters, especially the cycle period, could be optimized to better manage and mitigate injection-induced seismicity through controlling the phase shift and amplitude ratio. Although the pressure amplitude has little effect on phase shift and amplitude ratio, the absolute remote pressure increases with higher amplitude of the cyclic injection pressure. As revealed by the field data in the Denver earthquakes induced by fluid injection, there is a strong positive correlation between fluid pressure and seismicity rate (Healy et al. 1968). Thus, the pressure amplitude in cyclic injection is also a critical parameter, which should be adjusted according to the seismic response of the reservoir during injection, similar to the cyclic soft stimulation treatment performed in 2017 at the Pohang EGS site (Hofmann et al. 2018(Hofmann et al. , 2019. Our results also highlight the importance of initial fault permeability and stress sensitivity coefficient in controlling fluid pressure migration along the fault zone. Particularly for a fault zone with a low initial fault permeability and a small stress sensitivity coefficient, our simulations suggest that the cyclic pattern is significantly suppressed and delayed further away from the injection point (Figs 4a and 5b). This damping is more pronounced for higher cyclic injection frequencies. In these cases, the long-term cyclic injection modelled in this study could be accompanied by medium-term and short-term cyclic injection with higher peak pressures (Hofmann et al. 2018) to stimulate the fault rocks following the fatigue-hydraulic fracturing concept (Zang et al. 2013) to increase the fault permeability. Therefore, detailed reservoir characterization before the cyclic injection operation, and near-real time and high-resolution monitoring of the reservoir responses during the injection operation are important for cyclic injection parameter optimization.
The phase shift of remote pressure with respect to the injection pressure provides significant implications for injection-induced seismicity. We choose a typical case from Fig. 4(a) with an initial fault permeability of 5e-14 m 2 to show the spatial and temporal evolution of fluid pressure along the fault during a typical phase shift (Fig. 8a). In the later stage of this case, for example, from 10.5 day to 11.5 day, the injection pressure starts to decrease from 64 MPa to 56 MPa over 1 day (inset in Fig. 8a). The remote pressure at 500 m from the injection well keeps increasing while the injection pressure is reducing (Fig. 8a). Within 500 m from the injection  well, there is a critical zone with the local fluid pressure showing a tendency from rise to decline. As shown in Fig. 8(b), the fluid pressure at 200 m, which is within this critical zone, first increases to ∼61.5 MPa and then decreases to ∼60.9 MPa. According to the temporal evolution of fluid pressure gradient along the fault zone shown in Fig. 8(a), the first increase of the fluid pressure at 200 m from the injection well is due to the diffusion of fluid pressure from the zone near the injection well, suggesting that it serves as the fluid pressure sink before ∼940 000 s. Later, the local fluid pressure at 200 m from the injection well decreases because it starts to act as the fluid pressure source for flowback to the injection well and diffusion to the far-field after ∼940 000 s. Therefore, we refer to this zone as the fluid pressure pocket (Fig. 9). The presence of the fluid pressure pocket could possibly explain why the seismic hazard in the far field is less controllable by the wellhead pressure, especially in the later stage of an injection project (Healy et al. 1968;McGarr et al. 2015).
The amplitude ratio between remote pressure and injection pressure is important for understanding post shut-in seismicity. Our results confirm the previous finding that the closure of near-injection well fractures or faults could enhance the post shut-in seismic hazard (Ucar et al. 2017). Specifically, a higher stress sensitivity coefficient means a more significant fault normal closure in response to fluid pressure reduction near the injection well during the shut-in period. We show that the amplitude ratio increases with higher stress sensitivity coefficient (Fig. 5c), primarily because the more significant normal closure (i.e. reduction of aperture) of the near-injection well fault segment during injection pressure reduction further drives fluid pressure migration towards the remote region and increases the remote pressure (Fig. 9). The enhanced remote pressure may have induced the post shut-in seismic events near the periphery of the stimulated region, as observed at the Basel EGS in Switzerland (Häring et al. 2008;Mukuhira et al. 2016), and at the Paralana EGS in Australia (Albaric et al. 2014). One intention of the modelled long-term cyclic injection is to allow some time after a high-rate injection phase to observe seismicity (i.e. post shut-in events or events during low injection rate) to be able to react before another high-rate injection phase is started and introduced into the system (Hofmann et al. 2018). Therefore, our results may also confirm the feasibility of the diagnostic Seismogenic Fault Injection Tests (SFIT) using repeated injection proposed by Schultz et al (2022) to constrain the fault's seismogenic parameters in their statistical models for post shut-in seismicity. However, the diagnostic repeated injection in SFIT should be different from the simplified sine-wave cyclic injection in our models. As has already been pointed out by Schultz et al (2022), the injection in SFIT should be designed with non-periodic intervals and different injection and shut-in rates to better discern any time-dependent and rate-dependent features of post shut-in seismicity. This can be considered in future modelling studies.

Limitations of the numerical model
Although fluid pressure migration along a single planar fault subject to cyclic injection considering pressure-dependent permeability may be a primary process controlling the occurrence of injectioninduced seismicity (Shapiro et al. 1997;Shapiro & Dinske 2009), the linear flow along a pressure sensitive fault zone without diffusion into the surrounding host rock modelled in the present study cannot represent the case with more complex fault architecture. Fault damage zones and fault branches are the two main features of the complex fault architecture that are simplified in this study. The highly permeable zone of damaged host rock bounding the lowpermeability fault core facilitates fluid pressure diffusion and thus increases the extent of fluid pressure change along the fault length direction (e.g. Yehya et al. 2018;Yeo et al. 2020;Lei & Tsang 2022). The enhanced diffusion of fluid pressure may increase the likelihood of encountering critically stressed and favourably oriented fault segments and enhance the seismic hazard because small pore pressure changes could induce/trigger fault reactivation ). This has already been implicitly simulated in our cases with higher initial fault permeability values. However, the anisotropic permeability of the fault damage zone (Yehya et al. 2018) and the far-field stress orientation (Lei & Tsang 2022), which may exert important controls on the direction of pressure migration in the fault damage zone, are not explored in our model. Besides, the static stress changes induced by prior seismic events on the fault branch close to the injection site may trigger the sequential reactivation of adjacent fault branches (Sumy et al. 2014), such as in the case of the 2011 Prague, Oklahoma, earthquake sequence (Keranen & Weingarten 2018). Such complex stress changes associated with fault branches are not captured in the current model. Moreover, the shape of cyclic injection pressure, such as triangle, sine and rectangular waves, may also influence the fluid pressure migration and the associated seismicity. Particularly, the amplitude ratio produced in the case of rectangular wave is the largest followed by sine and triangular waves with the same pressure amplitude and cycle period, while the phase shift remains similar in the three cases (Zhu et al. 2022). But the differences in amplitude ratio caused by different shapes of cyclic injection pressure may not be significant (Zhu et al. 2022). In addition, other factors associated with fluid injection, mainly including shear slip (aseismic or seismic), poroelasticity and thermoelasticity, may also have significant effects on fluid pressure migration and seismic hazards along a fault zone (e.g. De Simone et al. 2017). The occurrence of injection-induced seismicity is also dependent on the fault frictional property (e.g. Larochelle et al. 2021) and stress state (e.g. Cappa et al. 2018). Taken together, these factors make the management of injection-induced earthquakes even more challenging and merit further studies in the future.

C O N C L U S I O N S
We numerically investigated the fluid pressure migration in a fault zone subject to cyclic injection pressure with pressure-dependent permeability. The parametric study demonstrates that a lower initial fault permeability, a smaller stress sensitivity coefficient, and a shorter cycle period (i.e. a higher frequency) of injection pressure collectively favour the attenuation of injection pressure, signified by a smaller amplitude ratio and a larger phase shift. We also show that the pressure attenuation seems to be independent of the amplitude of the pressure cycle, but the pressure amplitude is also important in controlling the seismic hazard and should be optimized. The spatial and temporal evolution of fluid pressure in the fault zone during a typical phase shift suggests the presence of a fluid pressure pocket, which provides the fluid pressure source for flowback and far-field diffusion in the later stage of shut-in. This explains the limited control that wellhead pressure changes at the injection well can have on the far-field seismic hazard, particularly for the later stage of fluid injection into a low-permeability fault. Our results also indicate that the post shut-in seismic hazard could be enhanced by the normal closure of the near-injection well fault segment. Overall, we contend that automatically adapted cyclic injection parameters based on near-real time and high-resolution monitoring of the reservoir responses could be promising in steering the flow rates required for EGS systems to operate economically and environmentally friendly.

A C K N O W L E D G M E N T S
This work has been supported by the Helmholtz Association's Initiative and Networking Fund for the Helmholtz Young Investigator Group ARES (contract number VH-NG-1516). We are grateful to the Editor Andrea Morelli, Ryan Schultz and an anonymous reviewer for their constructive comments and suggestions that improved the quality of the paper.

C O N F L I C T O F I N T E R E S T
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.

DATA AVA I L A B I L I T Y
The data supporting the analysis and conclusions are available from the corresponding author upon reasonable request.
Then the η solution considering pressure-dependent permeability with a variable inner boundary condition becomes: where, we denote the first term of Eq. (A18) as η 0 , which represents the solution with variable injection pressure but without considering pressure-dependent permeability; and we denote the second term as η 1 , which represents the contribution of pressure-dependent permeability to the solution (Garcez et al. 2020;Zhang & Mehrabian 2021b).
Therefore, theṗ(x, t) solution can be obtained aṡ Finally, the solution for absolute pressure is