Microseismicity-Based Modelling of Induced Fracture Networks in Unconventional Reservoirs

Abstract

A single planar hydraulic fracture is typically the primary component used to simulate hydraulic fracturing stimulation in conventional reservoirs. However, in ultra-low-permeability shale reservoirs, a large system of fracture networks must be generated to produce hydrocarbons economically. Therefore, traditional modeling approaches centered on single planar fractures are inadequate for accurately representing the intricate geometry and behavior of fractures in these reservoirs. In previous works, 2D fractal fracture networks (FFNs) have been used to generate sets of hydraulic and natural fractures based on microseismic event (MSE) data. Since microseismic data are retrieved in 3D space, the aforementioned model cannot accurately represent induced fracture properties. The objective of this paper is to study in detail the recently developed 2D FFN model and propose a novel solution by expanding the previous model to accommodate real 3D microseismic data. First, the definitions of the 2D FFN model are described, and associated calibration mechanisms are proposed. Next, the 3D FFN model and its calibration system are demonstrated. While the novel 3D calibration solution utilizes an identical matching concept to the 2D methodology, the residual distances between the nodes and the MSE are calculated in 3D spaces. Finally, a set of real microseismic data are used to calibrate the generation of 3D fractals using the proposed workflow. The interactions between microseismicity and fractured reservoir dynamics are represented through the integration of fractal fracture models and microseismic data. This work contributes to advancing the current understanding of hydraulic fracturing in unconventional reservoirs and provides a valuable framework for improving fracture modeling’s accuracy in reservoir engineering applications.

1. Introduction

In the past, organic-rich rocks with low permeability, such as shale, were classified by the petroleum industry as purely hydrocarbon source rocks. This classification stemmed from their uneconomical use in production. Hydraulic fracturing made it possible to economically produce oil and gas directly from these low-permeability source rocks. This relies on extensive horizontal wells, spanning from 3000 to 10,000 feet, combined with multistage hydraulic fracture strategies [1]. During hydraulic fracturing, highly pressurized fluids and proppants are injected into the source rock. This process continues until the pressure gradient exceeds the threshold at which the rock first fractures [2]. After the rock ruptures, tensile and shear failures propagate throughout the elastic matrix, which creates a permeable hydraulic fracture network [3,4]. This resulting hydraulic fracture network permits and enhances the flow of hydrocarbons from the reservoir through fractures in the rock that were not present prior to stimulation [5]. Some of these mechanical failures in rock induce MSEs, which are small-scale events with a low magnitude, typically less than zero [6].

The elastic rebound theory states that stress and strain around a fault accumulate gradually due to deformation, and are then released suddenly during an earthquake by relative motion along the fault [7]. Within a given stratum, stress is anisotropic and shear stress is concentrated along the naturally fractured surfaces. Before completion, the anisotropic ground stresses generated by natural fractures maintain an equilibrium state, stabilizing the fracture surfaces [8]. During hydraulic fracturing, the original stresses are disturbed. Stress concentrations expand into areas of existing fractures, generated fractures, and activated fractures, increasing strain energy. When stress increases to a certain degree, deformations such as fault slip occur and stresses relax. The associated energy is released in the form of elastic waves and induced microseismic events, as shown in Figure 1.

Microseismic monitoring has been utilized in previous studies to analyze and predict the general shape and spatial extent of the induced fracture network [5,9,10,11]. Additionally, numerical methods have been developed to predict and evaluate the relationship between hydraulic fractures and observed microseismicity. In the Bossier Sands of Texas, the distribution of microseismicity was found to depend on the regional stress state as well as the heterogeneity of the elastic media [3]. Urbancic et al. collected and studied 994 MSEs with 11 triaxial geophones during three separate phases of the hydraulic treatment plan [12]. Time-differentiated seismic waves (longitudinal and transverse waves) recorded by different sensors were utilized to determine the location, magnitude, and timing of microseismic events. The findings revealed that the occurrence of microseismic events can provide insight into the geometry of the hydraulic fracture. Moreover, the parameters relating to frequency and magnitude can elucidate differences in fracture propagation and failure mechanisms in the rock. Generally, the rate of microseismic events per unit of time varies with the stress state, the fracture distribution, the strength of the rock, and the geometric shape of the reservoir. It is important to note that microseismic events induced by fluid injections may not all be related to fractures that are connected with each other.

In most commercialized fractured reservoir modeling, models consisting of a single planar hydraulic fracture have been the primary models utilized to analyze stimulation performance [13]. These models assume a fixed fracture geometry and utilize a symmetric partially open fracture with a hyperbolic crack law to replace the tortuous fracture [11,14]. This simplification allows for the approximation of fracture properties using a modified Reynolds flow law. For ultra-low-permeability reservoirs, a large system of fracture networks must be generated by hydraulic fracturing to economically produce hydrocarbons [15]. Furthermore, in enhanced geothermal systems (EGSs) for geothermal energy production, fractured reservoirs must contain an interconnected system of predominantly shear fractures [16]. Thus, the simplification of a single, planar hydraulic fracture model used in industry software is inadequate for evaluating its performance in most unconventional and geothermal reservoirs.

While current 2D models are useful, they fall short in terms of accurately representing the nature of fracture networks due to their inability to capture the spatial complexity and interactions between fractures in three-dimensional space. This limitation results in less reliable predictions of fracture behavior and, ultimately, suboptimal hydraulic fracturing strategies. To address this gap, this study proposes an extension of the FFN model that incorporates 3D microseismic data to offer a more precise and detailed representation of fracture networks in unconventional reservoirs.

The proposed 3D fractal fracture network (FFN) model is an enhancement of the previously examined 2D model and is specifically designed to incorporate the complexities of 3D space. This advancement includes refined calibration mechanisms that extend the capabilities of the 2D FFN model to accommodate 3D microseismic data. Furthermore, MSE field data were employed to calibrate the generation of the 3D FFN model, demonstrating the potential of this modeling approach in terms of characterizing real-world fracture networks.

2. Background on Fracture Network Modeling

Accounting for fracture network complexity is critical to properly manage hydraulic fracture architecture and available monitoring technologies. The discrete fracture network (DFN) is utilized to model fracture propagation and fluid flow in hydraulically fractured reservoirs. In DFN, hydraulic and natural fractures are explicit, and all individual elements are contained within the modeling space. The geometric and physical properties of these discrete elements are statistically assigned to each specific fracture element based on measured and interpreted data. While fracture geometry behavior can be assessed by modeling the fluid flow for each fracture using DFN, localized fracture stress and heat transfer in reservoir fluid flow are often neglected [3,17]. Additionally, concurrent simulations must be performed as the fracture network is randomly generated based on statistical data [18]. Pham and Weijermars developed the time-stepped linear superposition method (TLSM) to observe the propagation behavior of multiple hydraulic fractures by modeling the induced localized stresses of pressurized fractures and boreholes [19,20,21]. Because the TLSM is purely based on elastic rock deformation, this approach neglects fluid flow and proppant transport. Diaz-Acosta et al. explored the use of shear wave splitting to investigate fracture properties in carbonate reservoirs in order to enhance the understanding of reservoir heterogeneity [22]. Similarly, Bouchaala et al. examine azimuthal compressional seismic wave attenuation to characterize fracture orientations and anisotropic properties in fractured reservoirs, aiding in the improvement of hydrocarbon recovery strategies [23]. Other prior works using boundary element methods (BEM) and finite element methods have successfully coupled fluid flow, proppant transport, and elastic rock deformation with model fracture propagation and stress interference [24,25,26]. Since BEM solves for the unknowns of the fracture boundary and FEM solves for the unknowns inside the fracture, their integration allowed for the number of unknowns on the fracture’s surface to be reduced. However, both BEM and FEM require one grid block to be in communication with another, and the accuracy of the calculation depends on grid properties. Consequently, the utilization of these methods heavily relied on the quality of the grid refinement process in order to deliver robust solutions [27].

In addition to the previously mentioned methods, fractal theory can be employed to model the complexity of fracture networks [28]. The Lindenmayer system (L-system) has previously been used to explore organic growth, simulate developmental rules, and represent topological structures by creating complex shapes with simple rules [29]. This system is a parallel rewriting system used to define complex structures by successively replacing the parts of a simple initial object using a set of rewriting rules. The fundamental principle of fractal geometry generation is recursion, which achieves successive results through the repeated application of rules. The rewriting concept used within the L-system is an effective method for applying recursion to achieve organic, complex geometries. Complex geometries can be observed in natural phenomena, such as in the fractal nature of drainage basins, plant growth, and lightning. In petroleum engineering, earlier fractal models have been utilized to accurately represent the pore space in sandstone cores [30]. In 2017, the concept of utilizing the L-system was expanded to capture the complex characterization of fracture networks [31].

The flexibility of developmental rules allows complex fracture networks to be characterized by L-system fractals. Therefore, the propagation and branching of a complex fracture network, resulting from the interaction of natural fractures and hydraulic fractures, can be represented as a flexible sequence of rule adjustments. Previous works based on the combination of fractal theory with MSE data allowed fracture network complexity to be captured [31,32]. In these modeling strategies, the rewriting sequence and its components are guided by the occurrence and location of detected MSEs during well stimulation. Zhang et al. modified the classical 2D L-system to represent fractures in natural rocks [33]. Figure 2 shows a computed tomography (CT) image of a fracture in natural rock, accompanied by a proposed fractal fracture illustration. Intuitively, the representative fractal closely resembles the image of the fracture. The larger main fracture bifurcates into a series of self-similar micro-fractures, which can be represented by fractals generated using the L-system. Modifications to the L-system model parameters allow for an adequate match between the fractal and the fracture to be achieved. The modeling of the 2D FFN is demonstrated and discussed in Section 3. Subsequently, Section 4 introduces the 3D calibration system, which is built upon the general concepts of the 2D FFN system. Additionally, the software used for development in this work is MATLAB.

3. Modeling of 2D FFN

L-system fractal fractures are generated by defining their variables and then calibrated to MSEs with two objective functions. Previous works have described and generated L-system fractals with three primary components: an alphabet, an axiom, and generating rules [29]. Appendix A provides detailed explanations of how the L-system’s parameters relate to physical fracture properties. As such, this work focuses on the calibration of the L-system fractal to MSE data to generate an induced fractal fracture network. When the generating rules are applied to the alphabet, each letter of the composed string can be replaced with an action, and the whole string is then updated. The number of possible rules is constrained by the number of valid letters within the alphabet. By extending the alphabet, axiom, or rules, any possible combination of fractals can be generated. Figure 3 shows examples of fractals generated at different complexities by altering the system’s rules. Evidently, this type of flexibility is useful for matching complex spatial datasets when the proper matching criteria are applied.

3.1. Calibration of MSE to L-System Fractal Models

The fractal fracture network (FFN) can be calibrated through the fractal nodes, as MSEs often occur at the point of bifurcation or changes in fracture propagation direction [32]. The MSEs can be considered to be data points, while the nearest fractal node serves as a query point. It is assumed that fracture bifurcations or changes in propagation direction may induce stress concentrations and strain redistributions, which could lead to microseismic activity. This assumption warrants further investigation to fully comprehend the relationship between microseismic events and fractal structures. A geometric representation of this process is shown in Figure 4, and the workflow is illustrated in Figure A1 of Appendix A. The matching problem can be solved using the following process:

• For each MSE, employ the nearest point search method to locate the nearest fractal node based on Euclidean distance. Then, couple the nearest MSE with its corresponding fractal node.
• Determine the total number of MSEs coupled to fractal nodes.
• Calculate the distance between each coupled MSE and its respective fractal node. Sum these distances for all pairs of coupled MSEs and fractal nodes.

Using this process, two important metrics are established for a set of MSEs and a fractal network:

• The total distances from MSEs to nodes.
• The number of matched MSEs to nodes.

Subsequently, the following criteria were employed as objective functions for the matching process:

$$\mathrm{O}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e} 1:\max \left(\mathrm{n}\right)={\displaystyle \frac{\mathrm{n}}{\mathrm{m}}}>\mathrm{T}\mathrm{h}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{d} 1,$$

(1)

$$\mathrm{O}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e} 2:\min \left(\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{L}}_{\mathrm{i}\mathrm{j}}\right)<\mathrm{T}\mathrm{h}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{d} 2.$$

(2)

where

L_ij: distance between the coupled ith MSE and jth fractal node;

n: total number of coupled MSEs and fractal node pairs;

m: total number of MSEs.

The first objective function involves calibrating the fractal geometry to MSEs using nodes, where the task is to calculate and minimize the total distance with each fractal generation. The second objective function involves maximizing the amount of available MSEs and constraining the number of generated nodes in the FFN. The threshold values were established through empirical calibration and statistical analysis, based on observed data, in order to ensure a satisfactory match using an optimization technique. In this study, the objective was to achieve at least 80% matching of microseismic events to fractal nodes and to keep the average residual distance to less than 3 times the average fractal length.

3.2. Demonstration of 2D FFN Calibration with Synthetic MSE Data

A synthetic case was generated to demonstrate the 2D L-system described in Section 3.1. A set of synthetic 2D MSE data was generated to resemble field MSEs, utilizing a point cloud algorithm, specifically, to initiate the calibration process in this study. This set of data is contained within boundaries defined by the well, which acts as the starting point for all hydraulic fractures. The starting positions of the fractals are determined by the perforation location in the well, with the initial orientation constrained to be perpendicular to the well. These fractals are then generated randomly away from the initiation point, representing induced microseismic events as the fracture network grows from the wellbore. The generated data only include spatial coordinates in the x, y, and z directions, while additional information such as uncertainties and magnitudes are disregarded. Additionally, it does not involve a fully physical model that considers isotropy, anisotropy, elasticity, or viscoelasticity. These MSEs started at the wellbore and extended laterally (up to 700 ft) away from the perforation cluster. The calibration workflow outlined in Section 3.1 was used to match the generated fractal to the synthetic MSE data. The resulting 2D fractal fracture network is shown in Figure 5b, and is denoted in blue. Additionally, the wellbore with associated bi-wing perforations is shown in bold black lines, and the MSEs are plotted in red. While an adequate match was achieved, not all the MSEs were coupled to fractal nodes. These events occurred in closer proximity to the wellbore, where the segment of the 2D fractal is least complex.

Previous works aimed to generate FFNs and calibrate them to 2D microseismic data [31,34,35,36]. These 2D calibrations of fractals to microseismic data restrict all MSEs to 2D space; hence, they are coplanar. However, this is an inaccurate representation when considering real MSEs, where the 3D spatial coordinates of events are obtained. Figure 6 shows a set of open-source field microseismic data obtained from the Halliburton Laredo Sugg 171-A HFTS Project. Visually, the data possess three dimensions: the x direction (eastings), y direction (northings), and z direction (depth). The front view denotes the x–z axes, the top view denotes the x–y axes, and the side view denotes the y–z axes. Additionally, the data were filtered with a signal-to-noise ratio (SNR) higher than 2 to ensure high certainty and accuracy.

Evidently, performing calibration in only two dimensions for a set of microseismic data would yield conflicting results. As indicated by Figure 6, the front, side, and top views revealed disparate perspectives. A 2D FFN result calibrated using data from the x–y axes would greatly differ from that of the y–z axes. It is apparent that, instead of having multiple disparate 2D FFNs, a single 3D FFN incorporating all spatial dimensions is more appropriate. Additionally, removing a third dimension would render any results counterintuitive as this method is data-driven. Therefore, we outline a 3D L-system in the next section that can accommodate the spatial quality of MSE data.

4. Modeling of 3D FFN

4.1. Three-Dimensional FFN System Description

In a similar manner to 2D L-system fractals, the same drawing framework is applied in 3D space [37,38]. The generated strings command the turtle to move in varying directions, enabling the fractal to grow. In addition to turning left and right, as in the 2D system, the 3D turtle can pitch up/down and roll left/right, as shown in Figure 7. The starting vector is composed of a rotation matrix, R. H represents the head of the turtle, L indicates the left direction, and U indicates the up direction. Then, the new coordinates of the turtle are obtained by multiplying the starting vectors with the rotation matrix. An example of a 3D FFN is shown in Figure 8, illustrating the x–y plane, x–z plane, and y–z plane.

4.2. Three-Dimensional FFN Conceptual Model

The 3D FFN model described in this study adopts a semi-infinite boundary. More precisely, the boundaries are defined by the well, which acts as the starting point for all hydraulic fractures. These fractures have the potential to extend deeper into the reservoirs, while generating microseismic events. A conceptual model for a section of a horizontal well is formulated with perforation clusters and the induced fracture network using a system of 3D FFNs. Figure 9a illustrates a plug-and-perf fracturing treatment schematic involving four perforation clusters. Fluid is injected through the inlet region and then flows through the outlet region and the perforations. When the fracturing fluid encounters the rock, the original stresses are disturbed, causing stress concentrations to expand into areas of existing fractures, generated fractures, and activated fractures. When the pressure reaches the fracture pressure, the associated energy is released in the form of elastic waves, inducing microseismic events [17]. The 3D spatial coordinates of the MSEs produced from the hydraulic fracturing process are used to calibrate the FFN. Additionally, the location of the cluster spacing for each well constrains the starting position of each fractal. Figure 9b depicts a conceptual model consisting of a wellbore (black) and induced fracture networks (blue). Different perforation schemes are also shown to demonstrate the flexibility of the 3D FFN model.

4.3. A Case Study of 3D FFN

A set of field microearthquake records is used to calibrate the generation of fractals using the calibration workflow proposed in Section 3.1. The field data selected for this process are taken from the publicly available Laredo Sugg 171-A HFTS Project [39]. The purpose of the hydraulic fracturing test site (HFTS) project was to investigate fracture network interactions within several multi-well completion schemes. Microseismic event (MSEs) data were obtained by monitoring arrays during the hydraulic fracturing process for all wells. The monitoring arrays consisted of 12 to 18 DDS-250 electrodynamic digital downhole shuttles. The DDS-250 is a four-channel, digital downhole shuttle, specializing in multi-level, real-time, continuous, and high-definition borehole microseismic data [40,41]. Additionally, the ultra-low-noise digitizer enables the recording of low-amplitude microseismic events. These microseismic monitoring arrays were constantly repositioned to minimize the listening distance, thereby increasing the quality and quantity of high-certainty microseismic data. The tomographic inversion process assessed velocity variations within the solid media of the study area. Moreover, the algorithm used for determining the coordinates of the microearthquakes foci incorporated the parameterization of discontinuous vertical layering with a simple two-dimensional Chebyshev polynomial representation for surfaces and velocities [42]. Both P waves and S waves were used to determine the location and magnitudes of the microearthquakes. The errors seen in determining the coordinates and depth of earthquake foci are vertical positioning errors, coming from seismic sources and receivers; horizontal positioning errors, coming from well spacing; deviation survey errors, based on accelerometer measurements or multi-shot surveys; and anisotropy, induced by positioning errors. Figure 10 illustrates the well pad configuration and completion sequence, while the monitoring array configuration used in the well pad is shown in Figure 11. This HFTS comprises 11 wells completed with a zipper fracturing sequence, multiple perforation schemes (3 and 5 perforation clusters per stage), multiple ranges of fracturing stages (37, 47, 45, 49), 3 horizontal observation wells for microseismic detection (Wells 3SU, 5SM, 6SM), and 1 vertical observation well for use in a tiltmeter.

Following the zipper fracturing sequence, detailed in Figure 10, the initial completion was conducted in well 7SU. To analyze high-certainty microseismic (MS) data without interference from other fractures and wells, data from the first stage of well 7SU were selected. In well 7SU, the associated MS events were detected and recorded using the monitoring array in well 5SM. The distinction between microseismic events (MSEs) and noise is pivotal for ensuring the reliability of any analysis, particularly in noisy environments like carbonates. Therefore, the MSEs were filtered to have a signal-to-noise ratio (SNR) greater than 2 in order to differentiate between MSEs and noise events. An SNR threshold value of 2 would result in a false alarm probability of approximately 2.5% [43]. The total number of MS events after applying the filter was 147, compared to the unfiltered count of 222 events. The size of the circles in Figure 12a does not indicate the magnitude of the microearthquakes. Additionally, the magnitude type is local, with a maximum magnitude of −1.61 and a minimum of −2.57. Using the calibration workflow depicted and the criteria detailed in Section 3.1, the FFNs calibrated with the 3D microseismic data results are shown in Figure 12 and Figure 13. Additionally, key parameters are shown in

Table 1. Parameters used in the calibration of 3D FFN to field microseismic data.

Fractal Property	Value
Starting position (ft)	(−1500, −10,337, −7697)
Starting orientation (°)	0° (Parallel to x-axis)
Number of iterations	5
Length of segment (ft)	35 to 70 ft
Deviation angle	30° to 60o

. While most MS events could be matched to a fracture, there were many events that could not be matched. Out of 142 microseismic events (MSEs), 121 were matched to fractal nodes, resulting in a matching ratio of 0.82. The average residual distance was calculated to be 104.63 ft and the standard deviation was 44.50 ft.

Table 2. Summary of resulting calibration parameters.

Quality of Match Property	Value
Total number of MSEs	147
Total number of matched MSEs	121
Percentage of MSEs matched	82.31%
Sum of residual distance	12,660 ft
Average residual distance	104.63 ft
Standard deviation of residual distance	44.50 ft

provides quantitative metrics on the quality of fit between the 3D and field microseismic data. Although the residual distance is relatively high, it accounts for all MSEs, including observed outliers. Empirically, the unmatched events are located far away from areas where a high number of MS events occurred in close proximity to each other. Additionally, the unmatched MS events are spatially located away from the main cluster, where the event density is higher.

4.4. Four-Dimensional Microseismic Events Observations

The proposed methodology allowed for the calibration of 3D fractals to 3D MSEs, which was an extension of current 2D methods [31,34,35]. Figure 14 shows the 20 min interval time slices from different perspectives of the real MSE dataset from the Halliburton Laredo Sugg 171-A HFTS Project. Empirical observations indicate that MSEs are spatially constrained near the wellbore relative to their occurrence time. This suggests that a longer stimulation period does not necessarily lead to a larger stimulated fracture network, as the spatial extent of the MSE remains unchanged with time. However, the density of MSEs increases over time in the vicinity of the fracture cluster, indicating that the complexity of the fracture network is growing.

The observation that MSEs are spatially constrained near the wellbore, even as stimulation time increases, can be explained by stress concentration and energy dissipation. Hydraulic fractures tend to be localized near the wellbore where the stress is applied during fluid injection [44]. While shales are generally considered brittle, they can exhibit ductile behavior under certain conditions. Under elevated temperatures and pressures, shale exhibits ductile behavior. This ductility allows the rock to absorb more energy before cracking, which contributes to the development of a more intricate and interconnected fracture system over time [45].

The highest stress concentrations are typically seen near the wellbore where hydraulic pressure is applied, leading to the most intense fracturing activity occurring in this region [46]. As fractures form and propagate, they dissipate the energy created by hydraulic fracturing. Beyond a certain distance from the wellbore, the energy is insufficient to initiate new fractures, limiting the spatial extent of MSEs. The increasing density of MSEs near the fracture cluster over a longer stimulation period suggests that while new fractures are not necessarily propagating further from the wellbore, the existing fracture network is becoming more complex. Initial fractures create new stress fields that can induce secondary fractures, increasing the density of the fracture network [47]. As induced fractures interact and collide, more pathways and branches are created, enhancing the network’s complexity.

Understanding these spatiotemporal patterns can significantly inform fracturing operations. The insight that MSE density increases near the fracture cluster over time allows operators to focus on maximizing fracture network complexity rather than merely extending stimulation periods. Additionally, the concentration of MSEs near the wellbore due to stress concentration informs better management of hydraulic pressure and fluid injection rates to minimize the risk of excessive fracturing or wellbore instability. Additionally, understanding the role primary fractures in inducing secondary fractures aids in the design of treatments to promote secondary fracture development. By leveraging these insights, the effectiveness of fracturing operations can be enhanced, leading to improved hydrocarbon recovery and operational efficiency.

4.5. Three-Dimensional FFN Model Discussions

The 3D FFN model is a significant advancement compared to the 2D FFN model, providing a more complete and nuanced understanding of fracture networks and their behavior within the reservoir. Unlike 2D models, which overlook propagation in a single direction, d, 3D modeling offers insights into fracture growth in all directions, leading to a comprehensive view of how fractures develop and connect. This enables more precise stress and strain analyses, enhancing the understanding of how applied stresses influence fracture formation and interactions. The improved accuracy in fracture network representation and reservoir characterization also leads to better predictions of fluid flow and resource recovery. By refining stimulation and production strategies based on the detailed 3D model, fracturing operations and overall reservoir management can be significantly improved.

However, the calibration of generated fractal fracture networks using microseismic events (MSEs) has several limitations. The first limitation is that fractal models assume self-similarity as they grow, which may not fully capture the heterogeneity of natural fracture networks. Another major limitation is the quality and quantity of microseismic data. The availability of MSE data can be limited, particularly in low-permeability reservoirs and high-temperature reservoirs, limiting the resolution and accuracy of the fracture network model. The deployment and maintenance of microseismic monitoring equipment can be challenging and costly, limiting continuous data collection [45]. Moreover, MSE data can be noisy, making it difficult to distinguish between real events and noise [48].

Temporal constraints and resolution are also significant limitations. Fracture networks evolve over time, and static models based on MSEs do not account for these dynamic changes. Additionally, the temporal resolution of MSE data might not capture rapid changes in the fracture network during hydraulic fracturing [46]. The heterogeneous physical and mechanical properties of rocks, such as variations in brittleness and ductility, are not always accurately represented in MSE-based models. The existence of natural fractures may significantly influence the propagation of new fractures, which further complicates the calibration of fractals to MSE data.

Spatial resolution is also a concern. MSE data generally offer insights into fractures near the wellbore, potentially facing reduced accuracy or overlooking distant fractures. Additionally, the spatial resolution of MSE data may not be sufficient to capture micro-fractures, which leads to incomplete network models. As such, calibrating fractal models to match MSE data can be complex and requires robust datasets and validation techniques. Addressing these limitations requires the integration of MSE data with other geophysical, geological, and operational data to create a more comprehensive fracture network model.

Lastly, while the computational demands of the 3D FFN model in this study are relatively manageable due to its focus on a single fracturing stage, extending the model to a larger multi-well scenario would significantly increase its computational requirements. This heightened demand could restrict the model’s applicability to large-scale projects and hinder the real-time analysis and calibration of microseismic events. Moreover, as the complexity of the model grows with more L-system iterations, it may face scalability challenges, making it more difficult to apply in extensive or highly detailed reservoir studies. Addressing these challenges will require advanced computational resources and optimization techniques to ensure effective application and scalability in more complex scenarios.

5. Conclusions

While the 2D FFN approach is valuable, it fails to accurately depict fracture networks because it cannot capture the spatial complexity and interactions of fractures in three dimensions. This limitation results in less reliable predictions of fracture behavior and suboptimal hydraulic fracturing strategies. To address this issue, this study introduces an advanced 3D fractal fracture network (FFN) model that integrates 3D microseismic data to achieve a more precise and detailed representation of fracture networks in unconventional reservoirs. This novel model enhances the traditional 2D FFN by incorporating the complexities of three-dimensional space and provides a calibration method designed to accommodate three-dimensional microseismic data. Field microseismic event data were used to validate the 3D FFN model’s effectiveness and compatibility with real-world fracture networks.

The results generated from this study reveal the following insights:

• The proposed novel solution successfully calibrated the 3D FFN using field microseismic event data for a single fracturing stage.
• The areal extent of MSE occurrence did not increase with time, suggesting that longer injection periods do not enhance the general size of the fracture network.
• The density of MSE occurrence increased with time in conditions of proximity to the wellbore, indicating that fracture network complexity and connectivity increase with time.

Future research should focus on several key areas to further advance the current 3D FFN model. These include integrating physical phenomena such as stress shadow effects, volume distribution, and proppant placement in order to provide a more comprehensive understanding of fracture growth and network complexity. Additionally, incorporating data from geological, geophysical, and production sources, along with high-resolution and real-time microseismic data, will enhance model accuracy and reliability. Field testing and validation through long-term studies and comparative analyses with other advanced models will further validate predictions and performance. Applying insights from the 3D FFN model to develop and optimize hydraulic fracturing techniques and reservoir management practices will lead to more efficient fracturing operations.