Fracture Propagation and Morphology Due to Non-Aqueous Fracturing: Competing Roles between Fluid Characteristics and In Situ Stress State

: Non-aqueous or gaseous stimulants are alternative working ﬂuids to water for hydraulic fracturing in shale reservoirs, which o ﬀ er advantages including conserving water, avoiding clay swelling and decreasing formation damage. Hence, it is crucial to understand ﬂuid-driven fracture propagation and morphology in shale formations. In this research, we conduct fracturing experiments on shale samples with water, liquid carbon dioxide, and supercritical carbon dioxide to explore the e ﬀ ect of ﬂuid characteristics and in situ stress on fracture propagation and morphology. Moreover, a numerical model that couples rock property heterogeneity, micro-scale damage and ﬂuid ﬂow was built to compare with experimental observations. Our results indicate that the competing roles between ﬂuid viscosity and in situ stress determine ﬂuid-driven fracture propagation and morphology during the fracturing process. From the macroscopic aspect, ﬂuid-driven fractures propagate to the direction of maximum horizontal stress direction. From the microscopic aspect, low viscosity ﬂuid easily penetrates into pore throats and creates branches and secondary fractures, which may deﬂect the main fracture and eventually form the fracture networks. Our results provide a new understanding of ﬂuid-driven fracture propagation, which is beneﬁcial to fracturing ﬂuid selection and fracturing strategy optimization for shale gas hydraulic fracturing operations.


Introduction
Shale gas production has substantial growth in North America due to the advances in horizontal drilling and multi-stage fracturing technology [1][2][3]. However, in some shale formations, impairment of productivity has been observed due to the use of water-based fracturing fluid: (1) water is trapped in the near-wellbore area and hinder the mass transport of oil and gas from the fractures to the wellbores; (2) clay minerals swell after water invades into formations which results to the formation permeability decrease [4][5][6][7]. Moreover, environmental issues, such as groundwater contamination and potentially induced seismicity by wastewater re-injection have been reported [8][9][10]. As an alternative to water-based fracturing fluids, non-aqueous or gaseous fluids, such as nitrogen (N 2 ), liquid carbon

Experimental Research
In this part, we firstly reported the apparatus and sample preparation procedures for fracturing experiments. Subsequently, we described how we performed the fracturing experiments with H 2 O, L-CO 2 , and Sc-CO 2 under various combinations of the stress state. Finally, the fracture morphology is extracted by the CT scanning method.

Apparatus and Sample Preparation
The fracturing experiments are performed with a triaxial loading system that can independently apply axial stress, confining stress, and perform fluid injection for fracturing ( Figure 1). Hand pumps are used to deploy confining stress and axial stress separately, which includes the maximum horizontal stress (σ H ), minimum horizontal stress (σ h ) and vertical stress (σ V ). An ISCO 260D plunger pump (manufactured by Teledyne ISCO, Lincoln, NE, USA) is used to inject fracturing fluid (Sc-CO 2 , L-CO 2 or distilled water) at a constant flow rate for each test.          Figure 2a. A blind central borehole is drilled whose length is 110 mm. The boreholes are perpendicular to bedding planes and analog as fracturing boreholes.
Since the experiments aim to investigate the effect of fluid characteristics and in situ stress on fracture propagation and morphology, the conditions for experiments are summarized in Table 3. It should be noted that even though the stress state cannot fully represent the in situ stress, the variation in horizontal stress difference is enough to investigate the role of fluid characteristics on the fracturing process and morphology [28][29][30][31][32]. We use three types of fracturing fluid, namely water, liquid carbon dioxide (L-CO 2 ), and supercritical carbon dioxide (Sc-CO 2 ) to perform fracturing experiments firstly. A heating unit is used to heat the injection fluid with a temperature controller. In the cases of water and liquid carbon dioxide fracturing tests, the system temperature is set at 25 • C. For supercritical carbon dioxide fracturing tests, the system temperature is set as 60 • C to take the carbon dioxide into the supercritical state. Before each test, the sample will be placed in the tank for 24 h to make the sample has a stable temperature. We inject the particular fluid with a constant flow rate of 30.0 mL/min into the shale sample until the injection pressure markedly decreases due to the initiation of fractures. Experiments 1-3 are used to identify the effect of fracturing fluid properties on the fracture propagation morphology. Moreover, we change the stress state while still using supercritical carbon dioxide in experiments 4-7, to identify the role of in situ stress on the fracture distribution and morphology in shale.
After each fracturing test, the samples are removed from thermotank for the CT scanning experiments. The scanning plane is perpendicular to the borehole direction, which is shown in Figure 2b.

Experimental Results
The relationship between fluid pressure and the elapsed time for seven fracturing experiments are summarized in Figure 3. Under the same stress state, H 2 O returns the highest breakdown pressure (29.25 MPa), followed by L-CO 2 (18.70 MPa), with Sc-CO 2 exhibiting the lowest breakdown pressure  The CT scanning images for different fracturing fluid under the same stress state were extracted. The vertical stress is 12.0 MPa, the maximum horizontal stress is 10.0 MPa and the minimum horizontal stress is 8.0 MPa. Then, to make the fractures clearly observable, we developed MATLAB code to adjust the threshold value of CT images and got the gray images in Figure 4. The images indicate that at the same stress state, H2O creates a smooth main fracture which propagates parallel to the maximum horizontal stress direction. This similar phenomenon has been commonly observed by other researchers [29]. The L-CO2 fluid also creates the main fracture. However, the induced fractures are not as smooth as the fractures driven by H2O. Also, the fracture deflection phenomenon is observed in the direction of minimum horizontal stress and some secondary fractures and branches are also observed for Sc-CO2 fracturing, as shown in Figure 4c. We believe more fractures and a complex fracture pattern could be formed by the Sc-CO2 fracturing method [26]. For L-CO2 fracturing, the fracture morphology between traditional hydraulic fracturing and Sc-CO2 fracturing, the main fracture could be formed with some branches and secondary fractures, but not as much as Sc-CO2 fracturing. Hence, the fluid characteristics play an important role in fracture propagation morphology. The CT scanning images for different fracturing fluid under the same stress state were extracted. The vertical stress is 12.0 MPa, the maximum horizontal stress is 10.0 MPa and the minimum horizontal stress is 8.0 MPa. Then, to make the fractures clearly observable, we developed MATLAB code to adjust the threshold value of CT images and got the gray images in Figure 4. The images indicate that at the same stress state, H 2 O creates a smooth main fracture which propagates parallel to the maximum horizontal stress direction. This similar phenomenon has been commonly observed by other researchers [29]. The L-CO 2 fluid also creates the main fracture. However, the induced fractures are not as smooth as the fractures driven by H 2 O. Also, the fracture deflection phenomenon is observed in the direction of minimum horizontal stress and some secondary fractures and branches are also observed for Sc-CO 2 fracturing, as shown in Figure 4c. We believe more fractures and a complex fracture pattern could be formed by the Sc-CO 2 fracturing method [26]. For L-CO 2 fracturing, the fracture morphology between traditional hydraulic fracturing and Sc-CO 2 fracturing, the main fracture could be formed with some branches and secondary fractures, but not as much as Sc-CO 2 fracturing. Hence, the fluid characteristics play an important role in fracture propagation morphology.
Minerals 2020, 10, 428 5 of 17 pressure (29.25 MPa), followed by L-CO2 (18.70 MPa), with Sc-CO2 exhibiting the lowest breakdown pressure at 14.11 MPa. When injecting Sc-CO2 as the fracturing fluid, the results indicate that with the increase of horizontal stress difference, a decreasing trend of breakdown pressure is observed. The CT scanning images for different fracturing fluid under the same stress state were extracted. The vertical stress is 12.0 MPa, the maximum horizontal stress is 10.0 MPa and the minimum horizontal stress is 8.0 MPa. Then, to make the fractures clearly observable, we developed MATLAB code to adjust the threshold value of CT images and got the gray images in Figure 4. The images indicate that at the same stress state, H2O creates a smooth main fracture which propagates parallel to the maximum horizontal stress direction. This similar phenomenon has been commonly observed by other researchers [29]. The L-CO2 fluid also creates the main fracture. However, the induced fractures are not as smooth as the fractures driven by H2O. Also, the fracture deflection phenomenon is observed in the direction of minimum horizontal stress and some secondary fractures and branches are also observed for Sc-CO2 fracturing, as shown in Figure 4c. We believe more fractures and a complex fracture pattern could be formed by the Sc-CO2 fracturing method [26]. For L-CO2 fracturing, the fracture morphology between traditional hydraulic fracturing and Sc-CO2 fracturing, the main fracture could be formed with some branches and secondary fractures, but not as much as Sc-CO2 fracturing. Hence, the fluid characteristics play an important role in fracture propagation morphology.   Figure 5 shows the role of in situ stress on the fracture propagation morphology. When injecting the same fluid (Sc-CO 2 ), the horizontal stress difference (maximum horizontal stress minus minimum horizontal stress) largely affects the propagation morphology for fluid-driven fracture propagation. When the horizontal stress difference is small, fracture networks are more likely formed, as shown in Figure 5a,b. However, with the increase of horizontal stress difference, the fracture propagates parallel to the direction of maximum horizontal stress direction, as shown in Figure 5c-e. This observation is similar to the previous research that the horizontal stress difference is the most important factor in determining the fluid-driven fracture propagation direction [29].   Figure 5 shows the role of in situ stress on the fracture propagation morphology. When injecting the same fluid (Sc-CO2), the horizontal stress difference (maximum horizontal stress minus minimum horizontal stress) largely affects the propagation morphology for fluid-driven fracture propagation. When the horizontal stress difference is small, fracture networks are more likely formed, as shown in Figure 5a,b. However, with the increase of horizontal stress difference, the fracture propagates parallel to the direction of maximum horizontal stress direction, as shown in Figure 5c-e. This observation is similar to the previous research that the horizontal stress difference is the most important factor in determining the fluid-driven fracture propagation direction [29]. It should be also noted that even though the horizontal stress is 4.0 MPa, the main fracture is observed in Figure 5e, rather than the fracture networks. We can still observe some small branches and secondary fractures. This is the result of natural fractures or micro-fractures which exist before the hydraulic fracturing experiments. Based on the sedimentary characteristics, the natural fractures are well-developed. One limitation we must admit to that it is difficult to locate all natural fractures before the experiment since the CT images are extracted layer by layer.

Numerical Simulation Research
We cannot conduct experiments under all combinations of stress state for different kinds of fracturing fluids. Hence, it is essential to develop a numerical model to further characterize the role of fluid characteristics and stress state on the fracture propagation morphology. In this part, a With the decrease of stress difference, a fracture network with more secondary fractures is more easily induced.
It should be also noted that even though the horizontal stress is 4.0 MPa, the main fracture is observed in Figure 5e, rather than the fracture networks. We can still observe some small branches and secondary fractures. This is the result of natural fractures or micro-fractures which exist before the hydraulic fracturing experiments. Based on the sedimentary characteristics, the natural fractures are well-developed. One limitation we must admit to that it is difficult to locate all natural fractures before the experiment since the CT images are extracted layer by layer.

Numerical Simulation Research
We cannot conduct experiments under all combinations of stress state for different kinds of fracturing fluids. Hence, it is essential to develop a numerical model to further characterize the role of fluid characteristics and stress state on the fracture propagation morphology. In this part, a numerical model coupling rock property heterogeneity, micro-scale damage, and fluid flow is built to run the simulation works and compare with experimental observations.

Governing Equations for Rock Deformation
Under fluid pressure conditions, based on poroelastic theory, the rock deformation could be expressed as the Navier-type equation [23]: where u is the displacement; i, j, and k denote the space coordinates; G is the shear modulus of rock; ν is Poisson's ratio of rock; α is the Biot coefficient; p is pore pressure and f is body force.

Governing Equations for Fluid Flow
The main difference between gas and fluid is compressibility. For conventional hydraulic fluid, like water and oil, the fluid has low compressibility. However, for gas like CO 2 and N 2 , the large compressibility and phase change cannot be neglected.
For slightly compressible fluid flow in the rock, the governing equation could be built based on the mass balance equation of fluid and neglect effect of gravity, which is [22]: where K f is the bulk modulus of fluid; K s is the bulk modulus of rock; Φ is porosity; k is permeability; µ is fluid viscosity; p is fluid pressure; Q s is the gas source. The governing equations for CO 2 (compressible fluid) flow in rocks could be written as: where N = RT(δ 2 ϕ r δδ + 2δϕ r δ + 1) which is a variable related with residual part of Helmholtz free energy, ϕ r , and reduced density of CO 2 , δ, ϕ r δ and ϕ r δδ are first and second derivatives of ϕ r concerning δ. A detailed explanation of this term is given in [22].

Damage Law for Rock
The damage law is defined at the microscopic scale which is detailed in [33]. Tensile damage and shear damage occur when REVs (representative elemental volume) satisfied the maximum tensile stress criteria or shear failure criterion, as defined by and shown in Figure 6 [33,34].
where σ H and σ h are maximum horizontal stress and minimum horizontal stress, f t and f c are axial tensile strength and axial compressive strength; θ is internal friction angle. F 1 and F 2 are threshold functions for tensile failure and shear failure occurs. The elastic modulus of rock decreases with the evolution of damage, which could be expressed as: where E 0 is the initial elastic modulus and D is the damage variable. A constitutive law defining D is given in reference [22,33] and, therefore, is omitted here for brevity. The elastic modulus of rock decreases with the evolution of damage, which could be expressed as: where E0 is the initial elastic modulus and D is the damage variable. A constitutive law defining D is given in reference [22] and [33] and, therefore, is omitted here for brevity.

Numerical Simulation Model and Parameters
For a better representation of the realistic heterogeneity of shale samples, the elastic modulus, Poisson ratio, and strength parameters of each element in the numerical model follow the Weibull distribution. The probability density function of Weibull distribution is expressed as: where u is the elastic modulus, Poisson ratio, and strength for each element, u0 is the average value, m is the homogeneity index parameter which describes the shape of the distribution.
Mathematically, the average value u0 is fixed in the Weibull distribution. Hence, the parameter m defines the degree of material homogeneity. The larger the value of m, the more homogeneous the material properties. Hence, if the sample is regarded as homogenous, the m value is always set as larger than 10, where the inhomogeneity properties will not largely affect the result. In our model, the shale rock samples are always inhomogeneous and the m is set as 5. Other physical, mechanical, and hydrological parameters of rock and fluids are obtained from previous experiments and summarized in Table 4. The geometry of the model is shown in Figure 7.

Numerical Simulation Model and Parameters
For a better representation of the realistic heterogeneity of shale samples, the elastic modulus, Poisson ratio, and strength parameters of each element in the numerical model follow the Weibull distribution. The probability density function of Weibull distribution is expressed as: where u is the elastic modulus, Poisson ratio, and strength for each element, u 0 is the average value, m is the homogeneity index parameter which describes the shape of the distribution. Mathematically, the average value u 0 is fixed in the Weibull distribution. Hence, the parameter m defines the degree of material homogeneity. The larger the value of m, the more homogeneous the material properties. Hence, if the sample is regarded as homogenous, the m value is always set as larger than 10, where the inhomogeneity properties will not largely affect the result. In our model, the shale rock samples are always inhomogeneous and the m is set as 5. Other physical, mechanical, and hydrological parameters of rock and fluids are obtained from previous experiments and summarized in Table 4. The geometry of the model is shown in Figure 7.   The original code is written with MATLAB. After that, the parameters related to the damage are also calculated by MATLAB to gain a better calculation speed. Finally, the related parameters are implemented by COMSOL Multiphysics software (Version 5.2a, COMSOL, Inc., Burlington, MA, USA) to complete the FEM (Finite Element Method) analysis.

Numerical Simulation Results
The numerical simulation results of fracture propagation and morphology for different fracturing fluid under the same stress state is shown in Figure 8. To extend the experimental content, we selected viscous oil, water, liquid carbon dioxide, supercritical carbon dioxide and nitrogen as the fracturing fluids. For all five kinds of fracturing fluid, the main fractures are formed, which propagates to the direction of maximum horizontal stress. However, when we use viscous oil and conventional water as the fracturing fluids, the smooth main fracture is formed without branches and secondary fractures, as indicated in Figure 8a,b. When the injection fluids change to liquid CO2, SC-CO2 and nitrogen gas, fractures initiate from the injection wellbore from several points, then the fractures propagate to the direction of maximum horizontal stress, as shown from Figure 8c-e. Moreover, the propagation morphology has longer pathways for supercritical carbon dioxide and nitrogen gas, which illustrate better propagation ability due to low fluid viscosity.

Numerical Simulation Results
The numerical simulation results of fracture propagation and morphology for different fracturing fluid under the same stress state is shown in Figure 8. To extend the experimental content, we selected viscous oil, water, liquid carbon dioxide, supercritical carbon dioxide and nitrogen as the fracturing fluids. For all five kinds of fracturing fluid, the main fractures are formed, which propagates to the direction of maximum horizontal stress. However, when we use viscous oil and conventional water as the fracturing fluids, the smooth main fracture is formed without branches and secondary fractures, as indicated in Figure 8a,b. When the injection fluids change to liquid CO 2 , SC-CO 2 and nitrogen gas, fractures initiate from the injection wellbore from several points, then the fractures propagate to the direction of maximum horizontal stress, as shown from Figure 8c-e. Moreover, the propagation morphology has longer pathways for supercritical carbon dioxide and nitrogen gas, which illustrate better propagation ability due to low fluid viscosity.  Figure 9 illustrates how horizontal stress difference affects fracture propagation morphology. When the horizontal stress is 0 and 1.0 MPa, fractures initiate from several points near the fluid injection wellbore and fracture networks are formed, as shown in Figure 9a,b. Moreover, the rock was damaged before the fracture across the whole sample. When the horizontal stress difference increases to 2.0 MPa, two fractures initiate near the injection wellbore; one propagates parallels to  Figure 9 illustrates how horizontal stress difference affects fracture propagation morphology. When the horizontal stress is 0 and 1.0 MPa, fractures initiate from several points near the fluid injection wellbore and fracture networks are formed, as shown in Figure 9a,b. Moreover, the rock was damaged before the fracture across the whole sample. When the horizontal stress difference increases to 2.0 MPa, two fractures initiate near the injection wellbore; one propagates parallels to the direction of maximum horizontal stress and the other one's propagation direction is relatively random, as shown in Figure 9c. One thing should be noted is that the fracture which propagates to the direction of maximum horizontal stress has longer pathways than the other fracture. This implies the role of stress on fracture propagation. When the horizontal stress difference changes to 3.0 MPa, a fracture perpendicular to the maximum horizontal stress direction is imitated near the injection wellbore. However, the fracture only propagates for a small distance, due to the stress restriction. When the horizontal difference is large at 4.0 MPa, a relatively smooth fracture is formed and tends to propagate parallel to the maximum horizontal stress direction.

Discussion
To further compare the experimental observations and numerical simulations results, we summarized them in combing figures, as presented in Figures 10 and 11. Figure 10 compares the fracture morphology based on experimental observations and numerical simulations where water, liquid CO2, and supercritical CO2 are injected into shale samples under vertical stress at 12.0 MPa, maximum horizontal stress at 10.0 MPa and minimum horizontal stress at 8.0 MPa. The fracture morphology images are similar between experimental results and simulation results when injecting water into shale samples, where a smooth main fracture is formed and fracture propagates to the direction of maximum horizontal stress. When the injection fluid is liquid CO2 and supercritical CO2, we can also observe the main fractures which propagate to the direction of the maximum horizontal stress direction. However, a fracture initiates from several weak points near the fluid injection borehole when we inject liquid carbon dioxide and supercritical carbon dioxide. As shown in Figure 10b, the fracture initiates from three points, which show an inverted "Y" shape fracture type. In Figure 10c, four (two groups) fractures are formed near the fluid injection borehole, which shows an "X" shape transverse fracture pattern. Furthermore, the apparent fracture width of supercritical carbon dioxide fracturing is larger than those of the fractures induced by water.

Discussion
To further compare the experimental observations and numerical simulations results, we summarized them in combing figures, as presented in Figures 10 and 11. Figure 10 compares the fracture morphology based on experimental observations and numerical simulations where water, liquid CO 2 , and supercritical CO 2 are injected into shale samples under vertical stress at 12.0 MPa, maximum horizontal stress at 10.0 MPa and minimum horizontal stress at 8.0 MPa. The fracture morphology images are similar between experimental results and simulation results when injecting water into shale samples, where a smooth main fracture is formed and fracture propagates to the direction of maximum horizontal stress. When the injection fluid is liquid CO 2 and supercritical CO 2 , we can also observe the main fractures which propagate to the direction of the maximum horizontal stress direction. However, a fracture initiates from several weak points near the fluid injection borehole when we inject liquid carbon dioxide and supercritical carbon dioxide. As shown in Figure 10b, the fracture initiates from three points, which show an inverted "Y" shape fracture type. In Figure 10c, four (two groups) fractures are formed near the fluid injection borehole, which shows an "X" shape transverse fracture pattern. Furthermore, the apparent fracture width of supercritical carbon dioxide fracturing is larger than those of the fractures induced by water.
Besides the similarities between the experimental observations and numerical simulation results, some differences can also be found. In Figure 10c, the fracture deflection is much sharper than numerical results, which almost propagate to the direction of minimum horizontal stress. We contribute this phenomenon to the potential bedding planes. In our numerical simulation model, the shale heterogeneity is considered by introducing Gaussian distribution of element mechanical properties, namely elastic modulus and Poisson's ratio. However, the potential pre-existing natural fractures and bedding planes in transversely isotropic shale may affect the hydraulic fracture propagation. In future works, we will consider this parameter. Besides the similarities between the experimental observations and numerical simulation results, some differences can also be found. In Figure 10c, the fracture deflection is much sharper than numerical results, which almost propagate to the direction of minimum horizontal stress. We contribute this phenomenon to the potential bedding planes. In our numerical simulation model, the shale heterogeneity is considered by introducing Gaussian distribution of element mechanical properties, namely elastic modulus and Poisson's ratio. However, the potential pre-existing natural fractures and bedding planes in transversely isotropic shale may affect the hydraulic fracture propagation. In future works, we will consider this parameter. Figure 11 compares the fracture morphology where the horizontal stress differences range from 0 to 4.0 MPa when injecting supercritical carbon dioxide into shale samples. There is also good consistency between the experimental observations and numerical results. When the horizontal stress difference is low (0 or 1.0 MPa), the fracture initiate from weak points near the fluid injection borehole and the fracture networks are formed easily, as shown in Figure 11a,b. In those circumstances, in situ stress is not the determining factor in controlling the fracture propagation. With an increase in horizontal stress difference, the main fractures are observed which propagate to the direction of maximum horizontal stress. Moreover, the main fracture is smoother and there is less deflection when the horizontal stress differences increase. Also, the effect of bedding planes and natural fractures on the fracture propagation is observed in the experimental research, while not being taken into consideration in our numerical simulation model. The comparisons between experimental observations and numerical simulation results show that the competing roles between fluid properties and horizontal stress difference in determining the fracture morphology, even though the natural fractures and pre-existing micro-fractures are not taken into account.  Figure 11 compares the fracture morphology where the horizontal stress differences range from 0 to 4.0 MPa when injecting supercritical carbon dioxide into shale samples. There is also good consistency between the experimental observations and numerical results. When the horizontal stress difference is low (0 or 1.0 MPa), the fracture initiate from weak points near the fluid injection borehole and the fracture networks are formed easily, as shown in Figure 11a,b. In those circumstances, in situ stress is not the determining factor in controlling the fracture propagation. With an increase in horizontal stress difference, the main fractures are observed which propagate to the direction of maximum horizontal stress. Moreover, the main fracture is smoother and there is less deflection when the horizontal stress differences increase. Also, the effect of bedding planes and natural fractures on the fracture propagation is observed in the experimental research, while not being taken into consideration in our numerical simulation model.
The comparisons between experimental observations and numerical simulation results show that the competing roles between fluid properties and horizontal stress difference in determining the fracture morphology, even though the natural fractures and pre-existing micro-fractures are not taken into account.
To quantitatively analyze the fracture propagation morphology, the tortuosity parameter is calculated as [22,25]: where C is fracture tortuosity, L is the total fracture length of a flow pathway; L 0 is the linear length between the two ends of the fracture pathway.
The CT scanning images are digitized by self-developed MATLAB code to calculate the tortuosity parameter. Moreover, the numerical results are also used to compare with the experimental observations. The fluid viscosity of H 2 O, L-CO 2 , Sc-CO 2, and N 2 will change with the change of fluid pressure. Hence, we plot the viscosity evolution of four kinds of fluid in Figure 12. We can find the four kinds of fluid show distinct fluid viscosity properties. Moreover, due to the state change from gaseous to liquid, we can also observe a jump of fluid viscosity for CO 2 under the temperature of 25 • C. Furthermore, we summarize the fracture tortuosity against the fluid viscosity in Figure 13. It should be noted that the fluid viscosity here in the X-axis is the fluid viscosity when the breakdown happens. We can observe the experimental results and numerical simulation results coincide considerably, which further validates our numerical model effectiveness. Moreover, with the increase of fluid viscosity, the fracture tortuosity decreases dramatically. We may conclude the effect of fluid viscosity on the fracture tortuosity without considering the role of the in situ stress state. Furthermore, we summarize the fracture tortuosity against the fluid viscosity in Figure 13. It should be noted that the fluid viscosity here in the X-axis is the fluid viscosity when the breakdown happens. We can observe the experimental results and numerical simulation results coincide considerably, which further validates our numerical model effectiveness. Moreover, with the increase of fluid viscosity, the fracture tortuosity decreases dramatically. We may conclude the effect of fluid viscosity on the fracture tortuosity without considering the role of the in situ stress state. Furthermore, we summarize the fracture tortuosity against the fluid viscosity in Figure 13. It should be noted that the fluid viscosity here in the X-axis is the fluid viscosity when the breakdown happens. We can observe the experimental results and numerical simulation results coincide considerably, which further validates our numerical model effectiveness. Moreover, with the increase of fluid viscosity, the fracture tortuosity decreases dramatically. We may conclude the effect of fluid viscosity on the fracture tortuosity without considering the role of the in situ stress state. However, combining these results with previous research results, we believe the in situ stress is still the dominant factor in determining the fluid-driven fracture propagation from the macroscopic aspect, as indicated in Figures 5 and 9 [29]. The fractures easily propagate parallel to the direction of maximum horizontal stress. With the increase of horizontal stress difference, a smooth fracture without secondary fractures or branches is formed. If the horizontal stress difference is small, fracture networks form easily, which may benefit shale gas production.
Hence, the competing roles between fluid viscosity and in situ stress determine the fracture propagation and morphology. As shown in Figure 14, low viscosity fluid is more likely to penetrate smaller pore throats in which large viscosity fluid will be excluded. By using the fluid entry pressure in shale rocks, the interfacial tension, and the contact angle between fluid and rock, the estimated However, combining these results with previous research results, we believe the in situ stress is still the dominant factor in determining the fluid-driven fracture propagation from the macroscopic aspect, as indicated in Figures 5 and 9 [29]. The fractures easily propagate parallel to the direction of maximum horizontal stress. With the increase of horizontal stress difference, a smooth fracture without secondary fractures or branches is formed. If the horizontal stress difference is small, fracture networks form easily, which may benefit shale gas production.
Hence, the competing roles between fluid viscosity and in situ stress determine the fracture propagation and morphology. As shown in Figure 14, low viscosity fluid is more likely to penetrate smaller pore throats in which large viscosity fluid will be excluded. By using the fluid entry pressure in shale rocks, the interfacial tension, and the contact angle between fluid and rock, the estimated result indicates that the Sc-CO 2 may enter pore throats as small as 5.7 nm while H 2 O cannot enter pore throats larger than 12.9 nm [25]. Hence, fracture propagation in shale formations is not only controlled by far-field stress state, but also by fluid viscosity. At the macroscopic scale, the in situ stress makes the fracture propagate parallel to the maximum horizontal stress direction. At the microscopic scale, low-viscosity fluid penetrates shale smaller pore throats and creates fractures with larger tortuosity; while high-viscosity fluid creates a relatively smooth fracture surface. result indicates that the Sc-CO2 may enter pore throats as small as 5.7 nm while H2O cannot enter pore throats larger than 12.9 nm [25]. Hence, fracture propagation in shale formations is not only controlled by far-field stress state, but also by fluid viscosity. At the macroscopic scale, the in situ stress makes the fracture propagate parallel to the maximum horizontal stress direction. At the microscopic scale, low-viscosity fluid penetrates shale smaller pore throats and creates fractures with larger tortuosity; while high-viscosity fluid creates a relatively smooth fracture surface. Low-viscosity fluid more easily penetrates smaller pore throats and creates secondary fracture or fractures with larger tortuosity.

Conclusions
We described a combined experimental-numerical investigation into the process of the non-aqueous fracturing process, with emphasis on the effect of fluid characteristics and in situ stress on the fluid-driven fracture propagation and morphology. A series of laboratory fracturing

Conclusions
We described a combined experimental-numerical investigation into the process of the non-aqueous fracturing process, with emphasis on the effect of fluid characteristics and in situ stress on the fluid-driven fracture propagation and morphology. A series of laboratory fracturing experiments have been conducted with water, liquid carbon dioxide, and supercritical carbon dioxide, and the induced fracture surfaces have been visualized by the CT scanning method. The numerical model coupling rock property heterogeneity, micro-scale damage and fluid flow is built and performed to gain further insights into the fracturing process. From the experimental and numerical results, the following conclusions are drawn: (1) The in situ stress state is the dominant parameter in controlling fracture propagation for fluid-driven fractures. Fracture networks are easily formed when horizontal stress difference is small and a smooth fracture parallel to the maximum horizontal stress direction will be formed when the horizontal stress difference is large. Moreover, the effect of natural fractures and bedding plans may deflect the fracture propagation direction. (2) The fracturing fluid with lower viscosity may induce secondary fractures and creates more branches during main fracture propagation. Due to the low-viscosity property, non-aqueous fluid is more likely to penetrate into smaller pore throats in which large viscosity fluid will be excluded in the shale matrix. More branches and secondary fractures may benefit the shale gas production through increasing reservoir permeability. (3) The competing effect of fluid viscosity and stress state decides the fracture propagation and morphology. The non-aqueous fluids like liquid carbon dioxide, supercritical carbon dioxide, and nitrogen are promising for working as fracturing fluid by creating more branches and secondary fractures and finally forming the fracture networks. However, economic, safety and filtration quality should be further analyzed before wide utilization of non-aqueous fluid for fracturing in shale reservoirs.