Effect of Joint Geometrical Parameters on Hydraulic Fracture Network Propagation in Naturally Jointed Shale Reservoirs

Chong, Zhaohui; Yao, Qiangling; Li, Xuehua

doi:https://doi.org/10.1155/2018/1852604

Geofluids

On this page

Abstract Introduction Results and Analysis Conclusions Data Availability Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2018 | Article ID 1852604 | https://doi.org/10.1155/2018/1852604

Effect of Joint Geometrical Parameters on Hydraulic Fracture Network Propagation in Naturally Jointed Shale Reservoirs

Zhaohui Chong,^1,2Qiangling Yao ,¹and Xuehua Li¹

Academic Editor: John A. Mavrogenes

Received18 Jul 2018

Accepted28 Aug 2018

Published25 Oct 2018

Abstract

The presence of a significant amount of discontinuous joints results in the inhomogeneous nature of the shale reservoirs. The geometrical parameters of these joints exert effects on the propagation of a hydraulic fracture network in the hydraulic fracturing process. Therefore, mechanisms of fluid injection-induced fracture initiation and propagation in jointed reservoirs should be well understood to unleash the full potential of hydraulic fracturing. In this paper, a coupled hydromechanical model based on the discrete element method is developed to explore the effect of the geometrical parameters of the joints on the breakdown pressure, the number and proportion of hydraulic fractures, and the hydraulic fracture network pattern generated in shale reservoirs. The microparameters of the matrix and joint used in the shale reservoir model are calibrated through the physical experiment. The hydraulic parameters used in the model are validated through comparing the breakdown pressure derived from numerical modeling against that calculated from the theoretical equation. Sensitivity analysis is performed on the geometrical parameters of the joints. Results demonstrate that the HFN pattern resulting from hydraulic fracturing can be roughly divided into four types, i.e., crossing mode, tip-to-tip mode, step path mode, and opening mode. As (joint orientation with respect to horizontal principal stress in plane) increases from 0° to 15° or 30°, the hydraulic fracture network pattern changes from tip-to-tip mode to crossing mode, followed by a gradual decrease in the breakdown pressure and the number of cracks. In this case, the hydraulic fracture network pattern is controlled by both (joint step angle) and . When is 45° or 60°, the crossing mode gains dominance, and the breakdown pressure and the number of cracks reach the lowest level. In this case, the HFN pattern is essentially dependent on and (joint spacing). As reaches 75° or 90°, the step path mode is ubiquitous in all shale reservoirs, and the breakdown pressure and the number of the cracks both increase. In this case, has a direct effect on the HFN pattern. In shale reservoirs with the same , either decrease in (joint persistency) and (joint aperture) or increase in leads to the increase in the breakdown pressure and the number of cracks. It is also found that changes in and result in the variation in the proportion of different types of hydraulic fractures. The opening mode of the hydraulic fracture network pattern is observed when increases to 1.2 × 10⁻² m.

1. Introduction

Hydraulic fracturing is widely used in the petroleum industry for boosting the yield of oil/gas wells [1–3]. It was later introduced to the mining industry to weaken the hard roof so as to avoid rock burst hazard [4]. Proved by research results in recent years, hydraulic fracturing is an efficient stimulation approach that can significantly enhance the permeability [5–7]. In the past studies, shale reservoirs are regarded as dual-porosity medium composed of shale matrix and fracture (joint) network [8, 9]. Through hydraulic fracturing, the high-pressure water is able to penetrate the shale matrix and joint, thus not only releasing the adsorbed gas but also creating a fracture network for gas to transport [10, 11]. However, the uncertainty regarding the composition of the shale reservoirs compounded by the complex geological setting makes it extremely difficult to control and predict fracture initiation and propagation in hydraulic fracturing [12]. Therefore, it is of great significance to gain a good understanding of fracture initiation and crack propagation, which can guide hydraulic fracturing operation in naturally jointed shale reservoirs.

A series of models were established by the predecessors to investigate the correlation between the hydraulic fracture (HF) and the naturally occurring discontinuities such as natural fracture (NF), joint, fissure, and cleat. Researchers and scholars also put forward some criterions regarding the relation between the HF and NF. Renshaw and Pollard [13] proposed a relatively simple criterion to determine the propagation direction of the HF after it traverses through the NF that is normal to it. Also, this criterion was proved effective by a following laboratory experiment. Gu and Weng [14] and Gu et al. [15] developed this criterion and employed an updated criterion to investigate the correlation between the nonorthogonal HF and NF. To date, the correlation between HF and NF has been divided mainly into 6 types [13, 16, 17]: (1) HF directly crossing NF, (2) HF intersecting NF, (3) HF crossing NF with an offset, (4) HF arrested by NF, (5) HF shear slipping along NF, and (6) HF branching or turning at the end of NF.

To simulate HF initiation and propagation in jointed rock mass, scholars made huge efforts. Xu et al. [18], Meyer and Bazan [19] put forth a wire-mesh model and discrete fracture network (DFN) model. In these pseudo-3D models, the hydraulic fracture network (HFN) is made up of two sets of cracks. One is the cracks that extend along the direction of the principal stress, and the other is the cracks that are in parallel with the direction of the principal stress by a fixed interval. Weng et al. [20] put forward an unconventional pseudo-3D model that considers several important characteristics of fracturing treatment simulation such as transport of proppant, behavior of non-Newtonian fluid, and permeability parameters of the wellbore. More importantly, the proposed model is able to simulate HFN resulting from the interaction of HF and joints. Olson and Taleghani [21] presented a 2D plane-strain extended finite element method (XFEM), which can model the formation of fracture wing and makes it possible for HF to change the propagation direction. The boundary element method [22] is also introduced to study the interaction between the joints and NF by analyzing the tensile and shear cracks. All the above methods are continuum-based, which can generate arbitrary fracturing paths without remeshing.

Recent years have witnessed the popularity of the discrete element method (DEM) in the study of HF propagation. In this method, the reservoir is generated by a great number of assembled bonded particles (in particle flow code (PFC )) [23] or by a great number of transformable blocks (in universal distinct element code (UDEC) and three-dimensional distinct element code (3DEC)) [24, 25]. The joints between the adjacent particles or blocks morphed into the fluid network which provides conditions for HF propagation [26, 27]. The advantage of DEM over the continuum-based method is that a great number of preexisting discontinuities can be added into the modeled reservoir and the mechanical behavior of these discontinuities can be modeled in the formation of the HFN in jointed rock mass.

In actuality, shale reservoirs contain a great amount of joints, exhibiting typical inhomogeneous behavior. Figure 1 is a site photo showing the Longmaxi shale reservoir in Chongqing, China, which contains a massive amount of joints. These joints occur in a laminar manner. The effect of these joints on the formation of HFN has already been captured by physical experiments [28]. In addition, the activation of these joints induced by HF initiation and propagation has been extensively studied [29, 30]. However, many challenges still lie ahead in the study of the interaction between HF and joints as well as the formation of HFN in naturally jointed shale reservoirs. In the process of performing horizontal well fracturing in shale reservoirs, a massive amount of HFs interact with the joints, damaging the integrity of the shale reservoirs further. In this case, the geometrical parameters of the joints have a pronounced effect on HFN propagation.

In this study, a DEM-based coupled hydromechanical model is developed to explore the effect of the geometrical parameters of the joints on HF initiation and propagation as well as the HFN pattern. The parameters of the matrix and joint used in the shale reservoir model are calibrated through physical experiment, and the hydraulic parameters used in the model are validated through comparing the breakdown pressure derived from numerical modeling against that calculated from the theoretical equation. In addition, the calibrated model is also employed to perform sensitivity analysis on the geometrical parameters of the joints such as joint orientation angle , joint persistency , joint spacing , joint step angle , and joint aperture .

2. Simulation Methodology

2.1. Parallel Bonded Model

The PFC, a simulation code based on the DEM is an effective numerical tool for simulating the behavior of jointed rock mass. In PFC, the rock matrix is modeled as an assembly of circular particles. Though PFC is based on DEM, it can still be employed to model the deformation characteristics and mechanical behavior of a continuous medium. The parallel bonded model (PBM) is used to define and simulate the particle contacts as it can deliver more realistic results in modeling rock fracture.

Figure 2(a) schematizes the micromechanical behavior of the particle contacts. At particle contacts, two sets of springs provide the normal and shear stiffness [31]. As shown in Figures 2(b) and 2(c), the first set consists of two springs that stand for normal and shear stiffness, which are denoted as and , respectively. In the model, the particles are idealized as a rigid body that is unbreakable and undeformable. However, the particles can overlap if the compressive force is great enough. The second set of springs provides the parallel bond normal and shear stiffness denoted as and , which are uniformly distributed over a disc-shaped cross section lying on the contact plane. In addition to the normal and shear stiffness, the parallel bonds also have normal and shear strength. A parallel bond is active over a finite area and can resist the rotation of particles. When the normal or shear stress exceeds the bond strength, the parallel bond is broken and bond stiffness is removed. Once the parallel bond is broken under shear stress, the shear strength falls to its residual value whereas upon breakage of a bond under tensile stress, the tensile strength is set to zero. The residual shear strength depends on the normal force and the friction coefficient of particles at contacts [32]. Microcracks are produced in PBM once the bonds fail due to external loading, and macrocracks are generated by the propagation and coalescence of these microcracks.

(a)

(b)

(c)

2.2. Smooth Joint Model (SJM)

Traditionally, joints are modeled in PFC by either removing the particles that fall within a specific range or assigning a low coefficient of friction to particles at certain contacts. However, this kind of treatment can result in an inherent microscale roughness of the joint surface and leads to an unrealistic increase in shear strength and dilation during shearing along the shear plane. To overcome this problem, a smooth joint model (SJM) was proposed [33]. A typical example of a two-dimensional SJM is shown in Figure 3. In this model, after the joint plane is defined, a smooth joint (SJ) is assigned to the contacts between the particles. The centers of these particles lie on the opposite side of the defined joint plane. At these contacts, the parallel bonds are removed and the smooth joints are defined in a direction parallel with the joint plane.

Forces and relative displacement can be expressed in local coordinates of the joint plane:

The strength of a SJ contact follows the Coulomb sliding model. The maximum possible magnitude of the shear force depends upon the coefficient of friction :

The value of the shear force is updated through recomputation at the end of each time step. If , then . The forces are updated as follows:

If , then sliding occurs at the contacts between the particles. The forces are updated through

For the SJ, the gap value between the particles can be defined. However, the SJ can only be activated when the gap value diminishes to zero. Therefore, by defining a specific gap value, joints with non-zero aperture at the initial can be modeled.

2.3. Coupled Hydromechanical Model

In PFC, the algorithm governing fluid flow is developed out of a network flow model. The fluid flow model is based on two important assumptions: (1) the flow domain is characterized by the polygons formed by straight lines linking the center point of the adjacent particles in a random plane, and (2) the flow domains are connected by the flow channels between the adjacent particles.

Figure 4 shows a schematic diagram of this fluid flow model. Due to the pressure difference in each flow domain, the fluid flows from one flow domain to another. This process is assumed as one in which the fluid passes in between two parallel plates in laminar flow. The flow process is governed by the cubic law. Therefore, the flow volume is expressed as [34]

In the fluid flow process, the pressure increment in the reservoir can be calculated with [35]

Figure 4(b) illustrates how pore pressure is formed within a rock matrix, where the size of the yellow spot is in direct proportion to the magnitude of the pore pressure. The changes to the stress state resulting from bond breakage at particle contacts are determined by the bond strength and Mohr-Coulomb failure criterion. For an open fracture, the permeability is assumed to be infinite. The fluid flow algorithm is based on two assumptions that fluid domains are connected by fluid flow channels (Figures 4(c) and 4(d)).

The fluid pressure applied to the contact between the particles results in the deformation of the reservoir. Therefore, the hydraulic aperture in the proposed model in this study is not the actual aperture of the flow path. It is an empirical value that characterizes the permeability of the rock under different confining pressures. The aperture of the flow path changes after the fluid is injected into the model (Figure 4(e)). In the fluid flow process, the aperture of the flow path is correlated with the normal stress at the contact. The relation is expressed as

Equation (10) characterizes the relation between and in hydraulic fracturing. As tends to infinity, gets close to . In (8), applies specifically to the fluid that flows from one flow domain to another. It is a microscopic quantity. Therefore, it cannot be used to calculate the permeability coefficient of the rock, which is a macroscopic quantity. It can be obtained as follows:

As the fluid pressure exceeds the tensile or shear strength of the bonded particles, cracks start to initiate between the adjacent particles. By this time, the aperture of the flow path is assumed to be infinite. This is an irreversible process, meaning the bonding between the particles loses its efficacy or is no longer effective. However, it should also be noted that is not left to be infinite. Instead, it is still set to a value to represent a different flow state.

Once the cracks are initiated, the infinite setting might undermine the stability of the simulation, delivering an adverse impact on the modeling results. This is because the sudden formation of the cracks results in the instantaneous flow of the fluids between the fluid passages and causes momentary changes in fluid pressure. In this paper, after the breakage of the particles, the fluid pressure is assumed to be the mean value of the pressures in the two flow domains prior to particle debonding [36, 37]:

3. Model Steps and Microparameter Calibration

3.1. Model Steps

Using DEM, Wang et al. [38] performed a hydraulic fracturing study on five typical coal mass models, i.e., intact coal, layered jointed coal, vertical jointed coal, orthogonal jointed coal, and synthetic jointed coal, and compared the HFN pattern produced in different models. A conceptual model for a jointed rock mass with multiple planes of nonpersistent joints is shown in Figure 5. The dimension of the model is 1.0 × 1.0 m, and the particle size falls within a range of 60 to 88 mm.

To create a coupled hydromechanical model in PFC, a confined square frame with four sets of walls is produced in the first place. Then, the shale reservoir model is established by the following steps: (1) Generate particles and arrange them randomly to fill the square frame. (2) Stop the moving particles until they are not embedded into each other. (3) Eliminate floating particles to ensure that there is no great dispersion in the magnitude of the hydraulic aperture. (4) Assemble SJM in specified zones and assemble PBM in the remaining zones. (5) Apply 10 MPa and 5 MPa to the model in the horizontal and vertical direction, respectively. (6) Excavate a borehole with a diameter of 50 mm in the center of the model and run the model until reaching equilibrium. (7) Form the flow domains in the model using an independent subprogram.

In this study, the SJM involves microparameters such as joint orientation angle , joint persistency , joint spacing , joint step angle , and joint aperture . Figure 5 clearly illustrates the above microparameters. The joint orientation angle is defined as the angle with respect to horizontal principal stress in plane , which changes from 0° to 90° at an interval of 15°. Joint persistency is defined by the joint length and the rock bridge length as

Joint spacing is defined as the straight-line distance between the centerline of the adjacent parallel joints. The joint step angle is defined as the angle of the line connecting the end of the present joint set with the starting point of the next joint set. Joint aperture is defined as the width of the SJ.

3.2. Microparameter Calibration

The microparameters of the jointed rock formations modeled in PFC cannot be directly obtained from the physical experiment. However, these microparameters can be calibrated by the physical tests performed on specimens with various joint orientation angles. Firstly, Brazilian disc tests were carried out on the samples. Secondly, sensitivity analyses were conducted on the mechanical parameters of the joint, including joint strength, ratio of cohesion to tensile strength, normal stiffness, shear stiffness, and angle range. The authors’ previous studies show that [39] (1) the strength of the SJ in the proposed model that represents the joints exerts a main effect on the strength of the reservoir; (2) the elastic modulus of the reservoir is related to the contact stiffness of the joints; and (3) the ratio of the SJ cohesion to its tensile strength mainly affects the number of the cracks generated in the reservoir, thus having further effects on the failure mode of the rock mass. Through trial and error, the authors finally obtain the microparameters that replicate the macromechanical behavior of the reservoir, which are tabulated in Tables 1 and 2.

To validate the accuracy of the microparameters, we carried out laboratory tests. The angle between the joint plane and the vertical direction is assumed as . Figure 6 summarizes the failure mode of the specimens tested at an angle ranging from 0° to 90° and presents the photos showing the distribution of the microcracks and macrocracks. Figure 7 presents the comparison results regarding failure load and splitting modulus in the physical experiment and numerical simulation. The yellow stands for the joint and the gray stands for shale matrix in figures of the second column. The comparison shows that the simulation results show good agreement with the experimental results. It should be noted that while performing the calibration on the microparameters of the shale reservoirs, we split the shale matrix into two types, i.e., brittle mineral and nonbrittle mineral, to have a model that delivers more precise results. For more particulars, readers can refer to the authors’ previous publication [39].

(a)

(b)

In Figure 6, the thickness of the core is 25 mm, which is rather small. In physical experiments, the specimens used in the Brazilian disc test are transversely isotropic shale, meaning that the specimens do not exhibit anisotropic characteristics along the -direction. The specimens in the Brazilian disc test fail under quasistatic conditions. In addition, the shale is very brittle. Therefore, from fracture initiation to failure, this process takes a very short time. The process for the fractures to evolve gradually is virtually nonexistent. In a nutshell, it is safe to say that the shale is isotropic along the -direction.

It should be noted that the 3D version of PFC is also able to model the above physical experiment. However, using PFC3D to model this experiment is very time-consuming and uneconomical. In addition, it is difficult to control the loading and boundary conditions if using PFC3D to perform the modeling. Instead, PFC2D is easy to use and can better illustrate the evolution of the fractures. Therefore, PFC2D is employed in this study.

A great number of scholars and researchers have validated the feasibility of modeling the flow dynamics in solid materials by using DEM [40–42]. However, to model hydraulic fracturing using DEM, not only should the microparameters of the matrix and joint in the shale reservoir be determined but also the microparameters governing fluid flow should be taken into account. Therefore, we performed calibration on the hydraulic fracturing breakdown pressure in a homogeneous shale reservoir. To date, no method has been available to calibrate the fluid parameters in a coupled hydromechanical model. However, Haimson and Fairhurst [43] and Fairhurst [44] put forward a regression equation for calculating the breakdown pressure where is the breakdown pressure of the reservoir, is the initial pore pressure, and is the tensile strength of the reservoir.

To validate the reliability of using the bonded particle method to model hydraulic fracturing, we set the initial pore pressure as 0 in this study. The tensile strength of the homogeneous shale reservoir is selected as 10.76 MPa. For the confining pressure in the -direction (), it is 20 and 10 MPa; the confining pressure in the -direction (). It varies from 5 to 10 MPa. The microparameters governing fluid flow in the model are listed in Table 3. The comparison between the numerical simulation and the values derived from the theoretical equation under different confining pressure ratios are shown in Figure 8. Results demonstrate that simulation results fit well with the theoretical values regardless of a certain disagreement. The disagreement can be reasonably explained. Equation (14) is derived based on elastic mechanics that assumes the rock mass as elastic, isotropic, and homogeneous. However, this assumption is against the reality.

4. Results and Analysis

To understand the geometrical parameters of the joints on the breakdown pressure, the number of hydraulic fractures, and the pattern of the generated HFN, we performed a series of sensitivity analyses using the calibrated model. In each study group, only the factor of interest is different and the other factors such as joint plane configuration and hydraulic parameters are kept the same. The fluid injection time is kept constant at 10 s. In this study, the factors of interest include joint orientation angle , joint persistency , joint spacing , joint step angle , and joint aperture .

The coupled hydromechanical model in this study can calculate the contact force between the particles, which is used to extrapolate the equivalent continuum stress. As the amount of the injection fluid increases, the bond between the particles breaks, resulting in HF. The cracks in the formation can be used for quantitative or qualitative analysis. To evaluate the macrocrack through monitoring the accumulative quantity of the microcracks has been proven effective by many scholars and researchers [45, 46]. Using an independent subprogram, various types of hydraulic cracks can be identified and recorded, i.e., tensile hydraulic fracture caused by shale matrix (THFSM), shear hydraulic fracture caused by shale matrix (SHFSM), tensile hydraulic fracture caused by joint plane (THFJP), and shear hydraulic fracture caused by joint plane (SHFJP).

In a homogeneous shale reservoir, HF in general propagates along the direction of the maximum principal stress. However, this is not applicable in jointed shale reservoirs. The modeling results demonstrate that the geometrical parameters of the joints have a significant effect on the HFN pattern formed in jointed shale reservoirs.

In general, after hydraulic fracturing, the HFN pattern in jointed shale reservoirs fall into four distinct types, i.e., crossing mode, tip-to-tip mode, step path mode and opening mode. Figure 9 illustrates the four types of HFN patterns characterized by crack distribution and stress concentration, where more dark arrows indicate greater stress concentration. Figure 10 presents the evolution of the four types of cracks with the injection time. (1)Crossing mode. This HFN pattern results from the stress concentration near the shale matrix between the joints and the slight stress between the rock bridges (Figure 9(a)). In the meantime, a small amount of shear HF initiates from the internal zones of the joints and propagates to the boundary along the joint orientation. In this scenario, HF no longer propagates in a direction normal to the principal stress. Instead, it propagates in a direction nearly vertical to the joints. In this HFN pattern, THFSM gains the dominance over other types of cracks, accounting for about 60% in number. As some joints are activated during HF propagation, the HF generated in the joints accounts for 30% in number (Figure 10(a)).(2)Tip-to-tip mode. This HFN pattern results from the stress concentration near the joint tips and the presence of a low stress zone between the rock bridges, as shown in Figure 9(b). The tensile HF initiates near the joint tips and coalesces with the HF initiated near other joint tips. The shear HF generates and propagates to the boundary along the joint. In this pattern, the number of THFSM further outweighs that of the other three HFs, which is about two times of the number of the cracks combined (Figure 10(b)).(3)Step path mode. At the beginning of fluid injection, a slight amount of HF initiates near the joint tips or within the internal zone. However, as more fluids are injected into the borehole, HF does not traverse the joints or propagate along the joint tip. Instead, along the joints, it initiates from the other end of the joints, resulting in joint slip (Figure 9(c)). In this pattern, the stress concentration zone forms in a direction parallel with the principal stress. Due to the propagation of HF along the joint, the number of SHFJP increases. Nonetheless, throughout the whole injection process, THFSM still maintains its role as the most dominant crack type (Figure 10(c)).(4)Opening mode. As joint aperture increases, the hydraulic pressure induced by fluid injection is completely absorbed by the joints. As a result, more fluid flows into the joints, generating THFJP. Stress concentration occurs only near the borehole, thus incapacitating the generation of HF in the shale matrix (Figure 9(d)). Due to the inherent low strength of the distant joints, the perturbation caused by the injection pressure generates HF. Therefore, SHFJP grows significantly in number, usurping the dominance of the THFSM (Figure 10(d)).

(a)

(b)

(c)

(d)

(a)

(b)

(c)

(d)

4.1. Effect of Joint Orientation Angle on Hydraulic Fracture Network Propagation

To investigate the effect of , we keep other microparameters constant, i.e., , , , and . Then, we performed numerical simulation with changing from 0° (horizontal joint) to 90° (vertical joint) at an interval of 15°. The modeling results are shown in Figures 11–13, which separately present the effect of on the breakdown pressure, the number of HF, and the HFN pattern. The four colors in Figure 13 represent the tensile and shear cracks generated in the shale matrix and the joint. To study to what extent the breakdown pressure in jointed shale reservoir decreases as opposed to that in homogeneous shale reservoir, the breakdown pressure in jointed shale reservoir is normalized to the breakdown pressure in a homogeneous shale reservoir, which is 13.22 MPa.

During hydraulic fracturing, reaches its maximum value when is 90° (vertical joint). While is 0°, the joint is in a direction parallel with the maximum principal stress. In this case, stress concentrates near the joint tip, and the tensile HF initiates near the joint tips and coalesces with the HF initiated near other joint tips. In the meantime, SHFJP is generated in the joints. Therefore, this HFN pattern is characterized by tip-to-tip mode. As can be seen from Figure 12(a), in this case, the number of HF reaches a peak, amounting to about 2000. When and 30°, the joint is angled with respect to the direction of the maximum principal stress. In this case, stress concentrates near the joint tip and in the contacting areas between the shale matrix. Moreover, the breakdown pressure declines as shown in Figure 11. The generated HF forms a pattern that is combined by tip-to-tip mode and crossing mode. As increases to 45° or 60°, stress concentrates within the shale matrix. Under the combined action of the joint orientation and the maximum principal stress, tensile HF propagates in a direction nearly vertical to the joint, forming a HFN pattern characterized by the crossing mode. In this pattern, HF firstly initiates in the internal zones of the joints, followed by the propagation and coalesce of HF between the joints. Finally, HF penetrates the joints. In this case, the breakdown pressure is roughly about 30% of the breakdown pressure . According to the Coulomb criterion, the jointed shale reservoir has its minimum tensile strength when varies from 45° to 60°. In addition, the previous studies have validated that the breakdown pressure is positively correlated with the tensile strength of the shale reservoir [47]. Therefore, theoretical solution and numerical simulation deliver a consistent result regarding the breakdown pressure, which further validates the correctness of the microparameters used in this study.

(a)

(b)

After reaches 75°, stress is less concentrated near the joint tip. Therefore, the breakdown pressure increases and is greater than the breakdown pressure while is 0°. When HF interacts with the joints, the joints are activated. As a result, HF initiates at the other end of the joint and propagates. However, when or 90°, cracks do not cover a wide area in the shale reservoir except for areas surrounding the injection borehole. Therefore, the number of HF is smaller than that when is 0° or 15°. However, it should be noted that the proportion of different types of HF remains unchanged despite the changes in the number of HF with (Figure 12(b)).

The modeling results demonstrate that exerts a significant influence on the breakdown pressure, the number of HF, and the HFN pattern. Therefore, will also be considered while investigating the effect of other microparameters of the joints.

4.2. Effect of Joint Persistency on Hydraulic Fracture Network Propagation

To study the effect of , we keep other microparameters constant except (, , , and changes from 0° to 90° at an interval of 15°). is changed by varying the length of the joint while the length of the rock bridge is kept unchanged (). Figures 14–16 present the effect of on the breakdown pressure, the number of HF, and the HFN pattern. Results show that as increases, has a more pronounced effect on the breakdown pressure and the HFN pattern, as shown in Figures 14 and 16. exerts a more apparent effect on the number of HF when varies from 45° to 60° as shown in Figure 15.

As shown in Figure 16, the HF pattern in all the shale reservoirs is characterized by tip-to-tip mode when is 0°. The number of HF varies slightly from one shale reservoir to another. In this scenario, an increase in leads to a decrease in breakdown pressure (Figure 14). This is due mainly to the aggravation of stress concentration near the joint tips. In fact, stress has already concentrated around the joint tips even when is 0.6. Therefore, the increase in further aggravates stress concentration near the joint tips. This is why the HFN pattern remains unchanged while the breakdown pressure diminishes. As increases to 15° and 30°, with a small ( or 0.7), the HFN pattern formed is a combination of the tip-to-tip mode and the crossing mode. However, the HFN pattern is still characterized by the tip-to-tip mode when changes to 0.8. As a larger results in a greater joint density, stress is more easily concentrated around the joint tips. Therefore, it is easier for HF to initiate at the joint tips and to propagate further (Figure 16). As a result, the breakdown pressure significantly drops (Figure 14). According to the theory of fracture mechanics, greater joint length results in greater stress intensity factor. Therefore, greater stress concentration occurs around the joint tips. This is why as increases the breakdown pressure drops.

As increases to 45° or 60°, the HFN pattern is characterized by the crossing mode. an increase in leads to an increase in the number of HF. In addition, the HF generated in joints accounts for a greater proportion, with SHFJP in particular. When reaches the extreme values ( and ), the difference in the number of HF amounts to over 400 and HF propagates in a direction nearly vertical to the joint, as shown in Figure 15. As a result, stress concentration gradually spreads to rock bridge areas, causing a difference in the breakdown pressure and the number of HF. Stress is less concentrated around the joint tips when is equal to or greater than 75°. Stress is more evenly distributed across the rock bridge area. Therefore, does not exert an apparent effect on the breakdown pressure. Its effect on the number of HF also weakens. The HFN pattern converts from the crossing mode to the step path mode.

Therefore, we can conclude that (1) when keeping the joint orientation angle and other parameters unchanged, an increase in leads to a decrease in the breakdown pressure and the number of HF; (2) an increase in causes no perceptible changes to the proportion of HF; and (3) though has an effect on the HFN pattern, the effect is not as significant as that of .

4.3. Effect of Joint Spacing on Hydraulic Fracture Network Propagation

To study the effect of , we keep other microparameters constant except (, , , and changes from 0° to 90° at an interval of 15°). In this study, is selected as 60, 120, and 180 mm, respectively. Figures 17–19 present the effect of on the breakdown pressure, the number of HF, and the HFN pattern. Results show that as increases, breakdown pressure grows as well and has a more pronounced effect on the number of HF. When increases to the largest value (), after hydraulic fracturing, the HFN pattern in all the reservoirs is characterized by the step path mode.

As shown in Figure 17, when is 0°, the increase in weakens the effect of joints on stress concentration in the shale reservoir. Therefore, the breakdown pressure grows gradually. However, due to the decline in the number of joints, it is more difficult for the HF to propagate and coalesce. As a result, the number of the HF in the reservoir as a whole diminishes. As increases to 15° and 30°, the HFN pattern varies from one reservoir to another. This indicates that exerts a more obvious effect than (compare Figures 16–19). When is small (), the stress concentration zones around the joint tips overlap. As a result, it is easy for HF to propagate along the joint tip, thus creating a HFN pattern characterized by the tip-to-tip mode. As reaches its maximum (), stress concentration zones around the joint tips are isolated from each other. As a result, after HF initiates at the joint tips, it propagates along the joints with a low strength. Finally, some zones show a HFN pattern characterized by the step path mode. This indicates that, with a small , an increase in results in changes to the HFN pattern. The HFN pattern converts from tip-to-tip mode to step path mode.

As increases to 45° or 60°, a different leads to different HFN patterns. This is a combined result of the changes in and . Due to the variation of , the breakdown pressure bottoms out at the same . In the meantime, crack propagation gains momentum and the joints collapse before they are activated (Figure 19). When gets larger ( or 90°), the effect of on the number of HF and the HFN pattern is attenuated. is the primary reason for the ubiquitous presence of the step path mode in all the shale reservoirs. However, changes in induce changes to the strength of the shale reservoir. As a result, the breakdown pressure still increases with increasing . When is 90° and is 180 mm, the breakdown pressure in the jointed shale reservoir is over 0.9 times of that in the homogeneous shale reservoir.

Therefore, we can conclude that (1) when keeping the joint orientation angle and other parameters unchanged, an increase in leads to a significant increase in the breakdown pressure and (2) when falls within the range of 45°~60°, variation in can cause changes to the number of HF and the HFN pattern. An increase in can lead to reduction of the proportion of HF induced by joint activation.

4.4. Effect of Joint Step Angle on Hydraulic Fracture Network Propagation

Joint step angle characterizes the relative position between the parallel joints (Figure 5). To study the effect of , we keep other microparameters constant except (, , and . varies from 0° to 90° at an interval of 15°). In this study, is selected as 60°, 90°, and 117°, respectively. Figures 20–22 present the effect of on the breakdown pressure, the number of HF, and the HFN pattern. As can be seen from Figure 20, has a negligible effect on breakdown pressure. The reason behind this imperceptible effect is that joint strength and joint density are not relevant to .

As shown in Figure 21, when is 0°, 15°, and 30°, there is a variation in the total number of HF. However, the proportion of different types of HF remains virtually unchanged. This is due to the change in HFN pattern in the reservoir. Following the injection of the fluids into the borehole, rocks surrounding the borehole start to fracture. Though joint strength and joint density remain unchanged, a different results in changes to the HF initiation locality and propagation direction. Therefore, different HFN patterns are generated (Figure 22).

As can be seen from Figure 22, when is equal to or greater than 45°, it has a dominant effect on the HFN pattern. At a constant , the HFN pattern does not vary with any more. Stress still concentrates around the joint tips and within the central section of the joints, irrespective of the changes of . The number of HF as well as the proportion of different types of HF are not sensitive to changes of . It should be noted that when is 60°, there is variation in the number of HF (Figure 21). This is due to the low strength of the reservoir when is 60°. In fact, the reservoir has minimum strength when is 60°. Therefore, cracks easily initiate, coalesce, and propagate in shale matrices and joints.

Therefore, we conclude that (1) when is not greater than 30°, and both have an effect on HFN pattern; (2) when is not smaller than 30°, has a dominant effect on the HFN pattern; and (3) at the same , has a negligible effect on breakdown pressure.

4.5. Effect of Joint Aperture on Hydraulic Fracture Network Propagation

To study the effect of , we keep other microparameters constant except (, , and , with changing from 0° to 90° at an interval of 15°). In this study, is selected as 1.2 × 10⁻⁴, 1.2 × 10⁻³, and 1.2 × 10⁻² m, respectively. Figures 23–25 present the effect of on the breakdown pressure, the number of HF, and the HFN pattern. Results show that as increases from 1.2 × 10⁻⁴ m to 1.2 × 10⁻³ m, breakdown pressure decreases, but not in an obvious manner, and there are no significant changes to the HFN pattern. Due to the increase in , the shear HF grows, with its proportion in HF increasing by about 10%.

However, as increases to 1.2 × 10⁻² m, significant changes take place in the breakdown pressure, the number of HF, and the HFN pattern. In this scenario, does not have a dominant effect on breakdown pressure any more. The normalized breakdown pressure falls within a range of 0.15~0.25. In addition, the proportion of tensile HF decreases by about 25%. The above changes are attributable to the opening mode generated after hydraulic fracturing. Due to the increase in , after injection, most of the fluids flow into the joints and only a slight amount of shale matrices surrounding the injection borehole are fractured. As the injection continues, distant joints start to slip and shift and HF is generated. As for joints that are more distant, it is less likely for them to be activated.

Therefore, we can conclude that (1) when keeping the joint orientation angle and other parameters unchanged, a decrease in leads to an increase in the breakdown pressure and the number of the cracks and (2) changes of result in the variation in the proportion of different types of hydraulic fractures.

5. Conclusions

In this study, a discrete element method-based coupled hydromechanical model is developed in particle flow code to explore the effect of the geometrical parameters of the joints on the breakdown pressure, the number of hydraulic fracture, and the hydraulic fracture network pattern in jointed shale reservoirs. The hydraulic fracture network pattern in the shale reservoir resulting from hydraulic fracturing can be roughly divided into four types, i.e., crossing mode, tip-to-tip mode, step path mode, and opening mode.

When is 0°, except for , most of the shale reservoirs have a hydraulic fracture network pattern that is characterized by tip-to-tip mode after hydraulic fracturing. When increases to 15° or 30°, the hydraulic fracture network pattern changes from tip-to-tip mode to crossing mode, as stress concentrates near the joint tip and in the contacting areas between the shale matrices. In this scenario, and both have an effect on the hydraulic fracture network pattern. Due to the change in the hydraulic fracture network pattern, the breakdown pressure and the number of cracks gradually decrease.

When is 45° or 60°, the crossing mode gains dominance, and the hydraulic fracture network pattern depends mainly on and . In addition, the magnitude of the breakdown pressure and the number of the cracks reach the lowest level. As reaches 75° or 90°, the step path mode is ubiquitous, and the magnitude of the breakdown pressure and the number of the cracks both increase. The breakdown pressure increases and is greater than the breakdown pressure when is 0°. In this scenario, the hydraulic fracture network pattern is largely dependent on .

Under the same , either decrease in and or increase in leads to the rise in the breakdown pressure and the number of the cracks. It is also found that changes in and result in the variation in the proportion of different types of hydraulic fractures. In particular, the opening mode of the hydraulic fracture network is observed when increases to 1.2 × 10⁻² m.

Nomenclature

Greek Letters

:	The aperture of the flow path
, :	The aperture values at infinite and zero
:	Joint orientation angle
:	Joint step angle
:	The decay speed of the aperture as increases, which usually equals −0.15 [40]
:	The fluid dynamic viscosity
:	The angle of dilation

Roman Symbols

:	The area of the smooth joint particle
:	Joint spacing
:	Joint aperture
, :	Normal force and relative displacement
, :	Shear force and relative displacement vectors
:	An updated value of shear force
:	Joint persistency
:	Bulk modulus of the fluid
, :	The normal and shear stiffness of smooth joint
, :	The pressure in different flow domains
:	The length of the flow path
, :	Fluid pressures in two domains
:	One time step
, :	The tangential and normal unit vectors
, :	The normal relative displacement and shear relative displacement vectors of smooth joint
:	The total volume of the shale reservoirs
:	Pore volume
:	Variation of the pore volume.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research is supported by the Fundamental Research Funds for the Central Universities (2017XKZD06).

References

K. B. Gregory, R. D. Vidic, and D. A. Dzombak, “Water management challenges associated with the production of shale gas by hydraulic fracturing,” Elements, vol. 7, no. 3, pp. 181–186, 2011.
View at: Publisher Site | Google Scholar
D. Rahm, “Regulating hydraulic fracturing in shale gas plays: the case of Texas,” Energy Policy, vol. 39, no. 5, pp. 2974–2981, 2011.
View at: Publisher Site | Google Scholar
Q. Wang, X. Chen, A. N. Jha, and H. Rogers, “Natural gas from shale formation – the evolution, evidences and challenges of shale gas revolution in United States,” Renewable and Sustainable Energy Reviews, vol. 30, pp. 1–28, 2014.
View at: Publisher Site | Google Scholar
B. Huang, Y. Wang, and S. Cao, “Cavability control by hydraulic fracturing for top coal caving in hard thick coal seams,” International Journal of Rock Mechanics and Mining Sciences, vol. 74, pp. 45–57, 2015.
View at: Publisher Site | Google Scholar
B. Yuan, D. A. Wood, and W. Yu, “Stimulation and hydraulic fracturing technology in natural gas reservoirs: theory and case studies (2012–2015),” Journal of Natural Gas Science and Engineering, vol. 26, pp. 1414–1421, 2015.
View at: Publisher Site | Google Scholar
C. Liu, F. Shi, Y. Zhang et al., “High injection rate stimulation for improving the fracture complexity in tight-oil sandstone reservoirs,” Journal of Natural Gas Science and Engineering, vol. 42, pp. 133–141, 2017.
View at: Publisher Site | Google Scholar
Y. Lu, Z. Ge, F. Yang, B. Xia, and J. Tang, “Progress on the hydraulic measures for grid slotting and fracking to enhance coal seam permeability,” International Journal of Mining Science and Technology, vol. 27, no. 5, pp. 867–871, 2017.
View at: Publisher Site | Google Scholar
Y. Su, Q. Zhang, W. Wang, and G. Sheng, “Performance analysis of a composite dual-porosity model in multi-scale fractured shale reservoir,” Journal of Natural Gas Science and Engineering, vol. 26, pp. 1107–1118, 2015.
View at: Publisher Site | Google Scholar
J. Hou, M. Qiu, X. He, C. Guo, M. Wei, and B. Bai, “A dual-porosity-stokes model and finite element method for coupling dual-porosity flow and free flow,” SIAM Journal on Scientific Computing, vol. 38, no. 5, pp. B710–B739, 2016.
View at: Publisher Site | Google Scholar
T. Fan and G. Zhang, “Laboratory investigation of hydraulic fracture networks in formations with continuous orthogonal fractures,” Energy, vol. 74, pp. 164–173, 2014.
View at: Publisher Site | Google Scholar
Y. Wang, X. Li, and B. Zhang, “Analysis of fracturing network evolution behaviors in random naturally fractured rock blocks,” Rock Mechanics and Rock Engineering, vol. 49, no. 11, pp. 4339–4347, 2016.
View at: Publisher Site | Google Scholar
Z. Chong, X. Li, X. Chen, J. Zhang, and J. Lu, “Numerical investigation into the effect of natural fracture density on hydraulic fracture network propagation,” Energies, vol. 10, no. 7, p. 914, 2017.
View at: Publisher Site | Google Scholar
C. E. Renshaw and D. D. Pollard, “An experimentally verified criterion for propagation across unbounded frictional interfaces in brittle, linear elastic materials,” International Journal of Rock Mechanics and Mining Sciences, vol. 32, no. 3, pp. 237–249, 1995.
View at: Publisher Site | Google Scholar
H. Gu and X. Weng, “Criterion for fractures crossing frictional interfaces at non-orthogonal angles,” in 44th U.S. Rock Mechanics Symposium and 5th U.S.-Canada Rock Mechanics Symposium, Salt Lake City, UT, USA, June 2010.
View at: Google Scholar
H. Gu, X. Weng, J. B. Lund, M. G. Mack, U. Ganguly, and R. Suarez-Rivera, “Hydraulic fracture crossing natural fracture at nonorthogonal angles: a criterion and its validation,” SPE Production & Operations, vol. 27, no. 1, pp. 20–26, 2012.
View at: Publisher Site | Google Scholar
G. D. Anderson, “Effects of friction on hydraulic fracture growth near unbonded interfaces in rocks,” Society of Petroleum Engineers Journal, vol. 21, no. 1, pp. 21–29, 1981.
View at: Publisher Site | Google Scholar
X. Weng, “Modeling of complex hydraulic fractures in naturally fractured formation,” Journal of Unconventional Oil and Gas Resources, vol. 9, pp. 114–135, 2015.
View at: Publisher Site | Google Scholar
W. Xu, M. J. Thiercelin, U. Ganguly et al., “Wiremesh: a novel shale fracturing simulator,” in International Oil and Gas Conference and Exhibition in China, Beijing, China, June 2010.
View at: Publisher Site | Google Scholar
B. R. Meyer and L. W. Bazan, “A discrete fracture network model for hydraulically induced fractures: theory, parametric and case studies,” in SPE Hydraulic Fracturing Technology Conference, The Woodlands, TX, USA, January 2011.
View at: Publisher Site | Google Scholar
X. Weng, O. Kresse, C. E. Cohen, R. Wu, and H. Gu, “Modeling of hydraulic-fracture-network propagation in a naturally fractured formation,” SPE Production & Operations, vol. 26, no. 4, pp. 368–380, 2011.
View at: Publisher Site | Google Scholar
J. E. Olson and A. D. Taleghani, “Modeling simultaneous growth of multiple hydraulic fractures and their interaction with natural fractures,” in SPE Hydraulic Fracturing Technology Conference, The Woodlands, TX, USA, January 2009.
View at: Publisher Site | Google Scholar
W. S. Dershowitz and C. Fidelibus, “Derivation of equivalent pipe network analogues for three-dimensional discrete fracture networks by the boundary element method,” Water Resources Research, vol. 35, no. 9, pp. 2685–2691, 1999.
View at: Publisher Site | Google Scholar
Z. Chong, S. Karekal, X. Li, P. Hou, G. Yang, and S. Liang, “Numerical investigation of hydraulic fracturing in transversely isotropic shale reservoirs based on the discrete element method,” Journal of Natural Gas Science and Engineering, vol. 46, pp. 398–420, 2017.
View at: Publisher Site | Google Scholar
F. Hamidi and A. Mortazavi, “A new three dimensional approach to numerically model hydraulic fracturing process,” Journal of Petroleum Science and Engineering, vol. 124, pp. 451–467, 2014.
View at: Publisher Site | Google Scholar
B. Zhang, X. Li, Z. Zhang, Y. Wu, Y. Wu, and Y. Wang, “Numerical investigation of influence of in-situ stress ratio, injection rate and fluid viscosity on hydraulic fracture propagation using a distinct element approach,” Energies, vol. 9, no. 3, p. 140, 2016.
View at: Publisher Site | Google Scholar
B. Hou, M. Chen, W. Cheng, and C. Diao, “Investigation of hydraulic fracture networks in shale gas reservoirs with random fractures,” Arabian Journal for Science and Engineering, vol. 41, no. 7, pp. 2681–2691, 2016.
View at: Publisher Site | Google Scholar
E. M. Llanos, R. G. Jeffrey, R. Hillis, and X. Zhang, “Hydraulic fracture propagation through an orthogonal discontinuity: a laboratory, analytical and numerical study,” Rock Mechanics and Rock Engineering, vol. 50, no. 8, pp. 2101–2118, 2017.
View at: Publisher Site | Google Scholar
S. Heng, Y. Guo, C. Yang, J. J. K. Daemen, and Z. Li, “Experimental and theoretical study of the anisotropic properties of shale,” International Journal of Rock Mechanics and Mining Sciences, vol. 74, pp. 58–68, 2015.
View at: Publisher Site | Google Scholar
C. Yan, H. Zheng, G. Sun, and X. Ge, “Combined finite-discrete element method for simulation of hydraulic fracturing,” Rock Mechanics and Rock Engineering, vol. 49, no. 4, pp. 1389–1410, 2016.
View at: Publisher Site | Google Scholar
Z. Yushi, M. Xinfang, Z. Shicheng, Z. Tong, and L. Han, “Numerical investigation into the influence of bedding plane on hydraulic fracture network propagation in shale formations,” Rock Mechanics and Rock Engineering, vol. 49, no. 9, pp. 3597–3614, 2016.
View at: Publisher Site | Google Scholar
W. Chen and H. Konietzky, “Simulation of heterogeneity, creep, damage and lifetime for loaded brittle rocks,” Tectonophysics, vol. 633, pp. 164–175, 2014.
View at: Publisher Site | Google Scholar
N. Cho, C. D. Martin, and D. C. Sego, “A clumped particle model for rock,” International Journal of Rock Mechanics and Mining Sciences, vol. 44, no. 7, pp. 997–1010, 2007.
View at: Publisher Site | Google Scholar
M. Pierce, R. Cundall, D. Potyondy, and D. M. Ivars, “A synthetic rock mass model for jointed rock,” in Proceedings of the 1st Canada-US Rock Mechanics Symposium - Rock Mechanics Meeting Society’s Challenges and Demands, pp. 341–349, Vancouver, Canada, May 2007.
View at: Publisher Site | Google Scholar
A. S. Berman, “Laminar flow in channels with porous walls,” Journal of Applied Physics, vol. 24, no. 9, pp. 1232–1235, 1953.
View at: Publisher Site | Google Scholar
Itasca Consulting Group, PFC^2D Manual, Version 5.0, Minneapolis, MN, USA, 2014.
A. Al-Busaidi, “Distinct element modeling of hydraulically fractured lac du bonnet granite,” Journal of Geophysical Research, vol. 110, B6, 2005.
View at: Publisher Site | Google Scholar
X. Zhao and R. Paul Young, “Numerical modeling of seismicity induced by fluid injection in naturally fractured reservoirs,” Geophysics, vol. 76, no. 6, pp. WC167–WC180, 2011.
View at: Publisher Site | Google Scholar
T. Wang, W. Hu, D. Elsworth et al., “The effect of natural fractures on hydraulic fracturing propagation in coal seams,” Journal of Petroleum Science and Engineering, vol. 150, pp. 180–190, 2017.
View at: Publisher Site | Google Scholar
Z. Chong, X. Li, P. Hou et al., “Numerical investigation of bedding plane parameters of transversely isotropic shale,” Rock Mechanics and Rock Engineering, vol. 50, no. 5, pp. 1183–1204, 2017.
View at: Publisher Site | Google Scholar
J. S. Yoon, A. Zang, and O. Stephansson, “Numerical investigation on optimized stimulation of intact and naturally fractured deep geothermal reservoirs using hydro-mechanical coupled discrete particles joints model,” Geothermics, vol. 52, pp. 165–184, 2014.
View at: Publisher Site | Google Scholar
A. Abdollahipour, M. Fatehi Marji, A. Yarahmadi Bafghi, and J. Gholamnejad, “Dem simulation of confining pressure effects on crack opening displacement in hydraulic fracturing,” International Journal of Mining Science and Technology, vol. 26, no. 4, pp. 557–561, 2016.
View at: Publisher Site | Google Scholar
J. Zhou, L. Zhang, Z. Pan, and Z. Han, “Numerical investigation of fluid-driven near-borehole fracture propagation in laminated reservoir rock using PFC^2D,” Journal of Natural Gas Science and Engineering, vol. 36, pp. 719–733, 2016.
View at: Publisher Site | Google Scholar
B. Haimson and C. Fairhurst, “Initiation and extension of hydraulic fractures in rocks,” Society of Petroleum Engineers Journal, vol. 7, no. 3, pp. 310–318, 1967.
View at: Publisher Site | Google Scholar
C. Fairhurst, “Stress estimation in rock: a brief history and review,” International Journal of Rock Mechanics and Mining Sciences, vol. 40, no. 7-8, pp. 957–973, 2003.
View at: Publisher Site | Google Scholar
J. Yoon, “Application of experimental design and optimization to pfc model calibration in uniaxial compression simulation,” International Journal of Rock Mechanics and Mining Sciences, vol. 44, no. 6, pp. 871–889, 2007.
View at: Publisher Site | Google Scholar
S.-Q. Yang, Y. H. Huang, H. W. Jing, and X. R. Liu, “Discrete element modeling on fracture coalescence behavior of red sandstone containing two unparallel fissures under uniaxial compression,” Engineering Geology, vol. 178, pp. 28–48, 2014.
View at: Publisher Site | Google Scholar
T. Ito and K. Hayashi, “Physical background to the breakdown pressure in hydraulic fracturing tectonic stress measurements,” International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, vol. 28, no. 4, pp. 285–293, 1991.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2018 Zhaohui Chong et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

772

Downloads

1021

Citations