Hydraulic fracturing is widely used to stimulate unconventional reservoirs, but a systematic and comprehensive investigation into the hydraulic fracturing process is rare. In this work, a discrete element-lattice Boltzmann method is implemented to simulate the hydro-mechanical behavior in a hydraulic fracturing process. Different influential factors, including injection rates, fluid viscosity, in-situ stress states, heterogeneity of rock strengths, and formation permeability, are considered and their impacts on the initiation and propagation of hydraulic fractures are evaluated. All factors have a significant impact on fracture initiation pressure. A higher injection rate, higher viscosity, and larger in-situ stress increase the initiation pressure, while a higher formation permeability and higher heterogeneity decrease the initiation pressure. Injection rates and heterogeneity degrees have significant impacts on the complexity of generated fractures. Fluid viscosity, in-situ stress states, and formation permeability do not change the geometrical complexity significantly. Hydraulic fractures are usually tensile fractures, but many tensile fractures also have shear displacement. Shear fractures are possible and the shear displacement can be significant under certain conditions, such as a high injection rate, and a high heterogeneity degree.

This paper introduces a comprehensive C++ software package, HatchFrac, for stochastic modelling of fracture networks in two and three dimensions. Two main methods, the inverse CDF method and acceptance-rejection method, are applied to generate random variables from the stochastic distributions commonly used in discrete fracture network (DFN) modelling. The multilayer per-ceptron (MLP) machine learning approach, combined with the inverse CDF method, is implemented to generate random variables following any sampling distribution. To make the code faster, we extend the Newman-Ziff to determine clusters in the fracture networks. When combined with the block method, the Ziff algorithm improves the coding efficiency significantly. The software generates the T-type fracture intersections in the network, which can be used in applications involving fracture growth or incorporating geomechanics. We introduce three applications of HatchFrac that demonstrate the versatility of our software: percolation analysis, fracture intensity analysis, and flow and connectivity analysis.

Fractures play an essential role in formations with low permeability; however, fracture sealing significantly reduces the permeability of fractures. The mechanism of how fracture sealing impacts the macro-scale fluid flow is rarely investigated. Here, we simulate sealing in two- and three-dimensional orthogonal fracture networks and investigate the impact of sealing on the percolation of these fracture networks. We find that a small amount of sealing can prevent the formation of spanning clusters, which suggests that global connectivity is rarely realized. Without significant stress perturbations, most fractures are partially sealed and non-critically stressed, and they usually do not contribute much to the fluid flow. However, under a significant stress perturbation, such as hydraulic fracturing, the well-connected and critically oriented fractures become critically stressed and slide because of the increased pore pressure. Partially sealed and non-critically stressed fractures can also contribute to the fluid flow by enlarging the stimulated reservoir volume (SRV). We estimate the stimulated reservoir volume in two dimensions by dividing the target distance (LSRV) into two parts. One is the distance limiting generation of hydraulic fractures (ΔLh), and the other is the limiting distance of making natural fractures slide (ΔLs). A rough estimation yields an elongated shape of the SRV, which is consistent with observations from microseismicity maps.

Fractures and their connectivity are essential for fluid flow in low perme-ability formations. Abundant outcrops can only provide two-dimensional (2D) information, but subsurface fractures are three-dimensional (3D). The percola-tion status of 3D fracture networks and their 2D cross-section maps are rarely investigated simultaneously. In this work, we construct 3D fracture networks with their geometries characterized by different stochastic distributions. Then, we take cross-section maps to mimic real outcrops and label clusters to check the percolation status of 3D fracture networks and their 2D cross-section maps. The properties, reflecting the connectivity of two essential phases, are summarized and analyzed. We find that clustering effects impact local intersections significantly but have negligible impacts on fracture intensities of 3D fracture networks. The number of intersections per fracture is not a proper percolation parameter for complex 2D and 3D fracture networks. Fracture intensities are scale-dependent and usually decrease with increasing scales. The real fracture networks in the subsurface should be geometrically well-connected and pervasive if their outcrop maps are well connected. In particular, the fracture intensity of the real fracture network can be several times (at least 3.6 times) larger than the intensity at percolation. However, if outcrop maps are not well-connected, but their intensities are large enough (at least 0.43 times as large as the intensity at percolation), corresponding 3D fracture networks can also form a spanning cluster and show good connectivity with a high possibility.

The fractal dimension and multifractal spectrum are widely used to characterize the complexity of natural fractures. However, a systematic investigation on the impact of different fracture properties (fracture lengths, orientations, center positions, system sizes) on the fractal and multifractal characterization of complex fracture networks is missing. We utilize an in-house developed DFN modeling software, HatchFrac, to construct stochastic fracture networks with prescribed distributions and systematically study the impact of four geometrical properties of fractures on the fractal and multifractal characterization. We calculate the single fractal dimension and multifractal spectrum with the box-counting method. The single fractal dimension, D, and the difference of singularity exponent, ∆α, are used to represent the fractal and multifractal patterns, respectively. We find that fracture lengths, orientations and system sizes have positive correlations with D and ∆α, while the system size has the most significant impact among the four parameters. D is uncorrelated with fracture positions (FD), which means that a single fractal dimension cannot capture the complexity caused by clustering effects. However, ∆α has a strong negative correlation with FD, which implies that clustering effects make fracture networks more complex, and ∆α can capture the difference. We also digitize 60 outcrop maps with a novel fracture detection algorithm and calculate their fractal dimension and multifractal spectrum. We find wide variations of D and ∆α on those outcrop maps, even for outcrops at similar scales. It means that a universal indicator for characterizing fracture networks at different scales or the same scale is almost impossible.

Stimulated reservoir volume (SRV), the high-permeable fracture network created by hydraulic fracturing, is essential for fluid production from low-permeable reservoirs. However, the configuration of SRV and its impacting factors are largely unknown. In this work, we adopt the stochastic discrete fracture network method to mimic natural fractures in subsurface formations and conduct a global sensitivity analysis with the Sobol method. The sensitivity of different fracture properties, including geometrical properties (fracture lengths, orientations and center positions), mechanical properties (fracture roughness and fracture strength), fracture sealing properties (probabilities of open fractures and segment lengths) and the fracture intensity, are investigated in two and three-dimensional fracture networks. JRC-JCS model is adopted to identify critically stressed fractures. We find that critically stressed fractures compose the backbone of SRV, while partially open fractures can significantly enlarge the size of SRV by connecting more critically orientated fractures. The fracture roughness is the most influential factor for the total length (area) of critically stressed fractures. For the relative increase of SRV (RI) in 2D/3D fracture networks, the probability of open fractures is the most significant factor. The fracture lengths and center positions are essential factors for RI in 2D fracture networks but insignificant in 3D fracture networks. This work provides a realistic scenario of the subsurface structure and systematically investigates the influential factors of SRV, which is useful for estimating the size of SRV and predicting shale gas reservoirs’ production in an accurate and physically meaningful way.

The fractal dimension and multifractal spectrum can characterize the complexity of fracture sets. However, studies of impacts of fracture geometries on their fractal and multifractal characteristics are largely insufficient, especially for three-dimensional (3-D) fracture networks (natural fractures are always 3-D instead of 2-D). In this work, we construct 3-D stochastic discrete fracture networks with an open-source DFN software, HatchFrac. Systematical investigations are then conducted to study the impact of geometrical fracture properties and system sizes on the fractal and multifractal characteristics. The box-counting method is adopted to calculate the fractal dimension and multi-fractal descriptors. The fractal dimension, D, and the difference of the singularity exponent, ∆α, represent the fractal and multifractal patterns, respectively. Two critical (percolative and over-percolative) stages of fracture networks are considered. 3-D fracture networks share similar characteristics with 2-D fracture networks at percolation. However, results at an over-percolative stage are systematically different. At the first stage, fracture orientations (κ), lengths (a) and system sizes (L) have positive correlations with D and ∆α. D is weakly correlated with fracture positions (FD), meaning that the fractal dimension is insensitive to clustering effects. However, ∆α is strongly correlated with FD, implying that ∆α can characterize the heterogeneity caused by clustering effects. a and L are positively correlated with ∆α, and κ and FD have negative correlations. At stage two, the sensitivity results on D are similar to stage one, but a and L become negatively correlated with ∆α. Impacts of κ and FD become more significant.