Hydraulic fracturing has been adopted globally in the resource sectors, for instance, to improve gas pre-drainage performance, promote hard roof falls, and alleviate abutment stress (Li et al. 2023a; Si et al. 2019; Xie et al. 2020; Yang et al. 2020). Nevertheless, achieving the expected stimulation performance at desired depths or in a fluid-pressure regions of a reservoir is challenging due to the geological complexity of rock formations. Rock in the field is known as a complex structure formed by long periods of geological activities (Li et al. 2023b; Peng 2023; Xiang et al. 2021). There are many randomly distributed discontinuous fracture planes of different sizes resulting in the deformable features of rock masses (Li et al. 2022). The properties of these distributed fracture planes, such as faults and joint sets, constitute the main part of the geological complex and make rock masses highly heterogeneous. They can directly dominate the rocks’ fracturing behaviours and cause a high level of permeability anisotropy (Liu 2005; Hou et al. 2023; Zhang 2013). Consequently, the fluid flow in the hydraulic fracturing operation becomes uncontrollable and cannot be guided directly to the target rocks as expected (Ren et al. 2015). To capture the geometric complexity of natural rocks, the discrete fracture network (DFN) has been widely used to simulate the fracture system in oil and gas reservoirs (Liang et al. 2019; Nejadi et al. 2017), which can also intuitively show the distribution of fractures in underground mining.

A discrete fracture network (DFN) can be generated by different methods including stochastic realisation, geological mapping and geomechanical simulation. It can characterize various sizes of rock fractures from joints, faults, to bedding planes. Geological mapping, which is an outcrop based DFN model, is often restricted to planar analysis and may not capture the geological complexities (Lei et al. 2017a, b). Therefore, the stochastic realisation combined with the field survey (normally 1D scanline survey) and geomechanical simulation becomes a better option. Since geological media typically exhibits discontinuous features under complicated boundary conditions, numerical simulation by the conventional continuum methods may not adequately capture the crucial mechanical behaviours of fractured rocks with fluid injection (Jing 2003), such as cracking of intact rocks (Hoek and Martin 2014), and interaction between multiple fracture walls among pre-existing fractures and hydraulic fractures (Pollard and Segall 1987). Thus, it is essential to consider discontinuum methods, for example, the distinct element method (DEM), to numerically solve a rock mass system with complex fracture geometries and fluid injection cycles.

Furthermore, both the DFN model and induced seismicity can be simulated by DEM, which has been widely applied to study the mechanical behaviour of fractured rock. Fluid injection and moment tensors using PFC has been reproduced to investigate interaction behaviours between hydraulic fractures and simple pre-existing fractures of different dip angles under various stress conditions (Zhao and Young 2011). Subsequently, more advanced numerical models have been developed using this software. The cyclic hydraulic fracturing operation in naturally fractured reservoirs by different injection controls were numerically reproduced by Yoon et al. (2014) to optimise the stimulation based on the observed borehole pressure and induced seismicity. Zhang et al. (2023) has also realized the PFC fluid injection controlled by different injection methods in coal sample with single pre-fracture and meanwhile provided some new insights into the correlation of dynamic aperture and stress with induced seismicity. These previous studies justified the utilities of DEM method to simulate induced seismicity and fluid injection on rock masses with complex fracture networks of varying scales.

While these fluid-geomechanically coupled models aim to investigate the impact of stress, injection rates, and fluid viscosity on the rock failure behaviours (Zhou et al. 2016), there have been limited efforts to further investigate how the geometries of fracture networks influence rock mechanical properties. Le Goc et al. (2014) conducted continuous research on the elastic properties of fractured rocks to examine the influence of fracture density and orientations under compressive loading using 3DEC simulations. However, it is not yet fully understood whether changes in elastic properties due to DFNs are related to hydraulic fracturing behaviours. Additionally, application of hydraulic fracturing in coal seams is getting increasingly popular as injection is becoming a necessary approach to stimulate gas drainage production or manage hard roofs (Fang, et al. 2023). Coal is fragile and heterogeneous, with extensive natural fracture networks (also known as cleats), which dominates its failure behaviour. However, the monitoring and evaluation of hydraulic fracturing operations in coal are still at an initial stage (Kang et al. 2023). The injection induced seismicity in coal seams remains largely unknown. There are limited numerical studies on hydraulic fracturing of naturally fractured coal samples, particularly with the application of DEM models.

This research aims to investigate the influence of fracture densities and orientations on hydraulic fracturing behaviours, induced seismic responses, and relative modulus changes during fluid injection in laboratory-scale coal samples using particle flow code (PFC). The evolution of computed acoustic emission (AE) events in coal samples is analysed with the installation of DFN systems governed by the power law distribution. The results are then used to determine the correlations between hydraulic fracturing behaviours and the elastic property of coal samples.

2.1 Algorithms of fluid injection in particle flow code (PFC)

Particle flow code (PFC) 2D based on distinct element method employs various contact models to design and bond individual particles which have an overall strength as required (Itasca 2008). To simulate fluid flow in PFC, three crucial steps need to be followed. Firstly, as shown in Fig. 1a, bonded particles and contacts of any contact models such as parallel bond or smooth joint contact model, can construct the domain skeleton with flow paths. Impermeable solid particles represent rock grains, voids enclosed by particles represent pores as shown in Fig. 1b, and broken bonds between contacts represent cracks that act as flow channels. Secondly, the equations of fluid flow and fluid pressure update should be defined using Darcy’s law. In the final step, the fluid pressure is coupled with the total stress surrounding the particles and contacts. Afterwards, pore scale fluid transport, which is forced by the pressure difference between neighbouring domains seen in Fig. 1b, can be simulated. The flow rate, defined by Eq. (1), is in proportion to the cubic of the aperture of flow channels and pressure variations:

where e denotes the hydraulic aperture; Δ \({P}_{\text{f}}\),the differential pressure between two neighbouring domains; η, the dynamic viscosity of injection fluid; and L, the length of a flow channel. The hydraulic aperture in a flow channel is dependent on the contact stress among particles according to Eq. (2).

where \({e}_{\text{inf}}\) is called the infinite hydraulic aperture with infinitely large normal force; e₀, the initial hydraulic aperture without any normal force applied; \({\sigma }_{\text{n}}\), the effective normal force on a bonded contact. Finally, the flow rate from Eq. (1) gives the temporal fluid volume variation throughout all channels, which is later plug into Eq. (3) to obtain the fluid pressure variation (ΔP) per computational time interval (Δt) in a fluid domain.

where K_f denotes fluid bulk modulus; V, the volume of individual flow domain; Q, flow rate; and ΔV, the real-time apparent volume change of a flow domain.

a Flow domain network (after Itasca, 2008): bounded particles (grey), the flow domain centroid (blue) and flow channel (black); b Individual fluid flow area (blue) outlined by grey particles

The update equations are only effective for unbroken particle bonds. When a crack is induced, the flow path joining the adjacent domains will be open fully and the two domains will be merged into one domain. The new fluid pressure will be the average of values of the two domains in the previous time step. Thereafter the updated fluid pressure will be added onto the elementary particles. All the channels will be detected and audited for individual elementary particles in every simulation timestep. At the same time, Eqs. (1)–(3) will refresh the fluid pressure for each fluid domain. All the channels will be detected and audited for elementary particles in every simulation timestep. At the same time, Eqs. (1)–(3) will refresh the fluid pressure for each fluid domain.

2.2 Algorithms of computing seismic moment tensor

The interpretation of rock failure mechanisms is widely facilitated using advanced seismic moment tensor analysis, as described by Aki and Richards (1980) and Wang et al. (2021). By examining the geometry of the resulting crack and its force–displacement relationship, the moment tensor of the crack failure can be constructed, thereby enabling the simulation of AE events. For this study, the simulation of AE events involves the recording of contact’s failure and its breakage information for every step once domain pressure is updated (seen Fig. 2). The moment tensor (\({M}_{ij}\)) is then given as follow:

where, ΔF_i is the i-th component of the space vector denoting force change, and \({R}_{j}\) is the j-th component of the displacement vector of moved particles. Then the product will be summed over the surface of the seismic event, S. The calculated seismic moment is a 2 × 2 symmetric matrix of the recorded crack (bond failure) in 2D model for this research. Moreover, the moment tensor can also describe the release energy of an AE event from that crack. One of empirical methods to calculate the magnitude of AE events (\({M}_{\text{w}}\)) is employed by Hanks and Kanamori (1979).

where \({M}_{0}\) is the scalar moment given by:

where \({m}_{j}\) is the j-th eigenvalue of the moment tensor.

Modelling procedure of fluid injection and seismicity

2.3 Design of discrete fracture networks

As previously mentioned, the challenge of obtaining a comprehensive estimate of natural fracture systems has resulted in the development and widespread use of stochastic approaches (Dershowitz and Einstein 1988; Wang and Cai 2020). The stochastic method to generate DFN evolved with the aim to investigate the percolation among a set of fracture populations (Balberg and Binenbaum 1983; Robinson 1984) and meanwhile the fluid flow through complicated fracture networks (Long et al. 1985). Assumed by the general stochastic DFN generation methods, fractures are usually lines in 2D planar modelling or round discs in 3D spatial modelling. The other geometrical parameters such as positions, sizes, orientations are treated as independent variables that follow certain probability distributions. They are usually obtained from in situ measurements, including outcrop traces scanline and borehole imaging (Zhang and Einstein 2000). The geometry of a DFN model is characterised solely by independent statistical distributions, namely the types of distribution and constrained parameters, of its geometrical properties, which are supported by the basic DFN generation functions are fracture size (length in 2D or diameter in 3D), orientation, and position distribution. A range of DFN realisations, depending on a random seed, can be acquired from a simple stochastic DFN model. In this model, fracture density and orientation are independently controlled as separate variables.

The fracture size follows a density distribution function, which specifies the fracture numbers in unit volume, n(l), with size range from [l, l + dl]. It is determined by a probability distribution function governing fracture size. The power law distribution function, which is merged by fracture maps with various scales and mechanical fracture growth models, has been proved to be a good candidate to characterise various sizes of discontinuities in rocks (Bonnet et al. 2001; Davy et al. 2018). For a power law distribution, the relationship is (Lei et al. 2017a, b):

where \(a\) is the scaling exponent and \(\alpha\) is the density term of the DFN system. The scaling exponent \(a\) is positive because the power law exponent would be always negative. The power law distribution is bounded by l_min and l_max, which represent the lower and upper limits of the fracture sizes, respectively. The value of \(a\) determines the ratio between small and large fracture sizes at any scale. The power law function includes both end member models, and infinite number models. The DFN will become constant size model with l_min when \(a\to \infty\) while the fracture number will become infinite when \(a\le 2\). With \(a\) increasing, the proportion of smaller fractures in this DFN system also increases.

The term \(\alpha\) determines the total fracture density based on the range of fracture sizes. The total fracture density in a volume of characteristic size L is also dependent on the range of fracture sizes. The fracture number with sizes within the range [l, l + dl] is denoted by N (l, L) (e.g., \(N\left(l,L\right)=\alpha \bullet {l}^{-a}\bullet {L}^{3}\)). The number of fractures with sizes between \({l}_{1}\) and \({l}_{2}\) is given by:

The distribution of cumulative fracture size density, which determines the number of fractures with the size larger than a certain value, is defined by:

Equation (9) is used to clarify that if the value of scaling exponent becomes large enough (e.g., ≥ 4), the yielding fracture population is determined by the smallest fractures because the term \({l}^{1-a}\) decreases quickly if \(a\) ≥ 4. In this model, the length range is restricted from 2 to 10 percent of the sample length, and \(a\) is defined within the range [1, 3] to represent real rocks (Maillot et al. 2016).

Fracture density is a scale-dependent property, whose distribution depends on the location and dimensions of the measurement area. It features local variability when the measurement area tends to be zero or define macroscopic properties if the measurement dimension becomes infinite. One of the most popular ways to define DFN density is referred to the cumulative surface area of fractures within unit volume. In this 2D model, it is the cumulative fracture length per unit area (P₂₁) which is employed to calculate the fracture density based on Eq. (9).

Integrating all the simulation algorithms discussed in Sects. 2.1, 2.1 and 2.3, the workflow in Fig. 2 illustrates the iteration process of the self-developed DFN fluid injection module in PFC2D. Initially, a simulated rock sample is calibrated and then the discrete fracture network with designed properties is integrated with the sample. Following this, a cubic sample with a flow channel structure is built based on the input values derived from a calibration process. The next step is to apply boundary conditions before iteration onset of the fluid flow algorithm and application of the tectonic stress. Upon initiating the injection process, the equations detailed in Sect. 2.1 are utilised to calculate the pore pressure, which is subsequently combined with principal boundary stresses and particle bonded force so as to determine the final state of effective stress (Potyondy and Cundall 2004). For each timestep, the AE computation module is executed to record bond contact failures and compute related seismic moment tensors. Notably, the algorithms governing fluid injection and AE monitoring are independent of each other without interference until the injection is terminated.

3.1 Input parameters calibration

To ensure desirable mechanical responses that align with lab observations, it is crucial to carefully determine certain essential input parameters. Previous lab test results can work as calibration data to parameterise the model. Empirical equations between macro mechanical properties of natural rock and values of input micro parameters in our modelling have been established by repeating lab tests of rock’s uniaxial compression strength (UCS) (Zhang et al. 2023). They can provide a preliminary range of input values before fine-tuning parameters are calibrated until the simulation results are consistent enough with lab tests.

Table 1 presents the proposed parameters for calibration. To validate the results, UCS test results from both simulations and laboratory experiments are compared in Fig. 3. The dashed rectangle highlights a strong correlation in the elasticity stage, as evidenced by the parallel curves. The red curve illustrates evident nonlinear features between stress and strain, which are prominent in the fragile and brittle nature of fractured sedimentary formations. Table 2 provides the calibrated strength values along with the average values from lab samples.

Comparison between calibrated UCS simulation and lab results

3.2 Assignment of discrete fracture network (DFN) and fluid injection

As previously mentioned, in situ rock masses are complicated systems where rock matrix may have much heterogeneity and are generally filled with different kinds of discontinuities such as cracks, joints, and faults. There may occur the failure in the matrix and discontinuities, or at both locations. To study the mechanical behaviour of these systems, the smooth joint contact model was developed in PFC and has proven effective (Mars 2010). When added to the inherent running logic of PFC the modelling approach can superimpose cracking with its geometries and mechanical properties onto a bonded-particle model. Due to the inherently discrete nature, failure can be deduced in not only intact bonded regions but also potentially along fracture planes. In this model, the parallel bond contact model, which is good at replicating many features of rock behaviours (Potyondy 2011), is used to model intact rock, while the smooth joint contact model is introduced to simulate the mechanical behaviour of pre-existing fractures at contacts intercepting fractures.

Based on the fracture size distribution function discussed in Sect. 2.3, different modes of the fracture network can be generated by designing various fracture densities and orientations (dip angles) as summarised in Table 3. DFN is only assigned within the flaw area outlined by a red square in Fig. 4 and the tensile strength and the cohesion, determined by groups of trials to have reasonable fracture length, are set to be nearly half value of the intact part of the sample (seen in Tables 1 and 4). This design aims to prevent the quick breakage of the entire block sample due to the weakening effect of fractures. For the same purpose, the DFN density range is determined from 16.22 m⁻¹ to 75.89 m⁻¹. If the fracture density is larger than 75.89 m⁻¹, the sample, calibrated based on Table 1, will start cracking and collapse shortly after the onset of injection. This is also observed if the fracture dip angle exceeds 60 degrees. Meanwhile, the designed fracture density range can ensure that the fracture numbers increase evenly, roughly from 10 to 50, to achieve better comparison between these modelling results. In this model, the total injection duration is determined based on the following criteria:1. Hydraulic fractures has crossed the boundary of flaw area; 2. Injection is terminated at the time when the sample is fully cracked but the sample integrity is maintained for easy observation.

Geometry of the simulated block coal sample

Taking the sample G1D1 in Table 3 as an example, the sample constructed with bonded particles and flow channels is shown in Fig. 5. The sample consists of 15296 bonded particles and 40277 active contacts. Black dots represent the bonded particles in the fractures, and particles in the rest part are marked in grey. Blue meshes depict the flow channel of individual domain structures is shown in Fig. 5.

Simulated block coal sample with bonded particles and flow channels in PFC

In our modelling, water is the injected fluid into the centroid borehole at a fixed flow rate. In our models, water is the injected fluid into the centroid borehole at a fixed flow rate. Although it is larger than the grain size in natural rock samples, the simulated particle size has met the accuracy requirements defined by the index of dimension resolution (RES) (Deisman et al. 2008). It is not appropriate to directly compare the injection flow rate in numerical models to lab tests due to the size difference. Additionally, a quasi-static state, which is usually affected by computation time interval, is vital to fracture’s smooth initiation and reasonable propagation (Tomac and Gutierrez 2020). Notably, the initial pore pressure would also influence rock’s failure by altering the effective stress (Terzaghi and Peck 1948), whereas it is not the focus of this model. Thereafter, there is no initial pore pressure as done in earlier models (Al-Busaidi et al. 2005; Zhang et al. 2023). By repeated trials and errors, some important input parameters for fluid injection with DFNs are ultimately determined in Table 5.

4.1 Distributed AE events and recorded borehole injection pressure

In this research, a total of 25 samples were simulated as listed in Table 3. One of these samples, G1D5, is presented here as an example to illustrate the distribution of AE events and the evolution of borehole injection pressure. With the most DFNs (and the largest fracture density), G1D5 is expected to display a more distinct interaction between hydraulic fractures and the pre-existing fractures than the other samples. To enhance visualisation, only the flaw area (X = -30 mm ~ + 30 mm) with the define DFNs in the sample is depicted in Fig. 6, and no AE events were recorded outside of this region. The colour-coded seismic events with time sequence in Fig. 6a are plot on together with the DFNs. The event size is normalised according to the energy magnitudes from −6.8 to −5.57 calculated by Eq. (5).

a The spatial distribution of induced AE events by hydraulic fracturing; b The spatial distribution of fluid pressure and AE events at the final stage of injection; c The displacement of particles after injection

The seismic events indicate a clear impact of pre-existing fractures on the propagation of hydraulic fractures, particularly for fractures propagating downward from the sample centre. This results in a horizontal deviation trend of hydraulic fractures, although they still extend along the maximum principal stress once departing from the front of pre-existing fractures. Moreover, most of the strong seismic events, labelled as Event 1 to Event 6 in Fig. 6a, are located near the pre-existing fractures due to higher level of stress concentration near the tips of pre-existing fractures (Zhang et al. 2023). The correlation between seismic activities and pre-existing fractures is supported by the fluid pressure distribution observed at the completion of the injection, as shown in Fig. 6b. In the areas with strong AE events marked in Fig. 6a, the corresponding fluid pressure is relatively higher either in the borehole or in the pre-existing fractures. The larger fluid pressure present in pre-existing fractures may cause higher energy release for stronger seismic events. However, despite that the fractures’ propagation generally follows the trace of fluid flow, there are some off-track points marked as fluid induced AE events in Fig. 6b. Although fluid cannot access these areas to directly trigger AE events, the stress perturbation caused by fluid injection may result in AE events in or adjacent to pre-existing fractures. Conversely, there are also some pre-existing fractures containing pressurised fluid flow but without hydraulic fractures indicated by AE events, as shown by the bottom of DFNs, which can be noted as well in Fig. 6b.

In this DEM model, different behaviours of discrete fracture networks can be simulated as representative occurrences in natural fractures in the field. They not only can convince the validation of injection simulation but also are very meaningful for future research on the failure mechanism of fractured reservoirs under fluid injection. Firstly, pre-existing fractures can alter the dominant propagation trajectory of hydraulic fractures. Secondly, they can be induced indirectly to crack but not being exposed to fluid flow. Thirdly, the DFNs can also serve solely as channels for pressurised fluid but without being reactivated.

To examine the effect of the distributed fluid on the sample’s deformation, the displacements of individual particles are displayed in Fig. 6c. Overall, the displacement is distributed symmetrically relative to the trace of hydraulic fractures, with the central part near the borehole exhibiting the largest displacement, which diminishes further away from the borehole. Specifically, the area with pressurised fluid creates a distributed zone as indicated by the displacement pattern. Above this zone, the displacement exceeds 0.2 mm, while below it, the displacement abruptly decreases to less than 0.15 mm.

In addition to the final state of fluid injection, the temporal borehole injection pressure of sample G1D5 is recorded with the travelling distance of AE events away from the borehole centre, as shown in Fig. 7a. During the first half stage of injection, strong events occur after dramatic jumps in the borehole injection pressure, for example at 5.5 min and 12 min. During the second half stage, more strong events are produced, but their occurrence has little correspondence with the pressure change after the flow has travelled far enough away from the borehole centre. Theoretically, the typical breakdown pressure (\({P}_{\text{wf}}\)) for a tension crack, as expressed by the following equation (Eq. (10)) proposed by Hubbert and Willis (1957), should be 6 MPa based on the input parameters in Table 1, assuming the initial pore pressure (\({P}_{0}\)) is 0 MPa. After the theoretical breakdown pressure, the hydraulic fracture should propagate at constant pressure until the sample is fully cracked.

whereas in this model, the fracture initiation pressure or breakdown pressure (P_B) is defined as the corresponding local pressure peak at the timepoint of AE events that were first observed in Fig. 6a, and after that, it is defined as the stage of the fracture propagation. The resulting breakdown pressure for the simulated sample is 3.77 MPa, which is less than the theoretically calculated breakdown pressure of 6 MPa. This difference could be due to the assumption of homogeneity and isotropy in the theoretical calculation, whereas in reality, the matrix near the borehole is weakened due to borehole drilling and DFN installation, resulting in an altered local stress distribution.

a The variation of borehole injection pressure against the onset time and travelling distance of AE events; The fluid pressure distributions at the time of b Breakdown pressure (P_B) and c peak pressure (P_P)

Furthermore, the observed borehole injection pressure curve is not linearly increasing before arriving at the breakdown pressure, unlike the idealised borehole pressure curve during fluid injection. This is due to the setup limitation of DEM models as reported before (Zhou et al. 2017), where the bonded particles, constructed domains and flow channels are larger than those existing in the natural rocks. When fluid is pumped into the borehole, according to the cubic law, more leakage can enter the neighbouring domains instantaneously through the flow channels, leading to fluctuated injection pressure near the borehole, which can also explain the non-linear propagation pressure after the breakdown pressure. As the cracks propagate further away from the borehole centre, particles in the sample becomes more homogeneous and better compacted so that more pressure is required for cracking, which may cause a larger global peak injection pressure (P_P = 4.61 MPa) compared to the initial breakdown pressure to fracture the borehole.

Accordingly, the contour plots of relevant fluid pressure distribution at P_B and P_P are illustrated in Figs. 7b and c, respectively. It can be observed that the first crack is initiated on the top of the borehole at P_B, which is reasonably considered as the breakdown point. Moreover, the fluid in Fig. 7b deviates bifurcately on the downside propagation to avoid crossing the pre-existing fracture, instead of following the theoretical propagation direction along the maximum principal stress. The branching flow confirms the influence of the pre-existing fracture on fluid distribution and crack initiation. In Fig. 7c, the isolated AE events (cracks) are indirectly induced during the period of 4.56 min to 7.96 min, which is earlier than the occurrence of the peak pressure (T = 8.87 min) and no more isolated cracks are indirectly induced by the fluid flow thereafter. This phenomenon applies to all the other 24 samples, which are not shown here to avoid repetition. In short, the breakdown pressure and peak pressure play different roles in this hydraulic fracturing simulation. The breakdown pressure indicates the onset of hydraulic fractures, while the peak pressure indicates the termination of indirectly induced cracks by fluid injection. They will be discussed separately in the following sections to investigate how they are impacted by DFNs.

4.2 Injection measurements and induced seismicity influenced by DFNs

To further quantify the differences among various samples, the peak injection pressure from all modelling groups is collected and shown in Fig. 8a. The peak injection pressure ranges from 3.6 to 5.1 MPa and decreases as fracture density increases. While fracture orientation does not have a significant impact on peak injection pressure, there is a slight difference between fracture dip angles less than or equal to 30 degrees and those greater than 30 degrees. For angles ≤ 30 degrees, the peak pressure initially decreases and then increases when the density is larger than 59.48 m⁻¹. For angles > 30 degrees, the peak pressure declines to a certain value and then remains unchanged when the density is larger than 59.48 m⁻¹.

The comparison of injection behaviour in different samples. a The peak value of borehole injection pressure; b The duration of the fluid injection process

Besides the peak injection pressure, the injection duration is also compared to see if it can describe the injection behaviours under different fracture densities. For each sample under fluid injection, it will eventually be fully cracked with enough time. The injection durations of all the groups have a wide variety ranging from 20 to 41.4 min, as seen in Fig. 8b. For a certain fracture orientation, the overall trend demonstrates that the injection duration will become shorter with the increase of fracture density but experiences a mild increase when the density exceeds 59.48 m⁻¹. This may be due to the combined influence of hydraulic fractures and natural fractures: the previous research on the stress shadow suggested that there is an optimal value between low density and high density DFNs, which can achieve the easiest crack propagation (Huang et al. 2023). For this point, a denser DFN may contribute to the production of complex fracture networks which are beneficial to permeability enhancement but do not necessarily speed up the propagation of hydraulic fractures and the injection completion. Moreover, at a given fracture density, the three upper dashed curves indicate that increasing fracture dip angles will reduce the injection time when the dip angles are in the range of 0 to 30 degrees.

The distribution of strong seismic events numbered in Fig. 6a provides evidence for the impact of DFNs on fluid-injection induced seismicity. This section also aims to establish the statistical correlation between seismic events and DFNs by analysing the ratio of strong AE events to the total number of AE events in each sample. Since all the samples share the same range of seismic magnitude from -6.8 to -5.57, the number of top 10% to 50% largest magnitude events is counted for every single sample trying to find a number to distinguish different samples. Consequently, the number of events with the top 25% largest magnitude are collected as strong events from every individual sample to calculate the proportion of strong events. Unlike the magnitude of the individual seismic event, the proportion of strong events is for comparison purposes among different samples to estimate the magnitude of the induced seismicity throughout the whole injection process. As seen in Fig. 9, the curves of strong event proportions generally follow a bell-shaped pattern. As the fracture density increases, the proportion of strong events increases and then decreases for a given angle, except for 30 degrees. At densities of 29.31 and 43.67 m⁻¹, the proportion reaches its peak value. For the angle of 30 degrees, the proportion trend is unclear, fluctuating from 0.02 to 0.17. When the angle is larger than 30 degrees, the peak proportions marked by blue dots are 0.25 for 45 degrees and 0.24 for 60 degrees, which are larger than the peak values of 0.11 for 0 degrees and 0.19 for 15 degrees.

Proportion of strong AE events during fluid injection

The comparison shows that the fracture density has a more significant impact on the peak injection pressure, injection duration and seismicity than fracture dips. Specifically, the peak injection pressure is negatively correlated with the fracture density. Moreover, denser DFNs can accelerate the hydraulic fracturing process by reducing the injection duration. Moreover, there are distinct effects observed between different dip angles, which will be discussed further in later sections.

4.3 Modulus influenced by DFN dip angles and rock’s friction angle

Section 4.2 explored the influence of fracture density on hydraulic fracturing, while this section focuses on fracture dip angle. By integrating key values that feature the initiation of hydraulic fractures, this section investigates how fracture orientations affect hydraulic fracturing behaviours from the perspective of elastic property, which is important to characterise a specific rock. According to the theory of linear elastic fracture mechanics (LEFM), the initiation of hydraulic fractures is assumed to follow Griffith’s law of energy balance at the breakdown pressure (Zeng and Roegiers 2002). Under plane strain conditions in a 2D model, along the single wing length of the fracture (c), the injection volume of a ‘line crack’ can be calculated from the fluid pressure p₀ (Valko and Economides 1995):

where L is the length of injection borehole in 3D, but it is assumed to be 1 in the 2D model. \({E}{\prime}\) is the plane strain modulus, which can be calculated theoretically by Young’s modulus (E) and the Poisson’s ratio (\(\upsilon\)):

Assuming no fluid penetration, the injection volume and the fracture volume are conservative:

Specifically in Eq. (13), t_p is called the initiation time when the fracture is initiated, q is the constant injection rate and p₀ is the borehole pressure (i.e., the breakdown pressure determined in Sect. 4.1) at time t_p. The length of single wing fracture (c) can be estimated using the seismic events at time t_p so the length can also be called breakdown fracture length. Therefore, the plane strain modulus can be expressed by the parameters obtained from the borehole injection data:

Since Eq. (14) is simplified for the 2D model, the calculated value of modulus may not be directly comparable to the real value obtained from UCS tests in the lab. However, it is still of great significance to compare the relative elastic properties between different samples with varying DFNs.

From hydraulic fracturing parameters, the relative modulus (E_rel) between one sample and another can be easily compared using Eq. (15). Following the same approach as in Sect. 4.1, G1D5 is taken as a reference sample. The plane strain modulus of other samples is divided by that of the sample G1D5 to calculate the relative plane strain modulus, which is referred to as relative modulus in the subsequent discussion.

The results of relative modulus with respect to fracture density are presented in Fig. 9a. The relative modulus shows an overall declining trend with an increase in fracture density when dip angles ≤ 30 degrees. In addition, the horizontal DFNs (angles = 0 degrees) exhibit the highest relative modulus among all angles at densities of 16.22 and 29.31 m⁻¹. In contrast, there is an opposite trend for fracture dip angles of 45 and 60 degrees. With an increase in fracture density, the relative modulus mildly increases. When the dip angle is larger than 30 degrees, the overall values of relative modulus become smaller than those of smaller angles. A similar trend for all angles (except for the horizontal DFNs) is that the values tend to remain unchanged when the fracture density is larger than 43.67 m⁻¹. There seems to be a critical angle between 30 and 45 degrees. To determine this critical value more accurately, supplementary modelling groups are conducted by adding dip angles at 35 degrees and 40 degrees. The supplementary groups shown in Table 6 have the same input parameters as other groups in Table 3, except for DFN orientations.

In Fig. 10a, the results from the supplemented groups are highlighted by the red curves, revealing a distinct transition in the relative modulus curves between 30 and 35 degrees. Another observation is that, for the seven samples with the smallest density (P₂₁ = 16.22 m⁻¹), the relative modulus is negatively proportional to the increase of fracture dip angles.

The calculated relative plane strain modulus. a All modelling cases in Tables 3 and 6; b Fracture dip angle = 35 degrees but with varying friction angles; c Fracture dip angle = 30 degrees and 40 degrees but with varying friction angles

To clarify if the distinguished change between 30 and 35 degrees of dip angles is related to the friction angle (ϕ) of the rock matrix, more groups of modelling with different friction angles are also conducted. Figure 10b shows the modelling results when the fracture’s dip angle is 35 degrees but with varying friction angles. Marked by blue dots and red diamonds, when the friction angle is 37 degrees, the relative modulus presents an overall increasing trend from the value of 1.057 to 1.111, except at the fracture density of 43.67 m⁻¹ (E_rel = 0.919). This trend is identified with the result highlighted by the red line with hollow diamonds in Fig. 10a where the friction angle is 33 degrees. Once the friction angle reaches 45 degrees, the trend becomes unclear and fluctuates with increasing density. When the friction angle continues to increase up to 50 degrees, the curve shifts to an overall decline trend. These results suggest that the trend of relative modulus is influenced by the friction angles. A distinguished change between 30 and 35 degrees can be justified only in conditions when the friction angle is no larger than 37 degrees. When the friction increases, the distinguished change may happen in other ranges of fracture dip angles. Further, when the friction angle is 45 degrees, the trend becomes fluctuated and is unpredictable. Therefore, in this case, the dip angle of 35 degrees sets a transition line for the relative modulus evolution.

To identify the transition line when the friction angle is 50 degrees, more models with the fracture dip angles at 30 and 40 degrees are simulated, whose results are shown in Fig. 10c. Similarly, the modulus evolution presents a fluctuation curve suggesting that 40 degrees can become a new critical value for trend transition when the friction angle is 50 degrees. No more friction angles larger than 50 degrees would be tested since this number is a reasonable upper limit of the friction angle for most natural rocks. On the other hand, a friction angle of 20 degrees, which is usually the lower limit of friction angle among natural rocks, is tested to check whether it can change the trend of the relative modulus. As Fig. 10c shows, the trend for the relative modulus of 30 degrees dip angle (Dip = 30^°) remains uncurbed with the change of DFN density. No further modelling with DFN dip angles less than 30 degrees will be tested since the results will remain the same trend according to the above findings.

Table 7 summarises the transition fracture (DFN) dip angles under different friction angles. The transition dip angle increases approximately with the increase of the sample’s friction angle. Although the resultant series of values cannot give a quantitative relationship between the samples’ friction angle and the transition dip angles, it is still useful to provide a range of transition DFN dip angles with its reflection on the friction angle. This information can aid in characterising the in-situ rock properties when combined with the abovementioned hydraulic fracturing parameters (fracture length, breakdown pressure and corresponding occurrence time). The characterisation given by Eq. (15) may not be as precise as that from conventional lab UCS tests, but it is not restricted to the location and time and can be easily obtained from seismic and injection data. This convenience privileges the relative modulus with a lot of practical meaning in in-situ scenarios.

In a word, the analysis conducted in Sects. 4.2 and 4.3 confirms the interconnectedness of hydraulic fracturing measurements, induced seismicity, discrete fracture networks, and the elastic properties of rock. Consequently, a flowchart illustrated in Fig. 11 is proposed to estimate the desired DFN density and geometry, as well as rock elastic properties, based on observed hydraulic fracturing results. First, the fracture density can be determined from injection duration, peak pressure, and seismic events. Subsequently, using the breakdown pressure, initiation time, and breakdown fracture length obtained from seismic events, the relative modulus can be calculated and the range of DFN orientation (dip angles) can be determined accordingly.

The inferring procedure of DFN properties by hydraulic fracturing measurements

Discrete fracture network (DFN) is essential in controlling the rock properties and fractures’ behaviours during fluid injection. Obeying the power law distribution function, fracture densities and sizes are designed to represent in-situ geometries of fractured rocks to simulate hydraulic fracturing and seismic/AE activities using the self-coded algorithm in PFC.

The modelling is effective in simulating three basic behaviours of DFNs that occur usually in naturally fractured rocks. The presence of DFNs results in lower breakdown pressure than the theoretical value and leads to a peak injection pressure higher than the breakdown pressure. The peak injection pressure and injection duration are influenced more obviously by the DFN density than by DFN dip angles. They would decrease comprehensively with the increase in DFN density. While the proportion of strong events, calculated from the induced seismicity, exhibits a bell-shape evolution with the increasing DFN density.

The relative modulus has been derived innovatively from the breakdown pressure, breakdown fracture length, and the related initiation time based on Griffith’s law of energy balance. It decreases with the increasing DFN density in larger dip angles (35–60 degrees) while would increase in smaller dip angles (0–30 degrees). The transition dip angle distinguishing the evolution trend of relative modulus will increase with an increase in the friction angle. These relationships hardly drawn from previous simulation research will contribute to inferring the distribution of discrete fracture networks in terms of density and dips by observed hydraulic fracturing measurements and induced seismicity as well as to obtain the elastic properties of rocks.