UBC Faculty Research and Publications

An enhanced pseudo-3D model for hydraulic fracturing accounting for viscous height growth, non-local… Dontsov, E. V.; Peirce, Anthony Feb 27, 2015

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


52383-Dontsov_Peirce_Enhance_P3D_model_27_Feb_2015.pdf [ 1004.5kB ]
JSON: 52383-1.0079338.json
JSON-LD: 52383-1.0079338-ld.json
RDF/XML (Pretty): 52383-1.0079338-rdf.xml
RDF/JSON: 52383-1.0079338-rdf.json
Turtle: 52383-1.0079338-turtle.txt
N-Triples: 52383-1.0079338-rdf-ntriples.txt
Original Record: 52383-1.0079338-source.json
Full Text

Full Text

An enhanced pseudo-3D model for hydraulic fracturing accounting forviscous height growth, non-local elasticity, and lateral toughnessE.V. Dontsov and A.P. PeirceDepartment of Mathematics, University of British Columbia, Vancouver, BC, V6T 1Z2, CanadaEmail address for correspondence: peirce@math.ubc.ca, phone: +1-604-822-2104AbstractThe goal of this paper is to develop an enhanced pseudo-3D (P3D) model for hydraulic fracturing (HF), whosepredictions are more accurate compared to that of the original P3D model, but which requires significantlyless computational resources than a fully planar HF simulator. We show that the lack of viscous resistancein the height growth and the local approximation in the computation of elastic interactions, which precludesthe incorporation of lateral toughness, are the primary weaknesses of the original P3D model that considerssymmetric stress barriers and no leak-off. To account for the viscous resistance, an apparent fracturetoughness is introduced. The apparent toughness is calibrated using a one-dimensional HF model resultingin an approximate expression that captures all regimes of propagation. To incorporate non-local elasticinteractions, the fracture opening in every vertical cross-section is approximated by a plane-strain solution,and then the 2D elasticity interaction integral is evaluated. To increase the computational efficiency, this2D integral is further approximated by two one-dimensional integrals. The use of non-local elasticity allowsus to include the asymptotic solution in the tip element, and, in particular, to include the effect of lateralfracture toughness. To further increase the accuracy of the P3D model, the flat fracture tip is replaced byits curved counterpart. This also permits us to capture radial behaviour at early times before the fracturehas reached the stress barriers. To evaluate the accuracy of the model we have developed, the results arecompared to the predictions calculated using a recently developed fully planar HF simulator, which is able tocapture viscous, toughness, and intermediate propagation regimes. It is shown that the enhanced P3D modelis able to approximate the propagation of hydraulic fractures accurately for various regimes of propagation,as well as for different fracture aspect ratios.Keywords: Hydraulic fracturing, pseudo-3D model.1. IntroductionHydraulic fracturing (HF) is a process in which a pressurized fluid is injected into rock to cause frac-ture initiation and propagation. Industrial applications of HF include accelerating the waste remediationprocess [1], waste disposal [2], preconditioning in rock mining [3], and oil and gas reservoir stimulation [4],where the latter stands out as one of the most common applications.To successfully apply HF, the treatment should be designed using an appropriate HF model or numericalsimulator. The earliest models started with the simplest one-dimensional Khristianovich-Zheltov-Geertsma-De Klerk (KGD) model [5], in which the fracture propagates under plane strain conditions and the couplingbetween viscous fluid flow and elasticity is included. Later, the Perkins-Kern-Nordgren (PKN) model [6, 7]was developed to simulate fracture propagation in a horizontally layered medium. The PKN model assumesthat the fracture propagates in the horizontal direction (i.e. along the reservoir layer) and has a constantheight that is equal to the thickness of the reservoir layer. Each vertical cross-section is assumed to beelliptic with a uniform fluid pressure, where the latter is calculated assuming that a state of plane strainprevails in each vertical plane. The primary advantage of the PKN model lies in its computational efficiency,since averaging over the vertical direction reduces the calculations to solving a one-dimensional problem.Preprint submitted to Elsevier February 27, 2015To allow for height growth, the pseudo-3D (P3D) model has been developed [8, 9, 10] as an extensionof the PKN model. There have been multiple variations of the P3D model. In particular, the geometryof the fracture was either approximated by ellipses or the fracture’s lateral dimension was divided intoelements, where each element had a specific height. The first case corresponds to the so-called lumped P3Dmodel [11, 12], while the second variation is called a cell-based P3D model [9, 10]. In addition, two types ofheight growth mechanisms were used: i) equilibrium height growth [9], where it is assumed that the pressureis uniform and a toughness propagation criterion is used to calculate the fracture height; ii) dynamic heightgrowth [10], where the height is calculated from the solution of a KGD problem (or its approximation)for each vertical cross-section. With the advent of more computational power, scientific effort has shiftedtowards developing more accurate planar 3D models (PL3D) [13, 14, 15]. In such models, the fracturefootprint is assumed to be planar and is discretized using a two-dimensional mesh. Two-dimensional fluidflow and elastic interactions between all elements are then considered. These modifications significantlyincrease the accuracy of HF simulators, but, at the same time, the CPU time increases dramatically. Thereason for this lies in the fact that both the PKN and P3D models require the solution of essentially aone-dimensional lubrication problem with a local relation describing elastic interactions, while the PL3Dmodels involve coupling between the two-dimensional lubrication operator and the two-dimensional elasticityoperator representing full non-local interactions. A more detailed review of various HF simulators can befound in [16]. Motivated by recent developments in the petroleum industry and the further developmentof computational power, researchers have started to investigate the interaction between multiple hydraulicfractures. As an example, the so-called phenomenon of stress shadowing [17] is investigated in [18]. Stressshadowing occurs when one attempts to grow multiple hydraulic fractures simultaneously. The mutual elasticinteractions produce compressive stresses that affect the evolution of neighbouring fractures and result indissimilar growth of fractures initiated from the different perforations in a stage. Finally, the complexity ofthe hydraulic fracturing problem can be increased even further by considering non-planar fractures [19].From the above discussion we observe that the complexity of hydraulic fracture modelling has increasedwith the advent of more computational power and the desire for models that represent real HF moreaccurately. Although more complicated models allow us to obtain accurate results, there are situationsin which the required computational time imposes severe limitations on the numerical experiments thatcan feasibly be performed. In this situation, a computationally efficient and accurate P3D model, whichallows us to obtain results rapidly, becomes an appropriate choice. As an example, the P3D model is usedin [20, 21] for production optimization purposes. Also, hydraulic fracture interaction is investigated usinga P3D modelling approach in [22]. In the latter case, the use of a P3D model permits the analysis ofcomplex phenomena without using excessive computational resources. Unfortunately, despite having theadvantage in computational time, P3D models are not accurate when the assumptions of the model arenot met. In particular, since viscous losses in the vertical direction are disregarded in a P3D model withequilibrium height growth, fracture height estimates are inaccurate when the value of fracture toughness issmall. Likewise, owing to the local elasticity assumption, the P3D model can be inaccurate for large valuesof fracture toughness, for which lateral fracture propagation occurs in the toughness regime. Recognizing theproblem, corrections to the P3D model, intended to increase the accuracy, have been made. Two-dimensionalfluid flow has been incorporated into a P3D model [23] making it possible to capture fracture height growthmore accurately. The work by [12] also considers two-dimensional fluid flow, but approximates the fracturegeometry by ellipses, which makes it possible to capture the effect of fracture toughness, but noticeablyrestricts the geometry of the hydraulic fracture. With regard to the PKN model, [24] mentions unpublishedwork, in which a toughness propagation criterion for PKN fractures is introduced. Also, a one-dimensionalnon-local elasticity equation for PKN fractures is derived and analyzed in [25]. The primary purpose of thatwork was to introduce a non-local relation between the fluid pressure and fracture opening, which can befurther used to develop PKN-based HF models that can capture fracture toughness in the lateral direction.Recognizing the necessity to develop a more accurate P3D model, this paper aims to introduce an en-hanced P3D (EP3D) model, which accounts for viscous height growth, non-local elasticity, and fracturetoughness in the lateral direction. The developments are based on the cell-based P3D model with equilib-rium height growth, analyzed in [26], and assume no leak-off (i.e. that the rock is impermeable). Correctionfor viscous height growth is implemented through the introduction of an apparent toughness, which ap-2proximately captures the viscous dissipation in the vertical direction. This provides a way of correctingthe vertical growth (an alternative to [23]), which is applicable to a P3D model with equilibrium heightgrowth. Following the developments in [25], non-local elasticity is incorporated in a P3D model and is thenimplemented into a HF simulator. This allows us to capture the effect of lateral fracture toughness with aP3D model. Moreover, since the EP3D model still relies on the solution of a one-dimensional problem, itpermits us to obtain results extremely rapidly.The paper is organized as follows. Section 2 briefly describes the classic P3D model with symmetric stressbarriers and equilibrium height growth, as well as its limitations. Section 3 incorporates the effect of viscousresistance in the vertical direction by calibrating the model using a one-dimensional KGD model. Then,Section 4 describes a procedure to implement non-local elasticity into a P3D model. Section 5 introduces acurved fracture tip, which also enables us to accurately capture radially symmetric behaviour at early stagesin the HF evolution. The details of the numerical implementation are presented in an appendix.2. Limitations of classical pseudo-3D model with equilibrium height growthTo motivate the necessity for improvements to the original pseudo-3D (P3D) model, this section aimsfirst, to briefly describe the P3D model with its assumptions, and then to compare the model predictionsto the output of a recently developed fully-planar hydraulic fracturing simulator that is able to capture themultiscale behaviour encountered when the processes of viscous dissipation and toughness energy releasecompete. This allows us to highlight the weak points of the P3D model that can lead to noticeable errorsfor some problem parameters.2.1. Pseudo-3D (P3D) modelThe P3D model under consideration assumes that the reservoir layer with height H is surrounded bytwo symmetric stress barriers that impose an additional confining stress ∆σ, as shown in Fig. 1. For thisproblem geometry, a fracture grows predominantly in the horizontal direction (i.e. along the x axis), butheight growth in the vertical direction is also allowed. The original P3D model formulated in [26] is based onthe following assumptions: i) the lateral fracture tip is vertical and its height is equal to H, ii) the verticalcomponent of the flux is negligible compared to the horizontal omcponent, which implies that the pressureis constant along the z axis, iii) plane strain elasticity conditions apply for any vertical (y, z) plane, iv)leak-off is limited only to the reservoir layer (i.e. no leak-off in the stress barriers) and behaves according tothe Carter model [27]. The assumptions of plane-strain elasticity and uniform pressure in each (y, z) planexzl(t)h(x, t) HFracture footprintStress barriers zyw(x, z)gFigure 1: Schematics of the P3D fracture.3allow us to find a solution for the fracture width in each cross-section asw(x, z) = 2E′ (p(x)−∆σ)χ+4∆σpiE′{χ arcsin(Hh)− zln∣∣∣Hχ+ 2zψHχ− 2zψ∣∣∣+ H2 ln∣∣∣χ+ ψχ− ψ∣∣∣}, (1)where w(x, z) is a fracture width, ∆σ is the amplitude of the stress barrier, h denotes fracture height (seeFig. 1), p(x) is the fluid pressure, χ =√h2 − 4z2, ψ =√h2 −H2, while E′ = E/(1 − ν2) is the planestrain Young’s modulus. The use of a constant pressure in the vertical direction implies that the fracturepropagates in the toughness regime. In this case, the fracture propagation criterion KI = KIc can be usedto findp = ∆σ[1 +√2piHKIc∆σ√Hh −2pi arcsin(Hh)], (2)where KI is the stress intensity factor and KIc is the fracture toughness. Equations (1) and (2) can becombined to yieldw(x, z) = 2E′√2pihKIcχ+4∆σpiE′{−zln∣∣∣Hχ+ 2zψHχ− 2zψ∣∣∣+ H2 ln∣∣∣χ+ ψχ− ψ∣∣∣}, (3)where the first term represents the fracture opening due to the fracture toughness, and the second termscapture the effect of the stress barriers. The aforementioned formulas (3) and (2) apply in the regionswhere h > H. In situations when h = H, an elliptic fracture width profile is used instead of (3), i.e.w = 2(E′)−1χp(x). The primary computational advantage of the P3D model lies in the fact that it can beformulated in terms of an effective width, defined asw¯ = 1H∫ h/2−h/2w dz. (4)Relation (3) can be integrated to obtainw¯ = HE′(√ pi2HKIc( hH)3/2+ ∆σ√h2H2−1), h > H, (5)which can be inverted to yield the function h(w¯). If one uses this result, then the fracture width becomes afunction of the effective width, i.e. w = w(w¯). In other words, the knowledge of w¯(x) is sufficient to obtainthe fracture width w(x, z). This makes the problem one-dimensional and computationally efficient. Thegoverning equation, that represents the vertically-integrated lubrication equation, can be written as∂w¯∂t +∂q¯x∂x +C ′√t− t0(x)= Q0H δ(x), (6)where the averaged flux is given byq¯x = −1Hµ′∂p∂x∫ h/2−h/2w3 dz, (7)Here C ′ = 2CL is Carter’s leak-off coefficient, Q0 is the total fluid volume pumped into the fracture per unittime, t0(x) is the inverse of the fracture length function x = l(t), and µ′ = 12µ, where µ is the fracturingfluid viscosity. Given (2)–(5) and (7), equation (6) can be solved numerically (see e.g. [26]) for w¯(x), whichcan then be used to find the fracture opening w(x, z), the height growth h(x), and the fluid pressure p(x).2.2. Drawbacks of the P3D modelWhile the computational time associated with the P3D model is not an issue, the accuracy and ap-plicability range have limitations. As was mentioned in the introduction, this paper deals with fracture4propagation in an impermeable medium (no leak-off), so, of the four P3D model assumptions, only the firstthree are relevant, namely: that the lateral fracture tip is vertical; that the pressure is constant in verticalcross sections; and that plane strain elastic conditions apply in any vertical plane (y, z).In the absence of leak-off, there are two dissipative mechanisms that control the fracture geometry: theviscous dissipation and the energy release due to formation of new fracture surface. If the viscous dissipationis dominant, then the fracture propagates in the viscous regime, otherwise it propagates in the toughnessregime. The assumptions of the P3D model imply that the fracture propagates in the toughness regime inthe vertical direction, while, since the fracture toughness in the lateral direction is not accounted for, thefracture is assumed to propagate in the viscous regime in the horizontal direction. This means that out ofthe two possible dissipative mechanisms, only one is considered in the vertical direction and, similarly, onlyone is considered in the horizontal direction. In situations when the fracture toughness is small, or evenzero, vertical height growth will be overestimated by the classic P3D model, since viscous resistance in thevertical direction, which should limit the fracture height, is ignored by the model. At the same time, for payzones with a large fracture toughness or at advanced times, the propagation in the horizontal direction maybecome toughness-dominated. In this case, the P3D model will overestimate the fracture length as the localelasticity assumption precludes accounting for the resistance due to fracture toughness when determining theprogress of the fracture front in the horizontal direction. Determining the bounds on the fracture toughness,for which the P3D model is valid, is a challenge, since hydraulic fractures are known [28] to change theirregime of propagation with time. Thus the bounds on the fracture toughness, within which the P3D modelis valid, are time-dependent. As a result, in spite of using the same problem parameters, the P3D modelmay be reasonably accurate over some time intervals and completely wrong over others.To illustrate the existence of a validity region for the P3D model and to highlight its limitations, Fig. 2compares the results of the P3D model with those obtained using the fully-planar Implicit Level Set Algo-rithm (ILSA) scheme [15, 29]. The parameters for the computation are chosen as H = 0.05 m, µ = 30.2 Pa·s,ν = 0.4, E = 3.3 GPa, Q0 = 1.7 mm3/s, ∆σ = 4.3 MPa. The rock is assumed impermeable in both models(i.e. no leak-off). Three values of fracture toughness are considered, KIc = 0, KIc = 0.16 MPa·m1/2, andKIc = 0.94 MPa·m1/2. The results of the simulations are compared at time instant t = 604 s. Note thatthese parameters (except for those with fracture toughness) correspond to the laboratory experiments in [30].Fig. 2 compares the fracture footprints (top picture), the variation of the fracture width and pressure versusx at z=0 (middle and bottom left pictures), and the variation of the fracture width and pressure versus zat x=0 (middle and bottom right pictures). The stress barrier is indicated by a thick black line. The ILSAresults are shown by solid lines, while the P3D data is shown by dashed lines. Blue lines correspond to theviscous regime, i.e. KIc = 0, red lines correspond to the toughness regime KIc = 0.94 MPa, while magentalines represent the intermediate case with KIc = 0.16 MPa·m1/2. As can be seen from the figure, the P3Dmodel is the most accurate for the intermediate case with KIc=0.16 MPa·m1/2. For this moderate fracturetoughness, the fracture length is accurately captured by the viscous regime lateral growth mode that is builtinto the P3D model, but the height growth is overestimated. The P3D solution for KIc = 0.94 MPa·m1/2provides very poor estimates of all characteristics of the fracture, namely, the fracture length, the height,the opening, and the pressure. With regard to the overestimated height growth, it is interesting to observethat the ILSA results are nearly identical for KIc = 0 and KIc = 0.16 MPa·m1/2. Since the P3D model ismore accurate for KIc= 0.16 MPa·m1/2, it might be useful to approximate the viscous ILSA solution by aP3D solution with some apparent toughness. In other words, the dashed magenta lines (i.e. the P3D resultfor KIc= 0.16 MPa·m1/2) are close to both the viscous solution and the corresponding ILSA solution withKIc=0.16 MPa·m1/2. For this reason, there is a possibility to introduce an effective or apparent toughnessin the P3D model, which can capture viscous resistance in the vertical direction. The next section describesa procedure to calculate this apparent toughness. It is also possible to upgrade the P3D model to capturethe solution for high values of toughness, e.g. KIc = 0.94 MPa, for which the original model works poorly.One option to limit the lateral growth is to introduce a fracture toughness at the corresponding fracturetip in the P3D model. However, this is not sufficient. The tip element can possibly satisfy the fracturetoughness growth criterion, but since local elasticity is used (i.e. plane strain conditions are assumed in eachvertical cross-section), there will be no elastic interaction with the rest of the fracture. Thus, the effect oftoughness will change the fracture shape only locally without affecting the global behaviour. To deal with5this problem, it is necessary to introduce a non-local or full elastic interaction, as is done for planar 3Dmodels. This inevitably leads to additional computational costs, but we demonstrate that this can be doneefficiently. The problem of including full elastic interactions into the P3D model is addressed in Section 4.3. Correction for viscous resistance in the vertical direction3.1. The conceptAs discussed in the previous section, one weakness of the P3D model is the lack of viscous resistance inthe vertical direction, which causes the fracture height growth to be overestimated, as shown in Fig. 2. Atthe same time, the comparison between the ILSA solutions for KIc=0 and KIc=0.16 MPa·m1/2 shows thata small value of fracture toughness does not significantly affect the fracture footprint, while the P3D solutionchanges significantly and is more accurate for KIc=0.16 MPa·m1/2. This establishes that a purely viscousILSA solution (i.e. with KIc=0) is approximated better by a P3D with a small toughness. The remainingquestion is how to find the appropriate value for this apparent toughness to achieve the best approximation.One possibility is to equate the fracture width calculated according to the viscous and toughness asymptoticsolutions at some distance d, in which casewK ≡√32pi∆KIcE′ d1/2 = βm(µ′VE′)1/3d2/3 ≡ wM , (8)where wK and wM denote respectively the toughness and viscous asymptotic solutions, ∆KIc is an apparenttoughness, βm=21/3 · 35/6, whileV = 12∂h∂tis the speed of the horizontal fracture fronts at the top or bottom of the blade-like fracture. Equation (8)yields the following formula for the apparent toughness∆KIc =√pi32βmE′2/3µ′1/3V 1/3d1/6. (9)One can also take an energetic approach and require a balance between the energy release rate for thetoughness regime and the viscous dissipation for the viscous regime. If the fracture propagates with avelocity V , then the energy released per unit time for a mode I crack is∂UK∂t = V∆K2IcE′ . (10)At the same time, since the fluid velocity profile inside the fracture is v(y) = 32V (1 − 4y2/w2), the totalviscous losses are∂UM∂t =∫ d0∫ w/2−w/2µ(∂v∂y)2dy dz = 3β−1m µ′2/3E′1/3V 5/3d1/3, (11)where the asymptotic solution for the fracture opening wM = βm(µ′V /E′)1/3z2/3 was used to evaluatethe integral. Note that the upper integration limit of the integral with respect to z is d, which is thecharacteristic length of the problem. By equating (10) and (11), one obtains∆KIc =√3β−1m µ′1/3E′2/3V 1/3d1/6. (12)Note that formulas (9) and (12) give almost identical results, especially given the fact that the ratio betweenthe corresponding numeric multipliers is β3/2m√pi/96 ≈ 1.01.60 50 100 1500204060x [mm]z[mm]0 50 100 150050100150x [mm]w[µm]0 20 40 60050100150z [mm]w[µm]0 50 100 1500246x [mm]p[MPa]0 20 40 600246z [mm]p[MPa]Figure 2: Comparison between the predictions of ILSA (solid lines) and the classic P3D model (dashed lines) for KIc = 0(blue lines), KIc = 0.16 MPa·m1/2 (magenta lines), and KIc = 0.94 MPa·m1/2 (red lines). The top picture shows the variousfootprints, the middle pictures show the variation of the fracture opening versus x at z=0, and versus z at x=0. The bottompictures show similar variations in the fluid pressure.73.2. Estimating the characteristic length dTo establish a framework for estimating the characteristic length d, first, without loss of generality, theconstant multiplier is omitted in (9) and (12), which effectively redefines d as∆KIc = µ′1/3E′2/3V 1/3d1/6 ⇒ d =∆K6Icµ′2E′4V 2 . (13)To find the value of the characteristic length d, note that the approach of introducing an apparent toughness isalso applicable to the Khristianovich-Zheltov-Geertsma-De Klerk (KGD) fracture geometry. For this reason,this subsection aims to utilize a KGD fracture geometry to find d. Also note that, consistent with the planestrain assumption, each vertical cross-section of a P3D fracture can be thought as a KGD fracture, seeFig. 1. Notation that is consistent with this statement is used throughout this subsection.A symmetric KGD fracture, that is subject to stress barriers and no leak-off, is governed by the lubricationequation∂w∂t −∂∂z(w3µ′∂p∂z)= Q(z), (14)and the elasticity equation, written asp−∆σ(z) = − E′4pi=∫ h/2−h/2w ds(s− z)2 , (15)where Q(z) is the source term, while ∆σ(z) denotes the stress resulting from the stress barriers (it is zerofor |z| < H/2 and ∆σ otherwise). The boundary and propagation conditions at the fracture tip arew3µ′∂p∂z∣∣∣z=±h/2= 0, w → 8KIc√2piE′(∓z ± 12h)1/2. (16)For KGD problems the source Q(z) is typically taken as Q0δ(z), i.e. as a point source. However, this is notthe case for the P3D fracture, where the fluid for each cross-section is supplied more uniformly. Away fromthe well-bore the source for each KGD-like cross-section is driven by the change in the horizontal fluxes qxin the x direction, i.e. ∂qx/∂x. Thus, since the pressure is assumed to be constant in each (y, z) plane, itfollows from Poiseuilles law that the source for each KGD-like cross section is proportional to w3(z). Forthis reason, the following distributed source is usedQ(z) = Q0w3∫ h/2−h/2 w3 dz. (17)It will be shown later that the difference between the results for this distributed source and a point sourceis small, so the source distribution does not significantly affect the characteristic length d.It is instructive to do a scaling analysis for this problem. The scales for pressure, width, length and timecan be introduced respectively asp∗ = ∆σ, w∗ =(Q0µ′E′∆σ2)1/2, h∗ =(Q0µ′E′3∆σ4)1/2, t∗ =µ′E′2∆σ3 . (18)In this case there are only two dimensionless parameters that determine the solution, namely, the normalizedfracture toughness and the normalized height of the reservoir layer, which can be calculated asK˜Ic =KIc(Q0µ′E′3)1/4, H˜ = H∆σ2√Q0µ′E′3. (19)The procedure for finding d, based on the numerical calculations for the KGD fracture, can be formulatedas follows. First, the numerical solution for the viscous regime (i.e. KIc = 0) is calculated. Then, for a8given fracture volume (Hw¯) and fracture length (h), equation (5) is used to find the apparent toughness.In other words, the fracture width profile for the viscous regime is approximated by a fracture width profilethat corresponds to the toughness regime by matching the fracture volume and length. Once the apparenttoughness is calculated, equation (13) is used to find d. This approach considers an extreme case where thefracture propagates in the viscous regime. If a good approximation is obtained for this extreme case, thenwe would expect such apparent toughness corrections to be even more accurate for less extreme cases wheresome fracture toughness is present. Since the viscous solution is used, only one dimensionless parameterin (19) affects the characteristic length d.We first consider a situation in which the fracture has not yet reached the stress barrier. In this situation,the parameter H˜ does not enter the solution, which implies that in this case d depends only on the currentfracture geometry, e.g. the length of the fracture h. Numerical calculations show that d is simply proportionalto h, i.e.d = Ch, (20)where C is a constant. Note that the above equation could also be deduced from the fact that h ∝ t2/3 for aKGD fracture propagating in the viscous regime (no stress barriers) together with (5) and (13). The valueof C obtained from numerical experiments is C = 0.175 for a distributed source (17), and C = 0.21 for apoint source. While this difference seems to be significant, keep in mind that the power 1/6 is applied to din (13). In this case, the discrepancy becomes less than 3%.The situation when the fracture has broken through the stress barrier is more complicated, since theparameter H˜ affects the solution and consequently d. It is useful to introduce a normalized apparenttoughness, which we denote by ∆K˜Ic, where the normalization is performed according to (19), and anormalized characteristic length d˜ = d/h. Equation (5) gives the relation between the normalized apparenttoughness, H˜ and time, which can be written as ∆K˜Ic(H˜, t). At the same time, equation (13) gives d˜(H˜, t).Combining these expressions and inverting the solution for the fracture length h(t), one may conclude thatd = h d˜(H2∆σ2h∆K2Ic, hH), (21)where the first parameter signifies the ratio between H˜2 and h˜∆K˜2Ic, while the second parameter is the ratiobetween the normalized height h˜ and H˜, where the latter defines h˜ = H˜h/H. The primary advantage of (21)lies in the fact that combining (20) with (13) yields an expression for ∆KIc, which does not directly dependon Q0. This makes the transition from KGD to P3D fractures easier. Given (21), there are two alternatives.Firstly, calculate d˜ numerically for a broad range of parameters, and then use interpolation to construct thecorresponding function. Secondly, find a sufficiently accurate approximation for d˜. We proceed here withthe second option, as it is simpler to implement and requires less computational resources, which makesthe overall program faster. Let’s consider the limiting case when ∆σ = 0, i.e. in which there are no stressbarriers. For this situation the answer is given by (20), so thatd˜(0, hH)= C. (22)In situations in which H2∆σ2/(h∆K2Ic)  1, the apparent toughness does not noticeably influence thefracture width (3), see (5) for the comparison of volumes. This can be interpreted as the “stress barrier”propagation regime, for which the resistance introduced by the stress barriers is larger than that of either theviscous and toughness resistances. The fracture geometry for this regime depends solely on the magnitude ofthe stress barrier and the volume of the fracture. An accurate evaluation of d˜ in this regime is not necessary,since the apparent toughness has little influence on the solution. For this reason, the constant approximationintroduced according to (22) may already be sufficiently accurate for the complete range of expected valuesof the parameter H2∆σ2/(h∆K2Ic) (since it works for both limiting cases). To introduce a more accurateapproximation, we try to find a characteristic distance from the fracture tip to the region where the influenceof toughness is suppressed by the stress barriers. To do this, it is necessary to equate the fracture widthterms in (3), one due to the toughness and the other due to ∆σ. Since d should be proportional to this9characteristic distance, one may use the asymptotic expansion in (3), also assuming that h/H  1, to findd˜ ∝ h1/2∆KIcH∆σ  1. (23)One of the simplest approximations that can be constructed from the above result and (22) isd˜(H2∆σ2h∆K2Ic, hH)≈ C1 + C2H∆σh1/2∆KIc, (24)where C2 is a constant.The combination of (13), (21) and (24) leads to∆KIc = C1/6µ′1/3E′2/3V 1/3h1/6(1 + C2H∆σh1/2∆KIc)−1/6, (25)which can be solved to find the apparent toughness ∆KIc. To check the error that is caused by theapproximation (24), the numerical solution of the KGD problem (without toughness) is compared to theanalytical solution (3) with the same fracture volume and the apparent toughness calculated accordingto (25). The numerical solution for the KGD fracture is calculated for the parameters E′ = 25 GPa,µ′ = 1.2 Pa, Q0 = 5×10−4 m2/s and H = 50 m. Firstly, simulations with a maximum pumping timet = 3×104 s and various magnitudes of stress barrier ∆σ = 0.25− 10.5 MPa are performed. The error, Eh,is defined as the maximum discrepancy for the fracture length for all time steps. To find the optimal valuefor C2, Fig. 3 shows the variation of Eh versus H˜ (as this is the only parameter that affects the viscoussolution of the KGD fracture with stress barriers) for different values of C2. One can observe that evenwith C2 = 0, the error is below 4%, while C2 = 1/2 reduces the maximum error to less than 1.5%. Notethat C2 = 1/2 does not exactly correspond the the minimum of the error, but is sufficiently close that itcan practically be considered a minimum. To illustrate the level of approximation of the fracture opening,Fig. 3 also shows the comparison between the widths stemming form the numerical solution (solid lines)and the approximation (3) (dashed lines) at t = {250, 500, 1000, 2000, 4000} s and C2 = 1/2. To coverdifferent regimes, three values of the stress barrier are considered, namely, ∆σ = 0.5 MPa, ∆σ = 1.5 MPaand ∆σ = 3.5 MPa. These values correspond to H˜ = 0.13 (a weak stress barrier), H˜ = 1.16 (a moderatestress barrier), and H˜ = 6.33 (a strong stress barrier). When the stress barrier is weak, the fracture widthresembles that without a stress barrier, but features a small bump at the centre. The difference between theviscous numerical solution and the approximation is visible, but not significant. When the stress barrier isstrong, neither the viscosity nor the apparent toughness affect the fracture opening, in which case there is asmall discrepancy between two solutions that occurs only in the near-tip regions. The case with H˜ = 1.16represents an intermediate case, for which the fracture length error reaches a maximum.3.3. Application to P3DTo apply the viscous height growth correction to the P3D model, first, the inclination of the fracturefootprint located within the stress barrier is neglected. This allows us to treat each vertical cross-section ofa P3D fracture as a KGD fracture and to apply the correction for the viscous resistance developed above. Interms of the numerical implementation of the P3D model with a viscous correction, the same equations (2)–(7) need to be solved, in addition to the value of the apparent fracture toughness for each vertical cross-section, which needs to be calculated independently. If there is some non-zero fracture toughness, then thesum of this material fracture toughness and the apparent toughness is used in the model, i.e. KIc+∆KIc.To find the apparent fracture toughness at each time step, (25) can be rewritten as∆KIc = C1/6µ′1/3E′2/3(12hj+1(KIc+∆KIc)− hjtj+1 − tj)1/3h1/6j+1(1 + C2H∆σh1/2j+1∆KIc)−1/6, (26)10−100 −50 0 50 10000.0050.010.0150.020.025z [m]w[m]10−2 100 10200.˜E h−100 0 10000.0050.010.015z [m]w[m]−150 −100 −50 0 50 100 15000.0020.0040.0060.0080.010.012z [m]w[m]H˜ = 0.13 H˜ = 1.16H˜ = 6.33 C2=0C2=1/4C2=1/2C2=3/4C2=1Figure 3: The top pictures and the bottom left picture show the comparison between fracture openings calculated nu-merically assuming no toughness (solid lines) and approximations computed using apparent toughness (dashed lines) fort = {250, 500, 1000, 2000, 4000} s and H˜ = 1.16, H˜ = 0.13, and H˜ = 6.33. The bottom right picture shows the fracturelength error versus H˜ for different values of C2 (markers), the solid lines represent the spline interpolant.where the dependence of h on fracture toughness, represented by hj+1(KIc+∆KIc), comes from (5). Here hjis fracture height at some point x at time ti, while hj+1 is an unknown fracture height at time instant tj+1.The solution of the nonlinear equation (26) allows us to find ∆KIc and consequently hj+1 at each pointin space, which are then used to solve the lubrication equation (6). The details of the numerical schemeare summarized in the Appendix. To check the accuracy of the viscous correction, Fig. 4 compares thefootprints calculated using the P3D model with the viscous correction with the footprints of the referencesolution given by ILSA. The parameters used for these calculations are the same as those for Fig. 2, exceptthat the fracture toughness is KIc = 0. To show the evolution of the fracture, the footprints are plottedat different time instants, namely, t = {37, 101, 200, 401, 604, 1048} s. Comparison between the ILSA andP3D footprints in Fig. 4 reveals that the corrected P3D model is able to capture the height growth of thefracture even in the absence of fracture toughness. This modification thus fixes one of the major drawbacksof the P3D model, namely its inability to capture the height growth for small (or zero) values of KIc. Atthe same time, for large values of the fracture toughness ∆KIc/KIc  1, the effect of the viscous correctionis effectively removed.4. Incorporating non-local elastic interactions into P3D model4.1. Non-local elasticityAs was mentioned in Section 2.2, there are two main drawbacks of the original P3D model [26], namely,its inability to capture the height growth for a small or vanishing fracture toughness, and the significantinaccuracies for large values of KIc, see Fig. 2. The first weakness is addressed in Section 3 through the110 50 100 150 2000204060x [mm]z[mm]Figure 4: Comparison of fracture footprints for KIc=0 generated by the P3D model with viscous resistance (dashed lines) andILSA scheme (solid lines) at different time instants t = {37, 101, 200, 401, 604, 1048} s.introduction of an apparent or effective fracture toughness that controls height growth, while this sectionaims to correct the second weak point of the P3D model by including non-local elastic interactions within thefracture. When full elastic interactions are used, formally, the pressure is no longer constant in each verticalcross-section and the fracture width cannot be calculated using the plane strain elasticity assumption (3).At the same time, Fig. 2 shows that the the fracture width variation in the vertical z direction, calculatedusing ILSA, resembles that obtained using the P3D model – even in the viscous regime, i.e. KIc = 0. Thisis further supported by Fig. 3, in which the viscous KGD solutions are compared to the apparent toughnesssolutions. Since the difference between the viscous and apparent toughness solutions is small (as soon asthe fracture height is correct), one may use the approximate analytical solution (3) to estimate the elasticinteractions. In particular, the pressure relation (2) can be replaced withp(x) = − E′8pi∫ l(t)−l(t)∫ 12h(x′,t)−12h(x′,t)w(x′, z′) dz′ dx′((x′−x)2 + z′2)3/2 , (27)where the function w is given in (3), while h and l are the height and the half-length of the fracturerespectively, see Fig. 1. Note that since the P3D model assumes that the pressure is uniform in the verticaldirection z, only one pressure value for each vertical cross-section is required. For this reason, (27) evaluatesthe pressure along the x axis, i.e. assumes z=0. For the purpose of numerical implementation, equation (27)is discretized and the integrals are evaluated numerically. This effectively gives the relation between p(x) andw¯(x), which is analogous to the elasticity equation for the KGD fracture (15). The details of the discretizationof (27) are enclosed in the Appendix. The use of (27) permits us to incorporate the asymptotic solutionfor the tip element, as was done in [15]. For the problem under consideration, either the viscous or thetoughness asymptotic solution is used, see the Appendix for details. To illustrate the effects produced bythe use of (27), Fig. 5 compares the fracture footprints generated by the P3D model with full elasticity (andthe correction for the viscous resistance) to those produced by the ILSA model. The parameters used forthe calculations are the same as those for Fig. 2, except that only the following two values of the fracturetoughness are considered, namely, KIc = 0 and KIc = 0.94 MPa·m1/2. One can see that the full elasticityinduces only slight changes for the viscous solution (see Fig. 4). However, the comparison between Fig. 2and Fig. 5 shows that the use of the non-local elasticity (27) (together with the asymptotic solution forthe tip element) allows us to approximate the ILSA solution extremely well for the case of large toughness.Thus, the use of the viscous correction and full elasticity allows us to overcome the two major limitationsof the original P3D model, namely, its applicability for small and large values of fracture toughness.In spite of the fact that Fig. 5 shows good agreement between the P3D model and the reference ILSAsolution, another consequence of using the full elasticity is a concomitant increase in computation time.120 50 100 150 20002040x [mm]z[mm]0 20 40 60 80 100 120 140 16002040x [mm]z[mm]Figure 5: Comparison of the fracture footprints for KIc = 0 (top) and KIc = 0.94 MPa·m1/2 (bottom) generated by the P3Dmodel with viscous resistance and full elasticity (dashed lines) and the ILSA scheme (solid lines) at different time instantst = {37, 101, 200, 401, 604, 1048} s.In particular, the numerical evaluation of the double integral in (27) requires O(N2xNz) operations (recallthat the pressure is evaluated at z = 0), where Nx and Nz denote the number of elements in the x and zdirection respectively. As a comparison, the original P3D model uses only O(Nx) operations to evaluatepressure, a KGD model with a similar mesh requires O(N2x) operations, while a fully planar model needsO(N2xN2z ) operations. This demonstrates that the introduction of non-local elasticity affects the efficiencyin the computation of the pressure appreciably, and poses a significant hurdle in developing an efficientcorrected P3D model. At the same time, the pressure evaluation is still much faster than that of a fullyplanar hydraulic fracturing model, which, together with the increased accuracy, still makes the correctedP3D model a good competitor.4.2. Approximate non-local elasticityFig. 5 indicates that the use of the non-local elasticity significantly improves the accuracy of the P3Dmodel, especially for large values of fracture toughness. At the same time, numerical evaluation of the doubleintegral in (27) requires extra computational resources and slows down the overall computation time. Thissubsection aims to introduce an approximation to the integral (27), which would decrease time needed tocalculate the pressure numerically.First, following [25], it is noted that if one substitutes an elliptical width profile w(x, z) = w0(x)√1− 4z2/h213in (27), the inner integral with respect to z′ can be evaluated analytically to obtainp(x) = −E′8pi∫ l(t)−l(t)w0(x′)∫ 12h(x′,t)−12h(x′,t)√1− 4z′2/h2 dz′((x′−x)2 + z′2)3/2 dx′= E′2pi∫ l(t)−l(t)w0(x′)h(x′, t)|x′−x|[K(− h(x′, t)24(x′−x)2)− E(− h(x′, t)24(x′−x)2)]dx′ (28)= E′2pi∫ l(t)−l(t)w0(x′)h(x′, t)− (x′−x)∂h(x′, t)∂x′dG(2(x′−x)/h(x′, t))dx′ dx′,where K(·) and E(·) represent the complete elliptic integrals of the first and the second kind. It should benoted here that the arguments of the complete elliptic integrals are defined through the relations K(m) =∫ pi/20 (1−m sin2 θ)−1/2 dθ and E(m) =∫ pi/20 (1−m sin2 θ)1/2 dθ, and the kernel G(·) is given byG(s) =√1 + s2s E( 11 + s2). (29)The asymptotic behaviour of the kernel G(·) isG(s) ≈ pi2 sign(s), |s|  1, G(s) ≈1s , s 1.Thus, at small distances it behaves similar to a KGD kernel, while for large distances it permits us to recoverthe local elasticity relation associated with the plane strain regime. Indeed, to obtain a plane strain limit,one may set h(x′) = const., w0(x′) = const. and take a limit l(t) → ∞ to obtain p = E′w0/(2h), which isthe local pressure-width relation that is consistent with the plain strain limit.The first term in (3) has an elliptic form and can be integrated using (28). Unfortunately, the rest ofthe solution, which consists of the logarithmic functions, cannot be integrated analytically. However, it canbe approximated by an ellipse with some suitable parameters. To this end, let’s approximate (3) byw(x, z) ≈ w1(x)√1−4z2/h2 + w2(x)√1− 4z2/h2∆σ, (30)w1(x) =2E′√2hpi KIc, w2 =2∆σHpiE′ ln(h+√h2−H2h−√h2−H2),where the first term, due to toughness, remains unchanged, while rest of the solution is approximated by anellipse, whose maximum opening is the same as that for the logarithmic terms in (3), while its height h∆σ istaken as an unknown. It is implicitly assumed that |z| < h/2 for the first term in (30), and |z| < h∆σ/2 forthe second term in (30). To find the corresponding pressure, one needs to substitute (30) into (28), whichyieldsp(x) = E′2pi∫ l(t)−l(t)w1(x′)h(x′, t)− (x′−x)∂h(x′, t)∂x′dG(2(x′−x)/h(x′, t))dx′ dx′,+ E′2pi∫ l(t)−l(t)w2(x′)h∆σ(x′, t)− (x′−x)∂h∆σ(x′, t)∂x′dG(2(x′−x)/h∆σ(x′, t))dx′ dx′ (31)To find the unknown h∆σ, one may consider a plane strain limit of (31) (h = const., h∆σ = const.,14w1 = const., w2 = const., and l→∞), which givespps =E′w12h +E′w22h∆σ=√2pihKIc +∆σHpih∆σln(h+√h2−H2h−√h2−H2), (32)where pps denotes the plane strain limit of p(x) calculated using (31). Since, for the case of plane strain,the pressure is known and is given by (2), one may find h∆σ by equating this plane strain pressure and (32)to obtainh∆σ =Hpi ln(h+√h2−H2h−√h2−H2)(1− 2pi arcsin(Hh))−1. (33)For completeness, it is noted that only one ellipse is used in situations when h=H, in which case w1(x) =4w¯(x)/pi and w2(x) = 0. The implementation of the approximation for the non-local elastic interactionsrequires us to replace the procedure for calculating the pressure using (27) with the approximation givenby (31). More details about numerical schemes for the P3D model with non-local elasticity and with thisnon-local elasticity approximation are given in the appendix.To illustrate the accuracy of the approximation, Fig. 6 compares fracture opening versus z at the borehole(i.e. x=0) and the pressure variation versus x calculated using the full elasticity integral (27) (dashed blueand red lines), the approximation given by (30) and (31) (solid black lines), and the reference ILSA solution(solid blue and red lines). The parameters of the problem are chosen to be the same as for the Fig. 2. Twovalues of fracture toughness KIc=0 (blue lines for ILSA and full elasticity P3D) and KIc=0.94 MPa·m1/2(red lines for ILSA and full elasticity P3D) are considered, and all results are presented for time instantt = 604 s. As can be seen from the figure, the pressure, calculated using the approximation (31), is very closeto that calculated using (27). Moreover, the difference between the two is much smaller than the differencebetween one of them and the reference ILSA solution. This shows that the use of an approximation inthe evaluation of elasticity integral does not introduce a noticeable error in the pressures. In terms of thefracture width, there is some visual difference between the two-ellipse approximation given by (30) and theP3D solution that utilizes full elasticity and whose width is calculated using (3). However, equation (3)can also be used for the purpose of calculating the fracture width of an approximate elasticity solutioninstead of (30). In this case, the two P3D solutions become nearly indistinguishable, similar to the pressurevariations.With regard to the numerical implementation, the pressure calculation using equation (31) requires onlyO(N2x) operations, where Nx is the number of elements in the x direction. This is much smaller than thenumber of operations required for the evaluation of (27), which is O(N2xNz), where Nz is the number ofelements in the vertical direction. Also, this calculation shows that the computation times for the P3Dmodel are comparable to that for a KGD fracture with the same number of elements. Note that althoughthe calculation of the pressure is one of the slowest parts of the algorithm, there are many other operationsthat affect the overall computational efficiency.5. Curved fracture tipFig. 5 shows that the use of non-local elastic interactions in a P3D model, described in Section 4.1, leadsto a more accurate approximation. Moreover, Section 4.2 describes how to avoid the numerical evaluationof the elasticity integral (27), which makes the calculations faster. The P3D model with these improvementsmay already be sufficient for an accurate and fast hydraulic fracturing simulator. However, it it possible toimprove the model even further. The modified P3D model has been improved to such an extent that thebiggest discrepancies between the P3D and the reference ILSA solutions, shown in Fig. 5, are due to the flattip assumed by the P3D model. To overcome this problem, this section aims to introduce a curved fracturetip to allow the P3D geometry to resemble the ILSA geometry more closely. Since full elasticity is usedin the P3D model, and a toughness propagation criterion can be used in both the vertical and horizontaldirections, it should be possible to match the solution to a radial fracture propagating in the toughnessdominated regime. This naturally introduces a curved fracture tip for early times when the fracture isradial, and then a similar approach can be used later when the fracture breaks through the stress barrier.150 20 40050100150z [mm]w[µm]0 50 100 1500246x [mm]p[MPa]0 20 40050100150z [mm]w[µm]0 50 100 1500246x [mm]p[MPa]10 12 14859095100105z [mm]w[µm]62 62.5 63 63.52.962.9833.023.04x [mm]p[MPa]14 16 18114116118120122z [mm]w[µm]36 38 40 42 444.84.95x [mm]p[MPa]10 12 14859095100105z [mm]w[µm]62 62.5 63 63.52. 62.9833.023.x [mm]p[MPa]14 16 18114116118120122z [mm]w[µm]36 38 40 42 444.84.95x [mm]p[MPa]0 12 14859095100105z [mm]w[µm]62 62.5 63 63.52.962.983.023.04x [mm]p[MPa]14 16 181116118120122z [mm]w[µm]36 38 40 42 444.84.95x [mm]p[MPa]10 2 14859095100105z [mm]w[µm]62 62.5 63 63.52.962.9833.023.04x [mm]p[MPa]14 16 18114116118120122z [mm]w[µm]36 38 40 42 444.84.95x [mm]p[MPa]Figure 6: Comparison between the reference ILSA solution (blue and red solid lines), the P3D solution with non-local elasticity(blue and red dashed lines) and the P3D solution with approximate non-local elasticity (black lines) in terms of fracture widthvariation versus z at x=0 (left column) and pressure variation versus x (right column). The parameters for the computationare the same as for Fig. 2. Two values of fracture toughness KIc = 0 (top row) and KIc = 0.94 MPa·m1/2 (bottom row) areconsidered, and all results are presented for time instant t = 604 s.To introduce a curved fracture tip in a P3D model, we try to match the solution for a radial fracturepropagating in the toughness regime. Using the same coordinate system as shown in Fig. 1, the fractureopening for the radial solution iswrad =4KIc√pilE′√l2 − x2 − z2, (34)where l is the half-length of the fracture and its radius at the same time. The corresponding fracture heightishrad = 2√l2 − x2. (35)Equation (34) can be used to calculate the effective width asw¯rad =1H∫ hrad/2−hrad/2wrad dz =piKIc2HE′√pilh2rad, (36)which can be inverted to findhrad = w¯1/2rad (2l)1/4√√2HE′√piKIc. (37)Equation (37) is the analog of (5), but for the radial fracture, since it relates the height of the fracture toits effective width. To find the fracture width as a function of the effective width, one may use (34), (35),16and (36) to findwrad =4Hpihradw¯rad√1−( 2zhrad)2, (38)which is the analog of (3).To be able to capture radial fracture propagation with the P3D model, equation (37) needs to beincluded into the P3D model. First, if the fracture toughness is not zero, equation (5) can only be invertedfor w¯ >√piH/2KIc/E′. For smaller w¯, an elliptic fracture opening with the height H is used in the originalP3D model. This can be clearly seen in Fig. 5, where there is a region of constant fracture height near thetip for the case of non-zero toughness. To capture a radial fracture and a curved fracture tip, equations (37)and (38) should be used in this near-tip region. In situations when 2l > H, however, equation (37) is nolonger valid since the fracture has reached the stress barrier. To ensure that the fracture height is continuouswhen 2l > H (given the continuity of the effective width w¯), equations (5) and (36) should give the sameresult for h = H. This implies that the fracture height and width in the tip region should be calculated ash = w¯1/2(min{2l,H})1/4√√2HE′√piKIc,w = 4Hpih w¯√1−(2zh)2, w¯ <√piH2KIcE′ , (39)while outside of the tip region, i.e. for w¯ >√piH/2KIc/E′, equations (5) and (3) are used. Usingequations (39) in the tip region, and (3) and (5) elsewhere, ensure that the P3D fracture propagating inthe toughness regime is radial at early times and features curved fracture tips later when the fracture hasbroken through the stress barrier.Unlike the toughness regime of propagation, the viscous regime has no near-tip region within which thefracture height is constant and equal to H. This makes the introduction of a curved tip more complicated.First, a zone with constant fracture height near the tip needs to be introduced, and only then one canincorporate a curved fracture tip. To produce a zone with a constant fracture height near the tip, someeffective fracture toughness, ∆KIc has to be introduced. Given the effective toughness, equation (39) can beused to produce a curved tip. In situations when both the fracture toughness and the apparent toughnessare relevant, (39) is replaced byh = w¯1/2(min{2l,H})1/4√ √2HE′√pi(KIc+∆KIc),w = 4Hpih w¯√1−(2zh)2, w¯ <√piH2KIc+∆KIcE′ , (40)which is able to capture both viscous, toughness, and intermediate regimes (approximately). Even though (40)comes from the solution for a radial fracture propagating in the toughness regime, given the level of ap-proximation of the P3D model, it can still be used for the fractures propagating in the viscous regime.The remaining question is how to find the value of the apparent toughness that would give a reasonableapproximation in the viscous regime. One possibility is to introduce an apparent toughness based on thelateral propagation. In particular, a formula that is similar to (25) with ∆σ = 0, can be used, i.e.∆KIc = C1/63 µ′1/3E′2/3(dldt)1/3d1/6, (41)where C3 is a constant. Equation (41) gives the value of the apparent toughness for the lateral (i.e. along xaxis) propagation. This apparent toughness is then used together with (40) to produce a curved fracture tip.The original formula (25) is obtained for a KGD fracture, however, it utilizes near-tip asymptotic solutionsthat also apply for other geometries with smooth boundaries (see [15] for a justification of this). For this17reason, (41) can be used to calculate the effective toughness in the lateral direction. The remaining taskis to determine the characteristic length d and the constant C3. At early times, when the fracture shouldbe radial, it is clear that d = l, i.e. equal to fracture radius. It is found through a series of numericalexperiments that C3 = 1.2 makes the aspect ratio of the P3D-modeled fracture equal to 1, even though thefracture is not exactly radial. Once the fracture has reached the stress barrier, ILSA simulations show thatthe size and shape of the region where h < H does not change as the fracture propagates. This impliesthat the characteristic length is constant and is related to the horizontal size of the region with h < H,which is a small fraction of H. Numerical experiments show that d = aH (a ≈ 1.1 × 10−2) allows us tomatch the ILSA solutions accurately. It is clear that the characteristic length d changes dramatically as thefracture breaks through the stress barrier, but it should be noted that the corresponding apparent toughnesschanges by less than a factor of two. To preclude a jump in the apparent toughness, a continuous transitionis implemented. In particular, the characteristic distance is calculated asd = H( 12H−l) l +H(l− 12bH) aH +H(l− 12H)H( 12bH−l)( 12bH − l) l + (l − 12H) aH12bH − 12H, (42)where H is a Heaviside step function and b = 1.5. Equations (40), (41) and (42) together with the values ofC3, a and b allow the P3D model to capture a radial solution at early times (it is actually not exactly radialfor the viscous and intermediate regimes) and introduce a curved fracture tip for fractures that have brokenthrough the stress barrier.The method described here has a rigorous justification for the toughness regime, while it is more empiricalfor the viscous and intermediate regimes, for which a number of fitting parameters are used. In addition,the introduction of a curved fracture tip is a correction, which improves the representation of the fracturegeometry, and increases the overall accuracy of the P3D model.To check the accuracy of the P3D model that is augmented with the curved tip, Fig. 7 compares theresults of the corresponding P3D model (dashed lines) to the reference ILSA solutions (solid lines). Theparameters for the computation are the same as for Fig. 6. Blue lines correspond to the viscous solution,i.e. with KIc = 0, while the red lines represent the results with KIc = 0.94 MPa·m1/2. The top twopictures compare fracture footprints at times t = {37, 101, 200, 401, 604, 1048} s, in which the footprintsthat correspond to t = 604 s are represented by thicker lines. The bottom four pictures are similar to thecorresponding bottom four pictures of Fig. 6, as they show the comparison in terms of fracture width, w,and pressure, p, versus x at z = 0, and fracture width, w, and pressure, p, versus z at x= 0 at t = 604 s.Fig. 7 shows that the P3D model with a curved fracture tip is more accurate than the equivalent P3Dmodel with a flat tip, see Fig. 5. The fracture footprint is represented accurately from early times, when thefracture is radial, to the end of the simulation, when the fracture shape is elongated due to the presence ofstress barriers. This is true for both viscous (blue lines) and toughness (red lines) regimes of propagation.The comparison of the fracture widths and pressures at t = 604 s further confirms the accuracy of the newP3D model. One can also observe that the constant pressure approximation in each vertical cross-section(i.e. along z direction) is a reasonable assumption even for the viscous regime, since noticeable pressuredeviations are localized near the source and near the fracture tip. It should be noted, that the numericalscheme for the P3D model switches between the toughness and viscous asymptotic solutions (for the tipelement), so it is able to accurately represent the viscous or toughness regimes, but not the intermediatecase. However, numerical simulations show that the error due to such a switch in the asymptotic solutionat the tip does not exceed approximately 5% in terms of fracture length. In addition, this error can beremoved by using an appropriate asymptotic solution for the M-K edge, see [29].It is important to keep in mind the assumptions behind the modified P3D model. Clearly, since fullelastic interactions are used, the fracture does not always have to be mature (i.e. elongated) for the P3Dmodel to apply, which is confirmed in Fig. 7. At the same time, a plane-strain solution for the fracturewidth is used for every vertical cross-section. So, considerable discrepancies are expected in situations whenthe actual fracture width is not approximated accurately by a plane-strain solution. This can happen whenthere are both pronounced height growth through the stress barriers and the aspect ratio of the footprint issmall. When both the break-through and the aspect ratio are small, the width profile is close to an ellipse,180 50 100 150 20002040x [mm]z[mm]0 20 40 60 80 100 120 140 16002040x [mm]z[mm]0 50 100 150050100150x [mm]w[µm]0 20 40050100150z [mm]w[µm]0 50 100 1500246x [mm]p[MPa]0 20 400246z [mm]p[MPa]Figure 7: Comparison between the predictions of the P3D model including the approximate full elasticity and curved fracturetip (dashed lines) with the reference ILSA solutions (solid lines) for the same problem parameters as for the Fig. 6. The toptwo pictures compare fracture footprints at different time instants t = {37, 101, 200, 401, 604, 1048} s for KIc = 0 (blue lines)and KIc=0.94 MPa·m1/2 (red lines). The thick solid lines correspond to t = 604 s. The bottom four pictures compare fracturewidths and pressures versus x at z = 0 and fracture widths and pressures versus z at x= 0 at t = 604 s. Similarly, blue linescorrespond to KIc=0, while solutions with KIc=0.94 MPa·m1/2 are indicated by red lines.190 10 20 30 40 50051015x [m]z[m]0 5 10 15 20 25 30024681012x [m]z[m]Figure 8: Comparison between the fracture footprints generated by the P3D model using a curved tip (dashed lines) andILSA scheme (solid lines) for the set of reference parameters, Young’s modulus E = 1 GPa, KIc = 0 (top picture) andKIc = 1.5 MPa·m1/2 (bottom picture). Footprints are shown at time instants t = {51, 100, 204, 403, 599, 801} s for the toppicture and t = {50, 100, 200, 399, 601, 801} s for the bottom picture.and the P3D model is accurate, see footprints at early times in Fig. 7. When the the fracture height growthis substantial, but the aspect ratio is high (i.e. the fracture is elongated), the conditions in each cross-sectionare close to plane-strain, in which case the P3D model again performs well.To examine the accuracy of the new P3D model at the limits of its applicability, two extreme cases areconsidered next. The set of reference parameters for the extreme cases is H = 10 m for the height of thereservoir layer, µ = 0.1 Pa·s for the fluid viscosity, ν = 0.2 for the Poisson’s ratio, Q0 = 0.01 m3/s for theinlet flux and ∆σ = 0.25 MPa for the magnitude of the stress barriers. Two values of Young’s modulus areconsidered, namely E = 1 GPa and E = 2 GPa. In addition, two values of fracture toughness are considered,KIc = 0 and KIc = 1.5 MPa·m1/2. Fig. 8 shows the comparison between the fracture footprint predictionsof the P3D model (dashed lines) and the ILSA solutions (solid lines) at different time instants for the caseE = 1 GPa. The top picture corresponds to the case with no toughness (blue lines), while bottom picturecompares the results for KIc=1.5 MPa·m1/2 (red lines). As can be seen from the figure, the P3D model isaccurate at early times for both the viscous and toughness regimes, while it becomes less accurate for themore mature fractures.Moreover, the discrepancy is larger for the toughness solution, for which the aspect ratio is smaller. Sincethe toughness solution becomes significantly less accurate, it is clear that it is not the viscous correction20that affects the discrepancy for the viscous solution, but the elasticity. The cause of the inaccurate elasticinteractions is the use of a plane strain solution for each vertical cross-section. The errors for the fracturelengths for t = 801 s are 2% and 5% respectively for the viscous and toughness solutions. The errors for thefracture heights at the origin for t = 801 s are 5.4% and 8.5% respectively for the viscous and toughnesssolutions. To highlight the degree of improvement of the P3D model, this error needs to be compared withthe accuracy of the original P3D model for t = 801 s, which predicts a fracture half-length of 48 m and half-height of 43 m for the viscous case, and a fracture half-length exceeding 80 m and almost no height growthfor the case KIc= 1.5 MPa·m1/2. The even more extreme case with E = 2 GPa is considered next. Fig. 9shows a similar comparison between the predictions of the new P3D model (dashed lines) and the referenceILSA solutions (solid lines) for the viscous case (blue lines) and KIc = 1.5 MPa·m1/2 (red lines). Thedisagreement between the predictions of the P3D model and ILSA solutions is more pronounced comparedto that in Fig. 8. This is related to the smaller aspect ratio of the footprints. In this case the plane-strainsolution for the fracture width in each vertical cross-section deviates significantly from the correspondingILSA width profile. The errors of the fracture lengths for the last considered time instant are 4% for theviscous solution and 7% for the toughness solution. The errors of the fracture heights at the origin fort = 800 s are 8% and 13% for the viscous and toughness solutions respectively. The magnitude of the errorshould be judged based on the fact that Fig. 9 represents an extreme case, which challenges the applicabilityof the P3D model. To establish a reference point, it is instructive to compare these results to the predictionsof the original P3D model. For the viscous case, i.e. the top picture in Fig. 9, the original P3D modelpredicts a fracture half-length of 31.5 m and half-height of 293 m for t = 800 s. These numbers show thatthe predictions of the original P3D model are very inaccurate compared to the ILSA solution, and cannotbe used for any practical applications. The situation with the toughness regime, i.e. KIc = 1.5 MPa·m1/2,is even worse, since the original P3D model becomes unstable as soon as the height of the fracture reachesa critical value. As follows from [26], the value hu beyond which the classic P3D model exhibits runawayheight growth, is given byhu =8 +√K4 + 64K2 H, K =√2piHKIc∆σ ,which is approximately equal to hu ≈ 14 m for the considered problem parameters. As can be seen fromFig. 9, the fracture height (the ILSA solution) exceeds the critical value at early times (note that the criticalhalf-height is 7 m). Thus, due to runaway height growth, the classic P3D model fails to give a solution,while the corrected P3D model continues to give a solution albeit somewhat less accurate.6. SummaryThis paper aims to improve the original P3D model for hydraulic fracturing by making the predictionsmore accurate, but not compromising the computational efficiency. In the case of no leak-off, the primaryweaknesses of the classic P3D model are its inability to capture the viscous resistance in the vertical directionand its inability to account for non-local elastic interactions. The latter, in particular, precludes capturingthe toughness propagation regime in the lateral direction. To include the viscous resistance in the verticaldirection, a toughness solution with an apparent toughness is used. To obtain an expression for the ap-propriate apparent toughness, a similar approach is applied to a one-dimensional KGD hydraulic fracturingmodel, where a fracture propagating in the viscous regime is approximated by a solution that corresponds tothe toughness regime and features an apparent toughness. A comparison with the results obtained using thefully planar HF simulator ILSA show that the introduction of an apparent toughness allows us to capturefracture height growth accurately using the P3D model even when there is no toughness. To account for thetoughness propagation criterion in the lateral direction, non-local elasticity is introduced. First, the fracturewidth in every cross-section is approximated by a plane strain solution, and then the appropriate elasticityintegral is evaluated numerically. This approach also allowed us to include the asymptotic solution in thetip element, which makes it possible to capture fracture propagation in either the viscous or the toughnessregimes. Comparison with the ILSA results demonstrate that the use of non-local elasticity significantlyincreases the accuracy of the P3D model, especially when the lateral fracture propagation is in the toughness210 10 20 30 40 500510152025x [m]z[m]0 5 10 15 20 25 30 35051015x [m]z[m]Figure 9: Comparison between fracture footprints generated by the P3D model using a curved tip (dashed lines) and ILSAscheme (solid lines) for the set of reference parameters, Young’s modulus E = 2 GPa, KIc = 0 (top picture) and KIc =1.5 MPa·m1/2 (bottom picture). Footprints are shown for time instants t = {51, 112, 203, 401, 600, 800} s for the top pictureand t = {52, 99, 201, 401, 600, 801} s for the bottom picture.22regime.To make the computation of non-local elastic interactions faster, the fracture opening in the verticaldirection is approximated by two ellipses, for each of which the elasticity integral is calculated analyticallyover the vertical axis. This effectively reduces the computation of the 2D integral to two 1D integrals, whichmakes the resulting P3D model significantly more efficient. The last improvement of the P3D model lies inreplacing the flat fracture tip by a curved tip. To do this, the P3D model is tuned to match the analyticsolution for a radial fracture propagating in the toughness regime. In this case, mature fractures (i.e. thosethat have reached the stress barriers) feature a tip region, in which the properties of a radial fracture dictatethe shape of the fracture tip.A radial fracture propagating in the viscous regime is also accounted for approximately by introducingan apparent toughness in the lateral direction. The accuracy of the resulting enhanced P3D model is verifiedby comparing the results to the predictions of the fully planar HF simulator ILSA. It is shown that the P3Dmodel is able to capture the evolution of hydraulic fractures from the initial stages, when the fracture isradial, to mature stages when the fracture has broken through the stress barrier. Moreover, this evolutioncan be captured when fractures propagate in either the viscous or toughness regimes. To examine the limitsof the applicability, extreme cases are considered. Even for such extreme parameters, the accuracy of theenhanced P3D model is on the order of 10%, while the error of the original P3D model is measured inhundreds of percent.AcknowledgementsThe authors gratefully acknowledge the support of the British Columbia Oil and Gas Commission andthe NSERC discovery grants program.Appendix: Numerical schemeThis appendix aims to describe the numerical solution for the P3D model with all its variations, i.e.starting from the original model with plane strain elasticity, and finishing with the model that captures fullelastic interactions together with the curved fracture tip. For all cases, it is assumed that there is no leak-offand that a point source is used to represent the wellbore. To effectively deal with the moving boundaryproblem, a scaled x coordinate is introduced, ξ = x/l(t) (0 6 ξ 6 1), where l(t) is a half-length of thefracture. In this case, equation (6) can be rewritten as∂w¯∂t −ξldldt∂w¯∂ξ +1l∂q¯x∂ξ =Q0Hlδ(ξ). (43)where the averaged flux is given byq¯x = −1Hlµ′∂p∂ξ∫ h/2−h/2w3 dz. (44)Equations (43) and (44) are common for all methods, while the differences come from the procedures forcalculating w, h and p based on the values of w¯. For the purpose of the numerical computations, boththe spatial variable ξ and the temporal variable t are discretized. Uniform discretization in ξ is used, i.e.ξi = i/Nξ, i = 0...Nξ, ∆ξ = 1/Nξ, where Nξ is number of spatial discretization points. At the same time,nonuniform temporal discretization is used to address the fact that fractures tend to grow rapidly at earlierstages requiring a smaller time step, and then slow down significantly later in the simulation allowing forcoarser time increments.To achieve the nonuniform discretization of the temporal variable that is consistent with the features ofthe problem, the time step is taken so that tj+1−tj ∝ tγj , where j = 0...Nt and γ is some power that isbetween 0 and 1. Here Nt is number of time steps and γ is chosen based on the regime of propagation. Thechoice of power law-type behaviour is based on the fact that the fracture length and width are known to obeypower-law time dependence in some significant regimes of propagation. By using w¯ji = w¯(ξi, tj) to denote23Computation of h Computation of w ConditionInner region Invert (5) Use (3) w¯ >√piH2KIcE′Tip region h = H w = 2E′√H2−4z2 p w¯ 6√piH2KIcE′Table 1: Calculation of fracture height h and fracture width w for different parameters in the original P3D model.the approximation for the effective fracture width and q¯jx,i±12= q¯x(ξi± ∆ξ2 , tj) to denote the approximationfor the fluxes, equation (43) can be discretized at point ξi asw¯j+1i − w¯jitj+1 − tj−ξi+ 12V (tj+1)l(tj+1)w¯j+1i+1 − w¯j+1i−12∆ξ +1l(tj+1)q¯j+1x,i+ 12− q¯j+1x,i−12∆ξ =Q0Hl(tj+1)δi0∆ξ , (45)where V (tj+1) = dl/dt is the velocity of propagation of the fracture, while δi0 is Kronecker delta that ensuresthat the source belongs exclusively to the element with index 0 located at the origin. Equation (45) featuresa backward time difference, in which case its solution utilizes an appropriate iterative scheme at each timestep to solve the nonlinear system of algebraic equations (45). In spite of this complexity, the L-stablebackward Euler scheme allows us to avoid stability issues and to use an arbitrary time step, see [15] inwhich a similar approach is used. Two boundary conditions are used, namely, w¯j+1Nξ = 0, and q¯x|ξ=1 = 0. Inaddition, a symmetry condition q¯j+1x,−12= −q¯j+1x,12is utilized. For completeness, it is noted that the fracturelength is updated using l(tj+1) = l(tj)+V (tj+1)(tj+1−tj). The difference in implementation of various P3Dmodels lies in computation of V (tj+1) and the fluxes q¯j+1x,i±12as a function of w¯.Original P3D model.. To complete the scheme for the numerical implementation of the original P3D model,the fluxes need to be prescribed. To do so, firstly the fracture height needs to be evaluated. There aretwo cases. First, when w¯ >√piH/2KIc/E′, equation (5) is inverted to find h. In situations when w¯ 6√piH/2KIc/E′, a constant height h = H is used. The fracture opening, appearing in (44) is taken from (3)for the first case, while if h = H, the following elliptic opening is chosen w = 2(E′)−1√H2−4z2 p. Thisprocedure is summarized in Table 1. Using this aforementioned procedure, the integral in (44) is evaluatednumerically using the trapezoidal rule for each spatial discretization point. To evaluate the pressure gradient,it is noted that in the case h > H, equations (2) and (5) can be used to compute∂p∂ξ =dpdhdhdw¯∂w¯∂ξ = Y (h)∂w¯∂ξ , (46)whereY (h) = HE′pih24∆σh1/2 −√2piKIc√h2H2 − 13√pi2KIc√h2H2 − 1 + 2∆σh1/2.This formula automatically reduces to ∂p/∂ξ = 2E′/(piH)∂w¯/∂ξ for the case h = H. The fluxes arecomputed using (46) asq¯j+1x,i+ 12= − 1Hl(tj+1)µ′Y(h(ξi+ 12)) w¯j+1i+1 − w¯j+1i∆ξNz∑k=0w3(l(tj+1)ξi+ 12, zk) ∆z, (47)where the sum appears from the numerical integration. Here the z coordinate is discretized as zk = −h2 +∆zk,k = 0...Nz, where ∆z = h/Nz, so that different element sizes are used for each cross-section (since h varies24with ξ). The quantities that need to be evaluated at the mid points (i.e. at ξi+ 12) are treated using linearinterpolation. To find the fracture velocity, equation (43) is integrated over the tip element, which, togetherwith the zero flux boundary condition, allows us to findV (tj+1) =q¯j+1x,Nξ−12w¯j+1i+ 12(1− 12∆ξ). (48)Using (47) and (48), the resulting system of nonlinear algebraic equations (45) is solved iteratively for eachtime step for w¯j+1i , i = 0...Nξ. Once the effective width, w¯, is calculated, formulas (5) and (3) allow us tofind fracture footprint and opening.Original P3D model with the correction for viscous height growth.. The implementation of the correctionfor the viscous resistance does not require many modifications. First, equation (26) can be rewritten as∆Kj+1Ic,i = C1/6µ′1/3E′2/3(12hj+1i − h˜jitj+1 − tj)1/3(hj+1i )1/6(1 + C2H∆σ(hj+1i )1/2∆Kj+1Ic,i)−1/6. (49)Here ∆Kj+1Ic,i denotes the correction to the fracture toughness at point ξi, while hj+1i are calculated from (5),where the fracture toughness is replaced by KIc+∆Kj+1Ic,i . Since the moving mesh is used, h˜ji is not justhji , because the calculated velocity is not exactly vertical in this case. Since the physical position x hasto be fixed during the calculation of the vertical velocity, one needs to take h˜ji = I[hji ](l(tj+1)ξi), wherethe latter means the interpolation of hji (i = 0...Nξ) at the point l(tj+1)ξi. The values of the constantsin (49), which are used for the computations, are C = 0.175 and C2 = 0.5. The solution of the nonlinearequation (49) together with hj+1i coming from (5), obtained using Newton’s method, produces both hj+1iand ∆Kj+1Ic,i . Finally, the fracture width, w, is calculated by replacing KIc with KIc+∆Kj+1Ic,i in (3). Withthese modifications in the calculation of h and w, the solution is obtained by solving equation (45) togetherwith (47) and (48).P3D with non-local elasticity.. The implementation of non-local elasticity requires more significant modi-fications. First, since the pressure is evaluated at z= 0, there should be an element located at this point,which requires Nz to be even (so that the total number of elements Nz+1 is odd). With this meshing,equation (27) can be discretized aspji =(C˜jikm + C˜ji(−k)m)W jkmnw¯jn = Cjinw¯jn, (50)where index “−k” means evaluation at negative ξk, which allows us to build symmetry into the influencecoefficients Cjin (to avoid capturing element with k = 0 twice, C˜ji(−k)m is added only for k > 1),C˜jikm = −E′l(tj)8pi∫ ξk+ 12 ∆ξξk−12 ∆ξ∫ zm+ 12 ∆zzm−12 ∆zdz′ dξ′(l(tj)2(ξ′−ξi)2 + z′2)3/2= −E′8pi[√l(tj)2(ξi − ξ′)2 + z′2l(tj)(ξi − ξ′)z′]ξ′=ξk+ 12 ∆ξ,z′=zm+ 12 ∆zξ′=ξk− 12 ∆ξ,z′=zm− 12 ∆z, (51)andW jkmn =wj(l(tj)ξk, zm)w¯jkδkn. (52)25Here wj(lξk, zm) in (52) denotes the evaluation of the fracture opening according to (3) (or ellipse) atx= l(tj)ξk and z= zm and δkn is the Kronecker delta. Relations (51) and (52) permit the computation ofthe elasticity matrix Cjin featured in (50). Note, however, that this elasticity matrix is a nonlinear functionof w¯jn, where this nonlinear dependence comes implicitly through the variation of h and explicitly from (52).The fluxes are calculated using the discretization of (44) asq¯j+1x,i+ 12= − 1Hl(tj+1)µ′pj+1i+1 − pj+1i∆ξNz∑k=0w3(l(tj+1)ξi+ 12, zk) ∆z, (53)where the values of the pressure are now calculated using (50).Unfortunately, calculating the pressure using the hypersingular integral (27) leads to poor accuracy forthe tip element. To deal with this, the asymptotic solution can be built into the tip element, as was donein [15]. To do so, it is assumed that the tip element follows the asymptotic solution (obtained for KGDfracture), and the fracture propagation velocity is calculated accordingly. In particular, if the fracturepropagates in the viscous regime, the velocity of propagation is calculated asV (tj+1) =E′µ′w3(l(tj+1)(1−∆ξ), 0)β3m(l(tj+1)∆ξ)2, (54)where βm = 21/3 ·35/6 and w(l(tj+1)(1−∆ξ), 0) is the fracture opening at z = 0 for the penultimate node. Iffracture propagates in the toughness regime, the velocity is adjusted in a way that the toughness propagationcriterion holds for the last element, i.e.w(l(tj+1)(1−∆ξ), 0) =√32piKIcE′√l(tj+1)∆ξ. (55)In situations when the fracture propagates in the intermediate regime, one of the asymptotic solutions isused. The criterion for switching between the two solutions is the following: at the beginning, the velocity ofpropagation is calculated for both viscous and toughness regimes, and then the minimum velocity is chosen.Physically, there are two resistance mechanisms - viscous and toughness. Taking the minimum velocityimplies that only the dominant resistance is considered. In situations when both mechanisms contribute,the velocity, and consequently the fracture length, are overestimated. Since fractures typically evolve fromthe viscous regime to the toughness regime (assuming no leak-off), the error introduced during the transitionzone is mitigated later when the toughness regime dominates. Moreover, comparison with ILSA simulationsshows that even during the intermediate regime, the length was overestimated only by approximately 5%.To compensate for the introduction of the extra condition comprising the formula for the velocity offracture propagation, and to ensure a zero flux boundary condition at the tip, the pressure at the tipelement needs to be treated as an unknown, as was also done in [15]. Knowing the velocity of propagation,equation (48) can be used to find the tip pressure, pj+1Nξ , aspj+1Nξ = pj+1Nξ−1 −Hl(tj+1)µ′∆ξV (tj+1)w¯j+1Nξ−12(1− 12∆ξ)∑Nzk=0 w3(l(tj+1)ξNξ− 12, zk) ∆z. (56)Finally, the algorithm with full elasticity consists of solving (45), where the fluxes are computed us-ing (53), the pressure (for all nodes, except the tip element) is calculated using (50)–(52), the fracturepropagation velocity is taken such that either (54) or (55) is satisfied depending on the regime of propaga-tion, while the pressure at the tip element is taken from (56). It should also be noted that the algorithmfor calculating the apparent toughness, based on (49), is included in its original form, as it is independentof the type of elasticity equation that is used for calculations.26Computation of h Computation of w Computation of ∆KIc ConditionInner region Invert (5) Use (3) Use (49) w¯ >√piH2KIc + ∆KIcE′Tip region Use (40a) Use (40b) Use (41) w¯ 6√piH2KIc + ∆KIcE′Table 2: Calculation of the fracture height h, the fracture width w, and the apparent toughness ∆KIc for different parametersfor the P3D model with curved tip.P3D with approximate non-local elasticity.. The only difference between the P3D algorithm with an ap-proximate non-local elasticity operator from the one described above, which includes full elasticity, lies incalculation of elasticity matrix Cjin in (50). By approximating both h and w¯ as piece-wise constant functions,equation (31) can be reduced topji =[(Cj1,ik + Cj1,i(−k))W j1,kn +(Cj2,ik + Cj2,i(−k))W j2,kn]w¯jn = Cˆjinw¯jn, (57)where Cˆjin is an approximation to Cjin, index “−k” denotes the evaluation of Cj1,ik at −ξk in order to capturesymmetry (as before, to avoid capturing element with k = 0 twice, terms with “−k” are added only fork > 1), andCj1,ik =E′2piG(2l(tj)(ξ′−ξi)h(l(tj)ξ′, tj))∣∣∣∣ξ′=ξk+ 12 ∆ξξ′=ξk− 12 ∆ξ, W j1,kn =w1(l(tj)ξk, tj)h(l(tj)ξk, tj)w¯kδkn,Cj2,ik =E′2piG( 2l(tj)(ξ′−ξi)h∆σ(l(tj)ξ′, tj))∣∣∣∣ξ′=ξk+ 12 ∆ξξ′=ξk− 12 ∆ξ, W j2,kn =w2(l(tj)ξk, tj)h∆σ(l(tj)ξk, tj)w¯kδkn. (58)The kernel G is defined in (29), w1 and w2 as a functions of h are given in (30), while h∆σ is defined in (33). Insituations when h=H, w1(l(tj)ξk, tj)=4w¯jk/pi and w2 =0. Evaluation of w1(l(tj)ξk, tj), w2(l(tj)ξk, tj), andh∆σ(l(tj)ξk, tj) all require the evaluation of h(l(tj)ξk, tj), which in turn depends on w¯(l(tj)ξk, tj) = w¯jk. Inthis case, similar to the full elasticity matrix, the approximate elasticity matrix, Cˆjin, is a nonlinear functionof w¯jk. Note that since a piece-wise constant approximation for h is used (because w¯ is also piece-wiseconstant), spatial derivatives of h and h∆σ in (31) do not influence the result.The P3D model with the approximate full elasticity requires the solution of the same equations as themodel for the full elasticity, except that the subroutine for evaluating the fluid pressure (50)–(52) is replacedby one corresponding to (57)–(58).Curved fracture tip.. To introduce a curved fracture tip, the algorithm needs to be updated in several places.The computation of the fracture height, h, and the width in the tip region is now performed using (40).Table 2 summarizes the procedure and replaces the corresponding procedure for the original P3D model thatis summarized in Table 1. In addition to the different formulas for h and w, the apparent toughness ∆KIc isevaluated differently. In particular, equation (41) is used to find the fracture toughness in the tip region. Theevaluation of the fluid pressure using (57) and (58) has two modifications. First, since h < H in the tip regionand only one ellipse is used to approximate the fracture width, w1(l(tj)ξk, tj) = 4H/(pih(l(tj)ξk, tj))w¯jk andw2 = 0. Unfortunately, since ∂h/∂x in (31) is singular at the tip element (if a curved fracture tip is used),the piece-wise constant approximation for h precludes an accurate evaluation of the elasticity integral. Todeal with this problem, we seek a correction that makes the elasticity integral accurate. One way to estimatethe accuracy of the elasticity equation is to consider a radial fracture propagating in the toughness regime.Since, in this case, every vertical cross-section is an ellipse the radial fracture can be accurately representedby the P3D model with a curved fracture tip. If one uses (35) and (36), and evaluates the fluid pressureusing (57) and (58), then the resultant pressure should be constant (since the fracture propagates in the270 0.2 0.4 0.6 0.8 101234x/lp/p0  No elasticity correctionWith elasticity correction0.98 0.99 10.9811.021.04x/lp/p 0  Fracture footprintExtra elementExtra fracturevolumeFigure 10: Variation of the normalized pressure versus the normalized distance for a radial fracture propagating in the toughnessregime, calculated numerically using (57) with and without the correction (left picture). The right picture schematically showsthe location of the extra element that is used for the correction.toughness regime). All deviations from the constant in this case are due to the inaccuracy caused by theevaluation of the elasticity equation. The results of implementing this procedure show that the calculatedpressure deviates from a constant near the tip and so a correction is required, see the dashed line in Fig. 10.Note that there is a fracture volume at the tip that is not accounted for in the pressure calculations. Inother words, since the fracture opening is zero at the tip, the opening of the corresponding element is zeroand so this element does not contribute to the pressure. To introduce a correction, one may account forthe extra volume by introducing a non-zero tip element for pressure calculations, see Fig. 10. By assumingthat the fracture opening follows the toughness asymptotic solution for the tip element, it can be shownthat the volume of the extra element is 1/8 of that of penultimate element. Then, to determine both h andw1 for the extra element, its height is taken as hNξ = hNξ−1/r, where r = 1.825 is a constant. The valueof r is found by minimizing the error between the numerically calculated pressure and the constant (whichcorresponds to the exact solution). Fig. 10 shows that the fluid pressure inside the radial fracture, that iscalculated using the extra tip element, is accurate everywhere to within 1% error.To summarize, the numerical scheme for the P3D model with a curved tip with the approximate non-localelasticity and correction for viscous height growth consists of solving (45), where the fluxes are computedusing (53), the pressure (for all nodes, except the tip element) is calculated using (57) and (58), the fracturepropagation velocity is taken such that either (54) or (55) is satisfied depending on the regime of propagation,while the pressure at the tip element is taken from (56). Procedures for calculating the fracture height, width,as well as the apparent toughness are summarized in Table 2.Initial condition.. Even though the initial condition does not affect the solution in the long term, it isimportant to specify an appropriate initial condition, since, otherwise, the numerical algorithm may notconverge at the first time step. For both viscous and toughness regimes of propagation it is assumed thatat a given small time instant t1 the hydraulic fracture is radial and the fracture opening corresponds to thesolution for “dry” crack, i.e. given bywrad =4(KIc+∆KIc)√pilE′√l2 − x2 − z2, (59)where the apparent toughness ∆KIc is calculated based on (41). By calculating the volume of the fracturefrom (59) and equating it to the pumped volume (no leak-off), we obtainQ0t =83√pi(KIc+C1/63 µ′1/3E′2/3l1/6dldt1/3)l5/2E′ . (60)28The latter differential equation (60) is solved numerically to find l(t1), dl(t1)/dt, and ∆KIc(t1). Note thatthe differential equation reduces to an algebraic equation in the toughness regime, in which casel(t1) =(3E′Q0t18√piKIc)2/5, dldt (t1) =25( 3E′Q08√piKIc)2/5t−3/51 , ∆KIc(t1) = 0.Having l(t1), dl(t1)/dt, and ∆KIc(t1), the effective width is calculated from (59) asw¯(t1) =2√pi(KIc+∆KIc(t1))l(t1)3/2HE′ (1− ξ2), (61)while, since the fracture is radial, its height is given byh(t1) = 2l(t1)√1− ξ2. (62)This procedure for calculating the initial condition is exact if the fracture propagates in the toughness regime,while it gives an approximate initial condition if the fracture propagates in the viscous regime. One mayalso consider using the solution for a radial fracture propagating in the viscous regime [31]. However, sincethe P3D model can represent such a fracture only approximately, it is sufficient to use an approximation forthe initial condition as well.References[1] U. Frank, N. Barkley, Remediation of low permeability subsurface formations by fracturing enhancements of soil vaporextraction, J. Hazard. Mater. 40 (2005) 191–201.[2] A. Abou-Sayed, D. Andrews, I. Buhidma, Evaluation of oily waste injection below the permafrost in prudhoe bay field,in: Proceedings of the California Regional Meetings, Bakersfield, CA, Society of Petroleum Engineers, Richardson, TX,1989, pp. 129–142.[3] R. Jeffrey, K. Mills, Hydraulic fracturing applied to inducing longwall coal mine goaf falls, in: Pacific Rocks 2000, Balkema,Rotterdam, 2000, pp. 423–430.[4] M. Economides, K. Nolte (Eds.), Reservoir Stimulation, 3rd Edition, John Wiley & Sons, Chichester, UK, 2000.[5] S. Khristianovic, Y. Zheltov, Formation of vertical fractures by means of highly viscous fluids, in: Proc. 4th WorldPetroleum Congress, Vol. 2, 1955, pp. 579 – 586.[6] T. Perkins, L. Kern, Widths of hydraulic fractures, in: J. Pet. Tech. Trans. AIME, 1961, pp. 937–949.[7] R. Nordgren, Propagation of vertical hydraulic fractures, in: J. Pet. Technol., 1972, pp. 306–314.[8] I. P. ID, H. Carroll, Three-dimensional hydraulic fracture propagation in the presence of stress variation, in: Proceedingsof the SPE/DOE/GRI unconventional gas recovery symposium, SPE/DOE 10849, 1982, pp. 870–878.[9] I. P. ID, H. Carroll, Numerical solution for height and elongated hydraulic fractures, in: Proceedings of the SPE/DOElow permeability reservoir symposium, Denver. SPE 11627, 1983, pp. 249–256.[10] A. Settari, M. Cleary, Development and testing of a pseudo-three-dimensional model of hydraulic fracture geometry (p3dh),in: Proceedings of the 6th SPE symposium on reservoir simulation of the Society of Petroleum Engineers, SPE 10505,1982, pp. 185–214.[11] J. McLennan, J. Picardy, Pseudo-three-dimensional fracture growth modeling, in: 26th US Symposium on Rock Mechanics,1985, pp. 323–331.[12] Y. Xiujuan, W. Tongtao, Y. Xiangzhen, W. Xin, A pseudo-3D model with 2D flow of hydraulic fracture propagation, in:Thin Interbedded Sandstone Reservoir. International Society for Rock Mechanics., 2010, pp. 429–433.[13] S. H. A. T. Lee, J. Lee, Three-dimensional modeling of hydraulic fractures in layered media: Part I–finite elementformulations, J. Energy Resour. Technol. 112 (1990) 1–9.[14] L. Vandamme, J. Curran, A three-dimensional hydraulic fracturing simulator, Int. J. Numer. Meth. Eng. 28 (1989) 909–927.[15] A. Peirce, E. Detournay, An implicit level set method for modeling hydraulically driven fractures, Comput. Methods Appl.Mech. Engrg. 197 (2008) 2858–2885.[16] J. Adachi, E. Siebrits, A. Peirce, J. Desroches, Computer simulation of hydraulic fractures, Int. J. Rock Mech. Min. Sci.44 (2007) 739–757.[17] L. Germanovich, L. R. D. Astakhov, J. Shlyopobersky, M. Mayerhofer, Hydraulic fracture with multiple segments II:Modeling, Int. J. Rock Mech. Min. Sci. 34 (1997) 472.[18] A. Peirce, A. Bunger, Interference fracturing: Non-uniform distributions of perforation clusters that promote simultaneousgrowth of multiple hydraulic fractures, SPE 172500.[19] J. A. L. Napier, E. Detournay, Propagation of non-planar pressurized cracks from a borehole, in: Research and Applicationsin Structural Engineering, Mechanics and Computation, 2013, pp. 597–602.29[20] M. Yang, Hydraulic fracture production optimization with a pseudo-3D model in multilayered lithology, in: Society ofPetroleum Engineers, SPE 149833, 2012.[21] W. Kordziel, W. Rowe, V. Dolan, S. Ritger, A case study of integrating well-logs and a pseudo 3D multi-layer frac modelto optimize exploitation of tight lenticular gas sands, in: Society of Petroleum Engineers, SPE 36866, 1996.[22] R. Wu, O. Kresse, X. Weng, C.-E. Cohen, H. Gu, Modeling of interaction of hydraulic fractures in complex fracturenetworks, in: Society of Petroleum Engineers, SPE 152052, 2012.[23] X. Weng, Incorporation of 2D fluid flow into a pseudo-3D hydraulic fracturing simulator, in: Society of Petroleum Engi-neers, 1992, pp. 331–337.[24] E. Sarvaramini, D. Garagash, Pressurization of a PKN Fracture in Permeable Rock during Injection of a Low ViscosityFluid, Tech Publishers, 2013, Ch. 30, pp. 629–640.[25] J. Adachi, A. Peirce, Asymptotic analysis of an elasticity equation for a finger-like hydraulic fracture, J. Elasticity 90(2008) 43–69.[26] J. I. Adachi, E. Detournay, A. P. Peirce, An analysis of classical pseudo-3D model for hydraulic fracture with equilibriumheight growth across stress barriers, Int. J. of Rock Mech. and Min. Sci. 47 (2010) 625–639.[27] E. Carter, Optimum fluid characteristics for fracture extension, in: Howard GC, Fast CR, editors. Drilling and productionpractices, 1957, pp. 261–270.[28] E. Detournay, Propagation regimes of fluid-driven fractures in impermeable rocks, Int. J. Geomech. 4 (2004) 35–45.[29] A. Peirce, Modeling multi-scale processes in hydraulic fracture propagation using the implicit level set algorithm, (accepted)Comp. Meth. in Appl. Mech. and Eng.[30] R. Jeffrey, A. Bunger, A detailed comparison of experimental and numerical data on hydraulic fracture height growththrough stress contrasts, in: Society of Petroleum Engineers, SPE 106030, 2007.[31] A. Savitski, E. Detournay, Similarity solution of a penny-shaped fluid-driven fracture in a zero-toughness linear elasticsolid, CR. Acad. Sci. II B 329 (2001) 255–262.30


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items