Bing Q.Li, Michel Csnov, Herbert H.Einstein

a Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, Cambridge, USA

b Department of Civil and Environmental Engineering, Western University, London, Canada

c School of Architecture, Civil and Environmental Engineering, école Polytechnique Fédérale de Lausanne, Lausanne, Switzerland

Keywords:Acoustic emission (AE)Hydraulic fracturing Enhanced geothermal system (EGS)Microseismicity

ABSTRACT Hydraulic fracturing is frequently used to increase the permeability of rock formations in energy extraction scenarios such as unconventional oil and gas extraction and enhanced geothermal systems(EGSs).The present study addresses uncertainties in the hydraulic fracturing process pertaining to EGSs in crystalline rock such as granite.Specifically, there is debate in the literature on the mechanisms (i.e.tensile and/or shear)by which these fractures initiate,propagate,and coalesce.We present experiments on Barre granite with pre-cut flaws where the material is loaded to high far-field stresses close to shear failure, and then the fluid pressure in the flaws is increased to move the Mohr’s circle to the left and observe the initiation and propagation of fractures using high-speed imaging and acoustic emissions(AEs).We find that the hydraulic fractures initiate as tensile microcracks at the flaw tips, and then propagate as a combination of tensile and shear microcracks.AE focal mechanisms also show elevated levels of tensional microfracturing near the flaw tips during pressurization and final failure.We then consider a numerical model of the experimental setup, where we find that fractures are indeed likely to initiate at flaw tips in tension even at relatively high far-field stresses of 40 MPa where shear failure is generally expected.

1.Introduction

Hydraulic fracturing is frequently used to increase the permeability of rock formations.This can be realized by creating new fractures as usually done for hydrocarbon extraction, or extending and opening existing fractures as usually done in enhanced geothermal systems (EGSs).For EGS, one also speaks about hydraulic stimulation since the idea is to extend and increase the existing fracture network.The main objective of EGS is to extract heat by circulating water through the fractures, usually at depths greater than 5 km,and transform it into electric energy(e.g.Tester et al., 2006; Manzella, 2019).Due to the generally great depths in EGS, only indirect information regarding the fracturing process such as in situ stress,pumping records and microseisms is known.Hence, detailed knowledge on hydraulic stimulation remains limited,as some of the fundamental aspects such as the fracturing processes are still unknown.

Specifically, it is not well established if hydraulic stimulation causes tensile and/or shear fracturing.In this context, one often thinks that fracturing in shearing(hydro-shearing)produces larger fracture apertures through dilatancy (Jaeger et al., 2007; Scholz,2018) but may also be problematic because of induced seismicity.Induced or triggered seismicity in the context of EGS has been addressed extensively in the Basel project (H?ring et al., 2008;Mukuhira et al., 2013), where fault reactivation occurred.These authors as well as Jung(2013)have proposed possible mechanisms such as shearing within cataclastic zones based on observed earthquake nodal plane orientations and multiplet cluster orientations as well as shearing in an en-echelon pattern involving tensile wing fractures of a strike slip fault (Jung, 2013; Cornet,2016).The recent ST1 geothermal pilot project in Helsinki,Finland (Kwiatek et al., 2019) has also shed some light on the geomechanics of deep EGS projects.The stimulation was conducted with a single injection well in locally jointed basement gneiss and granite at 6 km depth, and was instrumented with a seismometerarray to detect smaller-magnitude earthquakes compared to those in the Basel project.The hypocenter locations were predominantly located along ～100 m length scale joints oriented parallel and perpendicular to the maximum principal stress, suggesting again that failure occurred in mixed shear and tensile modes.These mechanisms,where the fluid pressure moves the Mohr circle such that shear failure is induced either in pre-existing fractures or in the intact rock, have been generally described as hydro-shear mechanisms (e.g.Schmittbuhl et al., 2014).

When reviewing recent studies on the topic of hydro-shearing,one finds that some authors suggest that water pressure causes slippage of pre-existing faults (Preisig et al., 2015; De Barros et al.,2016; Ishibashi et al., 2016; Cappa et al., 2019; Li et al., 2021).Others such as H?ring et al.(2008), Jung (2013) and Cornet (2016)believe that large-scale slip occurs through water pressure-induced connection of distinct but initially unconnected en-echelon fractures.

In any of these contexts, it is not quite clear how hydraulically stimulated fractures initiate and propagate.In this study,we intend to take an initial step to clarify the process by which new fractures initiate and propagate between existing non-persistent en-echelon discontinuities under high in situ stresses and elevated pore pressures.We do so by investigating two closely located but unconnected pre-existing fractures, as shown in Fig.1.

The concept of this study is to externally load Barre granite specimens until they reach a state close to shear failure(referred to as critical load(CL)in this study)of the rock matrix(intact rock).By hydraulically pressurizing the lower flaws at constant flow rate,hydraulic pressure will then induce the failure of the material by shifting the Mohr circle to the left (Fig.2a).As shown in Fig.2b,different types of failure are possible depending on the shape of the failure envelope and the location of the Mohr circle(e.g.Aydin et al.,2006).

In this study,we conduct an experiment where axial load,fluid pressure, and injected fluid volume as well as acoustic emissions(AEs) and visual images (using high-resolution and high-speed cameras) are recorded throughout the experiment.These are then compared to a finite element simulation to determine the extent to which the spatial-temporal distribution of fracturing may be explained by stress distribution.Based on this, it is possible to better understand the mechanisms of fracture initiation, propagation, and coalescence.

2.Experimental study

This section provides the necessary background information on the tested material, loading setup, testing parameters and results.

2.1.Material and specimen dimensions

In this study,we use Barre granite from Vermont,USA,which has a grain size of 2-3 mm,and approximate mineralogy of 36.5%plagioclase,31.9%quartz,17.8%K-feldspar,8%biotite,3%muscovite,and2.8%granophyre (Goldsmith et al.,1976).The material is purchased as a 25.4 mm thick slab,and is cut with a wet-saw to a rectangular prism with dimensions of 25.4 mm × 50.8 mm × 101.6 mm (Fig.3).A waterjet is then used to cut two flaws with width of 0.032 inch (1 inch = 2.54 cm) (Fig.3) in the specimen, which are intended to simulate two unconnected fractures,as illustrated in Fig.3.The mechanical properties are given in Table 1.They are used to determine the initial axial stress,which should be close to failure stress,in the experiment and numerical simulations.We assume a Coulomb strengthcriterionforshearfractures,andaGriffithcriterionfor tensile fractures,as presented in Fig.4.Note that the in situ state of stress in the subsurface is generally triaxial, where the differential principal stress σ1-σ3is considered as the key factor governing shear failure(Hoek and Brown,1980).The uniaxial loading case in our experiment is intended to represent the differential stress corresponding to deep enhanced geothermal reservoirs,although it does not include effects of the intermediate principalstress.Importantly,we alsonote thatour experimental setup isintended to shed insighton the phenomenonof fracturing through intact rock bridges rather than shearing along fully persistent discontinuities.

Table 1Mechanical properties of the Barre granite used in this study.

Table 2External axial load and flow rates applied during each of the pressurization cycles for the experiment.All stages are conducted with an injection rate of 0.019 mL/s.

2.2.Experimental setup

The experimental setup is shown in Fig.5a, and schematically presented in Fig.5b.The granite specimen is instrumented with AE sensors (red) placed within inserts in the steel platens (green),which transfer load from the Baldwin 200 Kip loading machine to the specimen.The bottom flaw is sealed (purple square) using a hydraulic fracturing enclosure developed at Massachusetts Institute of Technology (MIT) (see Fig.5c).Hydraulic oil (viscosity μ=3.89 cP)is injected into the bottom flaw at a selected rate using a pressure volume actuator (PVA) with a maximum pressure of 8.5 MPa.We assume no flow of the oil into the granite matrix,given the low permeability of Barre granite (10-6-10-7D from Kranzz et al., 1979).This is supported by discrete element modeling(DEM) of hydraulic fracturing in granite (Kong et al., 2021, 2022),which shows minimal pressure change in the intact rock matrix.These DEM results also show that microcracks directly connected to the fluid source carry a similar fluid pressure as the borehole,while distal en-echelon shear microcracks, which are triggered only by remote stress changes, exhibit minimal pressure increase.The experimental setup (Fig.5a) also includes a high-speed camera(Phantom v2512 with Tamron SP AF 90 mm f/2.8 Di MACRO 1:1 lens)operating at 3000 fps and a high-resolution camera at 0.5 fps(Canon 6D with Tokina AT X Pro D Macro 100 mm f/2.8 lens) to capture the short- and long-term fracture development, respectively.The imaging allows us to determine whether individual microcracks initiate and propagate in shear or tension.This is true even when the macrocrack may appear to be mixed-mode at the specimen scale (Wong and Einstein,2009).

Fig.1.Specimen with pre-existing fractures,so-called flaws,subjected to different external stresses and flaw pressurization conditions.Fluid pressure acts on the bottom flaw in the present study.

Fig.2.(a)Illustration of the experiment showing the Mohr circles at the initial state of stress,after application of axial load to a critically loaded state close to the failure envelope,and at failure induced by a minor increase in hydraulic pressure;and(b)Different types of failure depending on the shape of the failure envelope and the location of the Mohr circle.

2.2.1.AE setup

Fig.3.(a)Flaw pair geometry notation used in this paper;and(b)Dimensions and Gra-60-0 flaw pair geometry of the Barre granite specimens used for the experiments in this study.Abbreviation Gr refers to granite,a is the flaw half-length,L=a is the bridge length between flaws, β is the flaw inclination angle, and α is the bridging angle.

As shown in Fig.5b,the experiment is instrumented with eight Micro30S sensors from Physical Acoustic Corporation, USA, which sample at 5 MHz using four PCI-2 data acquisition cards.The sensors are coupled with honey to the specimen, and sensors in the load path are emplaced within an inset in the loading platens.The pre-amplification is set at 40 dB using the 2/4/6 pre-amplifier from Physical Acoustic Corporation.The sensors exhibit a frequency range of 125-400 kHz,with a resonant frequency of 225 kHz.Event locations are calculated in granite using a P-wave velocity of 4500 m/s and error tolerance of 1 mm.In our algorithm, the location error is defined,per channel,as the P-wave velocity multiplied by the difference between the observed arrival time and the predicted travel time between the inverted source location and the sensor.The location error must be less than the error tolerance on all detected channels, and be detected by a minimum of fourchannels to be an event.The P-wave velocities are measured under experimental conditions.Moment magnitudes are calibrated using the ball drop method developed by McLaskey et al.(2015).Twodimensional (2D) moment tensor inversion is implemented according to the SiGMA method (Shigeishi and Otsu,1996).Focal mechanisms are decomposed as specified in Li and Einstein (2017,2019).The overall AE acquisition system is modified to record close to continuously as necessary (Li et al., 2015).

Fig.4.Illustration of the combined Griffith-Coulomb failure envelope using the parameters described in Table 1.The envelope is used to determine the initial loads in the experimental program and the numerical results in Section 3.2.

Fig.5.Schematic diagram of experimental setup: (a) Photograph of experimental setup.(b) Schematic representation of experimental setup.Black lines denote specimen boundaries, red the AE sensors, blue the fluid injection device and flaw into which fluid is injected, and green the loading platens with inserts for AE sensors.The orange dashed inset specifies the dimensions of the pre-cut flaws.The specimen is 25.4 mm thick.(c)Illustration of the single flaw pressurization device used for hydraulic fracturing(taken from AlDajani et al., 2019).

Fig.6.Schematic diagram illustrating the concept of the experimental loading scheme.The rock is initially loaded to the starting axial load of 40 MPa, at which the pressure is increased up to the maximum allowable PVA pressure of 8.5 MPa.If no failure is observed,the axial load is increased by 5 MPa or 10 MPa,and the fluid pressure is increased again up to 8.5 MPa.The axial load and fluid pressure are alternately increased until we observe failure by hydraulic stimulation.Dashed red or green lines indicate that multiple load/pressure cycles are required until the rock reaches failure.

2.2.2.Load path

The experiment is a uniaxial test where we bring the rock to a stress state close to failure via external stresses and induce failure using fluid pressure, as shown schematically in Fig.6a.

Given the failure envelope shown in Fig.4, we estimate the uniaxial compressive strength at 65 MPa under dry conditions.To bring the rock to failure using fluid pressure, we set the initial uniaxial load to 40 MPa and inject fluid at 0.02 mL/s until we reach the maximum fluid pressure of 8.5 MPa allowed by the equipment,and if no failure is observed we increase the axial load by an increment of 5 MPa or 10 MPa at a rate of 5000 N/min.We continue to alternate the increase in fluid pressure and axial load until failure is observed(Fig.6)or the pressure cannot be increased anymore in the pressurized flaw due to leakage.

2.3.Experimental results

2.3.1.Load data

In this experiment, starting load is set to 40 MPa to limit microfracturing during the initial uniaxial loading, given the material’s UCS of 120 MPa.Table 2 and Fig.7 show the different pressurization cycles with the corresponding loads and maximum pressures reached as well as the data acquired during the experiment.Coalescence between the flaws occurs at uniaxial load of 65 MPa and fluid pressure of 1.68 MPa.Note that the maximum fluid pressure described in Table 2 decreases with increasing axial load, likely due to an increase in the number of microfractures surrounding the pressurized flaw.These microfractures increase the permeability around the flaw, resulting in leakage of the injected fluid.

Table 3Overview of the uniaxial loading cases discussed in this section.

2.3.2.High-speed imaging data

Fig.8(additional explanation given in Fig.A1 in the Appendix A)shows the visual observations obtained using the high-speed video to identify various fracture initiations and propagations during a 1.2 s window before final failure, which is defined as the time of Sketch 7 shown in Fig.7.Note that these high-speed video frames correspond to the very end of the experiment, as shown by the black X in Fig.7.We can see that wing fractures have initiated from the outer flaw tips by the time the high-speed video begins(Sketch 1),likely due to the application of uniaxial load.We first see white patching in the bridge between the flaws, indicating multiple microfractures corresponding to the process zone.White patching is the visually observed whitening of the specimen surface caused by multiple microcracks,as explained in Wong and Einstein(2009).The first fracture on the pressurized (left) inner flaw tip occurs in tension as predicted by the numerical simulation (discussed in Section 3.2)at approximately 1.2 s before the final failure.We then see a second tensile fracture from the non-pressurized inner flaw tip, followed by a complex series of tensile and shear fractures shown in Sketch 3 onwards, culminating in coalescence between the flaws in both tensile and shear modes in Sketches 5 and 6,respectively.Note that we find that the front and back faces of the specimen show very similar fracture patterns,which indicates that the fractures do not change significantly in the out-of-plane direction, as seen in other studies such as that by Morgan and Einstein (2017).

This section discusses the rate, hypocenter locations, and focal mechanisms of the AE detected throughout the experiment.The Gutenberg-Richter frequency-magnitude plot for the entire experiment can be found in Fig.A2.Fig.9 shows the locations of AE hypocenters in the external loading phase, which are clustered at all four flaw tips,while the high-resolution imaging showed no visible fracturing at the end of the loading phase.This suggests that while the fracture initiation and propagation shown in Fig.8 are generated via hydraulic pressure, the preceding external loads nevertheless induce zones of weakness in the rock surrounding the flaw tips.

In the pressurization cycles,it is interesting to look at the hit rate vs.time plot of the AE data(Fig.10).Note that the flow rate is kept at 0.019 mL/s throughout all the pressurization cycles.The number of hits over time curve presented in Fig.10a shows that during the four different pressurization cycles, higher numbers of hits are captured in the first and the last cycles.At the beginning of the first pressurization cycle, pressure increases up to 5.02 MPa (at t = 1355 s) where the pressure curve starts to be nonlinear and eventually decreases at 5.55 MPa.As in previous experiments (Li et al., 2019), AE data confirm that hydraulic stimulation occursduring this first pressurization cycle at the point where the pressure curve starts its nonlinear trend (1350 s in Fig.7a).The hypocenter locations (Fig.10g) show a dense cluster around top of the pressurized flaw tip, which is reasonably aligned with the fluid pressure-driven fracture observed visually (see Sketch 1 in Fig.8).In addition to the AE events that can be related to fluid pressuredriven fracture propagation, further AE events are observed at the middle right boundary of the specimen.

Fig.7.(a)Pressure and injected volume over time,with green boxes highlighting individual pressure cycles;(b)External axial load over time;and(c)Close-up of pressure cycle 4,where black crosses denote sketches shown in Fig.8.

Fig.8.Sketches of the visual observations made using the high-speed images during the last 2 s before failure.The sketches illustrate the different fracture mechanisms and their time of occurrence.A-F represent the order of initiation of the fractures.T and S denote the tensile and shear fractures, respectively (Morgan et al., 2013).The same reference also explains“white-patching” as the occurrence of multiple microcracks corresponding to the process zone.The commented version of this figure is included as Fig.A1 in the Appendix A.

During the second and third pressurization cycles (Fig.10a), no events could be related to hydraulic stimulation since no clusters of events are detected around the pressurized flaw (Fig.10d and e).The fact that the pressure cannot be increased to higher values and the absence of AE events mean that the fluid is able to pass through the existing fracture network without the creation of further microfractures.Nevertheless, in the second and third pressurization cycles,a dense cluster of events develops at the right boundary of the specimen (Fig.10e, circled in red).

The drop in axial load during the last external load step(purple time span) as well as the final failure (magenta time span) shows high magnitudes of AE events(marker size).During the short drop(～2284 s,Fig,10b)in the axial load during the last load increment from 60 MPa to 65 MPa,a large fracture forms(top right flaw tip in Sketch 1 in Fig.8)exactly in the area where high-magnitude events were observed (shown in Fig.10f).

The AE rate shown in Fig.10a raises possible questions about the effect of brittle creep on the failure observed during the experiment.As discussed in Brantut et al.(2013), brittle creep can lead to timedependent cracking of rocks at stresses close to the strength of the material.The time-to-failure of these tertiary creep processes can vary by many orders of magnitude from 10 s to 106s, affected byvarious factors such as humidity, pore fluid chemistry, confining stress, and temperature.As discussed by Heap et al.(2009) and Lockner and Byerlee (1978), the secondary creep process is accompanied by a constant rate of AE events,followed by accelerating AE rate during tertiary creep.These AE rates reflect the creep strain rate of the rock specimen.Fig.10a shows that the AE rate decays during all periods of constant σ1, and only accelerates when the pressure reaches a critical value of 1.68 MPa as shown in Fig.8,which suggests that brittle creep is not responsible for the observed cracking.

Fig.9.Illustration of the external loading phase divided into two time spans (first and second halves of the loading phase).The AE data of the hypocenter locations show dense clusters at locations where first fractures were observed visually.Pressurized flaw outlined in blue.

Fig.10.AE event locations throughout the experiment: (a) AE data showing the number of hits over time.(b) Illustration of the last pressurization cycle (during which failure occurred).(c-e) Hypocenter locations and focal mechanisms shown for the first three pressurization cycles.Red circles show that for each of the pressurization cycles, the magnitude of focal mechanisms becomes larger (marker size) and that the cluster grows.Blue circles indicate the fractures that initiate from the pressurized flaw.(f, g) AE event hypocenter locations and focal mechanisms.The color code for (c-g) indicates the relative time at which the events occur within their observation period.

It is additionally interesting to look at the ratios of the focal mechanisms for the loading phase(up to 1000 s of the experiment)before pressurization as well as for the pressurization phase(after 1000 s of the experiment; phases shown in Fig.7).Given that different sequences of failure are detected throughout the experiment,i.e.hydraulic stimulation in the first pressurization cycle,prefailure inducing a short drop in the axial load, and final failure,these three different sequences are also presented with respect to their relative ratios of focal mechanisms.

Fig.11.Relative proportions of AE focal mechanisms throughout the experiments:(a)Distribution of focal mechanisms throughout the entire experiment,with the loading stage up to 1000 s, compared with the pressure cycles after 1000 s(see Fig.7 for details on loading and fluid pressure histories).(b)Distribution of focal mechanisms for the first pressure cycle (Fig.10a), the decrease in axial load during the fourth pressure cycle, and during final failure of the rock.DC: Double couple, representing AE events where the dominant kinematic mechanism is in shear;NDC+:Non-double couple events where the dominant kinematic mechanism involves volume expansion;NDC-:Non-double couple events where the dominant kinematic mechanism involves volume contraction.

Fig.12.(a)Illustration of the mesh used for Gr-a-60-0 flaw pair;(b)Zoom of the inner bottom flaw tip with the black circle denoting the location where stresses are measured;and(c)Illustration of the notation related to the black circle indicating where stresses are measured showing the angle θ used to describe the location of the stresses measured and the colored circles indicating the location at which the stresses will be discussed in more detail.

Fig.13.Example of Gr-a-60-0 flaw pair geometry subjected to a vertical load of 10 MPa and a bottom flaw pressurization of 5 MPa:(a)Stress field of the minor principal stresses where the two squares denote the locations analyzed and discussed in terms of local vs.global failure;(b) Illustration of Mohr circles in locations shown in (a) and (c) where the local failure around the flaw tip boundary(one radius away from the tip)and the global“matrix”stresses without failure are represented by the two Mohr circles;and(c)Stress field of the major principal stresses where the two squares denote the locations analyzed and discussed in terms of local vs.global failure.

Fig.11 illustrates that for the relative proportions of focal mechanisms throughout the entire experiment and the different phases studied, the shear events (DC) are dominant.Even though the ratio of tensile events (NDC+) increased during the pressurization phase, notably during the drop in axial load, which occurs during the last load cycle, the analysis shows that tensile events account for less than 20% of all events.This discrepancy between the visually identified tensile fractures and large proportion of double couple contribution in the AE has been previously explored by researchers such as Nasseri et al.(2006) and Li and Einstein(2017, 2019), who showed that visually observed tensile fractures in crystalline rocks such as granite very often consist of en-echelon shear fractures at the grain scale,which are more representative of the length scales of the magnitude-6 to-8.5 AE events.However,we do note that the increase in the number of NDC+events towards failure suggests coalescence of microfractures into the visually observed tensile fractures.

3.Comparison with numerical simulations

Here,we use a simplified finite element code (linear-elastic) to better understand fluid pressure-driven fracturing mechanisms,in general, and for the laboratory results of this study, in particular.Specifically, the numerical simulation aims to obtain the stress distribution (just before fracture initiation) in Barre granite, when using a flaw pressurization device in laboratory studies similar to what was done by Gon?alves da Silva and Einstein(2014,2018).The specific conditions simulated and also used in the experiments are uniaxial external loading and pressurization of the bottom flaw.

3.1.Numerical model setup

The numerical modeling is performed using the Partial Differential Equation(PDE)Toolbox based on MATLAB.The model chosen in the PDE Toolbox is the “structural, static-plane-stress” analysis,in which the material is considered linearly elastic.It is important to be aware that this analysis is only valid before any fractures initiate and does not consider any plastic behavior of the material.This is particularly true near the flaw tips, where tensile stresses can readily exceed the tensile strength.We note that our numericalmodel is employed only to identify locations of microcrack initiation (e.g.bottom flaw), rather than correctly model the rock’s inelastic stress-strain behavior (Casanova, 2019).The material parameters described in Table 1 are estimated according to observations made by Hawkes et al.(1973), Miller and Einstein (2008)and Goldsmith et al.(1976).With these parameters, the combined Griffith-Coulomb failure envelope as shown in Fig.4 is constructed.

Fig.14.(a)Vertical load of 10 MPa and no flaw pressurization in the numerical model showing the locations studied in terms of the first failure mechanisms around the flaw tip in different colors; (b) Numerical outputs of the major principal stresses and the maximum shear stress represented in the same colors to match the locations shown in (a); and (c)Mohr circles corresponding to the principal stresses and the maximum shear stresses shown in (b).Colors correspond to the angular location around the flaw tip between 60°and -100°.

Fig.15.(a) Vertical load of 10 MPa and bottom flaw pressurization of 5 MPa simulated in the numerical model showing the locations studied in terms of first failure mechanisms around the flaw tip in different colors; (b,c)Numerical outputs of the major principal stress and the maximum shear stress represented in the same colors to match the locations shown in(a)for the non-pressurized and pressurized flaws,respectively;and(d,e)Mohr circles corresponding to the major principal stress and the maximum shear stress for the non-pressurized and pressurized flaws shown in (b) and (c), respectively.Colors correspond to the angular location around the flaw tip between 60° and -100°.

For the uniaxial analysis presented in this paper, the following boundary conditions are applied in the MATLAB PDE Toolbox model: The elements at the bottom horizontal boundary are fixed and the elements along the vertical boundaries are limited to move only in the vertical direction.

A zoom of the mesh is shown in Fig.12a.The mesh has been refined around the flaw tips allowing us to observe the stresses around the flaw tips more in detail.

The black circle in Fig.12b around the flaw tips(at a distance of 2r from the center of the flaw tip and equal to a distance of r from the flaw tip) illustrates where the stresses are determined and compared.The locations on the black circle are defined by the angles of the unit circle starting at the flaw tip where θ = 0°.Extensive studies by Gon?alves da Silva (2009) and Gon?alves da Silva and Einstein (2014) showed that using stresses at this location to determine where fractures initiate fit the experimentally observed fracturing very well.

Again, it is important to be aware of the limitations of such a numerical analysis.For instance,refining the mesh size for a planestress analysis is limited due to computer storage constraints.Applying high values of external stresses in the model could also produce misleading conclusions since under high load,the material could undergo plastic deformation, which is not considered in a linear-elastic model.

3.2.Numerical results

For comparison with the experimental results presented in Section 2.3, three loading cases are simulated, and the maximum stresses around both the pressurized and non-pressurized flaw tips are elaborated.This can provide substantial information regarding the mechanisms of first fracture initiation as well as on the loading that might favor shear or tensile failure.

The linear-elastic finite element study is used to determine that fracture initiation and propagation can occur.In doing so, one needs to differentiate what happens locally at the flaw tip (one radius away from the tip as shown in Fig.12)and what happens at a distance from the flaw tip in the matrix called “global” in the following.As to be expected, the local stress conditions at the tip are more extreme than those in the global matrix.This is illustrated in Fig.13.

Fig.16.(a) Vertical load of 40 MPa and bottom flaw pressurization of 5 MPa simulated in the numerical model showing the locations studied in terms of first failure mechanisms around the flaw tip in different colors;(b, c)Numerical outputs of the major principal stress and the maximum shear stress represented in the same colors to match the locations shown in(a)for the non-pressurized and pressurized flaws,respectively;(d,e)Mohr circles corresponding to the major principal stress and the maximum shear stress for the nonpressurized and pressurized flaws shown in (b) and (c), respectively.Colors correspond to the angular location around the flaw tip between 60° and -100°.

These local stress conditions determine the potential locations of fracture initiation and can be compared with what is observed in laboratory experiments.

In this paper, only uniaxial experiments are considered; specifically, the three loading cases listed in Table 3 are discussed in detail.

To analyze different points around the flaw tips,major principal stresses and maximum shear stresses around the flaw tips are extracted from the numerical model and are represented using different colors (the locations, i.e.60°, 0°, -60°and -100°, are shown in Fig.12c).The Mohr circles can then be represented using these numerical outputs and be compared to the combined Griffith-Coulomb failure envelope (Fig.4).

When looking at the numerical simulation of Case 1 where no flaw is pressurized and the vertical load is set to 10 MPa,we observe that two of the locations studied reach the tensile failure boundary(Fig.14a).Note that only one flaw is presented in Fig.14 because the stress field is considered symmetric around both flaws with no pressure inside the flaws.These results correspond well with the AE locations prior to pressurization as shown in Fig.9.

In Case 2 in which the bottom flaw is pressurized with 5 MPa and the vertical load is kept at 10 MPa,both locations that reached the tensile failure boundary for the previous loading case (Case 1,Fig.14)now cross the tensile failure boundary,as shown in Fig.15c,for the pressurized bottom flaw.We observe that both of the Mohr circles (at 0°and 60°) have shifted towards the left.Moreover, the change in the local tensile field around the pressurized flaw tip seems to extend up to the non-pressurized flaw, inducing tensile failure at 0°and 60°of the non-pressurized flaw tip.We can thus assume that local tensile fractures form under these loading conditions in these locations and that in this linear-elastic model, the first tensile fracture initiation can be expected at either flaw tip.

Overall, the stresses are not high enough to induce local shear failure,which is why a third loading case is discussed.By increasing the axial load to 40 MPa and keeping the bottom flaw pressure at 5 MPa(Fig.16),we can observe that the locations at-100°for both flaw tips cross the Griffith-Coulomb failure envelope at the shear failure boundary.The location at -60°around both flaw tips still remains intact.

Comparing the three cases provides interesting information:

(1) Locations-60°and-100°in all three cases and for both the non-pressurized and pressurized flaws have positive major and minor principal stresses.This might potentially cause shear failure but only does so in Case 3 with high external(uniaxial)stresses.Note also that in all cases the shear stress is the highest at the flaw tip (0°) (square or x symbols in Figs.14-16).At the locations above/below the tip (60°), the minor principal stresses are always tensile, and the major principal stresses are very often low and thus indicate tensile failure.Considering in detail the evolution of the Mohr circles for Cases 1-3, one can state that tensile failure at locations 0°and 60°will occur before shear failure at-60°and-100°.

(2) The reader will now ask if it would not be possible to study the stress behavior (regarding both external stresses and applied hydraulic pressure) in which we would clearly see under what conditions the Mohr circle just touches the envelope and if this is in tension or shear (or both).Whiledoable,this might be trying to get too much from a relatively simple stress distribution.Essentially, we can state that initial failure under the investigated conditions is most likely to be in tension at or near the flaw tip, and this can be compared to the experimental result.

(3) The numerical model can provide important information on the expected location and failure mode of the first failure around the flaw tips.For the uniaxial laboratory studies conducted on Gr-a-60-0, we can thus expect to observe initial tensile failure starting around the tip of the pressurized or non-pressurized flaw.Also, it seems to be highly unlikely to observe shear failure mechanisms as first failure mode.

Here we emphasize that the numerical study only provides information astowhere a fracture may initiate atthe flaw boundaryand if this is a tensile or shear fracture.The laboratoryexperiments,on the other hand,provide information on a sequence of fracturing events from initiation to coalescence.This means that the comparison between the two studies has to concentrate on initiation essentially comparing Figs.15 and 16 with Sketches 1 and 2 in Fig.8.For the pressurized(lower)flaw,a tensile fracture appears in the experiment but at a somewhat different location than in the numerical study;the white patching, however, matches the location where tensile fractures are simulated.This discrepancy may be attributed to the limitations of the linear-elastic model that does not account for damage accumulation in the rock.For the non-pressurized(upper)flaw,the numerical and experimental results are the same both regarding location and tensile type.It is also not surprising that during final coalescence,the tips of the flaws are experiencing shear failure,given that in the numerical study,the highest shear stresses were observed at the tip of the flaws(square or×symbols in Figs.14-16).

4.Conclusions

In this paper,we try to address the problem of different types of failures that might occur when a specimen containing flaws (preexisting fractures) is subjected to external stresses and hydraulic pressurization of one of the flaws.Of particular interest is the mechanism (type) in which fractures initiate from the flaws, specifically if this occurs in tension or shear.We did this with both experiments and numerical studies.The experimentally observed initial fracturing corresponds reasonably well to the numerical results.The simple elastic numerical study can provide the information on mechanisms and location of the initial fracture in relation to applied external stresses, flaw pressurization and failure criteria.In addition, the experiments provide information on fracture propagation and coalescence,which are observed both visually and with AEs again differentiating between tensile and shear events.As could be seen, these mechanisms are a reasonably complex combination of tension and shearing.The visual and AE observations show results that are generally similar with differences that can be related to the scale of the observations.To conclude, we can state that the elastically based numerical model provides good information on fracture initiation and that the experiments provide detailed information on the overall mechanisms.

Data availability

The data are available on Mendeley Data at https://data.mendeley.com/datasets/wtt9mbzs4k/1.The AE processing code is available on the MATLAB file exchange at https://www.mathworks.com/matlabcentral/fileexchange/72339-mit-rock-mechanics

group-ae-code.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This research was financially supported by the Total SA in France and the Abu Dhabi National Oil Company in Emirate of Abu Dhabi.The authors would like to thank Omar Aldajani for assistance in the laboratory, Professor John Germaine for advice, and Victoire Denis du Péage and Arabelle De Saussure for laying the groundwork for this publication.

Appendix A.Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2022.06.011.