Modelling Rock Fracture Induced By Hydraulic Pulses

Soft cyclic hydraulic fracturing has become an effective technology used in subsurface energy extraction which utilises cyclic hydraulic flow pressure to fracture rock. This new technique induces fatigue of rock to reduce the breakdown pressure and potentially the associated risk of seismicity. To control the fracturing process and achieve desirable fracture networks for enhanced permeability, the rock response under cyclic hydraulic stimulation needs to be understood. However, the mechanism for cyclic stimulation-induced fatigue of rock is rather unclear and to date there is no implementation of fatigue degradation in modelling the rock response under hydraulic cyclic loading. This makes accurate prediction of rock fracture under cyclic hydraulic pressure impossible. This paper develops a numerical method to model rock fracture induced by hydraulic pulses with consideration of rock fatigue. The fatigue degradation is based on S–N curves (S for cyclic stress and N for cycles to failure) and implemented into the constitutive relationship for fracture of rock using in-house FORTRAN scripts and ABAQUS solver. The cohesive crack model is used to simulate discrete crack propagation in the rock which is coupled with hydraulic flow and pore pressure capability. The developed numerical model is validated via experimental results of pulsating hydraulic fracturing of the rock. The effects of flow rate and frequency of cyclic injection on borehole pressure development are investigated. A new loading strategy for pulsating hydraulic fracturing is proposed. It has been found that hydraulic pulses can reduce the breakdown pressure of rock by 10–18% upon 10–4000 cycles. Using the new loading strategy, a slow and steady rock fracture process is obtained while the failure pressure is reduced.


Introduction
Hydraulic fracturing techniques have been widely used to extract energy in oil/gas and geothermal industries since the first successful engineering demonstration in late 1940s (Lecampion et al. 2018). This technology has witnessed significant progress in last few decades and made an important contribution to the development of subsurface energies (e.g., shale gas and enhanced geothermal). The ultimate aim of this technology is to achieve a fracture network with the highest possible surface area for a given energy input in order to drain the largest volume of rock of either gas or heat. Breakdown pressure is a key parameter in hydraulic fracturing and is dependent on in-situ stress, pore fluid pressure and rock fracture properties (Wu 2018). Due to the large in-situ stress and fracture toughness of rock, the breakdown pressure is normally very high. One consequence of the large pressures and volumes required to generate effective rock fracture networks is induced seismicity (Ellsworth 2013;Lee et al. 2019) which has resulted in significant public concerns about safety and environment protection (Altmann et al. 2010;Segall 1989;Segall and Lu 2015). In particular, the 2017 Mw5.4 earthquake in Pohang, South Korea, was confirmed by scientists as being related to the enhanced geothermal system (EGS) drilling site nearby (Grigoli et al. 2018;Kim et al. 2018). In the UK, seismicity induced by hydraulic fracturing operations at shale gas sites has resulted in moratoria being imposed by the UK and Scottish governments (UK's OGA 2019). For the widespread application of hydraulic fracturing in extracting subsurface energies, a new improved fracturing technology which overcomes these limitations is needed.
The conventional hydraulic fracturing system employs a high static hydraulic pressure which is rather limited in its capacity to generate complex fracture networks and hence optimise the enhanced permeability of the rock mass . In recent years, a new technique based on cyclic or pulsating injection of fluid to fracture the rock is being developed (Chen et al. 2020;Hofmann et al. 2018;Lee et al. 2019;Li et al. 2014;Patel et al. 2017;Stephansson et al. 2018;Zang et al. 2017;Zhou et al. 2017;Zhuang et al. 2019). Through cyclic or pulsating injection, rock fatigue fracture will occur and the strength of rock to cracking is reduced (i.e., rock softening). Li et al. (2014) carried out experiments to fracture rock-like materials under different frequencies of pulses and found that the fracture initiation pressure was reduced and more complex fracture networks were generated by the 6 Hz pulse. Patel et al. (2017) experimentally investigated the effect of cyclic injection on reducing breakdown pressure for wet and dry rock samples and found that in dry sandstone, cyclic injection could decrease the breakdown pressure by 16% whilst it did not have any effect on breakdown pressure in saturated sandstone. Zang et al. (2017) carried out field monitoring of hydraulic fracturing in a 410 m deep tunnel by cyclic injections with multiple-flow rates. Through acoustic emission (AE) monitoring, they found that the total number of AE events was lower when the continuous injection scheme was replaced by multiple flow rate cyclic injections. Further, Zang et al. (2019) proposed the concept of fatigue hydraulic fracturing as a combination of cyclic hydraulic fracturing and pulse hydraulic fracturing to drive hydraulic fractures in a more controlled way. Consequently, a lower breakdown pressure (i.e., 10% reduction at laboratory scale test and 15% at the mine scale test) and wider fracture process zone were obtained by a combination of the cyclic progressive and pulse pressurization. Hofmann et al. (2018) applied cyclic stimulation (injection) to an enhanced geothermal system with the aim of managing seismic risks and found that the maximum magnitude of induced seismic events was less than that from the conventional injection treatment in the same well. To understand strategies for mitigate coal seam gas risk, Chen et al. (2020) used pulses to fracture coal samples in the lab with pre-cracks and found that the frequency and amplitude of the pulsation could cause fatigue of coal and yield different failure mechanism.
The studies above suggest that hydraulic pulses could generate fatigue of rock which lead to reduced breakdown pressures and more complex fracture networks. The application of hydraulic pulses to create a desirable fracture system in rock may have the additional advantage of reducing induced seismicity during the injection process. However, the mechanisms for the response of rock under hydraulic pulses is still very unclear, posing significant challenges to the practical application of the pulsating technology for subsurface energy extraction.
To elucidate the pulse-induced fracture mechanism of rock, numerical modelling can provide unique insights given the difficulties in conducting pulsating hydraulic fracturing experiments. Zang et al. (2013) employed a particle flow code (PFC) to simulate the hydraulic fracture by cyclic injection and found that the total number of induced AE events as well as the occurrence of larger magnitude events were both lowered. Yoon et al. (2014) modelled the hydraulic fracture process for a 2 × 2 km reservoir using the particle flow method and found that the cyclic injection technique consumed more fluid but that the total number of events including the larger magnitude events was reduced. However, discrete-element based particle models (including the work by Ma et al. 2017 andXu et al. 2017) have not considered the fatigue effect on the contact/bonding properties between particles. Lu et al. (2014 compared the effective stress disturbance zone of coal seam under pulses by considering different material properties and boundary conditions using the elastic finite difference method. Fatigue is a key degradation form of rock under pulsed/cyclic loading, resulting in a reduction in the breakdown pressure during the hydraulic fracturing process . However, amongst the limited literature for numerical modelling of rock response under hydraulic pulses, fatigue of rock has not been included in the constitutive relationship. The cohesive crack model is a very useful approach for modelling the fracture process zone (FPZ) ahead of the crack tip in rock, concrete or other quasi-brittle materials (Lisjak et al. 2014;Mahabadi et al. 2012;Mergheim et al. 2005;Xi et al. 2018;Yang and Frank 2008). It is widely used for interface elements in the finite element method (FEM) as well as the interface/joint elements in the finite-discrete element method (FDEM). It accommodates the stress singularity in FPZ and constitutively represents the softening behavior of quasi-brittle materials that still exhibit load-bearing capacity after peak stress. In particular, the fluid cohesive element with a pore pressure degree of freedom has been developed (Chen et al. 2009). This new element is used for the coupled hydro-mechanical problems, especially for the fluid-driven hydraulic fracture of rock (Carrier and Granet 2012;Chen et al. 2009;Guo et al. 2017;Lisjak et al. 2017;Nguyen et al. 2017;Sarris and Papanastasiou 2010;Xiang et al. 2019;.  developed a fully coupled 2D/3D hydro-mechanical model in which the fluid flow in a crack was expressed by fracture seepage in the broken joint elements. Lisjak et al. (2017) developed a finitediscrete element code which discretized a mesh by 3-node triangular elements and 4-node interface elements. Their results showed the FDEM method with cohesive elements could be used to obtain unique geomechanical insights into the coupled hydro-mechanical phenomena. Guo et al. (2017) simulated the fracture initiation and propagation in a layered reservoir by cohesive elements representing the interaction between layers and found the crack could penetrate into the layer interface and extend along the interface to generate a branched fracture. The numerical results for conventional hydraulic fracturing by using the fluid cohesive crack model have been compared with analytical solutions of the KGD (i.e., Khristianovice-Geertsmae-de Klerk) model, which proved the accuracy and robustness of the fluid cohesive crack model (Carrier and Granet 2012;Chen et al. 2009;Nguyen et al. 2017). Although the fluid cohesive model extends the capability of the cohesive model in modelling the coupled hydro-mechanical behavior of the rock under static hydraulic pressure, it has never been used for pulsating hydraulic pressure. To date, as far as the authors are aware, there is no numerical model which simulates the crack propagation in rock subjected to hydraulic pulses with consideration of the pulse-induced rock fatigue and the injection fluid flow in fractures. This makes the accurate prediction of the response of rock under pulsating hydraulic fracturing impossible using the current state of the art numerical tools.
This paper develops a new numerical method to model rock fracture induced by hydraulic pulses. The fatigue of the rock under cyclic loading is considered in the constitutive model in form of S-N (cyclic stress vs cycles to failure) curves. The fluid cohesive crack model is then applied to simulate the discrete crack propagation which is coupled by the hydraulic flow. The fatigue damage model is then implemented into ABAQUS by in-house FORTRAN subroutines. Through a worked example, the developed numerical model is demonstrated with application in pulsating hydraulic fracking of rock. The breakdown pressure and the fracture initiation and propagation are calculated and simulated. The model is then verified by comparing the results with those from experimental work. Further, the effects of frequency and flow rate for cyclic injection on borehole pressure are investigated. Moreover, a new strategy for pulsating hydraulic fracture is proposed. The borehole pressure development, fracture propagation and fatigue crack characteristics under the new strategy are obtained.

Cohesive Crack Model
Quasi-brittle materials such as rocks exhibit tensile strain softening behaviour due to an inelastic zone being developed ahead of the crack tip, often referred to as fracture process zone (FPZ), as shown in Fig. 1a. When a crack propagates in rock, the cracked surfaces may be in contact and are rough (Shah et al. 1995), due to various toughening mechanisms such as bridging, void formation or microcrack shielding (Bazant and Le 2017). Therefore, the cracked surfaces may still be able to sustain the tensile stress which is characterized by the softening degradation curve as shown in Fig. 1b. Hillerborg et al. (1976) first proposed the cohesive model to simulate discrete cracking in the fracture process zone of concrete and since then the cohesive crack model has been employed to model cracking of rock given the same quasibrittle nature (Lisjak et al. 2014;Mahabadi et al. 2012;Xi et al. 2018;Xi et al. 2021;Yang et al. 2019). Normally, the cohesive element is of zero thickness before it is cracked or damaged. Because the distance between the nodes is used as the measure of crack opening rather than a change in strain (which depends on the element length), the mesh dependency is significantly reduced. A stress-displacement law is employed to constitutively control the cohesive elements. The stress of cohesive elements is a function of the corresponding relative displacements of the crack surface. In general, the stress-displacement relationship for the cracks can be expressed as follows: where and are the cohesive stress and crack opening displacement, respectively; f is the non-linear function defining the relationship between stress and displacement, which defines the linear elasticity and post-peak softening behaviour.
As shown in Fig. 1b, for a cohesive crack, the stress will linearly increase to the peak strength by a penalty stiffness and then gradually decrease to zero. It should be noted that in the cohesive crack model uses the engineering convention defining tensile stress as positive rather than compressive stress (the more common convention in geology). The stress-displacement relationship in the linear elasticity stage can be expressed as follows: where K p is the penalty stiffness, which should normally be significantly larger than the material stiffness to ensure that there is almost no deformation of the fictitious crack before damage initiation (Yang et al. 2018a).
After reaching the cohesive strength, the tensile stress decreases, following certain softening rules. The softening behaviour can be in a linear, bilinear or exponential curve. The area underneath the softening curve is known as the specific fracture energy G f which is also known as the critical energy release rate (ASTM 2007;Vasconcelos et al. 2008). The damage of the cohesive element at the softening stage under a static loading can be obtained as follows (Yang and Frank 2008): where 0 is the critical displacement when the stress reaches the peak value (i.e., D s = 0); m is the maximum displacement during the loading history; f is the failure displacement when the stress reduces to zero (i.e., D s =1); Further, the residual stress of the damaged cohesive elements can be obtained according to the damage factor D s as follows:

Fluid flow model
The cohesive crack model can include pore pressure by adding an additional degree of freedom Lisjak et al. 2017;. As shown in Fig. 2, this pore pressure cohesive crack model enables fluid pressure implementation in the cohesive element, to facilitate modelling of the hydraulic flow-driven fracture process. The nodes of the conventional cohesive element are added with a pore pressure degree of freedom (DoF) at each node, besides the translational DoFs. However, the middle nodes of the cohesive element only have the pore pressure DoF. The tangential flow within the cohesive element is expressed as a Newtonian fluid following Poiseulle's law, the relationship between the pressure gradient and fluid flux is formulated as follows: where q is the fluid flux of the tangential flow; ∇p is the fluid pressure gradient along the cohesive zone. μ is the fluid viscosity. d is the crack opening displacement, which is defined as follows: where d 0 is the initial gap opening for fluid entering; is the true crack opening displacement under mechanical response. d 0 should theoretically be zero before the crack is open because the cohesive element is the fictitious interface element with zero-thickness. d 0 is a non-physical parameter to ensure that the flow equations can be solved robustly when the physical separation is zero. The effect of d 0 with a very small value on the flow equations is diminished with the opening displacement δ increasing.
The normal flow is described by defining a fluid leakoff coefficient for the cohesive zone surface, presented as follows: where q t and q b are the flow rates into the top and bottom of the element surface, respectively; c t and c b are the corresponding leak-off coefficients, with the unit m/(Pa·s), for the top and bottom surface of the cohesive element, respectively. p f is the fluid pressure at the element gap. p t and p b are the pore pressures of the adjacent pore-containing elements on the top and bottom surfaces, respectively. The fluid leak-off describes the fluid flow from the cracked cohesive element into the surrounding porous rock. The fluid leak-off coefficient is related to fluid velocity, fluid viscosity, porosity and permeability of the formation, the reservoir fluid compressibility, etc. (Yarushina et al. 2013). The coefficients can be determined by experiments. In this paper, the leak-off coefficient is assumed as a constant parameter with reference to other studies (Ghaderi et al. 2019;Feng and Gray 2019;Carrier and Granet 2012).
Further, the continuity equation of mass conservation can be expressed by the lubrication model as follows (Detournay 2016): where t is the flow injection time; the singular term Q(t) (x, y) describes the boundary condition at the injection point. Q(t) is the flow rate at time t. (x, y) is the Dirac function which equals to zero everywhere except for the flow injection point.
The fluid pressure acts as a traction on the surfaces of the fracture and drive the opening of the fracture. The general stress field of the cohesive zones are balanced by the fluid pressure of the cohesive zones. Thus, a coupled fluid pressure-traction-separation relationship exists between the cohesive zone defined by the stress-displacement law (i.e., Eq. (1)) and the lubrication equation model (i.e., Eq. (9)).
The rock matrix beyond the cohesive crack is assumed to be an elastic porous medium fully saturated with a single-phase fluid. The fluid flow in the porous rock matrix is described by Darcy's law (Whitaker 1986): where k is the permeability coefficient in unit m 2 . It should be mentioned that the permeability in Abaqus is defined as hydraulic conductivity K in unit m/s, which can be expressed as follows: where γ is the specified weight of the fluid.
The mechanical behaviour of poroelastic rock is described by the effective stress principle as follows (Terzaghi 1943): where ′ and are the effective stress tensor and total stress tensor, respectively; I is the unit tensor; p is the pore pressure value; is the Biot's coefficient. Because the Biot's coefficient may vary non-linearly with the pressure and because the pressure in the matrix is not the dominant factor, the Biot's coefficient is simplified as a constant value 1 for calculating the effective stress (Alam et al. 2010).

Fatigue Damage Model
The conventional cohesive crack model only considers the crack behaviour under static and monotonic loading. As discussed, the cyclic pulsating pore pressure can degrade the mechanical property of rock during the cyclic loading/ unloading process and therefore the constitutive model for cohesive crack should be changed. The most significant difference between fatigue cracking and static cracking is the strength degradation induced by the cyclic load (Nojavan et al. 2016). In this paper, we introduce fatigue-controlled damage into the cohesive constitutive law. As shown in Fig. 3, for a cyclic loading profile, the peak stress (strength) is reduced to N1 0 , N2 0 and N3 0 , by different failure cycles, i.e., N1, N2 and N3. After the fatigue strength is reached, the crack is initiated. The stress then gradually decreases as the displacement of the cohesive element increases, following the linear softening law defined in Eq. (3). Whether the failure displacement f is the same for fatigue and static softening due to the fatigue-affected strain-softening mechanism for rock is unknown (Erarslan and Williams 2011;Fan et al. 2017;Liu et al. 2017). The fatigue crack strength after n cycles can be expressed as follows: where n 0 is the current strength of the materials after experiencing n loading cycles; 0 is the static strength; D f is the fatigue damage parameter.
Experimentally established stress-life curves (S-N curves) are the most widely adopted criteria for fatigue analysis in metals, composites, concrete and rocks (Cerfontaine and Collin 2017; Chen et al. 2017;Khoramishad et al. 2010;Nojavan et al. 2016). To quantitatively account for the fatigue-induced strength degradation, we implement the S-N curves for rock to model the pulsation effects on rock fracture. An empirical form of S-N curves for rock is expressed as follows (Cerfontaine and Collin 2017): where max 0 is the stress level as the ratio of the maximum stress during the fatigue loading cycles max to the static strength 0 . N is the maximum number of loading cycles (fatigue life). A and B are the fatigue parameter. Cerfontaine and Collin (2017) collected existing data on rock fatigue experiments from published research and fitted in an S-N curve with A = −0.0278 and B = 0.9455.
As shown in Fig. 4, the fatigue life N for a given stress level max ∕ 0 can be obtained from the S-N curve as the inverse function of the Eq. (11): It should be noted that once the fatigue life N is reached, we assume the cohesive strength is degraded from the static strength 0 to the fatigue crack initiation strength N 0 . The effect of cyclic loading on the softening behaviour of rock is completely unknown so in this paper we do not consider the fatigue effect on softening behaviour. All existing experimental results in the literature took the failure cycles as rock fatigue life without developing the softening-cycle relation (Cerfontaine and Collin 2017;Williams 2011, 2012). This is mostly because the tensile softening behaviour under both static and cyclic loading for rock is extremely hard to experimentally determine. Assuming linear strength softening (Nojavan et al. 2016), the fatigue damage upon Δn cycles with a maximum tensile stress max can be expressed as follows: If the rock experiences a varying amplitude loading as shown in Fig. 4(b), the strength degradation rate (i.e., fatigue damage evolution rate) will be changed from one slope to another as per the corresponding loading profile.
The fatigue damage model has been coupled with the fluid cohesive crack model through FORTRAN subroutines and implemented into ABAQUS standard finite element code. The cyclic load acting on the cohesive elements can be induced by cyclic hydraulic pressure and/or cyclic injection, in terms of boundary conditions. The degradation of strength under fatigue loading consists of two steps: first, find out the static loading capacity. In this step, the static fracture properties are used and the failure load is determined for further reference. In the second step, the recorded loading history is used to calculate the residual strength from the fatigue damage model and the static strength is replaced by the residual strength for the next increment. The steps run cycle by cycle until the maximum stress reaches the residual strength and softening damage occurs. Thus the fatigue damage and cracking damage will both accumulate during the cyclic loading. It should be mentioned that most experiments for rock fatigue fracture (e.g., Cerfontaine and Collin 2017) employed Brazilian disc tests to generate rock tensile cracking in which the pre-peak (before damage initiation) and post-peak (after crack damage initiation) cannot be distinguished. Therefore, the post-peak fatigue fracture behaviour of rock and the effect of cyclic load on fatigue damage accumulation are completely unclear to date. This model considers the effect of fatigue on the mechanical degradation of rock through reducing the peak (tensile) stress only.

Worked Example and Verification
As a demonstration of the developed numerical method and techniques in solving pulsation-induced rock fatigue fracture, a 2D plane strain worked example with the dimension 9 × 12 m for hydraulic fracturing KGD problem (Nguyen et al. 2017) is carried out. In-situ stress and symmetry boundary condition are applied to the model. Two kinds of elements are used, i.e., six-node pore pressure cohesive elements for the potential interface crack and four-node bulk elements. As shown in Fig. 5, there are 19,648 solid elements and 647 cohesive interface elements inserted along the predefined fracture path which is parallel to the direction of the maximum in-situ stress. A very fine mesh is generated for the area around the crack with the size 0.005 m. The values for all basic parameters are shown in Table 1, together with their sources.
The constant flow rate injection typically used in traditional hydraulic fracturing is first modelled as a reference. Figure 6 shows the injection pressure and fracture length development as a function of injection time under the injection flow rate 0.0005 m 3 /s. It can be seen that the borehole pressure rapidly increases up to the peak pressure (known as breakdown pressure), followed by a sharp drop to about half the maximum pressure, which then steadily decreases during the injection life. The breakdown pressure under monotonic or static loading is 13.0 MPa in these models. The speed of fracture propagation becomes slower because the amount of fluid leak-off increases as the fracture propagation increases the crack surface area. Finally, the fluid injection rate and leak-off rate are balanced at the residual pressure. Modelling the typical hydraulic fracturing problem under a constant fluid rate injection with a cohesive crack model has been intensively investigated and verified (Carrier and Granet 2012;Chen et al. 2009;Guo et al. 2017;Li et al. 2017;Lisjak et al. 2017;Nguyen et al. 2017;Xiang et al. 2019;. Therefore, detailed discussion and verification of this traditional problem are not presented in this paper. One of the main purposes of pulsating hydraulic fracturing is to reduce the breakdown pressure. A hydraulic pulse with a maximum pressure of 86% of the static breakdown pressure (i.e., 11.2 MPa) is used as the loading input (Fig. 7). The frequency of the pressure pulse is 10 Hz which is within the range shown to be effective by Chen et al. (2020). To verify the possible error that may be caused by the difference of static breakdown pressure for the boundary condition (flow rate or pressure) and loading rate, the monotonic fluid pressures with maximum values of 12.9 and 13.0 MPa are applied to the KGD model, respectively. And the loading Fig. 5 Mesh grid of the numerical example rate is the same as that for the cycle of pulses (i.e., the pressure increases to a maximum in 0.1 s). We found that the monotonic pressure 13.0 MPa can break the rock while the pressure 12.9 MPa cannot. However, after fracture initiation, pressure-based boundary condition will lead to different responses comparing with flow rate based boundary condition and hence different fracturing strategy should be developed accordingly. Therefore, the static breakdown pressure can be regarded as 13.0 MPa. The hollow points on the pressure line in Fig. 7 represent the hydraulic pressure at the calculated model increments. For fatigue modelling under cyclic loading, incremental points where the model calculates must include the maximum and minimum loads at the turning points. Thus the turning points for the pressure pulses are enforced by a predefined time points array. Figure 8 illustrates the pressure distribution and fracture process of the rock under the pressure pulse. It can be found that, during the first few cycles, the rock is not fractured, however, fatigue damage is being developed. As the loading cycle proceeds, crack initiation occurs at the peak pressure in the 88 th cycle (i.e., 87.5 s in Fig. 8b). The crack is then closed (Fig. 8c) due to the unloading and in-situ stress. As the cyclic loading is repeated, the crack re-opens and propagates further as shown in Fig. 8d and e. Figure 9 illustrates the tensile stress history, fatigue damage and static damage of the first cracked cohesive element. The crack element experiences cyclic tensile stresses smaller than the static strength, and Fig. 9 clearly demonstrates that the fatigue damage D f progressively increases  shows that hydraulic pulses can reduce the breakdown press for rock fracture by using cyclic loading. Figure 10 shows the fracture length development and fluid injection volume over injection time. It can be seen that the crack rapidly propagates after the crack initiation. The increase of crack length is in line with the growth of the fluid injection volume which confirms that the injection volume provides the driving force to crack propagation. The pressure driving the fracture propagation (See Fig. 6) is significantly smaller than the breakdown pressure. Therefore, if pressure-control persists after breakdown pressure is reached, the pressure applied will generate very fast crack propagation as demonstrated in Fig. 10. Moreover, the amount of fluid needed to maintain the pressure of the pulse is rapidly increased under the pressure-control, which at an actual site, could potentially exceed the capacity of the injection pump. Therefore, while pressure-controlled pulses can be effectively used at the beginning to soften (fatigue) the rock and generate crack initiation, the crack propagation stage would need additional and alternative attention. From the point of crack initiation, using pulses with a constant maximum pressure will generate runaway fractures, which may not necessarily be the optimal fracture networks for engineering applications.
The ratio of the maximum pulse pressure to the static or monotonic pressure (i.e., P c ∕P s ) is used to evaluate the effect of the hydraulic pulsation. Figure 11 illustrates the relationship between P c ∕P s and the number of cycles to cracking. The breakdown pressure is reduced to between 82 to 90% of the static breakdown pressure after between 10 and 4000 cycles. When higher maximum pressure of pulses is used, fewer cycles are needed to fracture the sample. The function for the pressure-cycles is obtained by linear fitting with the coefficient of determination R 2 = 0.89 as follows: Due to the lack of analytical solution on hydraulic fracture of rock under pulses, we cannot compare the developed numerical method with analytical solutions. To verify the proposed numerical model, the results are compared with experimental data from the literature (Zhuang et al. 2019). Zhuang et al. (2019) conducted tests on rock fracture by pressure pulses and obtained the relationship between the maximum pressure of pulses and the number of cycles to break down for granite samples. To obtain the fatigue parameters for granite in the S-N curve (Eq. 14), we carried out Brazilian disc splitting tests on Jinan granite with fine grains through applying monotonic and cyclic loads. First, granites were made to disc specimens with dimensions of 65 mm in diameter and 26 mm in thickness following the ISRM (18) P c P s = −0.026log 10 (N) + 0.91   (1978). Then the monotonic load was applied to the specimen by MTS rigid fatigue machine as shown in Fig. 12a. The average value of the peak loads from three tested specimens was used as the average peak load for monotonic peak load. Further, cyclic loads with the maximum load accounting for 88%, 84%, 73% and 70% of the monotonic peak were applied on the specimens, respectively. The minimum load for cyclic tests was set 2 kN to ensure the specimen did not move in the loading apparatus. Finally, the number of cycles to the failure of specimen were obtained. Figure 12b shows the relationship between the peak load and failure cycles from the cyclic Brazilian disc tests on 15 specimens. Therefore, the S-N curve can be obtained by fitting the experimental data. The fatigue parameters A and B for the granite were therefore calculated as −0.107 and 0.977, respectively, with the coefficient of determination R 2 = 0.834. A half-sample (due to symmetry) with the dimensions of 50 mm in diameter and 100 mm in height is modelled as illustrated in Fig. 13(a) using values for essential inputs from the tests by Zhuang et al. (2019) (i.e., Table 2). The measured fatigue parameters are used in the verification model. To accommodate the difference between Brazilian tensile strength values of Zhuang (2919) and direct tensile strength used in the model, the Brazilian tensile strength was reduced by about 15% in accordance with (Li and Wong 2012). From experimental and numerical results, the static breakdown pressure with a constant injection rate is 6.9 MPa. Figure 13b shows the results for the relationship between the breakdown pressure and cycles of hydraulic pulses from numerical results and experiments. It can be seen that the numerical results are generally in agreement with the experimental results, i.e., within the 95% prediction band. It can be seen that the predicted results well fit in the range of the experimental data from Zhuang et al. (2019). This validates this newly developed rock fatigue model under hydraulic pulses.

Controlling Fracture Growth with Cyclic Injection
As discussed, the pressure pulses with constant maximum pressure will result in a fast fracture propagation after initiation which is not viable for generating a complex fracture network. Some researchers have experimentally investigated pulsating hydraulic fracturing of rock by employing cyclic schedules of repeatedly starting and stopping the injection pump (Chen et al. 2020;Patel et al. 2017;Zhou et al. 2018;Zimmermann et al. 2018). The pressure at the injection point increases during fluid injection, while the fluid leakoff reduces the pressure when injection is paused. Using the values listed in Table 1, cyclic injections with different frequencies and injection flow rate are modelled. The period of the injection and pause in one cycle are the same and 'flow rate' refers to injection flow rate, not the average flow rate in one cycle. Figure 14 shows injection flow. Lower flow rate leads to slower force build up and, therefore, greater pressure variation across the time history. Figure 15 shows the effect of injection frequency on the borehole pressure. Larger frequency pulses result in smaller pressure variation during the crack propagation stage because the time for fluid leak-off in a cycle under a higher frequency is shorter given the same flow rate. These numerical results for the frequency effects on borehole pressure are qualitatively consistent with the experimental results from Chen et al. (2020). As the fracture propagates to a longer distance, a larger flow rate would be required to pump in more liquid to compensate for the leak-off through the increased crack surface area. To further investigate the effect of the increased flow rate, a cyclic injection with the flow rate increasing from 5.0 × 10 -6 m 3 /s to 5.0 × 10 -4 m 3 /s is modelled. The frequency of the cyclic injection is 0.25 Hz and the flow rate is increased by 5.0 × 10 -6 m 3 /s after every 20 cycles as the interval. Figure 16 shows the borehole pressure development as a function of time for 2000 cycles. With increasing injection flow rate, the variation of borehole pressure is initially significant and becomes smaller through time. The big variation in pressure which is caused by the relatively low injection flow rate at the start is favourable in terms of With an increasing flow rate for cyclic injection, the element experiences a cyclic tensile stress for several cycles and then increases to a higher cyclic stress profile. The higher cyclic stress profile is caused by the increased flow rate and such an intensified stress will favour the fatigue generation of rock. The fatigue damage with the cyclic injection is cumulative up to the point of crack initiation which is then followed by progressive crack propagation. The fatigue crack initiation occurs at a residual strength of 7.2 MPa which is 90% of the static rock strength. Therefore, cyclic injection with an increasing flow rate is suitable for fatigue fracture propagation in rock while the pressure variation becomes small with the fracture propagation and increasing flow rate.

A New Strategy for Pulsating Hydraulic Fracturing
Based on the modelling results, a new three-stage pulsating hydraulic fracturing strategy is proposed, shown schematically in Fig. 18. First, a high frequency pressure-controlled pulse is used to generate rock fatigue crack initiation at a reduced breakdown pressure. With the model parameters that we use, this fatigue results in a 10-18% lower breakdown pressure than the case where only a monotonic pressure is applied. In a real situation, the magnitude of the reduction will be dependent on the rock properties, in-situ stress and pressure, and the presence of heterogeneities such as sedimentary structures or existing fracture networks. Using higher frequency pulses prior to fracture initiation can reduce the time to generating the first crack. Crack initiation would be marked by a significant and sudden fluid injection volume increase. Once the crack has initiated, Stage II represents a change from high-frequency pressure-controlled pulses to lower frequency flow rate-controlled cyclic injection using increasing flow rates. Again in a real situation, the optimal frequency for generating the most complex fracture network possible will depend on the rock mass properties: higher frequencies may favour the generation of more fractures local to the borehole.
Stage III is at a lower frequency to incur more significant pressure variation within the lengthening crack. The high frequency in Stage II leads to a shorter time for injection in one cycle and thus less pressure build-up to prevent fast crack propagation after crack initiation. As the crack enters the steady development stage (Stage III), a lower frequency will pump in more liquid in one cycle to drive the crack propagation further into the rock mass. With the flow rate increasing and fracture propagation creating a larger volume crack, the cyclic injection frequency should be reduced in Stage III. The appropriate point to change from high frequency to low-frequency injection would be determined by a reduction in the variation of borehole pressure.
This new strategy outlined for pulsating hydraulic fracturing demonstrates its application in hydraulic rock fracture with the inputs shown in Table 3. The flow rate for cyclic injection increases by 5.0 × 10 -6 m 3 /s after every 20 cycles. Figure 19 illustrates the borehole pressure development over the new injection strategy period. Compared with the cyclic injection approach as illustrated in Fig. 16, the new approach produces an improved borehole pressure history. First, the breakdown pressure was dropped to 11 MPa compared to 13 MPa and achieved in much shorter time. Second, the long tail from Stage 3 in Fig. 19 has been enhanced with larger stress fluctuation. By changing the cyclic flow injection frequency from 0.25 to 0.125 Hz, the borehole pressure variation becomes greater which will generate more significant fatigue fracture of rock.
The fracture length development and injection fluid volume for the new strategy are shown in Fig. 20. The initiation of rock fracture occurs only after 88 cycles of   Fig. 9. At the same time, the injection fluid volume gradually increases. Figure 21 illustrates the amount of rock fatigue on the volumetric basis along the crack. It is important to understand how much fatigue it has actually been generated for every material point that is cracked under the given hydraulic pulse. This can allow identifying the most efficient rock fracturing and permeability enhancing method. It can be seen that the fatigue crack strength of the rock is generally reduced to the range of 7.2-7.8 MPa. Pulse is to generate fatigue of rock in line with the soft simulation technology. It should be noted that the magnitude reduced varies by the combination effects of the increasing flow rate, dynamic crack propagation and varying frequency. The new strategy for pulsating hydraulic fracturing that we propose can produce a steady and slow rock fracture with a reduced breaking stress. Delivering the required pulse amplitude, frequency and water volumes to the wellhead will require novel pump hardware which is the subject of ongoing research. Due to the assumption of homogeneity and the predefined fracture path, the current model cannot simulate the generation of a complex fracture network. Future work using modelling approaches that explicitly include rock mass heterogeneity (c.f. Xi et al. 2018) will be carried out to simulate the fracture of heterogeneous rock under pulses. Mitigation of the risk of seismicity is a key challenge for hydraulic fracking. Previous studies (Ellsworth 2013;Zang et al. 2019;Zhuang et al. 2019) have found that fracking-induced seismicity is closely related to highpressure injection and energy release during rock crack, etc. Through the developed numerical method and proposed injection strategy, the rock is softened by tensile strength degradation, the breakdown pressure is reduced and a steady and slow rock fracture is produced. Therefore, the hypothetical strategy provides a new approach for fracturing rock while reducing the risk of seismicity.

Conclusions
A fatigue damage model has been developed to investigate rock fracturing under hydraulic pulses, by introducing coupling of rock S-N curves with the fluid cohesive crack model. The fatigue cohesive crack model was implemented into ABAQUS by in-house FORTRAN subroutines. The progressive fatigue damage accumulation and fatigue crack processes of rocks under pressure pulses and cyclic injection were accurately simulated, showing a good agreement with experimental results. Hydraulic pressure pulses can reduce the breakdown pressure of rock by 10-18% of the static breakdown pressure upon 10-4000 cycles. The effects of flow rate and frequency for cyclic injection were investigated and led to a new strategy for pulsating hydraulic fracturing through three stages. Stage I involved a high-frequency pressure pulse to induce fatigue and lower the breakdown pressure required for rock fracture. After the crack was created, high-frequency hydraulic pressure pulses could result in a runaway single crack; to control the rapid crack growth, the hydraulic injection should be changed to a lower frequency flow rate-controlled cyclic injection as in Stage II. As the crack was propagated, a higher volume of fluid was needed to slowly drive the crack and as such lower frequency with increasing injection rates was proposed. Under this new loading strategy and protocol, a slow and steady rock fracturing process was obtained and the fatigue crack stress was reduced from 8 to 7.2-7.8 MPa. The modelling results in this paper showed that the trade-off between fatigue cracking, crack extension, and maintaining pressure within the crack must be understood to optimize the engineering of effectively, high surface area fracture networks in the subsurface.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.