Multifractal Analysis of Acoustic Emissions during Hydraulic Fracturing Experiments under Uniaxial Loading Conditions: A Niutitang Shale Example

Key Laboratory of Metallogenic Prediction of Nonferrous Metals and Geological Environment Monitoring (Central South University), Ministry of Education, Changsha 410083, China School of Geosciences and Info-Physics, Central South University, Changsha 410083, China Institute of Acoustics, Chinese Academy of Sciences, Beijing 100190, China Beijing Engineering Research Center for Offshore Drilling Exploration and Measurement, Beijing 100190, China School of Earth and Environmental Sciences, The University of Queensland, Brisbane QLD 4072, Australia


Introduction
The shale gas revolution in North America has affected the world energy market and greatly promoted the process of global shale gas exploration and production [1,2]. Shale gas is usually stored as adsorbed or free gas in organic-rich mud shale, which has (extremely) low porosity and permeability [3,4]. Hydraulic fracturing (HF) improves reservoir permeability by injecting high-pressure fluid into the ground to form a fracture network that is effectively connected to the wellbore [5][6][7]. Microseismic monitoring is generally implemented as a routine technology to monitor the development and expansion of fracturing networks in real time [8,9]. Additionally, microseismic monitoring technology has also been widely used for monitoring in other fields, including geothermal production, mining, carbon dioxide storage, and other industrial fields [10,11].
However, in-field microseismic monitoring, underground geological conditions are quite complex, such as the unknown distribution of natural/preexisting fractures and heterogeneous stress conditions. To effectively control the variables and focus when studying the evolution mechanism of HF, a feasible and important strategy is to simulate field-scale HF monitoring using laboratory-scale experiments with acoustic emission (AE) monitoring [12,13]. AE is primarily related to the initiation and expansion of cracks. The frequency of laboratory-and small-scale acoustic emissions is generally in a range between several kHz to MHz. By analyzing the AE characteristics of rock fractures, it is helpful to understand the rock fracturing mechanism and provide a theoretical and technical basis for field HF monitoring [14]. Numerous studies have reported laboratory-scale rock fracture analysis using AE monitoring, including the rupture process of brittle rock [15,16], the mechanical and acoustic properties of shale rocks [17,18], the characteristics of fracture propagation with bedding planes at different directions [19], and the focal mechanism of AEs during laboratory-scale HF experiments [20,21].
The processes of rock cracking and propagation have self-similarity and can be described by the fractal method (e.g., fractal dimension) [22]. There are many studies that have characterized rock fracture evolution using the fractal or multifractal analyses of the accompanying AE signals [23][24][25][26]. For example, Kong et al. [26] analyzed the fractal characteristics of AE parameters of methane-containing coal samples in triaxial compression experiments. The results of the correlation dimension revealed that the fractal dimension could describe the propagation of a crack. Compared with the conventional fractal method, the multifractal theory has the potential to describe the rock rupture evolution state comprehensively [27,28]. Kong et al. [29] used the multifractal method to study the AE characteristics of sandstone subjected to thermal treatment. The AE time series had multifractal features, and the fractal dimension followed an increasing-decreasing trend against the temperature. Kong et al. [30] applied the multifractal theory to analyze the AE characteristics of coal samples subjected to uniaxial compression. The results revealed that the time-varying multifractal characteristics could describe the AE mechanisms and the difference in strong and weak AE signals. Tan et al. [31] studied the multifractal characteristics of shale acoustic emission energy in uniaxial experiments under different immersion conditions. The results showed that the multifractal spectrum of the AE parameters was influenced by soaking times, which was directly related to changes in the pores and cracks of the shale samples. However, most of the above studies have focused on conventional compression experiments with intact samples. Therefore, quantitative multifractal analysis of AE parameters that results from shale hydraulic fracturing experiments requires further study.
In this study, the characteristics of AE events in hydraulic fracturing experiments under a uniaxial loading condition are analyzed to evaluate shale fracturing at different stages. The AE signals are quantified using multifractal analysis to study the rupture scale and energy characteristics. Finally, the reliability of the results is further verified using a timefrequency analysis and X-ray computed tomography (CT) scanning. The multifractal analysis can help provide quantitative guidance for fracture characterization during field hydraulic fracturing operations.

Sample Preparation.
The samples were collected from outcrop shale samples of Niutitang Formation in Yiyang, Hunan Province, China. The outcrop shales were cored into cylindrical samples (size of 50 × 100 mm) with errors of less than 0.5 mm, and small cylindrical holes (size of 3 × 55 mm) were drilled at the center of the samples. The size and other physical properties of the six samples are shown in Table 1. Three of the six samples had horizontal bedding characteristics (the direction of bedding plane is parallel to the axial direction), which are termed as H1, H2, and H3. The other three samples were characterized by vertical bedding (the direction of bedding plane is perpendicular to the axial direction), which are termed as T1, T2, and T3. There were no visible fractures at the surface of the samples. The mineral components of the shale samples were quartz, with an average value of 74.2%; feldspar, with an average value of 14.3%; and clay minerals, with an average value of 7.8%. The samples tested were brittle rock and were not influenced by weathering. During the experiment, six acoustic emission (AE) sensors were placed at heights of 15 mm, 50 mm, and 85 mm, as shown in Figure 1(a).

Experimental System and Scheme.
Compared with the true triaxial experiments in field applications, there is no confining pressure under uniaxial loading conditions. It is no doubt that the confining pressure affects the characteristics of both the propagation of the fractures and the accompanying AE events. [32][33][34]. However, hydraulic fracturing experiments under uniaxial loading conditions can better control associated variables and focus on the mechanism of fracture evolution. Besides, it can also clarify the fracture formation characteristic under different conditions, and reveal the fracture microscopic morphology in guiding field hydraulic fracturing [12]. Hydraulic fracturing experiments under uniaxial loading conditions are used to analyze AE characteristics associated with fracture propagation in various fields including geothermal production, salt well mining, and shale gas exploitation [35,36].
The uniaxial compressive loading and hydraulic fracturing experiments were conducted using a servo-hydraulically controlled deformation apparatus, the GCTS RTR-2000 rock mechanic system at the Institute of Acoustics, Chinese Academy of Sciences. The system can provide a maximum axial force of 4600 kN and a maximum confining pressure of 210 MPa. Figure 1 shows the experimental system and setup for the shale samples. The PAC AE monitoring system was used for the experiments, which is made by the Physical Acoustic Corporation Ltd., USA for the real-time monitoring of AE signals during hydraulic fracturing (HF) experiments. The components of the acoustic emission monitoring system include AE sensors, preamplifiers, and a central system with recording, processing, and display units (Figure 1(d)). The sampling rate of the AE acquisition was 1 MHz, and the bandpass filter was set to 125-750 kHz. The threshold value of the AE signals was set to 35 dB. The two primary parameters obtained using the AE monitoring system included the AE counts and energy, which can reflect the deformation and fracture process of rock samples. The number of times that the signal exceeded the threshold was the AE counts, and the 2 Geofluids total area of the envelope for the AE waveform was defined as the AE energy. Other AE parameters and a typical AE waveform along with the time-frequency characteristics are shown in Figure 2. The experimental scheme included the following five main steps: (1) Sealing Wellbore. A 50 mm long steel pipe was inserted into the injection hole ( Figure 1(a)). Epoxy resin was used to seal the injection section of the shale sample. The sealed sample was preserved for 24 hours to ensure a stable sealing.
(2) Sensor Installation. The heat shrinkage tube was closely bonded to the shale sample, and the AE sensors were placed on the tube surface using honey (Figures 1(a) and 1(c)).
(3) Axial Loading. The specimen was loaded in the axial direction, and the axial stress value was set to 35 MPa, referring to the measured minimum compressive strength. The axial loading rate was 2 MPa/min.
(4) Fluid Injection. Water was injected into the borehole at a rate of 2 MPa/min. When a sudden drop in pore pressure (injection pressure) occurred, this would indicate a break or damage of the sample and the system would shut down.
(5) CT Scanning. After the HF experiments, the fractured samples were removed for micro-CT scanning that produced both 3D images and 2D slices of the cross sections.

Multifractal Method.
The fractal method is used to describe the self-similarity of different local features in irregular phenomena in nature, and it is generally represented by the fractal dimension [37,38]. Fractal analysis has been successfully applied to many branches in petrology and rock physics, such as characterizing the pores and fracture networks in rocks [39][40][41]. To better quantitatively describe the spatial distribution characteristics of complex systems, the multifractal method has been proposed and utilized [27]. For example, the observed AE signals from rock samples could be multifractal. Hence, there is a need for more than one scaling exponent to characterize AE signals and related fractures [42]. The multifractal method has been used to describe the structural characteristics of rock samples using the generalized fractal dimension and multifractal spectrum [5]. In addition, many studies have adopted multifractal theory to study the AE activity of rock fracturing [26,30]. In this study, the box-counting method was used to analyze the multifractal characteristics of the time series of the AE parameters [31,43]. The detailed calculation process of the multifractal spectrum of the AE time series is as follows [31]:

Geofluids
First, a box with a width of L was used to divide the time series into N subsets. The normalized probability of each subset is fP i ðLÞ ; i = 1, 2, 3, ⋯Ng, and the singularity strength is defined as: When the time series has multifractal regularity, the following relationship can be obtained: where is the multifractal spectrum function, which can be also roughly defined as the fractal dimension with the same singularity subset. For calculating the multifractal spectrum of a time series, the following spectrum function was generally used: where τðqÞ is the quality index, q ∈ ð−∞, + ∞Þ; and the quality index, τðqÞ, is the slope of the relationship of ln ðXðq, LÞÞ and ln L under double logarithmic coordinates. Also, τðqÞ can determine the fractal characteristics of the time series: if the value of τðqÞ changes with q, the AE time series has multifractal characteristics. If it is a constant value for a certain q, the time series has a single regularity. Chhabra and Jensen [43] proposed the normalized single parameter measurement cluster as: The multifractal spectrum function, f ðaÞ, and singularity strength, aðqÞ, can be computed by: When the value of q is determined, a ðqÞ and f ðaÞ can be calculated, and we can obtain the multifractal spectrum f ðaÞ − aðqÞ of the time series. The heterogeneity of the time series can be indicated by the spectrum width Δα = α max − α min [30,44]. The larger Δα is, the greater the difference between the large and small signals is. In addition, the value of f ðaÞ represents the frequency at which the singularity strength, aðqÞ, occurs in the time series. Therefore, Δf = f ða max Þ − f ða min Þ indicates the frequency of occurrence of the large and small signals: Δf < 0 means the large-amplitude signals dominate, and the larger-scale fractures dominate the fracture system. In contrast, Δf < 0 means the amplitude/energy of most AE signals is small, and small-scale fractures dominate the fracture system.

Results and Discussion
3.1. Mechanical Behaviors. During the fluid injection process, the pore pressure and AE signals were recorded. Table 2 lists the vertical stress and pore pressure loading on the shale samples. Figure 3 shows the relation of the pore pressure and loading time. Since the mineral composition was basically the same and the samples had no visible preexisting fractures, the results suggest that the bedding structure played a controlling role on the process of rock deformation and fracturing. The results are consistent with Lin et al. [45], and it showed the effects of the anisotropy on fracture propagation in shale samples. As 4 Geofluids the pore pressure increased, the effective normal stress decreased continuously. Once the breakdown pressure was reached, the fracture initiated along the direction of the minimum principal stress and propagated along that of the maximum principal stress, which was in the vertical direction here. The bedding direction of the horizontal samples was parallel to the axial direction, and lower pore pressures were required to break the samples compared with the vertical samples. In addition, the bedding direction of which was nearly perpendicular to the direction of the fracture propagation [46][47][48][49]. The results are also consistent with the principle of the extensively-used horizontal well drilling technology in shale gas exploitation in which wells are drilled along the bedding direction (horizontal) to achieve a smaller breakdown pressure and a larger contacting volume with the target layer [50]. It is worth mentioning that shale anisotropy also has important effects on the field fracturing operation [51]. In the current study, we only focus on testing and validating the effectiveness of the multifractal method, without considering the shale anisotropy. We believe that further studies are required to address the anisotropic issues when generalizing the method to the field scale.

Multifractal Analysis of AE Signals.
There are several studies that have classified the fracturing process based on fracturing experiments. Damani et al. [52] reported that the fracture pressure response can be divided into six broad regions according to the pump pressure and cumulative AE events. Jiang et al. [53] divided the shale supercritical CO 2 fracturing process into the fracture initiation stage and the fracture propagation stage according to the AE parameter curves. Figure 4 shows the relationship of the AE counts, the AE energy, and the loading pore pressure for the horizontal samples (H1, H2, H3). Figure 5 shows the corresponding results of the vertical samples (T1, T2, T3). It is obvious that the AE counts and the AE energy shared a similar trend in all of the samples. In addition, the accumulative counts and energy steadily increased with the loading pore pressure. There is a difference between the detailed plots of the different samples, which was largely caused by the heterogeneity of the samples. As the loading pore pressure increased, the effective normal stress inside the rock sample decreased continuously. When the Mohr stress circle reached the rupture envelope surface, the rock broke down and generated a large amount of AE events [5,54,55]. Therefore, the hydraulic fracturing process can be roughly divided into three stages based on these plots [52,53,56]. First is the initial stage, which corresponds to the beginning of the water injection, where a small number of AE events occurred as the pore pressure increased instantaneously. Then, the quiet stage in which the fracturing system tended to be stable and the AE activity was inactive ("quiet"). Finally, the fracturing stage indicated the breakdown of the rocks and a sudden increase in the AE counts and energy. An important difference between the horizontal and vertical samples was the duration of the quiet stage. According to the dashed green lines in Figures 4 and 5, the duration of the quiet stage (stage II) for the horizontal and vertical samples was approximately 400 s and 700 s, respectively. This is consistent with the larger pore pressure required for the fracturing of vertical sample, which was largely caused by the direction of the bedding planes and discussed in the above section.
The multifractal method was used to further quantify the AE characteristics associated with the hydraulic fracturing experiments. It has been demonstrated that the time series of the AE parameters has multifractal characteristics [31]. The multifractal spectrum of the AE counts and energy time series throughout all stages is shown in Figure 6. The complete results for both the individual and overall stages are listed in Table 3. The results are closely related with the characteristics of fracture evolution. The values of the spectrum width, Δα, of the AE counts for the six rock samples vary between two and four, which means the dispersion of the AE signals was relatively obvious. The primary reason is that the fracturing process can produce both large and small AE events/signals. The multifractal spectrum of the AE energy showed that the Δ α values were between three and five, which shows a consistent result of the AE counts. The entire process of shale hydraulic fracturing can be divided into different stages, including the initiation stage, the quiet stage, and the fracturing stage, and these naturally produce AE events with varying characteristics (e.g., counts and energy).
The different Δα values of the samples with the same bedding direction partially revealed the heterogeneity of the AE parameters/signals and that of the shales. In sample T3       Figures 9 and 10) that resulted from the relatively large injection pressure and the connection between new fractures and weak bedding planes (preexisting fractures). More detailed discussions of the time-frequency analysis and CT scanning are presented in later sections. Figure 11 shows the variations in Δα and Δf values for the different samples at different stages. For the horizontal samples, Δα increased along with loading times and stages.

Geofluids
During the initial stage, the injection pressure was small and only excited small AE events. As the injection pressure increased, the strain energy continuously accumulated. At the fracturing stage, small fractures were extended and connected to form large-scale fractures. Therefore, both small and large AE events were generated, and the difference between small and large signals increased. Several previous studies also revealed that the Δα value rose gradually along with hydraulic fracturing [30,44]. However, for the vertical samples, the results were more complex. A possible reason was the inconsistency in the directions of the principal stress and the bedding planes. And there were few AE events at the quiet stage and led to low Δα values. For samples H2 and T2, the values of Δ α during different stages varied greatly, revealing a difference in the AE activity and/or fracture scales throughout all the stages. In addition, most large-scale fractures and AE events occurred during the final fracturing stage.

Time-Frequency
Analysis. AE waveforms and frequencies contain important information during hydraulic fracturing, including internal fracture size and fracture type [21,57,58]. The microscopic evolution characteristics can be inferred by statistically analyzing the time-frequency characteristics of the AE signals during the fracturing process [35,59]. First, the time domain AE signal was converted into the timefrequency spectrum using a wavelet transform to obtain the main frequency of its signal. AE events were classified to four types: low amplitude in low frequency band (LL), high amplitude in low frequency band (HL), low amplitude in high frequency band (LH), and high amplitude in high frequency band (HH) ( Table 4). The maximum amplitudes and main frequency contents are statistically summarized in Figure 12. It can be seen that a small number of LL-type and LH-type AE events were generated during the initial stage of hydraulic fracturing, and the amplitude was on the order of 10 -2 . As the     injection pressure gradually increased, the effective principal stress decreased continuously, and strain energy was accumulated. LL-type AE events were continuously generated during this stage. Once the breakdown pressure was reached, the crack initiated and propagated to form small-scale fractures, producing more LH-type AE events. Simultaneously, smallscale cracks might have been connected to each other, forming large-scale fractures and producing HL type AE signals. Therefore, large-scale and small-scale fractures coexisted simultaneously during the fracturing stage. This was consistent with the variations of the multifractal spectrum values Δ f at different stages and further verified the reliability of the results. And the HL-and LH-type AE signals were used as precursors before samples cracking [35]. Though the mineral composition of the rock samples is almost the same, the maximum pore pressures of the vertical and horizontal rock samples are also in the same range; the energy and the count values of AE events show a large difference. It is obvious that the anisotropy of shale seriously affects its physical properties. For the two horizontal rock samples H2 and H3, the number of events in the first and second stages is relatively small. How-ever, in the third stage, the HH-type and the HL-type AE events of rock sample H3 are much larger than those of H2.
3.4. CT Scanning. The internal structural characteristics of the horizontal sample H2 and the vertical sample T2 were studied using micro-CT scanning. Figures 7 and 9 show the overview of the micro-CT scanning results of rock samples H2 and T2. Figures 8 and 10 show four horizontal slices of the scanning results. The fractures initiated in the vicinity of the injection hole and continuously propagated along the axial direction (the vertical direction, which is also the direction of the maximum principal stress). The fractures were finally connected and formed a fracture plane, which was also nearly parallel to the axial direction. The morphological characteristics of the fractures are closely related to the confining pressure. When the axial pressure is much greater than the confining pressure, the fracture will propagate perpendicular to the direction of the local minimum principal stress, along with the borehole wall [60,61]. The different scales of the cracks of the two samples were also consistent with the results of the multifractal spectrum values.    14 Geofluids  For the horizontal sample H2, the scale of the fractures was small, and only a few fractures were visible. The fractures remained in the vicinity of the injection hole and did not propagate vertically far away from the injection hole (see Figure 8). The fracture widths corresponded to heights of 20 mm and 40 mm and were relatively large, being 0.1-0.15 mm and 0-0.1 mm, respectively [62]. There were no visible fractures in the horizontal slices corresponding to the heights of 60 mm and 80 mm. A possible reason is that the injection (pore) pressure was small, and only smallscale fractures were produced. To quantitatively study the morphological characteristics of the fractures, an image   16 Geofluids processing software (Image J) was used to calculate the length and curvature of the visible fractures in the horizontal slices.
The detailed results are listed in Table 5. The tortuosity was defined as the ratio of the total fracture length to the two direct lengths of the ends of the fracture [63]. The tortuosity of the visible fractures was small, and the fractures penetrated the rocks nearly straightly. For the vertical sample T2, the fracture width and scale were larger than those of horizontal sample H2. This was a result of the relatively larger maximum injection pressure. Fracture width gradually decreased from the bottom (0 mm) to the top (100 mm). The reason was that the fractures initiated at the bottom portion of the sample, where the injection hole was located. Another important and distinct characteristic was the nearly horizontal fracture plane at a height of approximately 80 mm, which was potentially caused by weakness of the bedding planes and/or nonvisible preexisting fractures. The fracture length, width, and tortuosity of samples H2 and T2 are quite different. The horizontal sample H2 is dominated by small-scale fractures, while the vertical sample T2 is dominated by large-scale fractures.

Conclusions
Hydraulic fracturing experiments were conducted under a uniaxial loading condition on shale samples from the Niutitang Formation in Yiyang, Hunan Province (China). The AE signals were recorded during the fracturing process, and the multifractal method was used to characterize fracture initiation and propagation. The study indicated that the multifractal method can be used to quantitatively characterize hydraulic fracturing and can contribute to optimize in situ hydraulic fracturing operations. The results were further verified using a timefrequency analysis of the AE signals and micro-CT scanning of the samples. The primary conclusions are listed below.
(1) The direction of the bedding planes largely affected the required injection pressure that fractured the rock. The required pressure needed for rocks with bedding planes parallel to the direction of maximum principal stress was less (nearly half) than those with vertical bedding planes under essentially the same conditions.
(2) The hydraulic fracturing process could be divided into three stages based on the characteristics of AE signals: an initial stage, the quiet stage, and the fracturing stage. During the initial stage, a small number of AE events were excited due to the instantaneous change in pore pressure. Also, the strain energy accumulated and a few AE events occurred during the quiet stage. During the fracturing stage, the fractures initiated and propagated promptly and produced a large amount of AE events. A multifractal analysis of the time series of the AE counts and energy quantitatively characterized the heterogeneity of the AE signals during hydraulic fracturing. The value of spectrum width, Δα, continued to increase as the energy accumulated until the fracturing stage began. The Δf value reflects the relationship between small and large signal frequencies and quantifies the fracture scale. Hence, the lower the Δf , the larger the fracture scale, and vice versa.
(3) The time-frequency analysis and CT scanning results further demonstrated the different magnitudes and scales of the AE signals and the fractures. Most AE events and large fractures were produced during the fracturing stage. The required fracturing pressure for the vertical samples was larger, while more largescale fractures and complex fracture planes might be generated due to the existence of bedding planes that vertically or obliquely cross the maximum principal stress. The time-frequency characteristics can verify the change of the multifractal spectrum parameter Δ α. And the LH-and HL-type AE signals can be used as precursors of rock failure. CT scanning results were consistent with the multifractal spectrum parameter Δf . These results further verified the stage classification and multifractal results of the AE signals.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.