Simulation of hydraulic fracture propagation in poroelastic rock with application to stress measurement techniques

Int. J. Rock Mech. Min. Sci. & Geomech. Abstr. Vol. 28, No. I, pp. 1-14,

1991

Printed in Great Britain.All rights reserved

0148-9062/91 $3.00+ 0.00 Copyright © 1991 PergamonPressplc

Simulation of Hydraulic Fracture Propagation in Poroelastic Rock with Application to Stress Measurement Techniques T. J. BOONEI" A. R. INGRAFFEA~ J.-C. ROEGIERS~

Hydraulic fracturing is a technique that is now employed relatively frequently for measuring in situ stresses at depth. The in situ stresses are typically interpreted from the pressure vs time logs. The interpretive methods are typically based on hypothetical models for the initiation, extension and closure of the fracture during the test. However, in porous media where coupled poroelastic effects may be important, the validity of these hypothetical models can be questioned. This paper presents the results from numerical simulations of the complete process of fracture initiation, propagation, closure and reopening. These analyses were performed using a specially developed numerical procedure that includes a nonlinear, Dugdale-Barenblatt fracture model, a poroelastic model for the rock mass, and full coupling between fluid flow in the fracture and the rock mass. The analyses provide insight into the influence of poroelastic effects on the techniques used for stress measurement and into the kinematics of fracture closure and reopening. The results should be useful in improving the hypothetical models employed in the interpretation of in situ stress state.

INTRODUCTION Hydraulic fracturing has become a recognized method for stress determination at depth. Numerous numerical models have been developed to simulate these processes [1], but they typically assume that the rock mass is an elastic medium. However, when fractures are created for the purpose of enhancing oil or gas production the rock, virtually by definition, is a porous and permeable medium. Commonly, fluid leak-off into the rock mass has been assumed to conform to relatively simple I-D models and coupling between deformation of the rock and fluid flow in the rock mass has been ignored. A distinctive aspect of the present work is that the rock mass is modelled as a fully-coupled poroelastic tEsso Resources Canada Limited, 3535 Research Road N.W., Calgary, Alberta, Canada T2L 2K8. ~Program of Computer Graphics, and School of Civil and Environmental Engineering, CorneU University, Ithaca, NY 14853-3501, U.S.A. §School of Petroleum and Geological Engineering, The University of Oklahoma, Energy Center, Suite T301, Norman, OK 73019, U.S.A.

medium. The numerical methods used herein permit the simulation of a hydraulically-driven, plane strain fracture with 2-D fluid flow in the rock mass. In addition, the fluid flow in the rock mass is coupled to the volumetric strains. Another distinctive feature is the use of a nonlinear, Dugdale-Barenblatt or equilibrium model for the fracturing process. Poroelastic effects have generally been ignored when modelling hydraulic fractures. The inherent assumption is that: (1) the time-scale of the problem is such that poroelastic effects have not had time to develop; or (2) that the limiting magnitude of these effects, given the time to develop, is small enough to be neglected. However, in most cases involving fractures in porous rock, poroelastic effects of significant magnitude can develop [2-4]. It is therefore important to establish the time scales over which these effects develop. It is reasonable to assume that a porous rock behaves as an undrained, elastic rock mass if the time T required to create a fracture is such that T<
2

BOONE

eta/.:

SIMULATION OF HYDRAULIC FRACTURE PROPAGATION

is relatively impermeable or the fracturing fluids and pore fluids are highly viscous, c tends to have a small value and fractures of the desired length can be generated without poroelastic effects developing. This is, generally, the case for classical hydraulic fracturing treatments in oilfield applications. However, as the selection of fracturing fluids becomes more diverse, this assumption may need to be reconsidered. A case where poroelastic effects can clearly be important arises during operations such as waterflooding where fluid injections occur over periods of years at pressures that are often sufficient to induce fractures [5], This paper focusses on a second case: stress testing where water is commonly used to generate micro-fractures which are typically only a few metres in length. A prime objective here is to illustrate how poroelastic effects may influence the determination of stresses from the field data. There have been considerable efforts directed towards practical development of techniques for in situ stress measurement [6]. However, most of these have been based on the tacit assumption that the fracture closes uniformly along its length after shut-in. In fact, there have been relatively few detailed simulations of these processes, and most stress measurement techniques are based on conceptual models. One example of a previous effort to model the details of a mini-frac test is work by Ishijima and Roegiers [7]. the purpose of the present study is to extend understanding developed in previous work by presenting a series of simulations of hydraulic fracturing tests, where it is assumed that: (1) the rock mass is a poroelastic material; and (2) fluid flow in the fracture, deformation of the rock mass and fluid flow in the rock are all fully-coupled. The numerical procedures used in these simulations are especially appropriate for this task since the fracture geometry is a natural product of the solution procedure. In contrast to previous simulations that have attempted to model fracture closure [8], no assumption is made or constraint imposed during the simulation with regard to fracture extension, closure, or whether the fracture closes first from the tip, at the borehole or uniformly along its length. The results have application to micro-frac stress tests, as well as mini-frac tests which have been used for determining fluid-loss parameters [9].1" The presentation of these results includes several images generated on a high-performance engineering workstation with the intent of providing the reader with insight into the transient development of the coupled phenomena. It is hoped that by understanding the processes more fully, through the use of computer graphics, information such as pressure logs can be interpreted more accurately. In the following sub-sections: (i) the theory of poroelasticity is reviewed; (ii) the numerical methods are tMini-frac tests require fracturing fluid volumes on the order of 10 m 3, whereas micro-frac tests require fluid volumes of one to two orders of magnitude less [22]. :~The relation between the intrinsic permeability k which has dimensions of length squared (Darcies), and the permeability coefficient is ~c = k/l~, where Iz is the viscosity of the saturating fluid.

briefly described; and (iii) two analytical solutions for fractures propagating in poroelastic materials are examined. These solutions are used to illustrate the extremes of slow and fast fracture propagation. In the next section, a series of simulations is presented that illustrates poroelastic effects related to fracture initiation, propagation, closure and reopening. These results are subsequently discussed in detail and the practical implications of the simulations are highlighted.

Theory of linear poroelasticity Five material constants are required to define a poroelastic system. Typically, these are: the shear modulus G, the drained Poisson's ratio v, the undrained Poisson's ratio v,, Skempton's pore pressure coefficient B, and the permeability (or mobility) coefficient r21 These parameters can be related to a set of micromechanical ones which include the porosity ~b, the bulk modulus of the fluid Kr, the bulk modulus of the solid grains Ks and the bulk modulus of the porous, solid skeleton K; as well as v and k [10, 11]. The governing equations have been recently presented succinctly by Detournay and Cheng [12]; hence, their notation has been adopted herein. Three,derived poroelastic constants figure commonly in both analytical solutions and numerical approximations. Biot's coefficient of effective stress ~t, which is independent of the fluid properties, is defined as: 3(v. - v)

K

= B(1 - 2v)(1 + v.) = 1 - ~.

(1)

The common range for at in rock is about 0.5-1 [13], where the upper bound represents the typical assumption made in soil mechanics. A related parameter r/is defined as:

3(v~- v) =~t(1 - 2 v ) ~/ = 2B(l - v ) ( l + vu) 2(1 - v)

(2)

The range for r/ is 0-0.5, and like ~t, is independent of the fluid properties. It is often appropriate to introduce a nondimensional time scale ~, defined as ct/L2; where the diffusivity coefficient c is given by: e=

2GB2r(I - v)(l + v,) 2 9(v. - v)(l - v.)

(3)

and L is a characteristic length. Finally, it should be noted that, by definition, the total stress is related to the effective stress (assuming tension is positive), as follows: a~ = a ~ - p ,

(4)

where p is the pore pressure and ~ is the Terzaghi effective stress and is assumed to govern failure of the material.

Numerical approximations The numerical procedure and approximations employed herein have been described in detail by Boone and Ingraffea [14] and can be summarized as follows:

BOONE et al.:

SIMULATION OF HYDRAULIC FRACTURE PROPAGATION

1. The finite element method is used to model the domain of the rock mass. The formulation incorporates all aspects of the fully-coupled theory of linear poroelasticity. The isoparametric, plane strain elements used in the simulations have quadratic displacement fields but linear approximations to the pore pressure. 2. The fracture is modelled by zero-thickness interface elements that are placed along a predetermined fracture path. These dements allow for a nonlinear, equilibrium fracture mechanics model to be incorporated. 3. Fluid flow along the fracture is modelled using a 1-D finite difference scheme that assumes laminar flow and the cubic law to be valid. 4. Fluid leak-off through the fracture faces is fullycoupled with 2-D flow in the rock mass. Both leak-off rates and fluid pressure along the fracture are matched as boundary conditions for the finite element mesh. 5. It is assumed that the fluid entering the fracture has the same viscosity as the fluid that saturates the rock mass.

Analytical solutions for fractures propagating in poroelastic materials There are two analytical solutions that provide insight into the nature of the pore pressure fields in the vicinity of hydraulic fractures. The solutions represent the extremes of very slow and very fast propagation. First, Ruina [15] has found a solution for a semi-infinite fracture moving at a constant rate, loaded as shown in Fig. l a by an applied stress P over a length L along the crack faces which are assumed to be impermeable. The velocity of the fracture v is subject to the condition Y

(a)

l =ry=p dp/dy=O

))))))

I

= v

:-X

(()))) k_ I--

L

_1 -7

vL/c >>1. The solution for the pore pressure field in the near-tip region is given as [15]:

P =~'~-"-K e 2B(13+ v " ) { l

Y I O'y=P .p=P

))))))

=--X

)))))¢ t__ '-

v"O

----

d L

Fig. 1. (a) The dimensions and boundary conditions for Ruina's analytical solution [15]; and (b) the dimensions and boundary conditions for an analytical solution of a static fracture.

-e-"~l+~/(2c)}c°s

(o)

2 ,

(5)

where r is the radial distance from the crack tip (normalized with respect to L), 0 is the angle measured about the crack tip from the X-axis and K o is the apparent fracture toughness of the rock. The crack tip is at the origin and the reference frame is assumed to translate with the fracture. Figure 2a shows a 3-D relief plot of Ruina's solution for the pore pressure in the close vicinity of the crack tip. Several important features are illustrated by this solution. First, along the crack faces the pore pressure is equal to the ambient pore pressure, which in this case is zero. This is indicative of the fact that there is no volumetric strain along the crack faces. Second, the gradient of the pore pressure field normal to the crack faces is zero everywhere along the crack faces. This complies with the no-flow boundary condition. The gradients of pore pressure field when plotted in relief are proportional to the flow rate since the theory of poroelasticity incorporates Darcy's model for fluid flow. Third, ahead of the crack tip is a negative pore pressure which is indicative of positive volumetric strain and reduced effective (tensile) stress in that region relative to the initial state. However, at the crack tip itself, the pore pressure equals the ambient value and the material is effectively drained. This effect is due to the fact that the volumetric strain gradient approaches infinity near the crack tip. Thus, if the material has a finite permeability and the crack is propagating, then the tip must be drained since the fluid would be subject to an infinite volumetric strain rate. It is apparent that, even for the case of impermeable crack faces, the flow pattern is away from the region near the crack faces and into the region ahead of the crack tip. Figure 2b is a relief plot of the pore pressure solution for a static fracture with the boundary conditions shown in Fig. lb. In this case, the crack walls are assumed to be permeable and steady state conditions have been reached (i.e. vL/c = 0). The solution for the pore pressure field is well known from potential theory and is given by: P

(b)

3

f 2r

---tan - ) ~ r - - - ~ c o s ~ ) ~ ,

0~0

~<~,

(6)

where r is the radial distance from the crack tip, normalized with respect to the crack length, 0 is the angle measured with respect to the crack path and P is the pressure applied at the crack face. Under these conditions, the pore pressure field is positive everywhere. The stress singularity at the crack tip is of the standard form where the stress intensity K~ is given by [16]:

K,=P

{l-r/}.

(7)

BOONE et al.:

4

SIMULATION OF HYDRAULIC FRACTURE PROPAGATION

K~ is found to be positive for all material properties. The upper limit to equation (7), corresponding to ~ --0, is the solution for the case with impermeable crack faces. The asymptotic limit for very fast fracture propagation is identical for both permeable and impermeable crack faces [15]. Therefore, the static solution and Ruina's solution represent the extremes of very slow and very fast propagation of hydraulically driven fractures, respectively. In practical situations, at finite crack speeds, the pore pressure field will lie between these extremes. A S I M U L A T I O N O F F R A C T U R E IN A PERMEABLE MATERIAL

In this section the results are presented for a detailed simulation of fracture propagation initiated from a borehole in a poroelastic material. The finite element mesh used in this simulation is shown in Fig. 3, and the material parameters are those listed in Table 1. The assumed stress state was al = 2.5 and ¢r2 = 2 MPa. The initial pore pressure P0, along with the tensile strength ¢rT, and the fracture energy G~¢, were assumed equal to zero. Two analyses were conducted. In the first it

was assumed that, on closing, the fracture "healed" completely. In the second it was assumed that the fracture was partially propped open, either by an added prappant or by existing asperities on the fractured surface. Plotted in Figs 4a and b are the crack-mouth pressure (CMP), crack-mouth opening (CMOD), crack length and the fracture volume as a fraction of time, for the two cases. During the analyses, the boundary conditions at the borehole were varied as follows: (i) from 0 to 0.23 sec the CMP was increased at a rate of 20 MPa/sec; (it) from 0.23 to 1 sec the flow rate at the crack-mouth was held approximately constant at 2 x 10 -4 m~/m. sec; and (iii) from 1 to 4 sec the flow was shut-in, so that the flow rate at the crack-mouth was approximately zero. The pressure along the borehole wall was imposed to be the same as that at the crack-mouth. The simulation for the fracture that was allowed to heal was terminated after 5 sec. For the second case, the flow rate at the crack-mouth was shut-in until 4 sec, at which point the CMP was again increased at a rate of 4 MPa/sec. Note that the two analyses differ only after shut-in once the fracture begins to close.

,Y

,'o E

,=

o,:

I I

1,o

:

/1\ /I ",-.

12

3t- I II . r i_ - r ~

SHUT-IN r

N•

12m

-

i

"G',,. s-

CMP

', " ' £ . . . . . . . .

2~

I /-CMOO

0 ~

(o)

I

2

3

4

5

TIME (s)

X Io

'

7

i

s

/ I

(b) SHUT-IN

CLOSURE

REOPENING

2.Om

Fig. 3. Dimensions for the finite element mesh used in the simulations of fracture propagation. The isoparametric elements have quadratic displacement fields and linear pore pressure field. Table I. Material properties for the poroelastic simulations

-

i

'14 tw5o o =,

I:

~E

..reX>

Shear modulus Poisson's ratio, drained Poiuon's ratio, undrained Skempton's coefficient Permeability Tensile strength Fracture energy Fluid viscosity (rock and fracture)

G v v, B u aT Glc #

6000 MPa 0.2 0.33 0.62 2 x 10-~ mi/MPa .sec 0 MPa 0 MPa" m 1 x 10-9MPa.sec

ilK.

3/I

zL!

"--~

;

:/

',"-.

/',

LENGT.: """--:=/ --

I/~"

, /

',

! ,-CMO0

I 0 I~ 0

',

FCMP

/'~

1

l VOLUME-" --~"t . . i. . . . i

i I

2

3

2' "~

i I 4

i 5

TIME (s)

Fig. 4. Crack length, crack volume, C M O D and C M P vs time for the simulation of a fracture in a poroelastic medium: (a) without proppant; and Co) with proppant.

Fig. 2. (a) A 3-D plot of Ruina’s analytical solution for a fracture propagating at a constant rate in a poroelastic medium; (b) a 3-D plot of an analytical solution for a static fracture in a poroelastic medium. In both cases, pore pressure is plotted in relief and highlighted by the colour contours, as well.

Fig. 5-Caption

overleaf.

Fig. 5. A series of 3-D plots showing results at particular times for the same numerical simulations that are plotted in Fig. 4. The scale of the colour contours is shown on the right-hand-side in units of MPa. Information about the displacement magnification and analysis time is located in the top left comer of each photo. The finite element mesh is located in the Z = 0 plane of the 3-D display. XY dimensions can be related to the mesh shown in Fig. 3. The scale in relief can be related to known boundary conditions such as the pressure at the borehole which can be determined from Fig. 4. Pore pressure is plotted in relief, and effective stress (u;) displayed by colour contours in parts (a-g). Image (a) is taken at 0.22 set of simulated time just prior to breakdown; (b) is at 0.28 set just subsequent to breakdown; (c) is at 1set just prior to shut-in; (d) is at 1.61 set, for the fracture with proppant, at the point when the fracture length is stationary; (e) and (f) comps.re “overhead views” of the fracture just prior to shut-in (1 set) and at the stationary point (1.61 set), for the case without proppant; and (g) is at 4.71 set, for the case without proppant, illustrating fracture closure. The display parameters are reversed in parts (h), (i) and (j), which illustrate results for the propped case. Effective stress (CT;)is plotted in relief and pore pressure is displayed by colour contours. Image (h) is taken at 2.55 set of simulated time as the fracture closes on proppant; (i) is also at 2.55 set, however, a scaled grid has been introduced to allow more quantitative estimation of the stresses. Finally, (j) shows the fracture reopening from the borehole.

BOONE eta/.: SIMULATIONOF HYDRAULIC FRACTURE PROPAGATION This series of events was designed to illustrate-a number of effects associated with techniques used to determine the state of stress at depth. Several important observations will now be noted and illustrated with 3-D plots. First, it is observed that the breakdown pressure is approx. 4.6 MPa, although the fracture was observed to initiate at a lower pressure. Figure 5a shows pore pressure in relief and effective stress 0"~ using colour contours. The surface plot is shown relative to the deformed mesh so that the extent of fracture propagation can be observed. The figure is taken at 0.22 sec, prior to the pressure reaching the breakdown pressure and after the fracture has initiated. Figure 5b is taken shortly after breakdown, at 0.28 sec. It is now observed that the flow into the fracture is significantly altering the pore pressure field, and there is no longer a steep gradient in pore pressure at the fracture mouth. Figure 5c corresponds to 1 sec into the simulation. Ahead of the crack tip, the pore pressure field is very similar in form to Ruina's analytical solution which is shown in Fig. 2a. Along the crack faces the effective stress 0"~is zero which is the correct boundary condition. At this point, flow into the fracture was shut-in, but the fracture continues to propagate at a progressively declining rate. At approx. 1.6 sec, as shown in Fig. 5d, the crack length is stationary. This figure shows notable similarities with the plot of the analytical solution for a static fracture in Fig. 2b. The rate of fracture propagation in the simulation has slowed to the point where the pore pressure has had sufficient time to diffuse ahead of the propagating tip. Figures 5e and f show the contrast between the crack-opening profile just after shut-in and at 1.6 sec when the fracture length is stationary. As the fracture begins to recede, the two cases diverge. It was assumed for the propped case that any point along the fracture that had opened more than 0.05 mm could only close to that width. It was also assumed that there was no appreciable flow resistance induced by the asperities or proppant, other than that due to the width of the channel. In Fig. 5g fracture closure for the unpropped case is seen to be progressing back towards the borehole. In an extensive region around the fracture, both the stress state and pore pressure are significantly altered. Clearly, a 1-D leak-off model would have been inappropriate at this stage. Figure 5h shows a plot at 2.5 sec for the propped case where the fracture has partially closed on the proppant. Effective stress 0"~ is plotted in relief and colour contours are used for the pore pressure. It can be seen from the relief plot that along the section of the fracture that has closed on proppant, the effective stress is negative. This is more clearly illustrated in Fig. 5i where a scaled grid has been introduced along the fracture path. The fracture finally closes completely on the proppant at a time of 2.7 sec. As might be expected the fracture closes more quickly for the propped case. The CMP at closure is about 2.5 MPa, which significantly exceeds the value of 0"2 (2 MPa). It was observed that after closure there was a sudden increase in the rate of pressure decline, and it appears that the pressure is then asymp-

9

totically approaching P0 (0 MPa), as shown in Fig. 4b. After closure, it was noted that there was a peak in the magnitude of 0"~.at the ends of the propped section and that stress decreased towards the centre of the fracture, with a second stress concentration at the borehole. The CMP was then increased in order to determine at what pressure the fracture would reopen. Figure 5j plots effective stress in relief and pore pressure using colour contours, at a simulated time of 4.3 sec, just after the fracture reopened. It was observed that the fracture reopened at the well bore and again at a CMP of about 2.5 MPa. INTERPRETATION OF PRESSURE VS TIME LOGS Breakdown pressure

The breakdown pressure Pb is normally defined as the maximum pressure observed during a hydraulic microfrac test. Prior to reaching breakdown, the pressure is ramped or progressively increased until there is a sudden increase in the flow rate into the borehole which identifies fracture opening. Typically, the pumping rate is limited so the fluid pressure then declines as the fracture propagates. It is commonly assumed that this peak corresponds to the instance of fracture initiation which can in turn be determined using a strength-of-materials approach [17]: Pb ----30"2-- 0"t -t- 0"T,

(8)

where aT is the tensile strength of the rock. Given an estimate of a2 and 0.r, 0"1 can be calculated from the measurement of Pb. However, in practice, as discussed in the following paragraph, equation (8) is often an unreliable method for estimating 0"1. For the simulations reported herein, 0"T--0, 0"1--2.5 and 0"2--2. Therefore, based on equation (8), a breakdown pressure of 3.5 MPa would be predicted. In contrast, Fig. 4a shows a simulated breakdown pressure of 4.5 MPa. Potentially, these simulations can provide an understanding of the shortcomings of a strength-of-materials approach for predicting Pb. Zoback et al. [18] has conducted a series of experiments that provide insight into the nature of breakdown at a borehole and the nature of this discrepancy: (1) breakdown pressures can be reduced by poroclastic effects; (2) breakdown is typically preceded by detectable acoustic emissions which proves to be a better indicator that the tensile strength has been reached; and (3) the breakdown pressure is clearly dependent on the rate of pressurization. The first and third points had also been shown previously by Haimson and Fairhurst [17]. It has been shown numerically for an impermeable material that the breakdown pressure can be influenced by the development of a nonlinear fracture process zone [19]. It was concluded that the initiation of the fracture process at a borehole can be predicted using equation (8) and that the pressure required for breakdown Pb may be measurably greater than the pressure at the instant of fracture initiation Pi. If it is reasonable to assume that fracture initiation would be accompanied by acoustic emissions, it can be concluded from these numerical

10

BOONE et al.: SIMULATION OF HYDRAULIC FRACTURE PROPAGATION

results and experimental observations that at breakdown there may have already been significant fracture development causing Pb to be greater than P~. This effect can be clearly associated with the development of a fracture process zone. In the analyses described herein, the fracture toughness and tensile strength are assumed to be zero which precludes the existence of a fracture process zone. Therefore, there must be other factors which contribute to the observed increase in Pb relative to that predicted by equation (8). The simulations presented herein show that prior to breakdown there may be a significant pressure gradient along the length of the developing fracture as is illustrated in Fig. 5a. Linear elastic fracture mechanics concepts can be used to show that for a crack to remain open when there is a significant gradient in the fluid pressure along the fracture, the pressure at the crackmouth may significantly exceed Pb, as predicted by equation (8). Since the pressure distribution along the fracture is a function of parameters such as the crackopening profile, fluid viscosity and the pressurization rate, it follows that Pb is related to these parameters as well. Figure 5b shows that the pressure gradient at the crack-mouth decreases rapidly after breakdown. In poroelastic materials, assuming a strength-ofmaterials approach and an effective stress failure criterion, P~ is found to fall within the following limits [3]: 30"2 -- 0.1 + 0.T -- 2~/P0

2(1 -,I)

~< Pi ~< 30.2 - a~ + 0.T-- P0, P0 ~< 0.2 ,,< 0.~, (9)

where P0 represents the initial pore pressure. When a steep pore pressure gradient prevails near the borehole, as seen in Fig. 5a, P~ should approach the upper bound in equation (9). The lower bound would be approached in cases where pore pressure approaches a uniform value in the vicinity of the borehole [12]. Results have been presented in Ingraffea and Boone [20] for similar poroelastic simulations which specifically investigated the experiments made by Zoback et al. [18]. The time scale on these simulations suggested that the lower bound to equation (9) should be approached. As predicted, poroelastic effects were found therein to lower PiThe upper limit for P~ given in equation (9) does not account for the observed increase in Pb resulting from higher rates of pressurization. This effect can be associated with two physical phenomena. The first is related to the previous definition of breakdown and its relation to the pressure gradient along a crack. Since the pressure profile in a fracture is a function of the crack-opening and fluid viscosity, the fluid's viscous flow resistance causes a rate dependent effect which implies that, as the rate of pressurization increases, Pb will increase. Furthermore, it must also be recognized that fracture initiation is a process whereby numerous micro-cracks are initiated over a relatively broad region. Subsequently, the strain localizes within this region and forms a single discrete fracture, as discussed by Boone et al. [19]. The point at

z

P'Po, =r~"o ' ~ - ' - ~ k

~f-

~"'~

r ~'~

te

ct

g

:= I

2

3

RADIUS (r/a)

Fig. 6. A plot taken from Detournay and Cheng [12] illustrating the transient development of o0e for a pressurized borehole (o, = 0, o00= P0 at the borchole wall). which fluid can flow into a fracture, or perhaps one of several micro-cracks, with a minimal pressure gradient at the fracture mouth is therefore a function of the fluid properties, the time one allows for the fluid to penetrate the fracture, and the process of strain localization. The second phenomenon is related to the diffusion of pore pressure from the borehole. Detournay and Cheng [12] have published an analytical, plane strain, solution for arbitrary loading of a borehole in an infinite poroelastic medium. They found that the loading due to fluid pressure only, 0., = 0, p = P0, at the borehole wall, results in transient development of the circumferential stress (0.00) as shown in Fig. 6. At early times, 0.0e varies from tension at the borehole wall to compression a short distance away from it. It has been observed in several simulations that this advancing compressive zone limits the rate of fracture propagation [20]. The fracture tip progressed at approximately the same rate as the boundary of the tensile and compressive zones, provided fluid could not enter freely into the fracture. Since the crack-mouth opening generally increases with the crack length this effect may influence the time at which fluid can enter freely into the fracture. At loading rates which are fast relative to the diffusion of the pore pressure this effect can be expected to cause an increase in Pb which is consistent with experimental observations [! 7, 18]. Shut-in pressure and fracture closure

Commonly, the instantaneous shut-in pressure Pi=p, is identified from a slope change on the pressure log after pumping has stopped. For the example in Fig. 4b, there is a distinct change of slope in the C M P curve at 2.7 sec from which the P~sipis determined to be 2.5 MPa. In fact, this is the point at which the fracture was observed to close completely on the proppant and is, therefore, equivalent to the so-called closure pressure Pf=. This observation leads to the conclusion that the physical process of closure can be clearly identified and that it may be reasonable to equate the P~,~p to the Pr=. However, in this case, the P~sapis found to be significantly greater than 0.2, as has been recognized in previous experimental work [21]. The challenge is therefore to determine the relationship between the P~,~pand 0.2.

B O O N E et al.:

SIMULATION

OF

~i~ULiC

IO A

=0--

8

i<

E

CLOSURE ( I )

CMP ( I )

u

"1

== 4 o =E o.

.

.

CMP 121 . . . .

"

CMOD (2

\

FRACTURE

I1

PROPAGATION

increase in the Pfoc is an upper bound. While the effect of the proppant may not be negligible, it cannot account for the observed difference between the Pfo~ and S2. Ar~upper bound to increase in the Pr~ due to poroelastic effects can be readily derived. First, it is known [26] that the long-term limit for the apparent change in confining stress A0.c, induced along a fracture due to pressure in the fracture P r - P 0 , is:

CLOSURE121 J . . . .

4-- I

P0).

AO"c ---- r / ( P f -

(10)

2

If Pf is the fluid pressure in the crack prior to closure and

O

that pressure is held constant for a relatively long period of time T prior to closure such that T>>L2/c, where L is the length of the crack at closure, then it follows that the pressure at closure Pf~, is:

-2

~

.

=

=

•

b

-I O I 2 3 TIME AFTER SHUT-IN (MIN.)

FiB. 7. A plot taken from M e d l i n and Masse [22] illustrating the C M P

and C M O D following shut-in for a small-scale laboratory experiment fractured with grease. Results for two experiments are illustrated: (I) where $2 = 55 MPa; and (2) where $ 2 = 4 1 MPa.

The pressure log shown in Fig. 4b is similar in form to the signature obtained for small-scale experiments by Medlin and Masse [22, Fig. 17], and for propped, minifrac tests by Smith [23, Fig. 1 I] and Nolte [24, Fig. 1). McLennan and Roegiers [21] have observed similar features in a series of laboratory experiments. Figure 7 is the figure from Medlin and Masse. It includes plots of both the CMP and CMOD vs time for two experiments. These authors did not introduce proppant in the experiments; however, the CMOD was observed to be effectively propped on closure, supposedly due to creep. Regardless of the mechanism for residual crack-opening, it is apparent that closure at the crack-mouth can be correlated with a characteristic increase in the slope of the pressure-time curve and that the CMP at this point exceeds the minimum in situ stress of 2 MPa. Plateaus in the pressure log, such as that between 1.5 and 2.7 sec, and distinct increases in the rate of pressure decline such as that at 2.7 sec in Fig. 4b, are not always observed in the micro-frac tests used for stress measurement [6]. This is most probably due to the fact that no wellbore storage capacity is included in the simulation. In practice the zero flow condition is commonly applied at the surface, so that the fluid and piping system from the surface to the crack-mouth also acts as reservoir of fluid, and there may be significant flow in or out of the fracture. Two factors which may contribute to the difference between the Pf= and $2 are: (1) stresses induced by the displacements due to the proppant, and (2) poroelasticallyinduced confining stresses. For these simulations, the elastic solution for a Griffith crack has been used to estimate that the effect of the proppant is to increase the Pro~ by about 0.06 MPa. Timoshenko [25] also provides the analytical solution for the stresses induced by a rigid circular punch. In this analogous case, the stress normal to the surface is found to be a minimum at the centre of the punch. Therefore, based on the observation that the fracture tip recedes back towards the borehole as it closes, it would seem that the estimate of a 0.06 MPa

Pfoc -'- 0"2 -I- A0. c ,

(11)

The upper limit to Pfoc c a n be determined by equating P f to Pfo~ which gives: 0.2 ~P0 1 -- r/ -

Pfoc.upperbouad - -

(12)

In practice, since closure is a transient event the poroelastic effects cannot fully develop, so that equation (12) represents an upper bound on Pfoc- For the example problem (02 -- 2 MPa, P0 = 0 and r/ffi 0.33), the upper bound is estimated to be 3 MPa, which appears reasonable. Potentially, situations could be contrived where this upper bound was exceeded, however, for practical purposes it should be acceptable. This simple estimate of an upper bound, in combination with a simulated shut-in process, leads to the conclusion that the difference between Pisip and $2 is predominantly a transient poroelastic effect. Smith came to the same conclusion based on the observation that the closure pressure increases with successive testing of the same fracture [23]. He also found the changes in closure stress to be comparable in magnitude to the difference between P~p and $2 in the simulation presented herein. It should be noted that Settari and Raisbeck [27] have also arrived at a result similar to equation (12) and that it can be ascertained from analogous work in thermoelasticity that this upper bound is applicable to the 3-D problem of a planar elliptical fracture [28]. Both the 2-D and 3-D solutions are limited to infinite domains of a single poroelastic material. For the case of a vertical fracture which is confined to a poroelastic layer between two impermeable horizontal layers as is assumed in the PKN fracture model, the upper bound comparable to that in equation (12) is [29]: Pfoc, upperbouad,PgN ----

S2 - 2r/P0 1 - 2r/

(13)

This upper bound can significantly exceed that in equation (12) and indicates that geometric effects and bounding conditions are an important consideration.

Closure stress along a propped fracture For the closed fracture, as seen in Fig. 5i, there is a sharp increase in the stresses at the end of the propped

BOONE et al.: SIMULATION OF HYDRAULIC FRACTURE PROPAGATION

12

section. These stress concentrations are a characteristic observed in problems with an uniform applied displacement applied over a specified area. For example, one can refer to the linear elastic solution for a rigid punch on a half space which was referenced previously [25]. Poroelastic effects can be expected to have a similar and compounding influence on the stress distribution. Cleary [30] has shown that for a stationary fracture the apparent increase in compressive stress resulting from leak-off and poroelastic effects, which he termed backstress, is maximum at the crack tips and minimum in the middle of the crack. It follows that if the stresses acting along a closed fracture are greatest near the crack tip, as in Fig. 5i, then during the process of closure the crack tip is likely to progressively recede back towards the borehole. This provides a simple explanation of the closure mechanisms observed in the simulation.

Reopening pressures The fracture reopening pressure has also been suggested as a measure of a2 and often its value has been observed to be approximately equal to the estimated pressure at closure [31,32]. Enever has suggested that this value is most reliably determined by the use of a flowmeter located within the packed-off section. In practice, the fluid pressure in the borehole is incremented until there is a distinct increase in the flow rate which identifies the fracture reopening pressure. For the case of a fracture which has closed on proppant, there is, however, a residual stress concentration at the borehole wall. This leads one to question whether the fracture would first reopen at the borehole or away from the borehole. It was observed during the simulation that as the borehole and propped fracture were re-pressurized, the stress concentration at the borehole decreased as the pressure increased and the fracture first reopened at the borehole. Reopening was associated with a distinct increase in the flow rate. The pressure at this point was approximately equal to the pressure at closure, 2.5 MPa, which is consistent with field observations. However, the reopening pressure, like the Pisip, did not equate to a2. The difference between a: and the reopening pressure can again be clearly identified as a poroelastic effect as discussed in the following section. SIMULATIONS OF FRACTURE REOPENING A final series of simulations investigated the influence o f poroelastic effects on fracture reopening. It was assumed that the fracture was propped along a known length, such that there is negligible flow resistance in the propped channel (i.e. pressure along the channel is constant). The initial conditions were P0 = 0, a~ = 2.5 MPa and a2 = 2 MPa. T h e pressure along the fracture was then increased at a constant rate, until the fracture reopened. This procedure was applied in a series of tests that investigated the effect of varying the rate of pressurization. The results are plotted in non-dimensional form

1.0 fir"

//

RFopIr NeNG CMP

o,

o u

/

i / , 1 o.oo

.

Wok,~.-~,;

z W 0.5

~

/

rlMtr

/

OZ ,,2.0 MPo

/

I~

//

8 N O Z

/ / 0

I -4

..4..

0

i

I -2

LOGARITHM

~

~'I'2.SMPo

~

""

~"

2L

I

I 0

I

OF NORMALIZED

I 2 TIME,

I

I 4

('% c T Lz

Fig. 8. Normalized plot of the fracture reopening factor f,. as a function of the pressurization rate t* = cT/L 2 where L is the propped length and T is the time from initial pressurization to fracture reopening. The pressurization rate is assumed to be constant and at the beginning of pressurization the pore pressure is constant, P0, throughout the field. in Fig. 8 with the reopening pressure P,, defined as follows:

P,=S2+f,(z)~

/

ze°

l-r/

}

$2 ,

"

r=~-~,

(14)

where L is the half length o f the propped fracture, t, is the length of time from commencement of pressurization to fracture reopening, andfr is the influence function for the poroelastic contribution to P , which can vary from 0 to 1. In order to obtain a more rigorous estimate of a,, reopening tests using various rates o f pressurization should be performed at each horizon of interest. Indeed, if one can produce reasonable estimates or measurements of c, r/and L, then it is possible that equation (14) could be used in conjunction with Fig. 8 to improve the estimate o f a2 that one would interpret from a reopening test.

DISCUSSION AND CONCLUSIONS The intent of this paper has been to illustrate clearly that poroclastic mechanisms can have a significant effect on hydraulically-driven fracture propagation. It has been shown that a fully-coupled, 2-D simulation o f fracture propagation, closure and reopening exhibits many significant events than can be identified in pressure logs from both laboratory and field tests. On the basis of these comparisons, the following conclusions can be drawn: I. Poroclastic mechanisms contribute to increases in the breakdown pressure for a borehole as a function of the rate of pressurization. 2. Fracture closure can be identified from a distinct increase in the rate o f pressure decline on a pressure log after shut-in. 3. Poroelastic effects can cause the pressure at closure to be significantly greater than the minimum in situ stress $2.

BOONE et al.: SIMULATION OF H~RA~L~LiC FRACTURE PROPAGATION 4. A practical u p p e r b o u n d to the shut-in pressure-has been derived.

13

infinite impermeable elastic layers. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr. 27, 189-197 (1990). 5. Koning E. J. L. Waterflooding under fracturing conditions. Ph.D. Dissertation, Delft Technical University, Netherlands (1988).

5. The reopening pressure for a p r o p p e d fracture was f o u n d to be approximately equivalent to the shut-in pressure and also influenced by poroelastic mechanisms.

6. Pr~. Second Int. Workshop on Hydraulic Fracturing Stress Measurements, Vol. 1 (B. C. Haimson, J.-C. Rocgiers and M. D.

Based on the numerical simulations it is also reasonable to conclude that poroelastic effects influence the kinematic process o f fracture closure and reopening. The fracture was observed to close progressively f r o m the tip back towards the borehole. The fracture reopened f r o m the borehole progressively along the length o f the propped section. These mechanisms can be attributed in part to stresses induced t h r o u g h poroelastic coupling which act to close the fracture. The stresses are m i n i m u m at the borehole and m a x i m u m at the ends o f the p r o p p e d section. It is a p p a r e n t that a better understanding o f these processes could lead to i m p r o v e d estimates o f the in situ stress state. The dependence o f the reopening pressure for a p r o p p e d fracture on the pressurization rate has been d e m o n s t r a t e d numerically. The results, which are plotted in Fig. 8, exhibit a form similar to plots developed by others [e.g. Ref. 3, Fig. 5] for other related poroelastic effects. Specifically, the time-scale over which poroelastic effects develop are similar. In general, poroelastic effects can be o f potential significance if the treatment time exceeds 0.01 .L2/c. Plots such as Fig. 8 can be used to determine if poroelastic effects are potentially significant and they can be used to estimate appropriate corrections for estimates o f the in situ stress state. Finally, hydraulic stress measurement techniques are largely based on conceptual models for the behaviour o f fractures during closure. The a u t h o r s have f o u n d that applying visualization techniques to problems o f this nature can provide an extremely valuable insight into the complex, coupled, transient physical processes.

8. Settari A. and Cieary M. P. Three-dimensional simulation of hydraulic fracturing. JPT 36, 1177-1190 (1984). 9. Nolte K. G. Determination of fracture parameters from fracturing pressure decline. Paper SPE 8341, SPE-AIME 54th Annual Fall Technical Conf. and Exhibition, l.,as Vegas (1979). i0. Blot M. A. General solutions of the equations of elasticity and consolidation for a porous material. J. Appl. Mech. 23, 91-96 (1956). 1I. Nut A. and Byerlec J. D. An exact effective stress law for elastic deformation of rock with fluids. J. Geophys. Res. 76, 6414-6419 (1971). 12. Detournay E. and Chang A. H.-D. Poroelastic response of a borehole in a non-hydrostatic stress field. Int. J. Rock Mech. Mm. Sci. & Geomech. Abstr. 25, 171-182 (1988). 13. Ri~ J. R. and Cleary M. P. Some basic stress diffusion solutions for.fluid-saturated elastic porous media with compressible constituents. Rev. Geophys. Space Phys. 14, 227-241 (1976). 14. Boone T. J. and Ingraffea A. R. A numerical procedure for simulation of hydraulically-driven fracture propagation in poroelastic media. Int. J. Numer. Analyt. Meth. Geomech. 14, 27-47 (1990). 15. Ruina A. Influence of coupled deformation-diffusion effects on the retardation of hydraulic fracture. Proc. 19th U.S. Rock Mechanics Syrup., pp. 274-282 (1978), also M.S. Thesis, Division of Engineering, Brown University (1978). 16. Huang N. C. and Russdl S. G. Hydraulic fracturing of a saturated-porous medium--lI. Special cases. Theoret. Appl. Fract. Mech. 4, 215-222 (1985). 17. Haimson B. and Fairhurst C. Hydraulic fracturing in porous permeable materials. JPT 21, 811-817 (1969). 18. Zoback M. D., Rummel F., Jung R. and Raleigh C. B. Laboratory hydraulic fracturing experiments in intact and pre-fractured rock. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr. 14, 49-58 (1977). 19. Boone T. J., Wawryznek P. A. and lngraffea A. R. Simulation of the fracture processes in rock with application to hydrofracture. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr. 23, 255-265 (1986). 20. Ingraffea A. R. and Boone T. J. Simulation of hydraulic fracture in poroelastic rock. Nwnericai Methods in Gemnectamics (lnnsbruck 1988) (Swoboda, Ed.), pp. 95-105. Balkema, Rotterdam 0988). 21. McLennan J. D. and Roegiers J.-C. How instantaneous are instantaneous shut-in pressures? Paper SPE 11064. 1982 SPE Annual Technical Conf. and Exhibition, New Orleans (1982). 22. Medlin W. L. and Masse L. Laboratory experiments in fracture propagation. SPEJ 24, 256-268 (1984). 23. Smith M. B. Stimulation design for short, precise hydraulic fractures--MHF. SPF_J 25, 371-380 (1985). 24. Nolte K. G. Fracture design considerations based on pressure analysis. Paper SPE 10911, 1982 SPE Cotton Valley Syrup., Tyler (1982). 25. Timoshenko S. P. and Goodier J. N. Theory of Elasticity, 3rd Edn, pp. 408-409. McGraw-Hill, New York (1970). 26. Detournay E. and Cheng A. H.-D. Plane strain analysis of a stationary hydraulic fracture in a poroelastic medium. Int. J. Solids Structures. In press (1990). 27. Scttari A. and Raisbeck J. M. Analysis and numerical modeling of hydraulic fracturing during cyclic steam stimulation in oil sands. JPT 33, 2201-2212 (1981). 28. Kassir M. K. and Sih G. C. Threc-dimensional thermoclastic problems of planes of discontinuities or cracks in solids. Developments in Theoretical and Applied Mechanics, Vol. 3 (W. A. Shaw, Ed.), pp. 117-146. Pergamon Press, Oxford (1967). 29. Boone T. J. and Detournay E. Response of a vertical hydraulic fracture intersecting a poroelastic formation bounded by semiinfinite elastic layers. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr. 27, 189-197 (1990). 30. Cleary M. P., Crockett A. R., Martines J. I., Narendran V. M. and Slutsky S. Surface integral schemes for fluid flow and induced stresses around fractures in underground reservoirs. Paper SPE/IX)E 11632, 1983 SPE /DOE Syrup. on Low Permeability, Denver (1983).

Acknowledgements--The work has been sponsored in part by a grant from the National Science Foundation, No. PYI 8351914. Corporate sponsorship by DoweU-Schlumberger,Digital Equipment Corporation and the Hewlett-Packard Company is gratefully acknowledged. The research was conducted at the Program for Computer Graphics, Cornell University. The counsel, comments and suggestions provided by Dr Emmanuel Detournay have been exceptionally valuable. The advice of Mr Paul Wawryznek has also been greatly appreciated.

Accepted for publication 27 February 1990. REFERENCES 1. Mendelsohn D. A. A review of hydraulic fracture modeling--l.

General concepts, 2-D models, motivation for 3-D modeling. J. Energy Res. Technol. 106, 369-376 (1984).

2. Cleary M. P. Analysis of mechanisms and procedures for producing favorable shapes of hydraulic fractures. Report SPE 9260 (1980). 3. Detournay E., Cheng A., Rovgiers J.-C. and McClennan J. Poroclastic considerations in in situ stress determination by hydraulic fracturing. Proc. Second Int. Workshop on Hydraulic Fracturing Stress Measurements, Vol. I (B. C. Haimson, J.-C. Rocgicrs and M. D. Zohack, Eds), Minneapolis (1988). 4. Boone T. J. and Detournay E. Response of a vertical hydraulic fracture intersecting a poroelastic formation bounded by semi-

Zoback, Eds), Minneapolis (1988). 7. Ishijima Y. and Rocgiers J.-C. Fracture initiation and breakdown pressure---are they similar? Proc. 24th U.S. Symp. on Rock Mechanics, College Station (1983).

14

BOONE et al.: SIMULATION OF HYDRAULIC FRACTURE PROPAGATION

31. Enever J. Tell years experience with hydraulic fracture stress measurements in Australia. Proc. Second Int. Workshop on Hydraulic Fracturing Stress Measurements, Vol. 1 (B. C. Haimson, J.-C. Roegiers and M. D. Zoback, Eds), Minneapolis, pp. 1-92 (1988).

32. Whitehead W., Gatens J. and Holditch S. Determination of/n situ stress profiles through hydraulic fracturing measurements in two distinct geologic areas. Prac. Second Int. Workshop on Hydraulic Fracturing Stress Measurements, Vol. 2 (B. C. Haimson, J.-C. I~oegiers and M. D. Zoback, Eds), Minneapolis (1988).