Hydrofracturing is the primary method of hydrocarbon reservoir stimulation that enhances unconventional gas recovery. Unconventional natural gas reservoirs usually are composed of different minerals and have a composite or heterogeneous microstructure (Fouche et al. 2004; He et al. 2015). As a result, the initiation and propagation of hydrofracturing cracks are usually influenced by the contrast between Young’s modulus and Poisson’s ratio of different components (Garcia et al. 2013; Sarmadivaleh and Rasouli 2015; Stanchits et al. 2014). Understanding the initiation and propagation of hydraulically induced cracks is key to predicting and measuring hydrofracturing behavior of low-permeability reservoirs.

Numerical models have been developed over recent decades to investigate the crack behavior in rock media. However, due to the heterogeneity of unconventional reservoirs, classical fracture models do not accurately predict the fracture pattern in these rocks. In recent years, different techniques have been used to incorporate heterogeneity into numerical models, such as FEM (Wang 2013), RFPA (Tang et al. 2002; Yang et al. 2004; Li et al. 2013; Wang et al. 2013), and PFC (Al-Busaidi et al. 2005; Hiyama et al. 2013; Shimizu et al. 2011, 2014) methods. In addition, there are other methods, such as unconventional fracture models (UFM) (Kresse et al. 2011, 2013), BEM (Zhang et al. 2007), XFEM (Wang et al. 2015), DDM (Gu et al. 2008). These numerical studies are useful for understanding and characterizing the initiation and propagation of hydrofracturing cracks. However, these numerical models are generally applied to simple media and become extremely complicated when random discontinuities or material heterogeneities are involved. Thus, the effectiveness of these models usually needs to be verified. This requires studying actual hydrofracturing crack patterns in the field or laboratory. However, it is challenging to observe these patterns in the rock media. In general, the acoustic emission method is widely used to monitor hydrofracturing crack patterns (Cipolla et al. 2008; Maxwell and Cipolla 2011; Maxwell et al. 2013; van der Baan et al. 2013; Warpinski 2014; Maxwell et al. 2015). However, uncertainties in the microseismic data, limitations of azimuthal coverage, and shortcomings of velocity models make it difficult to determine complex fracturing processes (Castano et al. 2010; Cipolla et al. 2012; Johnston and Shrallow 2011). In addition, the CT identification method has been employed to characterize fracture geometry in traditional hydrofracturing experiments carried out on rock or cement samples (Renard et al. 2009; Guo et al. 2014; Hampton et al. 2014; Zou et al. 2016). This method depends on the resolution and capabilities of the CT system. As a result, the size of the sample is usually limited and microfractures are difficult to identify. Therefore, a visual representation of the geometry of hydrofracturing cracks is of particular use to understanding the hydraulic fracturing pattern. Usually, Polymethyl methacrylate (PMMA), glass, and resin are widely used because of cost effectiveness and ease of processing. Wu et al. (2004) investigated hydrofracturing crack propagation in a layered formation which was composed of rigid and soft layers, and drew a map of fracture behavior in the vicinity of an interface. Wu (2006) revealed mixed mode fracturing by applying mode III loading to the cylindrical PMMA samples with a circle internal fracture. Wu et al. (2008) carried out a series of hydraulic fracturing test to observe fracture initiation and propagation in well orientations based on transparent PMMA samples. It was observed that far-field stresses had little influence on fracture initiation, but these stresses did influence the final fracture orientation. Kaswiyanto et al. (2014) used resin specimens to validate fracture mechanisms in the laboratory. Bunger et al. (2013) utilized both glass and PMMA samples to verify the efficiency of a numerical model by observing the crack path, and the evolution of the fluid and fracture fronts. The results demonstrated that considerable coupling among fluid flow, elastic deformation, and radially symmetric crack growth captures enough of the relevant physical processes. However, traditional transparent materials (PMMA, glass, resin) are homogeneous and are of limited use for understanding heterogeneous samples. Currently, 3D printing techniques applied to manufacturing, military use, the medical and food industry, and the aerospace industry have shown to be a promising tool (Berman 2012; Vaezi et al. 2013). Despite limitations of the 3D printing technique applied to rock mechanics (e.g., low strength and poor brittleness of printing material), this technique shows obvious advantages for rock mechanics studies. First, it is precise and allows flexibility for controlling the geometry of the specimen. This is especially useful in producing samples with complex components. Second, the printing process can be repeated, making it is possible to produce multiple identical specimens. These advantages in specimen preparation make 3D printing suitable for manufacturing complicated transparent 3D solids (Miller et al. 2011; Campbell and Ivanova 2013; Dal Ferro and Morari 2015).

The purpose of the research reported in this study was to observe hydrofracturing crack distribution in heterogeneous media experimentally. The transparent specimens containing heterogeneous gravel were constructed by means of the 3D printing technique. The effect of the geostress condition was tested through different confining pressures. The remainder of this paper is organized as follows. In Sect. 2, we outline the preparation of the 3D printing specimens, including the test of materials, reconstruction of the model, conditions and specifications. In Sect. 3, we present morphologies of the hydrofracturing cracks in the transparent specimens. The study’s conclusions are presented in Sect. 4.

2.1 CT identification of heterogeneous structure

To obtain the distribution of gravel in the rock materials, the inner structures of two concrete specimens that contain gravel were probed by means of CT scanning. The concrete specimens were constructed by mixing cement and aedelforsite gravel, in which the amount and grading of the gravel was probed by detecting the gravel distribution of natural glutenite cores drilled from the Shengli SINOPEC oil field at a depth of 4000 m. The reason for using the inner structure of cement specimens for reference is that the matrix and gravels can be easily distinguished in the artificial specimens. In order to construct the 3D numerical model, a multi-thresholding segmentation method and a self-developed computer program (Ju et al. 2013, 2014b) were used to enhance and digitize the original images into binary pixels of gray and white that represented matrix and gravel, respectively. Furthermore, square images were cut into circles in order to construct cylindroid models. The whole process is shown in Fig. 1a. Moreover, we used the MIMICS^© (http://biomedical.materialise.com/mimics) software as a platform to build the 3D model as shown in Fig. 1b. The model ensures the accurate characterization of the spatial distribution of the gravel because of the high precision of the CT images.

The construction of 3D models; a the process of image processing; b the construction of the heterogeneous model in MIMICS software

2.2 Preparation of the 3D printing samples

We employed the Object Connex 500 3D printer (Mueller et al. 2015) to manufacture the cylindroid samples (see Fig. 2a). The device used the Polyjet 3D printing technique (Barclift and Williams 2012; Lifton et al. 2014) for the heterogeneous sample construction. This technique involved spraying liquid photopolymer through the print head on the build tray, and then solidifying the liquid simultaneously utilizing ultraviolet light, in layer by layer fashion. This process is sketched in Fig. 2b. This technique is available for multiple printing by spraying various types of printing materials from different print heads. This approach makes it possible to construct heterogeneous objects. In addition, the 3D printer has advantages of high precision and automation (Ju et al. 2014a), with print resolution up to 600 × 600 × 1600 dpi, a dot accuracy of 10–50 μm, and molding thickness ranging 16–30 μm.

a Object Connex 500 3D printer; b schematic diagram of Polyject printing (Barclift and Williams 2012); c 3D printing samples. Samples are denoted as type i and ii

The general procedure for preparing a printing sample is as follows. First, the 3D numerical model was transited to the STL files. Second, the STL files were put into the Objet Studio software to create a 3D model. To distinguish different components of the heterogeneous model, each component, such as matrix and gravel, was put in independent STL files. From this, the Objet Studio software automatically assembled the whole model in accordance with the location coordinates of the parts. Third, the material properties were accordingly assigned to the matrix and particles in line. Using the Objet Studio, we were able to set the additional printing parameters, such as the material properties of contact surfaces and the thickness of supporting materials, in order to finalize the printing setup. In this study, we adopted the transparent photopolymer Vero Clear to print the matrix, and white RGD to print the gravel. In the hydrofracturing experiment, an open hole was required to inject water that would fracture the sample. Therefore, the light type lattice supporting material Fullcure 705 was employed to form the open hole without fillings. The supporting material features low-transparency, loose structure, low intensity, and can be physically removed or washed through high-pressure water easily after the sample is prepared. In the end, all the data of the model were transmitted to the 3D printer. Consequently, two types of cylindroid samples were printed, which are shown in Fig. 2c and denoted as type (i) and (ii). The samples were 54 mm in diameter and 108 mm in length. A 4-mm-diameter vertical hole was left at the center of the sample with a length of 54 mm. In addition, all the surfaces of the samples were polished to make the inner structure available for observation.

It is noteworthy that the printing model is structurally complicated, Thus it is possible that parts of the resinous materials did not polymerize completely or even remained in a liquid state once the printing was completed. In addition, solidification and polymerization of the resinous printing materials are sensitive to ultraviolet light and temperature. As a result, the rigidity and strength of the printing sample may change over time. To avoid this problem, the printed samples were stored at 60 °C for 48 h in order to ensure mechanical stability of the resinous printing material. Furthermore, the physical and mechanical properties of the printing materials were tested before implementing the hydraulic fracturing tests. The measured results of Vero Clear and RGD are listed in Table 1. Accordingly, the compressive strength of Vero Clear is similar to rock material, but there is a difference between the tensile strength of the printing materials and actual rock, while RGD was even weaker than Vero Clear. In fact, there are other printing materials that have greater brittleness than those used in our experiments, such as plastic and gypsum (Jiang and Zhao 2015; Jiang et al. 2016). However, these materials were not used for two reasons. First, these materials are poorly suited for printing multiple types of materials in one sample. Second, these materials are usually opaque, while our samples had better transparency.

2.3 Triaxial hydrofracturing experiment

The experiments were performed using the triaxial hydrofracturing device at Monash University which was developed based on the triaxial apparatus (Wasantha et al. 2013). Figure 3a illustrates the assembly of the hydrofracturing system, in which the fluid was injected into the sample through a syringe pump that is capable of applying a maximum injecting pressure of 50 MPa. The triaxial pressure was applied by the loading cell which is shown in Fig. 3b. The cell was built from high tensile steel and was designed to operate at a 30 MPa maximum confining pressure. The top of the cell was placed on the base of the apparatus and the bolts were inserted and tightened. A steel brace fixed the position of the loading ram such that it stayed centered and vertical during loading. The loading frame was capable of applying an axial load of up to 245 kN. Because horizontal geostress greatly influent crack propagation (Koceir and Tiab 2000; Nasehi and Mortazavi 2013), three magnitudes of confining pressures were applied while implementing the hydrofracturing tests. The three magnitudes of confining pressure applied to the samples were 1, 3 and 5 MPa. The vertical pressures were always kept at 20 MPa. The fracturing fluid was deionized water, and it was injected at a rate of 3 mL/min in all experiments. Because 3D printing techniques make it possible to construct numerous samples with the same structure, three samples for type (i) models were prepared for all confining pressure conditions. Another two samples for type (ii) models were prepared for confining pressure of 1 and 3 MPa.

Triaxial hydrofracturing device in Monash University. The photo of the device is shown on the left, while the schematic diagram of the loading cell which is indicated with a yellow square is shown on the right

3.1 Experimental observation

Because of the transparency of the matrix, the morphology of the hydrofracturing cracks can be identified by the naked eye observing the trajectories of the tracer. However, it is still challenging to record the crack morphology clearly with camera. To achieve this goal, the samples were laid on a rotatable platform while two fluorescent lamps were fixed at two sides of the fractured sample. The best observation position was found by rotating the platform. Thus, the cracks inside the sample were clearly viewable with a digital camera from only one side.

Figure 4 illustrates the top-view and side-view photographs of the fractured samples and the injection pressure curves, respectively. It shows that the cracks emerged vertically without transverse cracks appearing and that the gravels along the crack trajectory were almost broken apart at weak strengths as shown in Table 1. When the confining pressure was 5 MPa (Fig. 4a), a bi-winged crack formed with an included angle smaller than 180°. However, once confining pressure was <5 MPa as shown in Figs. 4b, c, only a radial crack formed. Radial-crack height was confined to the region near the simulated wellbore, even in different types of samples. In addition, the initial positions of the cracks are generally the same, which is also the initial position of one branch at confining pressure of 5 MPa. It seems that there was a preferential azimuth for crack initiation near the simulated wellbore.

Crack morphologies under different confining pressure. The blue tracer indicates the crack plane. The rows from a to c refer to the results of various confining pressures (5, 3, and 1 MPa)

Note that the original crack arrested at the interface of different materials as shown in Fig. 4a and b, but a secondary crack emerged. Previous studies have implied that there were several scenarios for crack behavior as it approaches the interface between soft regions and rigid regions (Wu et al. 2004). Fracture pattern depends on the energy release rate, \(G_{1} \, = \,\frac{{K_{lC}^{2} }}{{E^{\prime}}}\), where, \(E^{\prime}\, = \,E\) and \(E^{\prime}\, = \,E\,\left( {1\, - \,U^{2} } \right)\) for plane stress and plane strain problems, respectively. Unstable crack propagation was onset as K _l  = K _lC . The stress intensity was \(K_{l} \, = \,\sigma \sqrt {\pi a}\), where σ was the normal stress along the crack faces, and a was the length of the crack. Once the cracks exited the weaker gravel the energy release rate (G _l ) tended to be zero (Nuller et al. 2001). This meant that the crack would arrest according to the Griffith criterion (Griffith 1921). However, σ became discontinuous across the interface and jumped to a value larger than the tensile strength of the rigid matrix. Thus, the high value of σ fostered the generation of a secondary crack that was oriented parallel to the original crack.

3.2 CT identification

CT technology was employed to identify the distribution of the cracks in detail. Slices of the fractured samples used in all cases are shown in Fig. 5. It is evident that the cracks initiated not only from gravel as expected but also from the matrix. In addition, it is noteworthy that the gravel was not homogeneous. Flaws can be found easily inside, but they were insufficient to break the gravel. In fact, the number and development of the faults in the gravel would enhance crack growth. Consequently, few crack would propagate through the boundary of the gravel. To assess crack volumes, the width of the cracks were measured based on the CT images. The results show that the crack width is about 600 μm under a confining pressure of 5 MPa. However, this value decreased to 400 μm when the confining pressure was lowered to 3 MPa. This suggests that the crack volume decreases with the confining pressure.

CT slices of the fractured samples. The rows a–c refer to the results under confining pressure of 5, 3, and 1 MPa, respectively. The columns from left to right in b and c show results for samples with type i and ii structures

3.3 Fracture breakdown

Crack morphology is associated with the fracture breakdown pressure, which is the critical pressure during pressurization of a borehole. According to the injection-time curves in Fig. 4, the wellbore pressure drops instantaneously at a peak point, which is defined as the breakdown pressure and marks onset of unstable fracture propagation. In general, the rock is assumed homogeneous and isotropic, following Hooke’s law of linear elastic behavior. According to the Hubbert-Willis expression (Hubbert and Willis 1972), the condition for crack initiation is reached when,

The 3D printing samples are a desirable impermeable media with no pores contained. Therefore, the pore pressure is zero (P ₀ = 0). According to the experimental results, the tensile strength, T _o, of the matrix is 38 MPa, compared to 33 MPa for the gravel. The experimental and theoretical results plotted in Fig. 6 illustrate that the experimental breakdown pressures follow a variation trend similar to theoretical breakdown pressures, but that the values appear much lower in the experimental results. In the heterogeneous medium, the gravel has significant influence on stress distributions, which leads to weakened regions in the crack initiation. Mathematically, the linear variation of the breakdown pressure obtained in the experiments could be expressed as,

Comparison of theoretical and tested breakdown pressure

Comparing Eq. (2) with Eq. (1), the expression yields,

It seems that heterogeneity enhances the influence of confining pressure on the breakdown pressure of cracks. The smaller the confining pressure, the more difficult it is to achieve the theoretical breakdown value. Consequently, heterogeneity of the material would have greater effect on the initiation of the cracks.

This study reports a series of hydrofracturing tests designed to observe the hydrofracturing crack distribution in heterogeneous media based on transparent samples. Advanced 3D printing techniques were employed to construct heterogeneous samples that represented accurate distribution of gravel. To probe the effect of confining pressure on hydrofracturing crack pattern, three confining pressure states were utilized. The main conclusions are as follows:

1. (1)

   Geostress has significant influence on the characteristics of hydrofracturing cracks in the heterogeneous media. When the confining pressure is high (i.e., 5 MPa), a bi-winged crack formed. However, when the confining pressure was <5 MPa, only a radial crack formed. It seems that hydrofracturing cracks are reduced as the confining pressure increases in the heterogeneous media.

2. (2)

   The heterogeneous gravel led to segmentation of the cracks and significantly affected the initiation of cracks by changing the stress state around the simulated wellbore.

It should be noted that the mechanical properties of the 3D printing samples exhibits a certain degree of deviation, despite the preliminary attempts in this study that employed 3D display and visualization of crack patterns under the coupled effect of heterogeneous gravel and confining pressure. Thus, further research will aim at improving the 3D printing materials to make their properties approach that of actual rocks. With such improvement, we would return to the hydrofracturing tests to verify and characterize crack patterns.