TY - JOUR AU - Zhang,, Ruxin AB - Abstract Bedding structure and strong tectonic compression of deep shale reservoirs have a great influence on the vertical propagation of hydraulic fractures and proppant transportation, posing a severe challenge for the hydraulic fracturing design. A combined experimental and numerical study was carried out to investigate the vertical propagation of hydraulic fractures in laminated shales. As a supplement to the physical experiment method, a flow-stress-displacement-damage coupling method based on continuum mechanics was developed, which can effectively simulate complex hydraulic fracture geometries observed in mining experiments. The validity of the numerical method was verified by experiment results. The experiment and simulation results show that: (1) strong horizontal tectonic compression is not favorable for hydraulic fractures to pass through bedding planes; (2) the decrease of bedding interface strength and increase of bedding dip angle will cause hydraulic fracture diversion, offset and containment; (3) increasing the viscosity of the fracturing fluid can reduce the filtration loss of fracturing fluid into beddings, and improve the ability of hydraulic fractures to penetrate bedding planes in various unfavorable conditions; (4) under the same geo-stress condition, increasing the injection rate of fracturing fluid can promote the penetration of beddings with low dip angles, but the high injection rate may lose its role when low viscosity fluid is used under conditions of high tectonic compression and high bedding dip angles (over 30°). In such conditions, fracturing fluid with higher viscosity and low damage to reservoirs is recommended. laminated shale, vertical propagation of hydraulic fractures, physical experiment, FSDD method 1. Introduction Shale gas is natural gas produced from shale formations. Un-fractured shales have permeability on the order of 0.01–0.000 01 md. Because of the ultra-low permeability of the shale matrix, most shale gas today is produced through hydraulic fracturing, often in combine with horizontal drilling (Denney 2009, Laffin and Kariya 2011, Psarras et al2017). Mechanical parameters of reservoir rock, such as rock brittleness, have been adopted as key parameters for testing reservoir fracturing performance (Rickman et al2008, Mahanta et al2016). In southeastern Sichuan Basin, buried depth of gas shale reservoirs ranges from 2000 m to around 4500 m and field data show that most of the formation is in a stress state of the strike-slip fault because of significant tectonic extrusion. The brittleness of the reservoir shale is high based on the brittle evaluation method by Rickman et al (2008) and good fracturing performance can be obtained when the reservoirs lie shallow (with buried depths range from 2000 to 3000 m) by injecting slick water at a high rate. But, field experience showed that the same technique was not applicable for all, especially for deep buried reservoirs. Problems such as ultra-high breakdown pressure and proppant transport difficulty were encountered during fracturing operations, therefore, apart from brittleness, there are other factors that impact fracturing performance. Further investigations are urgently needed for improved hydraulic fracturing performance in these reservoirs. It has been proved that, with the existence of natural fractures in Shales, hydraulic fracture geometries are highly affected leading to high fracture tortuosity. Mineback experiments of formations hydraulically fractured have firstly provided results about how a hydraulic fracture interacted with natural fractures (Warpinski and Teufel 1987, Warpinski 1991, Warpinski et al1993). Four types of encounter modes between hydraulic fractures and natural fractures have been concluded as: penetration, diversion, offset and arrest (Thiercelin et al1987, Beugelsdijk et al2000). Especially, bedding planes sometimes act as weak interfaces, leading to the arrest or diversion of hydraulic fractures propagating in the vertical direction. The fracture height may be overestimated without considering the effect of bedding planes, resulting in suboptimal pumping schedules. In addition, the non-planarity of fracture system caused by fracture offsets will also hinder proppant transport and placement. Therefore, more reasonable prediction of fracture geometry in laminated shale reservoirs plays an important role in hydraulic fracturing design. The mechanisms of fracture growth in fractured reservoirs have been studied for decades using both experimental and modeling techniques. The effect of interface properties between adjacent formations on fracture containment was firstly simulated by small-scale laboratory experiments (Daneshy 1978, Anderson 1981, Hanson et al1981). Experiment results showed that friction coefficient, bonding strength and normal strength of interface were among the important factors that affect fracture containment. It was concluded that frictional slip along the interface reduced the chances of hydraulic fracture crossing the frictional interface. Large-scale hydraulic fracturing physical experiments have been used to study the interaction between hydraulically induced and pre-existing fractures with more external conditions considered. It was found that hydraulic fractures tend to cross pre-existing fractures under high differential stresses and high angles of approach (Blanton 1982, Zhou et al2008). An increase in flow rate or fluid viscosity will enhance propagation across the pre-existing fractures and reduce the hydraulic fracture tortuosity (Beugelsdijk et al2000). Some experiments were conducted on real laminated rock outcrops and similar results have been obtained (Hou et al2014, Tan et al2017). Effect of cemented natural fractures on hydraulic fracture propagation has also been examined through experiment tests, and three interaction forms were found: (1) bypass, (2) arrest then diversion, and (3) a combination of bypass and diversion (Olson et al2012). The results of physical experiments are direct while the factors that can be examined are less. Modeling techniques can make up for the deficiencies by simulating more complex conditions at relatively low cost. Former researches on how natural fractures affect hydraulic fracture propagation can be divided into two categories, one is for the impact of natural fractures with the high dip on hydraulic fractures propagating along the horizontal direction; the other one is the impact of bedding planes on hydraulic fractures propagating in the vertical direction. Both finite element method (FEM) and boundary element method (BEM) have been typically used in rock mechanics applications. BEM simplifies boundary value problem by only discretizing external and internal boundaries and is well suited for the investigations of fractures in homogeneous materials. Based on BEM, the influence of deformation along bedding contacts on fracture intersection with bed contacts has been explored and concluded three contact modes: (1) arrest, (2) penetration and (3) step-over (Cooke and Underwood 2001, Zhang and Jeffrey 2007). Hydraulic fracture reorientation affected by faults with different approach angles has also investigated based on displacement discontinuity method and results showed the fault evidently influences fracture reorientation. Based on continuum mechanics, a coupled fracture and flow modeling approach can be developed in modeling the problem of fracture development (Vishal et al2015). Coupled with damage evolution and fluid flow, FEM method can also model fracture propagation, so-called the flow-stress-damage (FSD) coupling method. Based on this method, not only propagation of a single fracture can be studied but also the influence of fractures on each other be easily considered (Tang et al2002, Li et al2012, Wang et al2016). Cohesive zone model is another one of widely used methods to simulate the fracture initiation and propagation and a finite element-discrete element method has been developed to simulate the complex fracture propagation in shale formations and fracture height growth in laminated formations in recent years (Li et al2016, Yan et al2016, Li et al2017, Li et al2017, Lisjak et al2017). Simulation results showed that without considering the effect of bedding planes, the fracture design may be failed by the inaccurate prediction of fracture height. In this paper, a group of physical experiments on shale outcrops was first conducted to investigate the influence of tectonic compression, fracturing fluid viscosity and injection rate on vertical propagation of hydraulic fractures; then numerical simulations based on flow-stress-displacement-damage (FSDD) method were carried out to make a further analysis of the effect of bedding strength and bedding dip angle. The research findings would give valuable suggestions to field fracturing practice. 2. Experimental study 2.1. Experiment apparatus and specimens Physical simulation experiments have been widely used in the study of hydraulic fracture propagation because they can offer visible results. In this research, experiments were performed using a tri-axial hydraulic fracturing test system, which was composed of a pressure-bearing system, a hydraulic intensifier system, a fluid injection system, a data acquisition system, and other auxiliary devices (figure 1). The specimens used in the experiments were collected from Longmaxi outcrops in southeastern Sichuan Basin. For better experimental control and result contrast, all samples that used in experiments were cut out from the same large shale outcrop. The samples were cut into 300 mm × 300 mm × 300 mm cubes (figure 2). To study the influence of bedding planes on the vertical propagation of hydraulic fractures, bedding planes of each block were maintained parallel to one group of the opposite block faces when cutting out rock samples from a large rock mass and the direction of the overlying stress was set perpendicular to the bedding planes (figure 2). Because the hydraulic fracture initiation process is not the focus of this study, and fracture initiation configuration, transverse or longitudinal, that decided by wellbore direction (Valko and Economides 1995) has little effect on the study of the influence of beddings on hydraulic fracture propagation, so the direction of the simulated wellbore (figure 2) in every experiment is not necessarily set along the minimum horizontal stress according to the field practice. In addition, experiment experience shows that setting the direction of the simulated wellbore along the maximum horizontal stress can significantly reduce the fracture initiation pressure and ensure the success rate and safety of experiments, therefore in this study, the direction of simulated wellbore was set along the maximum horizontal principal stress in each experiment to facilitate initiation of hydraulic fractures. The simulated wellbore was actually a centralized hole with a diameter of 20 mm and a depth of 170 mm drilled in each block (figure 2). A steel pipe with a length of 130 mm cemented in the simulated wellbore was used to pump in fracturing fluid, and an open hole with 40 mm length was left uncased in each block for the initiation of hydraulic fractures. Figure 1. View largeDownload slide Schematic of a tri-axial hydraulic fracturing test system. Figure 1. View largeDownload slide Schematic of a tri-axial hydraulic fracturing test system. Figure 2. View largeDownload slide Shale sample used in experiments. Figure 2. View largeDownload slide Shale sample used in experiments. 2.2. Experiment scheme To study the effects of tectonic compression, injection rate and viscosity of fracturing fluid on hydraulic fracture propagation across bedding planes, experiment parameters were designed as shown in table 1. Stress states of the strike-slip fault (overlying stress is the intermediate principal geo-stress) were set in experiments, similar to the stress condition of the shale formation in southeastern Sichuan Basin. The effect of tectonic compression on the hydraulic fracture height is one focus of the investigation. Two sets of tri-axial stress parameters with different minimum horizontal stress (6 and 15 MPa) were set in experiments. The fracturing fluid viscosities of 3 mPa s and 120 mPa s can represent the viscosities of slick water and glue liquid widely used in the field. Two different injection rates, 10 and 60 ml min-1, were selected based on similarity criterion (Pater et al1994) and experiment experience (Hou et al2014, Hou et al2016, Tan et al2017) to study the influence of different fracturing pump rates, such as low pump rates (1 ∼ 4 m3 min-1), and high pump rates (14 ∼ 16 m3 min-1), used in the field. Luminous green fluorescent dye was added to the fracturing fluid before every test to help observe hydraulic fracture geometries. The amount of luminous green fluorescent used in the experiment was strictly controlled and had little influence on the viscosity of fracturing fluid. After each experiment, the block was cut open along the main hydraulic fracture to observe the inside fracture geometry. Table 1. Experiment parameters. Experiment number Triaxial stress σσ/σH/σh(MPa) q (ml min-1) Fluid viscosity μ (mPa s) Fracture geometry after experiments 1 18/25/6 10 3 ±-shape 2 18/25/6 10 120 Single fracture 3 18/25/15 10 3 H-shape 4 18/25/15 60 3 ±-shape Experiment number Triaxial stress σσ/σH/σh(MPa) q (ml min-1) Fluid viscosity μ (mPa s) Fracture geometry after experiments 1 18/25/6 10 3 ±-shape 2 18/25/6 10 120 Single fracture 3 18/25/15 10 3 H-shape 4 18/25/15 60 3 ±-shape View Large Table 1. Experiment parameters. Experiment number Triaxial stress σσ/σH/σh(MPa) q (ml min-1) Fluid viscosity μ (mPa s) Fracture geometry after experiments 1 18/25/6 10 3 ±-shape 2 18/25/6 10 120 Single fracture 3 18/25/15 10 3 H-shape 4 18/25/15 60 3 ±-shape Experiment number Triaxial stress σσ/σH/σh(MPa) q (ml min-1) Fluid viscosity μ (mPa s) Fracture geometry after experiments 1 18/25/6 10 3 ±-shape 2 18/25/6 10 120 Single fracture 3 18/25/15 10 3 H-shape 4 18/25/15 60 3 ±-shape View Large 2.3. Experiment results and analysis The experiment results are shown in figure 3 and on the right side of each picture shows an illustration of fracture morphology in each block. The bedding planes (grey surfaces) in all illustration pictures are horizontal, making the propagation of hydraulic fractures (blue surfaces) in different situations easily compared. Some conclusions and discussions about the experiment results are summarized as follows: From the experiment results, three different morphologies of fracture system can be observed. They are the ±-shape fracture (figures 3(a) and (d)), the H-shape fracture (figure 3(c)) and the single main fracture (figure 3(b)), which represents the effect of beddings on the fracture propagation. Due to the filtration of fracturing fluid into beddings, the hydraulic fracture height was influenced by bedding planes to different degrees (figures 3(a), (c), and (d)) because of the loss of hydraulic power. Bedding planes can be penetrated more easily under the situation of lower tectonic compression. Compared with Experiment #1, the Experiment #3 had the same overlying stress and maximum horizontal stress but increased minimum horizontal stress (15 MPa), which resulted in high hydraulic pressure (or hydraulic power) needed for hydraulic fracture propagation to overcome the increased horizontal compression stress. Therefore, with limited total hydraulic power, hydraulic fractures tend to be arrested by beddings more easily in Experiment #3 than in Experiment #1, due to the loss of hydraulic power caused by the filtration of fracturing fluid into beddings. In Experiment #4, the total hydraulic power was increased by increasing the injection from 10 ml min-1 in Experiment #3 to 60 ml min-1 (table 1), and a fracture morphology similar with Experiment #1 was obtained (figures 3(a) and (d)). It can be found that under the situation of high tectonic compression, to get the equivalent fracturing performance in low tectonic compression, larger hydraulic power (or higher injection rate) is needed, reflecting the higher degree of difficulty in fracturing the formations under high tectonic compression. The increase of fracturing fluid viscosity can significantly reduce the loss of fracturing fluid into bedding. Limited loss of fracturing fluid into the bedding will reduce the hydraulically pressurized area of bedding plane and maintain enough interface shear strength, which is conducive for hydraulic fracture growth across the interface (Anderson 1981, Hanson etal1981). Therefore, In Experiment #2, a simple main hydraulic fracture was obtained at a low injection rate (10 ml min-1) (figure 3(b)). Figure 3. View largeDownload slide Fracture geometries after experiments. Figure 3. View largeDownload slide Fracture geometries after experiments. The physical experiments on real shale outcrops have several advantages over those conducted on the artificial samples due to the properties of real rock samples are much closer to reservoir rocks. However, the variables of the real rocks that can be controlled are limited. For example, the mechanical properties and dip angles of beddings are similar or the same among samples above, resulting in the incomplete display of the encounter modes between the hydraulic fractures and natural fractures that discovered in former investigations (Warpinski and Teufel 1987, Thiercelin et al1987, Beugelsdijk et al2000), such as the fracture offset and diversion are not presented in above tests. What’s more, lots of physical experiments will consume much effort, time and money. Therefore, proper numerical simulations are necessary to make up for the shortage of physical experiments. 3. Numerical study 3.1. Modeling method The continuum damage theory was employed in this study due to its flexibility in modeling complex fractures (Weinberger et al2000, Busetti et al2012, Busetti et al2012). A FSDD coupling model was developed based on the conventional FSD method (Tang et al2002, Li et al2012). Table 2 compares the typical conventional FSD method and the FSDD method used in this study. The main modification is the damage behavior implemented in FSDD model was characterized by a stress-displacement response rather than a stress–strain response, which could moderate unreasonable mesh sensitivity of the modeling results (Hillerborg et al1976, ABAQUS 2014), and also made it convenient to calculate the permeability of fractured elements (elements completely damaged under tensile stress and with ‘hypothetical fractures’ formed inside) based on the fracture aperture. The main feature of the FSDD model is that it takes into consideration the influence of fracturing fluid filtration on the bedding properties, and can effectively simulate the complex hydraulic fracture morphology under the influence of the bedding planes. The following is a brief description of the new numerical method. Table 2. Comparison between the conventional FSD method and the FSDD method. Comparing items Conventional FSD method FSDD method Description of damage Based on stress–strain relationship Based on stress–displacement relationship Description of permeability Function of stress and damage Function of stress, displacement and damage Comparing items Conventional FSD method FSDD method Description of damage Based on stress–strain relationship Based on stress–displacement relationship Description of permeability Function of stress and damage Function of stress, displacement and damage View Large Table 2. Comparison between the conventional FSD method and the FSDD method. Comparing items Conventional FSD method FSDD method Description of damage Based on stress–strain relationship Based on stress–displacement relationship Description of permeability Function of stress and damage Function of stress, displacement and damage Comparing items Conventional FSD method FSDD method Description of damage Based on stress–strain relationship Based on stress–displacement relationship Description of permeability Function of stress and damage Function of stress, displacement and damage View Large In the simulation, the shale was regarded as a porous medium with ultra-low permeability and effect of flow-stress coupling on the stress field of the model was considered. The control equation and continuity equation in the finite element calculation are as following (Detournay and Cheng 1993): ∇(σ′-αpI)+b=0∇kγl(∇p-bl)+mTε̇=0 1 where σ′ is the effective stress tensor, compression is taken as positive; α is the Biot coefficient; p is the pore pressure; ε is the strain tensor; k is the permeability; γl is the specific weight of liquid; b is the body force vector of rock; bl is the body force of liquid phase; I is the Kronecker delta; m=[1,1,1,0,0,0]. To predict damage initiation of rock material at both compression and tension stress, the Drucker–Prager criterion of hyperbolic form was employed (ABAQUS 2014): F=(d0-pt0tanθ)2+q′2-p′tanθ-d=0, 2 where F is the yield function; q′ is the deviatoric stress; p′ is the effective mean stress; θ and d are the friction angle and cohesion of linear Drucker–Prager criterion in p′-q′ space; d0 is the initial value of cohesion; pt0 is the initial hydrostatic tension strength. When the stress condition of one element satisfies the failure criterion showed in equation (2), the element damage occurs and the strength and modulus of the element degrade, which can be described as (Tang et al2002): E=(1-D)E0σf=(1-D)σf0, 3 where D represents the damage variable; E and E0 are the elastic moduli of damaged and undamaged elements; σf and σf0 are the element strengths of damaged and undamaged elements. When a tensile failure occurs in one element, i.e. ∣σ3∣>σt0, taking uniaxial tension condition as an example (figure 4) the damage variable D can be expressed as: D=0,u3≤ut0u3-ut0utr-ut0⋅1-σtrσt0,ut0utu, 4 where u3 is the change of an element length under tensile stress, which can be obtained by the product of character length and strain of the element; σt0 and σtr are the element peak strength and residual strength, respectively. Other parameters in the equation are defined in figure 4. When the element length exceeds utu, the element strength can be considered damaged completely and macroscopic fractures begin to produce. Here the fracture width was assumed approximately equal to the plastic deformation displacement of fractured elements (Zhou and Hou 2013). Figure 4. View largeDownload slide Stress–displacement relationship of an element subject to the uniaxial tensile stress (at lower left) and uniaxial compressive stress (at upper right). Figure 4. View largeDownload slide Stress–displacement relationship of an element subject to the uniaxial tensile stress (at lower left) and uniaxial compressive stress (at upper right). When a compressive failure occurs in one element, taking uniaxial compression condition (figure 4) as an example the damage variable D can be expressed as: D=0,u1≤uc0u1-uc0ucr-uc01-σcrσc0,uc0ucr, 5 where u1 is the change of an element length under compressive stress; σc0 and σcr are the element peak strength and residual strength, respectively. Other parameters in the equation are defined in figure 4. Different with the tensile condition, the element damage variable will not reach 1 because of the existence of residual stress due to the shear friction of failure surfaces under compression condition (figure 4). The rational description of fluid flow in rock is the key to accurately simulating the propagations of hydraulic fractures. In the conventional FSD model, the permeability evolution of damaged element was described by considering effects of stress and damage based on rock compression experiments (Zhang et al2000, Fredd et al2001, Tang et al2002, Wheaton 2016). However, the fractures formed during hydraulic fracturing are mainly caused by the tensile failure. Here the cubic law for fluid flow in a deformable rock fracture was implemented to represent permeability of fractured element under tensile stress (Witherspoon et al1980) as a supplement to the conventional method. Through the analysis above, the permeability of a rock element under different stress, deformation and damage conditions can be expressed as: k=k0e-β(σm-αp),D=0ξ(D)k0e-β(σm-αp),0σt0, 6 where k is the effective permeability of rock matrix under stress; k0 is the permeability of undamaged rock matrix; β represents the coupling parameter that reflects the influence of stress on the coefficient of permeability; α is the coefficient of pore-fluid pressure; σm is the mean stress; ξ(D) represents the permeability multiplier assumed to be a linear function of the damage factor, expressed as ξ(D)=1+D(ξu-1), ξu is the permeability multiplier when the rock element is completely damaged ( D=1 ⁠); l is the character length of the rock element; b is the width of ‘hypothetical fracture’, represented by maximum tensile principle plastic displacement of the element. When a model geometry was built, the whole geometry would be divided into several sections and the narrow sections were selected as bedding planes while the rest as rock matrix (figure 6). The bedding plane sections should be narrow enough to be composed of only one layer of elements after the whole model geometry meshed (figure 6(c)). Because only one constitutive relationship was used in the model, to distinguish the differences between the mechanical behavior of bedding and matrix, a nonzero initial state variable of damage was set to beddings while the initial damage variable of matrix was set 0, then the bedding element would have initial lower elastic modulus, lower strength and higher permeability than matrix element, according to equations (3) and (6). During the process of calculation, the damage of a bedding plane element would show a dynamic change with the element deformation so as the elastic parameter, strength, and permeability. According to the above method, the whole model was built using ABAQUS as a platform. A user subroutine was implemented in the model to relate rock element parameters (e.g. elastic modulus, strength, and permeability) to the damage (ABAQUS 2014). The calculation process in every analysis step is shown in figure 5. Figure 5. View largeDownload slide Flowchart of the FSDD model. Figure 5. View largeDownload slide Flowchart of the FSDD model. 3.2. Model verification The validity of the model needs to be verified by experiment results before application. Through the analysis of fracture geometries after the experiments (figure 3), it was found that there were no large deflections of fractures along the wellbore axial direction. Therefore, the numerical simulation can be simplified to a two-dimensional plane strain problem. Figure 6 shows the model geometry and loading method in the following simulation, where the model size and stress data was the same as that in the experiments. The model geometry was divided into sections of three rock matrices and two bedding planes (figure 6). The distances between the bedding planes and the wellbore were determined according to the real sample conditions as much as possible, but there was also another factor to be considered. During simulation, the injection rate cannot be applied on the borehole immediately or in a short time as in real experiments, which would cause non-convergence of calculation, therefore the injection rate was set gradually increased to that used in experiments. Then a proper distance between the wellbore and the bedding is needed to make sure the inject rate reaches to the experiment value before the hydraulic fracture meets the nearest bedding because hydraulic fracture will propagate when the injection rate is still increasing. After several simulation attempts, a distance of 45 mm was proved to satisfy all experimental conditions. Eventually, to make modeling results more general, the vertical distances of the two bedding planes to the center of the wellbore were set to two unequal values, which were 45 and 75 mm (figure 6(b)). Free quad meshes of an average edge length of 1.5 mm were used to divide the model geometry (figure 6(c)). The 4-node plane strain quadrilateral, bilinear displacement, bilinear pore pressure, reduced integration element type (CPE4RP in ABAQUS) was selected for the simulation purpose. Figure 6. View largeDownload slide Model geometry and loading method. Figure 6. View largeDownload slide Model geometry and loading method. Material parameters of shale used in the simulations can be checked in table 3. Parameters such as elastic parameters, strength parameters, and initial permeability were obtained by rock laboratory experiments. Some parameters were borrowed from formers’ studies on the similar rock samples, such as the influence coefficient of stress on permeability ( β ⁠) (Vishal et al2013, Chen et al2016). The mutation coefficient of permeability ( ξu ⁠) used in this model was obtained by comparing the permeability of fractured shale cores and intact shale cores under the same stress condition. The initial damage of bedding planes in the experiments was acquired by comparing tensile strengths of shale cores containing bedding planes and complete cores, and an average damage value 0.6 is taken as the initial damage of bedding planes. Other parameters such as the viscosity and injection rate of the fracturing fluid were set up regarding the experimental parameters (table 1). The simulation results are described as follows. Table 3. Material parameters of shale used in the simulations. Parameters Value Initial Young’s modulus ( E0 ⁠), MPa 28 000 Poisson ratio ( μ ⁠) 0.2 Initial tensile strength ( σt0 ⁠), MPa 14 Residual tensile strength ( σtr ⁠), MPa 1 Initial compressive strength ( σc0 ⁠), MPa 150 Residual compressive strength ( σcr ⁠), MPa 40 Friction angle ( θ ⁠), ° 58 Initial hydrostatic tension strength ( pt0 ⁠), MPa 5.4 Dilation angle ( φ ⁠), ° 50 Initial liquid saturation ( sl0 ⁠) 0.1 Initial permeability ( k0 ⁠), μm2 0.3 × 10–6 Influence coefficient of stress on permeability ( β ⁠) 0.15 Mutation coefficient of permeability when D=1 ( ξu ⁠) 2 × 103 Coefficient of pore-water pressure ( α ⁠) 0.8 Initial damage of bedding planes ( D ⁠) 0.6 Parameters Value Initial Young’s modulus ( E0 ⁠), MPa 28 000 Poisson ratio ( μ ⁠) 0.2 Initial tensile strength ( σt0 ⁠), MPa 14 Residual tensile strength ( σtr ⁠), MPa 1 Initial compressive strength ( σc0 ⁠), MPa 150 Residual compressive strength ( σcr ⁠), MPa 40 Friction angle ( θ ⁠), ° 58 Initial hydrostatic tension strength ( pt0 ⁠), MPa 5.4 Dilation angle ( φ ⁠), ° 50 Initial liquid saturation ( sl0 ⁠) 0.1 Initial permeability ( k0 ⁠), μm2 0.3 × 10–6 Influence coefficient of stress on permeability ( β ⁠) 0.15 Mutation coefficient of permeability when D=1 ( ξu ⁠) 2 × 103 Coefficient of pore-water pressure ( α ⁠) 0.8 Initial damage of bedding planes ( D ⁠) 0.6 View Large Table 3. Material parameters of shale used in the simulations. Parameters Value Initial Young’s modulus ( E0 ⁠), MPa 28 000 Poisson ratio ( μ ⁠) 0.2 Initial tensile strength ( σt0 ⁠), MPa 14 Residual tensile strength ( σtr ⁠), MPa 1 Initial compressive strength ( σc0 ⁠), MPa 150 Residual compressive strength ( σcr ⁠), MPa 40 Friction angle ( θ ⁠), ° 58 Initial hydrostatic tension strength ( pt0 ⁠), MPa 5.4 Dilation angle ( φ ⁠), ° 50 Initial liquid saturation ( sl0 ⁠) 0.1 Initial permeability ( k0 ⁠), μm2 0.3 × 10–6 Influence coefficient of stress on permeability ( β ⁠) 0.15 Mutation coefficient of permeability when D=1 ( ξu ⁠) 2 × 103 Coefficient of pore-water pressure ( α ⁠) 0.8 Initial damage of bedding planes ( D ⁠) 0.6 Parameters Value Initial Young’s modulus ( E0 ⁠), MPa 28 000 Poisson ratio ( μ ⁠) 0.2 Initial tensile strength ( σt0 ⁠), MPa 14 Residual tensile strength ( σtr ⁠), MPa 1 Initial compressive strength ( σc0 ⁠), MPa 150 Residual compressive strength ( σcr ⁠), MPa 40 Friction angle ( θ ⁠), ° 58 Initial hydrostatic tension strength ( pt0 ⁠), MPa 5.4 Dilation angle ( φ ⁠), ° 50 Initial liquid saturation ( sl0 ⁠) 0.1 Initial permeability ( k0 ⁠), μm2 0.3 × 10–6 Influence coefficient of stress on permeability ( β ⁠) 0.15 Mutation coefficient of permeability when D=1 ( ξu ⁠) 2 × 103 Coefficient of pore-water pressure ( α ⁠) 0.8 Initial damage of bedding planes ( D ⁠) 0.6 View Large Figure 7 shows the simulation results of Experiment #1. It can be seen the hydraulic fracture initiated from the wellbore firstly contacted and penetrated the nearer bedding plane I at 12.9 s, then the upper bedding plane II was also penetrated at 27.31 s. With the hydraulic fracture in the lower side contacted with the model boundary at 65.73 s, the calculation stopped because of calculation non-convergence. It can be found that the propagation of hydraulic fracture in the upper side was obviously influenced by the bedding plane due to the loss of hydraulic fluid into the beddings. Moreover, the distribution of pore pressure was almost the same as the distribution of damage because the fracturing fluid almost did not infiltrate into the shale matrix. Therefore, in the following parts, if there is no special need, pore pressure distributions will be displayed to represent the hydraulic fracture shape. Figure 7. View largeDownload slide Evolutions of pore pressure and damage in the simulation of Experiment #1 (the deformation scale is 50). Figure 7. View largeDownload slide Evolutions of pore pressure and damage in the simulation of Experiment #1 (the deformation scale is 50). The four simulation results are summed up in figure 8(a), which can reasonably simulate the influence of the tectonic compression, the viscosity, and injection rate of the fracturing fluid on the shape of hydraulic fractures showed in experiments (figure 3). The difference is that, in numerical simulations, beddings on both sides of wellbore would often be penetrated by hydraulic fractures (figure 8(a)-#1, #2 and #4), while in experiments there were situations that only beddings on one side were crossed (figure 8(b)-#1 and #4). This might be related to the strength non-uniformity of beddings in one sample. On the whole, the simulation results are in a good agreement with the experimental results and the effectiveness of the FSDD method is verified. Based on the simulation model, a sensitive study can be carried out to make up for the inadequacy of the factors considered in the physical experiment. In the following part, a further analysis of the effect of bedding strength and bedding dip angle is carried out. Figure 8. View largeDownload slide Comparison between the modeling and experiment results, (a) is the simulation results (the deformation scale is 50); (b) is the experiment results (the blue lines represent bedding planes and the red lines represent hydraulic fractures). Figure 8. View largeDownload slide Comparison between the modeling and experiment results, (a) is the simulation results (the deformation scale is 50); (b) is the experiment results (the blue lines represent bedding planes and the red lines represent hydraulic fractures). 3.3. Sensitive study based on the numerical model Usually, in order to cut out complete rock samples, the outcrops with bedding planes of relatively high bonding strength would be selected and the cutting direction was generally maintained perpendicular or parallel to beddings to prevent the specimen from breaking along the beddings due to mechanical vibration. The limitation of the properties of the samples leads to the lack of enough factor analysis in the physical experiments. Numerical simulation method can serve as a good supplement to make up the shortage of physical experiment. In this part, the influence of bedding strength and dip angle on the vertical propagation of hydraulic fractures was investigated by the FSDD model. Figure 9. View largeDownload slide Simulation results of experiments with increasing bedding damage (the deformation scale is 50). Figure 9. View largeDownload slide Simulation results of experiments with increasing bedding damage (the deformation scale is 50). 3.3.1. Effect of bedding strength In the numerical simulation, the bedding strength is evaluated by bedding damage. High degree of damage represents the low strength of beddings. When the damage of one bedding equals one, the bedding is completely damaged and the bedding strength is 0. Because bedding damage increases dynamically with the deformation of bedding elements, therefore the initial damage of beddings should be defined before every simulation. Four kinds of initial bedding damage were taken into consideration, which were 0.4, 0.6, 0.8 and 1.0. Parameters of the geo-stress, injection rate and viscosity of the fracturing fluid in the numerical simulations were consistent with the corresponding parameters in the four physical experiments (table 1). Figure 10. View largeDownload slide Fracture network obtained from numerical simulation (a) and visualization of fracturing treatment results in a jointed rock mass (b). (b) Reproduced with permission from Warpinski and Teufel (1987), copyright 1987. Figure 10. View largeDownload slide Fracture network obtained from numerical simulation (a) and visualization of fracturing treatment results in a jointed rock mass (b). (b) Reproduced with permission from Warpinski and Teufel (1987), copyright 1987. The simulation results of models with different bedding damage have been summarized in figure 9 and several findings can be obtained as follows: The smaller the bedding damage, the more easily the bedding is passed through. When initial bedding damage reduced to 0.4, all beddings are crossed by hydraulic fractures (figures 9(a)-#1, (b)-#2, (c)-#3 and (d)-#4), therefore, it can be inferred that beddings with initial damage of or less than 0.4 have no effect on the vertical propagation of hydraulic fractures under all current experiment conditions. When bedding damage increases to 1.0, a fracture network is obtained under conditions of low tectonic compression (6 MPa), low fluid viscosity (3 mPa s) and low injection rate (10 ml min-1) (figure 9(a)). In such conditions, after contact with bedding planes, hydraulic fractures offset for certain distances and then initiated fractures propagating vertically again, in good accordance with fracture network morphology depicted from mineback experiments (Warpinski and Teufel 1987) (figure 10), reflecting the effect of unbonded beddings on propagation of hydraulic fractures. Proppant transport and placement would be hindered because of the non-planarity of the fracture system (Warpinski and Teufel 1987). Under conditions of high tectonic compression (15 MPa), low fluid viscosity (3 mPa s) and low injection rate (10 ml min-1), hydraulic fractures will be arrested by bedding planes when bedding damage is greater than 0.6 (figure 9(c)). With varying initial bedding damage, high fluid injection rate (60 ml min-1) and high fluid viscosity (120 mPa s) are always beneficial to the penetration of beddings (figures 9(b) and (d)). 3.3.2. Effect of dip angle of beddings In high-steep structures, the increase of dip angle ( ϑ ⁠) of bedding plane will cause the decrease of normal stress on the bedding plane (figure 11). The results of previous studies have shown that the decrease of the normal stress is not conducive for a hydraulic fracture to cross the bedding plane (Anderson 1981). It is necessary to study the influence of the dip angle of bedding planes on the vertical propagation of hydraulic fractures and the effect of four different dip angles, 0°, 15°, 30° and 45° was considered in this part. The initial bedding damage in different cases was maintained 0.6. The vertical distance between the two beddings to the center of the wellbore in every model remained the same. Figure 11. View largeDownload slide Effect of overlying stress on high-steep structural bedding planes. Figure 11. View largeDownload slide Effect of overlying stress on high-steep structural bedding planes. Effect of bedding dip angle on hydraulic fracture propagation is summarized in figure 12 and some conclusions can be drawn as follows: Under conditions of low tectonic compression (6 MPa), low fluid injection rate (10 ml min-1) and low fracturing fluid viscosity (3 mPa s), increasing dip angle (greater than 30°) would lead to the offset of hydraulic fractures (figure 12(a)), which will hinder the transportation of proppant and result in invalid hydraulic fractures above the beddings during the production stage (Warpinski and Teufel 1987). High fracturing fluid viscosity (120 mPa s) makes the propagation of hydraulic fractures insensitive to the influence of variation of the bedding dip angle (figure 12(b)). Under condition of high tectonic compression (15 MPa), most hydraulic fractures driven by low viscosity (3 mPa s) fluid at a low injection rate (10 ml min-1) tend to be arrested at bedding s with different bedding dip angles, except at the dip angle of 15°, where the hydraulic fracture penetrates the lower bedding plane (figure 12(c)). Combination of low fluid viscosity (3 mPa s) and high injection rate (60 ml min-1) loses its role in promoting penetration of bedding planes when the dip angle increased over 30° under conditions of high tectonic compression (15 MPa) (figure 12(d)). Therefore, the bedding dip angle of 30° might be considered a key parameter for hydraulic fracturing in the high-steep structure with strong tectonic compression. In such a situation, only increasing the injection rate of fracturing fluid cannot guarantee sufficient height of hydraulic fractures and appropriately increasing the viscosity of fracturing fluid should be considered (Ribeiro et al2016). Figure 12. View largeDownload slide Simulation results of experiments with increasing dip angles (the deformation scale is 50). Figure 12. View largeDownload slide Simulation results of experiments with increasing dip angles (the deformation scale is 50). 4. Summary and conclusions Geologic discontinuities, such as joints and bedding planes, and stress condition can significantly influence hydraulic fracture treatments. In this study, a combination of experimental and numerical study was carried out to investigate the influence of tectonic compression, bedding planes, hydraulic fracturing fluid viscosity and injection rate on the vertical propagation of hydraulic fractures and several meaningful conclusions can be reached as follows: Physical experiments show that high tectonic compression has an adverse influence on the vertical propagation of hydraulic fractures in laminated shale. Under the situation of high tectonic compression, larger hydraulic power (or higher injection rate) is needed for shale reservoir treatments to get the equivalent fracturing performances of formations under low tectonic compression. The FSDD method was developed to provide a supplement to the physical experiments. The effectiveness of this method on simulating complex fracture morphology was verified by the experiment results and some meaningful results about the influence of bedding strength and dip angle on fracture propagation were obtained through numerical simulations. For the results of numerical simulation, When low injection rate (10 ml min-1) is used to inject low viscosity fracturing fluid (3 mPa s), hydraulic fractures prone to offset, divert or stop propagation at beddings, resulting in the formation of complex fracture morphology, which will affect the transportation and placement of proppant. The increase of bedding dip angle causes more offset and arrest of hydraulic fractures driven by low viscosity fluid (3 mPa s). The bedding dip angle of 30° might be considered a key parameter. As bedding dip angles increase over 30°, the combination of low fluid viscosity (3 mPa s) and high injection rate (60 ml min-1) loses its role in promoting the penetration of bedding planes under the situation of high tectonic compression (15 MPa). Both the experimental and numerical results show that the increase of fracturing fluid viscosity can reduce the loss of fracturing fluid into beddings, making the vertical propagation of hydraulic fractures less sensitive to various situations even at a low injection rate (10 ml min-1). Hence, when current pump power cannot satisfy the successful treatments of reservoirs in high-steep structures or under strong tectonic compressions, fracturing fluid with higher viscosity and low damage to reservoirs is recommended. Acknowledgments This work is financially supported by National Natural Science Foundation of China (No. 51234006, No. 51325402, No. 51490651). References ABAQUS . , 2014 Abaqus Theory Guide, version 6.14 (Dassault Systèmes Simulia Corp.) (https://www.sharcnet.ca/Software/Abaqus/6.14.2/v6.14/books/stm/default.htm?startat=ch03s09ath96.html) Anderson G D . , 1981 Effects of friction on hydraulic fracture growth near unbonded interfaces in rocks , Soc. Pet. Eng. J. , vol. 21 (pg. 21 - 29 ) https://doi.org/10.2118/8347-PA Google Scholar Crossref Search ADS Beugelsdijk L J L , Pater C J D , Sato K . , 2000 Experimental hydraulic fracture propagation in a multi-fractured medium SPE Asia Pacific Conf. on Integrated Modelling for Asset Management pg. SPE-59419-MS https://doi.org/10.2118/59419-MS Blanton T . , 1982 An experimental study of interaction between hydraulically induced and pre-existing fractures SPE Unconventional Gas Recovery Symp. pg. SPE-10847-MS https://doi.org/10.2118/10847-MS Busetti S , Mish K , Hennings P , Reches Z . , 2012 Damage and plastic deformation of reservoir rocks: II. Propagation of a hydraulic fracture , AAPG Bull. , vol. 96 (pg. 1711 - 1732 ) https://doi.org/10.1306/02011211011 Google Scholar Crossref Search ADS Busetti S , Mish K , Reches Z . , 2012 Damage and plastic deformation of reservoir rocks: I. Damage fracturing , AAPG BULL , vol. 96 (pg. 1687 - 1709 ) https://doi.org/10.1306/02011211010 Google Scholar Crossref Search ADS Chen D , Pan Z , Ye Z , Hou B , Wang D . , 2016 A unified permeability and effective stress relationship for porous and fractured reservoir rocks , J. Nat. Gas Sci. Eng. , vol. 29 (pg. 401 - 412 ) https://doi.org/10.1016/j.jngse.2016.01.034 Google Scholar Crossref Search ADS Cooke M L , Underwood C A . , 2001 Fracture termination and step-over at bedding interfaces due to frictional slip and interface opening , J. Struct. Geol. , vol. 23 (pg. 223 - 238 ) https://doi.org/10.1016/S0191-8141(00)00092-4 Google Scholar Crossref Search ADS Daneshy A A . , 1978 Hydraulic fracture propagation in layered formations , Soc. Pet. Eng. J. , vol. 18 (pg. 33 - 41 ) https://doi.org/10.2118/6088-PA Google Scholar Crossref Search ADS Denney D . , 2009 Evaluating implications of hydraulic fracturing in shale-gas reservoirs , J. Pet. Technol. , vol. 61 (pg. 53 - 54 ) https://doi.org/10.2118/0809-0053-JPT Detournay E , Cheng H . , 1993 Fundamentals of poroelasticity , Comprehensive Rock Engineering: Principles, Practice & Projects Oxford Pergamon Press (pg. 113 - 171 ) Fredd C N , McConnell S B , Boney C L , England K W . , 2001 Experimental study of fracture conductivity for water-fracturing and conventional fracturing applications , Soc. Pet. Eng. J. , vol. 6 (pg. 288 - 298 ) https://doi.org/10.2118/74138-PA Hanson M E , Shaffer R J , Anderson G D . , 1981 Effects of various parameters on hydraulic fracturing geometry , Soc. Pet. Eng. J. , vol. 21 (pg. 435 - 443 ) https://doi.org/10.2118/8942-PA Google Scholar Crossref Search ADS Hillerborg A , Modéer M , Petersson P E . , 1976 Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements , Cement Concr. Res. , vol. 6 (pg. 773 - 781 ) https://doi.org/10.1016/0008-8846(76)90007-7 Google Scholar Crossref Search ADS Hou B , Chen M , Cheng W , Diao C . , 2016 Investigation of hydraulic fracture networks in shale gas reservoirs with random fractures , Arab. J. Sci. Eng. , vol. 41 (pg. 2681 - 2691 ) https://doi.org/10.1007/s13369-015-1829-0 Google Scholar Crossref Search ADS Hou B , Chen M , Zhimeng L I , Wang Y , Diao C . , 2014 Propagation area evaluation of hydraulic fracture networks in shale gas reservoirs , Pet. Explor. Dev. , vol. 41 (pg. 833 - 838 ) https://doi.org/10.1016/S1876-3804(14)60101-4 Google Scholar Crossref Search ADS Laffin M , Kariya M . , 2011 Shale gas and hydraulic fracking 20th World Petroleum Congress Doha, Qatar pg. WPC-20-3262 Li H , Zou Y , Valko P P , Ehlig-Economides C . , 2016 Hydraulic fracture height predictions in laminated shale formations using finite element discrete element method SPE Hydraulic Fracturing Technology Conf. pg. SPE-179129-MS https://doi.org/10.2118/179129-MS Li L C , Tang C A , Li G , Wang S Y , Liang Z Z , Zhang Y B . , 2012 Numerical simulation of 3D hydraulic fracturing based on an improved flow-stress-damage model and a parallel fem technique , Rock Mech. Rock Eng. , vol. 45 (pg. 801 - 818 ) https://doi.org/10.1007/s00603-012-0252-z Li Y , Deng J G , Liu W , Feng Y . , 2017 Modeling hydraulic fracture propagation using cohesive zone model equipped with frictional contact capability , Comput. Geotech. , vol. 91 (pg. 58 - 70 ) https://doi.org/10.1016/j.compgeo.2017.07.001 Google Scholar Crossref Search ADS Li Y , Deng J , Liu W , Yan W , Feng Y , Cao W , Wang P , Hou Y . , 2017 Numerical simulation of limited-entry multi-cluster fracturing in horizontal well , J. Pet. Sci. Eng. , vol. 152 (pg. 443 - 455 ) https://doi.org/10.1016/j.petrol.2017.03.023 Google Scholar Crossref Search ADS Lisjak A , Kaifosh P , He L , Tatone B S A , Mahabadi O K , Grasselli G . , 2017 A 2D, fully-coupled, hydro-mechanical, FDEM formulation for modelling fracturing processes in discontinuous, porous rock masses , Comput. Geotech. , vol. 81 (pg. 1 - 18 ) https://doi.org/10.1016/j.compgeo.2016.07.009 Google Scholar Crossref Search ADS Mahanta B , Tripathy A , Vishal V , Singh T N , Ranjith P G . , 2016 Effects of strain rate on fracture toughness and energy release rate of gas shales , Eng. Geol. , vol. 218 (pg. 39 - 49 ) https://doi.org/10.1016/j.enggeo.2016.12.008 Google Scholar Crossref Search ADS Olson J E , Bahorich B , Holder J . , 2012 Examining hydraulic fracture: natural fracture interaction in hydrostone block experiments SPE Hydraulic Fracturing Technology Conf. pg. SPE-152618-MS https://doi.org/10.2118/152618-MS Pater C , Weijers L , Cleary M P , Quinn T S , Barr D T . , 1994 Experimental verification of dimensional analysis for hydraulic fracturing , SPE Prod. Facil. , vol. 9 (pg. 230 - 238 ) https://doi.org/10.2118/24994-PA Google Scholar Crossref Search ADS Psarras P , Holmes R , Vishal V , Wilcox J . , 2017 Methane and CO2 adsorption capacities of kerogen in the eagle ford shale from molecular simulation , Acc. Chem. Res. , vol. 50 (pg. 1818 - 1828 ) https://doi.org/10.1021/acs.accounts.7b00003 Google Scholar Crossref Search ADS PubMed Ribeiro L H , Li H , Bryant J E . , 2016 Use of a CO2-hybrid fracturing design to enhance production from unpropped-fracture networks , SPE Prod. Oper. , vol. 32 (pg. 28 - 40 ) https://doi.org/10.2118/173380-PA Rickman R , Mullen M J , Petre J E , Grieser W V , Kundert D . , 2008 A practical use of shale petrophysics for stimulation design optimization: all shale plays are not clones of the barnett shale. Society of petroleum engineers SPE Annual Technical Conf. and Exhibition pg. SPE-115258-MS https://doi.org/10.2118/115258-MS Tan P , Jin Y , Han K , Hou B , Chen M , Guo X , Gao J . , 2017 Analysis of hydraulic fracture initiation and vertical propagation behavior in laminated shale formation , Fuel , vol. 206 (pg. 482 - 493 ) https://doi.org/10.1016/j.fuel.2017.05.033 Google Scholar Crossref Search ADS Tang C A , Tham L G , Lee P , Yang T H , Li L C . , 2002 Coupled analysis of flow, stress and damage (FSD) in rock failure , Int. J. Rock Mech. Min. Sci. , vol. 39 (pg. 477 - 489 ) https://doi.org/10.1016/S1365-1609(02)00023-0 Google Scholar Crossref Search ADS Thiercelin M J , Roegiers J C , Boone T J , Ingraffea A R . , 1987 An investigation of the material parameters that govern the behavior of fractures approaching rock interfaces 6th ISRM Congress Montreal, Canada 30 August–3 September pg. ISRM-6CONGRESS-1987-049 Valko P , Economides M J . , 1995 , Hydraulic Fracture Mechanics Chichester Wiley Vishal V , Jain N , Singh T N . , 2015 Three dimensional modelling of propagation of hydraulic fractures in shale at different injection pressures , Sustain. Environ. Res. , vol. 25 (pg. 217 - 225 ) Vishal V , Ranjith P G , Singh T N . , 2013 CO2 permeability of Indian bituminous coals: Implications for carbon sequestration , Int. J. Coal Geol. , vol. 105 (pg. 36 - 47 ) https://doi.org/10.1016/j.coal.2012.11.003 Google Scholar Crossref Search ADS Wang Y , Li X , Zhang Y X , Wu Y S , Zheng B . , 2016 Gas shale hydraulic fracturing: a numerical investigation of the fracturing network evolution in the Silurian Longmaxi formation in the Southeast of Sichuan Basin, China, using a coupled FSD approach , Environ. Earth Sci. , vol. 75 pg. 1093 https://doi.org/10.1007/s12665-016-5696-0 Google Scholar Crossref Search ADS Warpinski N R . , 1991 Hydraulic fracturing in tight, fissured media , J. Pet. Technol. , vol. 43 pg. 146 https://doi.org/10.2118/20154-PA Google Scholar Crossref Search ADS Warpinski N R , Lorenz J C , Branagan P T , Myal F R , Gall B L . , 1993 Examination of a cored hydraulic fracture in a deep gas well (includes associated papers 26302 and 26946) , SPE Prod. Facil. , vol. 8 (pg. 150 - 158 ) https://doi.org/10.2118/22876-PA Google Scholar Crossref Search ADS Warpinski N R , Teufel L W . , 1987 Influence of geologic discontinuities on hydraulic fracture propagation (includes associated papers 17011 and 17074) , J. Pet. Technol. , vol. 39 (pg. 209 - 220 ) https://doi.org/10.2118/13224-PA Google Scholar Crossref Search ADS Weinberger R , Lyakhovsky V , Baer G , Agnon A . , 2000 Damage zones around en echelon dike segments in porous sandstone , J. Geophys. Res.: Solid Earth , vol. 105 (pg. 3115 - 3133 ) https://doi.org/10.1029/1999JB900361 Google Scholar Crossref Search ADS Wheaton R J . , 2016 Dependence of shale permeability on pressure , SPE Reservoir Eval. Eng. , vol. 20 (pg. 228 - 232 ) https://doi.org/10.2118/183629-PA Google Scholar Crossref Search ADS Witherspoon P A , Wang J S , Iwai K , Gale J E . , 1980 Validity of cubic law for fluid flow in a deformable rock fracture , Water Resour. Res. , vol. 16 (pg. 1016 - 1024 ) https://doi.org/10.1029/WR016i006p01016 Google Scholar Crossref Search ADS Yan C , Zheng H , Sun G , Ge X . , 2016 Combined finite-discrete element method for simulation of hydraulic fracturing , Rock Mech. Rock Eng. , vol. 49 (pg. 1389 - 1410 ) https://doi.org/10.1007/s00603-015-0816-9 Google Scholar Crossref Search ADS Zhang J , Bai M , Roegiers J , Wang J , Liu T . , 2000 Experimental determination of stress-permeability relationship 4th North American Rock Mechanics Symp. American Rock Mechanics Association pg. ARMA-2000-0817 Zhang X , Jeffrey R G . , 2007 Hydraulic fracture propagation across frictional interfaces 1st Canada –U.S. Rock Mechanics Symp. American Rock Mechanics Association pg. ARMA-07-204 Zhou J , Chen M , Jin Y , Zhang G Q . , 2008 Analysis of fracture propagation behavior and fracture geometry using a tri-axial fracturing system in naturally fractured reservoirs , Int. J. Rock Mech. Min. Sci. , vol. 45 (pg. 1143 - 1152 ) https://doi.org/10.1016/j.ijrmms.2008.01.001 Google Scholar Crossref Search ADS Zhou L , Hou M Z . , 2013 A new numerical 3D-model for simulation of hydraulic fracturing in consideration of hydro-mechanical coupling effects , Int. J. Rock Mech. Min. Sci. , vol. 60 (pg. 370 - 380 ) https://doi.org/10.1016/j.ijrmms.2013.01.006 Google Scholar Crossref Search ADS © 2018 Sinopec Geophysical Research Institute TI - Experimental and numerical investigations on the vertical propagation of hydraulic fractures in laminated shales JF - Journal of Geophysics and Engineering DO - 10.1088/1742-2140/aac12f DA - 2018-08-01 UR - https://www.deepdyve.com/lp/oxford-university-press/experimental-and-numerical-investigations-on-the-vertical-propagation-gTTYjykWf5 SP - 1729 VL - 15 IS - 4 DP - DeepDyve ER -