Pore pressure effects on fracture net pressure and hydraulic fracture containment: Insights from an empirical and simulation approach

Pore pressure and its relationship with fracture net pressure has been reported qualitatively from both field and experimental observations. From a modeling perspective, the ubiquitously used pseudo 3D (P3D) models that are based on linear elastic fracture mechanics (LEFM) do not include the effect of reservoir depletion (or overpressure). Models that utilize effective stress as propagation criteria with a cohesive zone description, introduce the pore pressure directly into the simulation and hence can potentially capture the effect of pore pressure on fracture propagation. This work investigates the effect of pore pressure on hydraulic fracturing net pressure and geometry using empirical and numerical simulation approaches. We carried out an analysis of more than 400 datafrac injections spanning a wide range of geological ages and depositional environments in order to investigate the relationship between observed net pressure and reservoir pore pressure. The net fracture propagation pressure from the fracture treatment analysis was seen to be correlated with the effective stress in the reservoir. Fracture propagation simulations were performed using a coupled finite element – finite difference fracture simulator. The code uses a cohesive zone model (CZM) to describe fracture propagation. Four different effective stress scenarios were used to study the effect of effective stress on net pressure. The simulation results closely match the empirical relation between net pressure and effective stress as obtained from the analysis of actual frac treatment data. It is observed from the simulations that the magnitude of the effective stress also has an effect on the fracture geometry with a high effective stress leading to wider, shorter and more radial fractures. The derived empirical correlation is hence useful as a fracture design parameter. The datafrac net pressure diagnostics workflow in the pseudo 3D models can incorporate local tip pore pressure as a new pressure matching parameter. The pore pressure effect can thus explain high net pressures routinely observed in frac operations and also as a containment mechanism.


Introduction
Pore pressure changes affect the minimum stress, which results in well known effects on fracture propagation pressure and fracture height growth. Field fracturing data had qualitatively shown that hydraulic fracturing net pressure observed from frac jobs performed in depleted reservoirs are much higher than was initially reported from treatments done in virgin reservoir condition (Shaoul et al., 2003). It was also observed that the closure stress in the pay has decreased on depletion. Schmitt and Zoback (1989) conducted laboratory tests on the effect of pore pressure on tensile failure. Bruno and Nakagawa (1991) presented theoretical analysis based on Griffith strain energy release criteria that suggested that tensile fracture is controlled by general effective stress.
Experiments performed by them in rocks with induced pore pressure gradients showed that fracture orientation is biased towards regions of higher pore pressure or lower effective stress and the fracture initiation pressure is less when fracture tip is in the vicinity of these regions. It was argued that the observed differences in pore pressure magnitude around the crack tip (local pore pressure) and its effect on fracture orientation and direction was supportive of an effective stress law for tensile failure in porous rocks. Visser (1998) also performed a series of experiments on the pore pressure effect on extensile fracturing of saturated porous materials. Reconciling high observed net pressure with model net pressure has been an issue in fracture modeling for decades. The magnitude of net pressure during hydraulic fracturing is influenced by multiple factors and this is depicted in Fig. 1. A high net pressure can be explained by linking its dependence to simple physical parameters like injection rate, rock properties, fluid viscosity etc. However more complex mechanisms may also be at play like hydraulic fracture interaction with natural fractures, multiple fracture propagation and growth through layer interfaces. The commonly used pseudo 3D fracture models that use stress intensity factor based criteria for the fracture propagation problem do not explicitly include pore pressure. Recent experimental work indicates that the LEFM approach for fracture propagation is not always valid and a cohesive zone description is necessary (Van Dam and de Pater, 2001;2002). Such cohesive zone models (CZM) would then use effective stress at the tip exceeding tensile strength as a propagation criterion which introduces the effect of pore pressure. It is worth investigating this relatively simple effect in order to explain high observed net pressure and its influence on fracture geometry. This paper attempts to quantify the effects of this simple mechanism by comparing net pressure from actual fracture treatment data to that predicted using a coupled cohesive zone model.

Fracture treatment analysis
To investigate the pore pressure effect, the criteria for treatment selection and analysis was such that the pore pressure effect would be dominant over other factors such as tortuosity, fluid friction changes by proppant and fracture containment. Only those wells that have an inclination of 10 or less are chosen for the analysis. This is to eliminate the effect of multiple fracture growth and near-wellbore tortuosity that occurs when fractures initiate from the wellbore perf orientation and align to preferred fracture plane orientation, which is not aligned with the wellbore. The basis for this threshold were experimental studies on differences between perforation orientation and minimum stress direction (Behrmann andElbel, 1991Abass et al., 1994;Weijers, 1995). In order to eliminate the effects of fluid friction owing to fluid viscosity, proppant transport and tip screen out, only breakdown and minifrac injections were considered for the analysis.
A total of 421 datafrac treatments spanning 13 fields were studied to investigate the pore pressure effect (Table 1). The treatments include a wide range of reservoir depths, pore pressures, geological settings and geological age. All the fields however were clastics from different depositional environments. Since most of the cases involved two datafrac injections (breakdown and minifrac) followed by the main propped treatment, the main treatment pressure data (additional 215 jobs) was also collected in order to compare with the level of net pressure from the datafracs. The net pressure is defined as the difference between pressure inside the fracture and the closure stress.
The effective closure stress is defined as the difference between the closure stress and the pore pressure.
Using the Instantaneous Shut-In Pressure (ISIP) from the measured pressure data and the closure stress obtained using the graphical method from G-Function plots (Nolte, 1979) (Barree and Mukherjee, 1996), the net pressure was obtained for each treatment, at the end of each injection. Using this value and the pore pressure data, the Terzaghi minimum principal effective stress (or the fracture closure stress) was also calculated. The net pressure variation in time is an important diagnostic quantity that is recorded during the hydraulic fracturing process using either surface or bottomhole gauges. In the absence of microseismic data the net pressure is usually the only parameter that can give information on fracture geometry. An idealized schematic of fracture pressure variation with time for a hydraulic fracturing injection is depicted in Fig. 2. The breakdown pressure is depicted as the peak of the fracture pressure plot which is then followed by fracture propagation till pumping is  stopped. The ISIP point and the closure stress points are also shown in the pressure decline curve after pumping is terminated. The graphical closure stress determination is based on plotting three parameters: the bottomhole pressure, the derivative of the pressure with respect to the G-function, dP/dG and the superposition derivative, G dP/ dG. This type of plot is indicative of the leakoff mechanism as well. Example plots for determination of closure stress for four different types of leakoff behavior is depicted in Appendix (Figs. A1-A4). In the case of normal leakoff behavior, the hydraulic fracture area is constant and leakoff occurs through a homogenous rock matrix. Pressure derivative is constant and the superposition derivative lies on a straight line through the origin. The fracture closure is identified as the point when the superposition derivative data deviates downward from the straight line ( Fig. A1). In the case of pressure dependent leakoff, the superposition derivative exhibits a hump owing to the presence of fissures. Subsequent to fissure opening, the superposition derivative exhibits similar behavior to normal leakoff and the fracture closure stress is again identified as the point when the curve starts to deviate downward from the straight line (Fig. A2). In the case of fracture height recession, the superposition derivative exhibits a concave upward trend. The point where the concave upward trend of the G dP/dG curve ends and intersects the normal leakoff line is the fracture closure pressure (Fig. A3). In the case of fracture tip extension, when the fracture continues to grow even after shut-in, a straight line fit for the superposition derivative exhibits an intercept with the y-axis (Fig. A4). From the database of treatments, we depict three treatments from Field G, D and L, which corresponds to high effective stress (Fig. A5), intermediate effective stress (Fig. A6) and low effective stress (Fig. A7), respectively.
Plotting the treatments for closure stress gradient against pore pressure gradient as in Fig. 3, it can be observed that the selection covers a wide range of pore pressure scenarios. It is evident that the high pressure reservoirs are also high stress. All the treatments were performed in the past 20 years. Some of the fields were extremely faulted & compartmentalized and hence had considerable pore pressure variation over relatively small areas. Some data points correspond to fields under secondary recovery and in which injection fractures were the goal of the treatments. The pore pressure and closure stress variation with depth is depicted in Fig. 4 alongside the hydrostatic gradient (0.45 psi/ft) and the normal lithostatic gradient (1.1 psi/ft). Fields A, E and H are subhydrostatic. Fields B and G are overpressured. Fig. 4 is illustrative of the difference between the pore pressure and the closure stress. A very high difference is observed in the case of Field G, K and M. In the case of Fields B, C and G, the closure stress is very close to a value of 1.1 psi/ft. As can be seen in Fig. 4, the closure stress varies within the treatments from the same field.

Synthetic test cases
Using the range of parameters from the treatment analysis, four synthetic reservoir test cases were formulated. The reservoir is approximated as three layer systems with a payzone bounded by overlying and underlying shales. It may be mentioned beforehand that these synthetic cases may not be representative of both conventional reservoirs and unconventional reservoirs where much of the current fracturing activity in the world is going on. The rationale behind the selection of simplified  three layer systems is to exclude net pressure increase due to factors like composite layering, PKN type growth and focus on just the influence of effective stress. The first three synthetic cases represent under pressured, hydrostatic and over pressured resevoirs. The fourth synthetic case represents a high effective stress reservoir in which there is a large difference between the pore pressure and closure stress. The four synthetic cases are shown in Fig. 5 to highlight the relative difference in pore pressure and effective stress. The first three synthetic cases utilize a common reservoir geometry and geomechanical property descriptions. This was obtained using the statistical mean from the reservoir geometry for all treatments (see Appendix Fig. B). For the fourth high effective stress, high pressure synthetic case a different depth is used (Appendix: Table A and B). The fluid parameters and injection schedule are maintained constant (Appendix: Table C). The pore pressure and closure stresses are varied for the four synthetic cases so as to represent a variation in pore pressure and hence the effective stress (Appendix: Table D).

Frac3D coupled simulation
Frac3D is used for the fracture propagation simulations. Frac3D is a 3D planar hydraulic fracture model developed as an advanced model within the fracture modeling suite of Fracpro fracture design and analysis software. The model is based on a finite element solution for the fracture opening and a finite difference solution for slurry transport in the fracture. In previous work detailed in Hsu et al. (2012), the model had been benchmarked against data from physical laboratory tests done by Shokir and Al-Quraishi (2007). The Frac3D model was used by De Pater (2015) to match microseismic event clouds from the M-Site B-Sand injection. The MWX project was carried out in the Mesaverde sands of the Piceance Basin in Colorado (Warpinski et al., 1999). The M-Site fracture treatment was a test case for an example of contained fracture growth in an environment where the stress contrast was insufficient to explain the observed fracture containment.
The Frac3D model uses an effective stress propagation criterion rather than LEFM based stress intensity computations for fracture propagation. The crack opening is related to fluid pressure by deriving a relation between nodal fluid pressure on a specified region and the crack width of open nodes. This relation is obtained before the simulation run and allows fast computation of crack width for a given fluid pressure distribution. The crack width relation is only dependent on the geomechanical parameters and layering so that it can be re-used for different pumping schedules. Frac3D requires a potential crack area (PCA) which is just the  limits of the planar mesh. This is the maximum fracture height and length which can be input manually or the FracPro end of job fracture dimensions may be directly used. The initial solution to the fracture initiation problem yields a fracture width profile, initial fracture dimensions (half-length and height) which corresponds to the wellbore pressure at perforation depth as calculated by the wellbore model in FracPro. The flow of control between the FracPro lumped model and Frac3D is depicted in Fig. 6. A quarter symmetry is used for simulation and hence the control volume is half of a symmetrical bi-wing fracture which is again divided into symmetrical half width volumes.
In the general theory of poroelasticity in rock mechanics and which is used in many conventional geomechanical simulators, the Biot definition holds. However for rock tensile failure, the experimental work by Schmitt and Zoback. (1989), Bruno and Nakagawa (1991) and Visser (1998) indicate that the failure criterion is based on Terzaghi effective stress. The Frac3D code uses single component elasticity and in the propagation criterion effective stress is used. In future development with the Frac3D code, for problems of computation of backstresses which assume significance for multistage fracturing simulations, Biot poroelasticity would be made use of.
The Frac3D code uses a cohesive zone model. A simple linear traction separation relationship that corresponds to a rigid softening behavior characteristic of ductile rocks is used. The loading branch of the constitutive relationship is not considered, the implication being that the behavior of soft, ductile rocks is perhaps not adequately captured. In future work we would envisage to modify the Frac3D code so that it can accommodate bilinear or linear-exponential traction separation relationships. The importance of using a complete traction displacement relationship, its effect on cohesive zone size and its influence on net pressure was highlighted by Sarris and Papanastasiou (2011). Since our endeavor was to numerically investigate the effect of pore pressure keeping other variables constant, we deliberately chose a simple linear traction separation law. Fig. 7 depicts this simple linear constitutive relationship along with the cohesive zone where TS is the rock tensile strength and uc is the critical fracture width.

Treatment analysis
Plotting observed net pressure obtained by picking graphical ISIP versus the closure stress obtained from G-Function plots, we see that there is a strong positive correlation. G-Function is a dimensionless time function relating shut-in time after fracture creation to total pumping time and is used to linearize pressure behavior during normal leakoff from a fracture. Barree and Mukherjee (1996) developed idealized curve plots to aid interpretation and Barree et al. (2009) summarized methods for closure stress interpretation. Fig. 8 shows the treatment analysis results for the breakdown injections and Fig. 9 for the minifrac injections. The linear fit equations are as follows: Breakdown:y ¼ 0:44x þ 139:36Minifrac:y ¼ 0:43x þ 257:14 The presence of an intercept for observed net pressure implies a positive net pressure for an effective stress of zero, assuming a linear relation. In shallow, unconsolidated formations which have low values of effective stress the net pressure is observed to be high. Hence it is hazardous to assume a linear intercept and it is more likely that the slope of the correlation is steep for very small values of effective stress. There are not enough treatments for this range in our analysis to discern this clearly. It is also noted from these equations that the ratio of the net pressure to effective stress is less than 1. It is seen that the median net pressureeffective stress ratio is higher in the case of minifracs than for breakdown injections and there are more treatments for a ratio greater than 1.

Coupled simulations
Using the simple three layer geomechanical models, coupled simulations are run for the four synthetic reservoir test cases. The net pressure and fracture geometry results are depicted in Fig. 10 and Fig. 11. The discontinuous part of the curves indicate the fracture initiation pressure and geometry calculated by the pseudo 3D model. This is used as an initial input for the Frac3D simulator. From the numerical results, it is observed that a higher effective stress elevates net pressure. This is true for an effective stress increase owing to low pore pressure and also in the case when there is a large difference between pore pressure and closure stress. The fracture widths are also higher with higher effective stress. In the case of fracture dimensions, it is seen that for the lowest effective stress case i.e. the overpressured case, the fracture is longer or possesses a higher aspect ratio (fracture length to height). Increase in effective stress results in a more radial fracture. In all the simulations, the volume of fluid pumped was just enough to cover the payzone. This ensures that there is no effect of fracture containment on net pressure.

Comparison with fracture treatment data
By directly comparing the simulation results to the pressure data from the fracture treatments (see Fig. 12), it can be observed that the Frac3D model net pressures match closely with the ISIP net pressure pertaining to the breakdown and minifrac injections. However the mainfrac ISIP net pressure is considerably higher and the scatter is so high that there is no discernable correlation. Assuming a linear fit for the range of effective stress used for the four synthetic cases, it can be observed that the breakdown injection data fits more closely and the minifrac datapoints exhibits more scatter.
An obvious explanation is that the viscosity of the fluid used for minifracs are much higher in general than for breakdown injections which are mostly performed using water or brine. There are multiple reasons by which the scatter in the breakdown and minifrac data can be explained. Firstly, there is uncertainity in picking the ISIP net pressure and the closure stress values. Since graphical tangent methods are used there is an element of interpretation involved. The value of effective stress corresponding to each treatment is a function of the ISIP and the closure stress and errors can be accumulated. Secondly, the database of treatments involve those in which extreme depletion have yielded inaccurate measurements of bottomhole pressure since the wells are on vacuum after pumping stops. Some treatments used surface gauges instead of bottomhole gauges adding to inaccuracy in determination of the bottomhole pressure. A third reason could be complex interplay of other factors (see Fig. 1) like multiple fracture propagation, shear decoupling and tip plasticity that dominates and hence mask the effective stress effect. A fourth reason, could simply be the fact that the treatment volumes and pumping times were different for each treatment. Hence the ISIP net pressure picked was not the end of radial growth when the fracture interacts with the stress barrier but simply the time when injection was terminated. Sensitivity analysis on the hydrostatic case revealed that the tensile strength of the rock was an important parameter that influenced the net pressure. Though tensile strength is generally presumed to be low in rock masses, it can possibly have higher values with depth. An analysis of the slope of a linear fit line through the points obtained from the simulations indicate that the slope is the ratio of numerical net pressure to the sum of effective stress and the tensile strength.
There is hence a component of net pressure utilized to maintain fracture opening and one that maintains fluid flow within the fracture. The fracture volumes created in the simulations are small compared to viscosity dominated propped fracture treatments and the width difference is nearly four times when the overpressured and high effective stress cases are compared. This indicates that the numerical correlation of net pressure as a function of effective stress is a function of the material properties of the payzone defined and is strongly related to the input tensile strength and the critical cohesive zone width. It has been observed that experimental hydraulic fractures are toughness dominated and follow a different scaling law as compared to field scale fractures which are viscosity dominated (Detournay et al., 2007). This viscosity dominance is evident from the mainfrac data.

Conclusions
The results from the simulations and analysis of treatments suggest  From the strong empirical correlation between net pressure and effective stress, it is evident that the effective stress law for fracture propagation is valid over a wide range of effective stress. It is remarkable that such a correlation is observed for a wide range of reservoir permeabilities, elastic modulus and treatment fluids. The scatter in the treatment data or deviations from the net pressure effective stress relationship can be attributed to effects of fracture complexity, fracture tip plasticity, near wellbore tortuosity and complex interactions between fluid driven fractures and pre-existing fractures.
From the numerical synthetic simulations the net pressure is around 48 percent of the effective stress and this matches closely with a fitted empirical correlation from the database of treatments. This fraction depends on material properties like tensile strength and the size of the cohesive zone and we ascertained this from sensitivity analysis performed by varying parameters with respect to the hydrostatic synthetic test case. The relatively simple physics of pore pressure (and high effective stress) can be used to explain high net pressures observed in fracturing operations in depleted and high stress regions. The correlation obtained can be used for fracture treatment design. The local tip pore pressure can be used as a pressure matching parameter and should be considered before applying very large values of modulus or fracture toughness to match high net pressures, as is commonly done using many commercial fracture models. Accurate values of closure and pore pressures are hence useful as a predictive tool for fracture design and pressure diagnostics. Since measured fracture dimensions from the treatment data were not available, we could not make any direct comparisons between the containment effects of pore pressure from the treatments to that predicted by the simulations. However the simulations indicate that the magnitude of effective stress influences the fracture geometry. A  higher effective stress would result in a larger net pressure, giving a wider, shorter and more radial fracture. Lowering the effective stress magnitude would give a lower net pressure, and therefore a smaller width and a longer length for the same fracture height, which would increase the fracture aspect ratio. It has been observed that fractures in depleted reservoirs remain contained within payzones and rarely exhibit fracture height growth. It can be thus be extrapolated that in addition to the effect of higher closure stresses in bounding shales (higher stress contrast) the higher effective stress owing to depletion is also playing a major role for fracture containment.
Comparison between net pressure profiles as a function of height growth and fracture dimension evolution as a function of injection time for the pseudo 3D and Frac3D coupled simulations, indicate that the lumped models in general underestimate net pressure and overestimate fracture size. Such an outcome is due to the fact that in the pseudo 3D models pore pressure does not explicitly influence the propagation criterion but just the fluid leakoff calculation. Thus fracture models need to be based on effective stress propagation criteria with a cohesive zone description in order to represent the physics of pore pressure change. Since the current Frac3D model is very computationally intensive, the pore pressure effect can also be incorporated in an approximate way for pseudo 3D models.