UBC Faculty Research and Publications

Proppant transport in hydraulic fracturing : crack tip screen-out in KGD and P3D models Dontsov, E. V.; Peirce, Anthony Feb 27, 2014

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

Item Metadata


52383-Dontsov_Peirce_Proppant_Transport_2014_08_13.pdf [ 3.88MB ]
JSON: 52383-1.0079344.json
JSON-LD: 52383-1.0079344-ld.json
RDF/XML (Pretty): 52383-1.0079344-rdf.xml
RDF/JSON: 52383-1.0079344-rdf.json
Turtle: 52383-1.0079344-turtle.txt
N-Triples: 52383-1.0079344-rdf-ntriples.txt
Original Record: 52383-1.0079344-source.json
Full Text

Full Text

Proppant transport in hydraulic fracturing: Crack tip screen-out inKGD and P3D modelsE.V. Dontsov and A.P. PeirceAugust 13, 2014AbstractThe aim of this study is to develop a model for proppant transport in hydraulic fractures capable ofcapturing both gravitational settling and tip screen-out effects, while prohibiting the particles fromreaching the crack tips by imposing a width restriction based on the particle size. First, the equationsthat govern the propagation of hydraulic fractures and the proppant transport inside them are formulated.They are based on the solution for the steady flow of a viscous fluid, mixed with spherical particles, in achannel, which is obtained assuming an empirical constitutive model. This proppant transport model isapplied to two fracture geometries – Khristianovich-Zheltov-Geertsma-De Klerk (KGD) and pseudo-3D(P3D). Numerical simulations show that the proposed method makes it possible to capture proppantplug formation and growth, as well as the gravitational settling for both geometries. A dimensionlessparameter, whose magnitude reflects the intensity of the settling, is introduced for the P3D fracture.1 IntroductionHydraulic fracturing is a process whereby the fluid pressure due to fluid injection into a crack is the drivingforce for the fracture opening and propagation. Among the multiple uses of hydraulic fracturing, such asaccelerating the waste remediation process [1], waste disposal [2], or preconditioning in rock mining [3], oil andgas reservoir stimulation [4] stands out as one of the most common applications. Recognizing the significanceof hydraulic fracturing, many studies have been devoted to the modelling and numerical simulation of thisphenomenon. Starting from the work of [5], further examples of the analytical modelling can be foundin [6, 7, 8, 9], where the near tip solutions and regimes of propagation are studied, while reviews of theexisting numerical approaches aiming to predict hydraulic fracture propagation are given in [10, 11].The problem of hydraulic fracturing is challenging to analyze due to a variety of physical processes that areinvolved in the problem, such as fluid flow inside the fracture, fluid leak-off to the surrounding rock, the rockfracturing due to crack propagation, and, in some cases, elastic interaction with natural fractures or otherhydraulic fractures. Moreover, the fracturing fluid can be non-Newtonian, and its properties may vary withtime and temperature. To effectively model the process, however, many assumptions are typically made. Forinstance, the fluid is assumed to be Newtonian, the flow is assumed to be laminar, the behaviour of the rockis taken as linear elastic, poroelastic effects are typically neglected, the geometry of the fracture is greatlysimplified to one-dimensional, radial, or planar etc. Even with these simplifications, the phenomenon ofhydraulic fracturing is difficult to model, as it requires the solution of a nonlinear problem with a singularity,in which the nonlinearity comes form the lubrication equation and the singularity typically appears atthe crack tip. This study aims to add an additional aspect to the problem, namely, the movement ofproppant within the fracture. Typically, proppant is used to prevent the fracture from closing once thewell is depressurized. In this case, modelling the fracture propagation driven only by a viscous fluid isnot sufficient, since the proppant, blended with the fracturing fluid alters the properties of the fracturingfluid. Incorporating the effects associated with the presence of particles poses an additional challengingproblem, which is addressed in this study. As mentioned in [10], in hydraulic fracturing problems, the slurryis typically modelled as a Newtonian fluid with the effective viscosity given by an empirical function ofproppant content. In addition, a uniform particle distribution across the fracture is assumed and the slip1velocity only due to gravity is considered. In contrast, the current study utilizes an approach developedin [12], where the governing equations for the slurry and the proppant transport are obtained based on theempirical constitutive model for the slurry introduced by [13]. This model accounts for the non-uniformparticle distribution across the channel, slip velocity induced by the slurry flow, and captures the transitionfrom Poiseuille’s flow for small particle concentrations to Darcy’s law for nearly maximum proppant content.It is also important to highlight a two phase model proposed by [14], and work by [15], where proppantplug formation near the crack tip is studied. In the latter study, the one-dimensional problem of KGDfracture propagation is considered, and the problem is tackled using a double moving coordinate system.One coordinate is scaled by the length of the crack and another by the distance from the inlet to theproppant plug, in which case the boundaries of the plug are tracked automatically based on the calculateddistances. While this approach works well for a 1D geometry, its generalization to 2D fractures seems tobe tremendously difficult since two moving boundaries cannot easily be resolved using scaling. Some of thestudies, however, investigate just the flow of the slurry, i.e. the mixture of fluid with the particles, and donot apply it to the fracture propagation problem. One example is a study by [16], in which the granulartemperature is used to account for micro-level particle movements. To address the problem, this study aimsto develop a model for hydraulic fracturing by a slurry, which accounts for the mechanics of the slurry, while,at the same time, is sufficiently simple that it can be implemented into a hydraulic fracturing simulator.The paper is organized as follows. First, the governing equations for the slurry flow and proppant trans-port inside hydraulic fractures, obtained in [12], are in summarized in Section 2. Then, in Sections 3 and 4,the governing equations for the slurry and proppant transport are embedded into the fracture propagationproblems (respectively KGD and P3D) and the complete problems are solved numerically.2 BackgroundThis section aims to summarize background information that is necessary to develop a computational schemefor proppant transport inside hydraulic fractures. The approach is based on the slurry flow solution in thechannel, developed in [12]. This solution is, in turn, based on the empirical constitutive model for themixture of a Newtonian fluid and spherical particles introduced in [13]. Fig. 1 shows the schematics of thefracture and the associated coordinate system, where x is the coordinate along the fracture in the horizontaldirection, z is the vertical axis (it is assumed that the fracture is contained in the vertical plane), while y isthe coordinate across the fracture. As shown in [12], the balance equations for the slurry and proppant can2axyz gezvxzwFigure 1: Schematics of the hydraulic fracture (left) and the slurry flow inside it (right).be written as∂w∂t +∇·qs + g = 0,∂wφ¯∂t +∇·qp = 0, (1)2where w is the width of the fracture, φ¯ = 〈φ〉/φm is the normalized proppant concentration averaged overfracture width, i.e. in the y direction, φm = 0.585 is the maximum allowed proppant concentration, grepresents leak-off, while qs and qp denote respectively the slurry and proppant fluxes. Note that the fluxeshave two components, namely x and z, and consequently, ∇ = (∂/∂x, ∂/∂z) in (1). The expressions for thefluxes areqs = − w312µf Qˆs(φ¯, wa)∇pˆ,qp = B(wa)Qˆp(φ¯, wa)qs −B(wa) a2w12µf (ρp−ρf)gezGˆp(φ¯, wa), (2)where µf is the clear fluid viscosity, pˆ is the fluid pressure corrected by hydrostatic pressure, ρp−ρf is thedifference between proppant and fluid mass densities, g is the gravitational acceleration, a is the particleradius, B is a, so-called, blocking function, while the functions Qˆs, Qˆp and Gˆp come from the slurry flowsolution [12].The blocking function B is introduced to capture proppant bridging that occurs when the fracture widthis on the order of several particle diameters. For the purpose of calculations, the blocking function is takenasB(wa)= 12H( w2a−N)H(wB − w2a)[1 + cos(piwB − w2a)]+H(w − wB2a), (3)where N represents “several” particle diameters, H denotes Heaviside step function, while wB = 2a(N+1),which provides a continuous vanishing of the function and helps in the numerical implementation. N =3 ischosen for all examples considered in this paper.0 0.5φ¯Gˆp0 0.5 100.51φ¯Qˆp 0.9 0.95φ¯Gˆpwa = 6wa = 10wa = 100wa = 6wa = 10wa = 1000 0.5 100.51φ¯Qswa = 6wa = 10wa = 100QˆsFigure 2: Variation of the functions Qˆs, Qˆp and Gˆp versus normalized proppant concentration φ¯ for differentvalues of the parameter w/a.Functions Qˆs, Qˆp and Gˆp can be expressed in a simpler form asQˆs(φ¯, wa)= Qs(φ¯) + a2w2 φ¯D,Qˆp(φ¯, wa)= w2Qp(φ¯)w2Qs(φ¯) + a2φ¯D , (4)Gˆp(φ¯, wa)= Gp(φ¯)− w2Gs(φ¯)Qp(φ¯)w2Qs(φ¯) + a2φ¯D ,where Qs, Qp, Gs and Gp are functions of φ¯ only and are calculated numerically, D = 8(1−φm)α¯/3φm isrelated to the permeability of the packed particles, α¯ = 4.1, see [12]. As an illustration, Fig. 2 plots thefunctions Qˆs, Qˆp and Gˆp versus normalized proppant concentration φ¯ for different values of the parameterw/a. Function Qs represents the reciprocal of the effective slurry velocity, while the term with D in Qˆsaccounts for a Darcy’s law. So, the slurry flux in (2) is able to capture the transition from Poiseuille’s flow(with effective viscosity) to Darcy’s filtration flow as the concentration reaches nearly the maximum value.The function Qˆp in the proppant flux in (2) describes the proppant motion caused by the slurry flow, whileGˆp captures gravitational settling of the particles.3It is the difference between the fluid pressure and hydrostatic pressure, denoted by pˆ, that drives heslurry, see (2). According to [12], pˆ is given bypˆ = p− ph = p+ ρfgz + (ρp−ρf )gφmφ¯z + (ρp−ρf )g∫ z0a2Gs(φ¯)w2Qs(φ¯) + a2φ¯D dz′, (5)where p is the fluid pressure and ph is the hydrostatic pressure. The hydrostatic pressure is the drivingforce for buoyancy-driven fractures [17], which is typically neglected for other hydraulic fracturing problemsfor simplicity. For this reason, all calculations in this paper neglect hydrostatic pressure, which effectivelyreplaces pˆ by the fluid pressure p in (2).3 Numerical solution for a KGD fracture3.1 Problem formulationTo better understand the effects associated with the presence of particles (or proppant) on hydraulic fracturepropagation, it is instructive to consider the simplest one-dimensional case of a KGD fracture in the presenceof stress barriers. To this end it is assumed that the fracture lies along z axis and occupies interval (−l1, l2)(see Fig. 3), where l1 and l2 are the extensions of the fracture respectively in the negative and positive zdirections. By using Carter’s leak-off model [18] and adding the source terms, the governing equations forfluid flow within the fracture including proppant can be deduced from (1) as∂w∂t +∂qsz∂z +C ′√t− t0(z)= Q0δ(z),∂wφ¯∂t +∂qpz∂z = φ¯0Q0δ(z), (6)where −l1 6 z 6 l2, andqsz = −w3µ′ Qˆs(φ¯)∂p∂z ,qpz = B(wa)Qˆp(φ¯, wa)qsz −B(wa)a2wµ′ (ρp−ρf)gGˆp(φ¯, wa). (7)Here µ′ = 12µf is the scaled viscosity, C ′=2CL (CL - Carter’s leak-off coefficient), Q0 is the injection volumeof the slurry per unit time, while φ¯0 is the normalized volume fraction of proppant at the source. To closethe system of equations (6), one needs to add the elasticity equation (see e.g. [19]), which can be written asp−∆σ(z) = − E′4pi=∫ l2−l1w ds(s− z)2 , (8)where E′ = E/(1−ν2) is the plane-strain Young’s modulus, the integral is understood in the sense of aHadamard finite part, while ∆σ(z) is an additional confining stress coming from the stress barriers. Theboundary and propagation conditions at the tips areqsz|z=−l1 = 0, qpz |z=−l1 = 0, w →K ′E′ (z + l1)1/2, z → −l1,qsz|z=l2 = 0, qpz |z=l2 = 0, w →K ′E′ (l2 − z)1/2, z → l2, (9)where K ′ = 8KIc/√2pi is the scaled fracture toughness.4Q0Plane strain¯0Q0l2l11 0 1V1 V2Total fluxProppant fluxz⇣1 ⇣2Figure 3: Schematics of asymmetric KGD fracture.3.2 Numerical algorithmThe problem under consideration is split into two steps: i) solve for the propagation of the fracture, i.e.coupling (6a), (8) and (9) and ii) solve for the proppant transport, i.e. (6b) with the boundary conditionsgiven in (9). In other words, at each time step, first (6a), (8) and (9) are solved to update the fracturewidth profile and length, and then (6b) is solved to find the new proppant concentration distribution overthe fracture length. This subsection is aimed to cover both steps in the procedure.To facilitate the numerical solution of the moving-boundary problem under consideration, a double mov-ing mesh is introduced. In this case, negative and positive components of the x coordinate are normalizedrespectively by l1 and l2, so thatζ1 =zl1, −1 6 ζ1 6 0, ζ2 =zl2, 0 6 ζ2 6 1, (10)and∂(·)∂t∣∣∣z= ∂(·)∂t∣∣∣ζ1− V1ζ1l1∂(·)∂ζ1∣∣∣t, ∂(·)∂z∣∣∣t= 1l1∂(·)∂ζ1∣∣∣t, z 6 0,∂(·)∂t∣∣∣z= ∂(·)∂t∣∣∣ζ2− V2ζ2l2∂(·)∂ζ2∣∣∣t, ∂(·)∂z∣∣∣t= 1l2∂(·)∂ζ2∣∣∣t, z > 0, (11)where V1 = dl1/dt and V2 = dl2/dt are the velocities of crack propagation in the negative and positive zdirections, as indicated in the Fig. 3. By substituting (11) into (6), one can write∂w∂t −V1,2l1,2ζ1,2∂w∂ζ1,2+ 1l1,2∂qsz∂ζ1,2+ C′√t− t0(l1,2ζ1,2)= Q0l1,2δ(ζ1,2),∂wφ¯∂t +V1,2l1,2wφ¯+ 1l1,2∂q˜pz∂ζ1,2= φ¯0Q0l1,2δ(ζ1,2), (12)whereq˜pz = qpz − V1,2ζ1,2wφ¯,and the indices 1, 2 correspond respectively to the regions z 6 0 and z > 0 (see (11)).Fracture propagation. To formulate the numerical scheme, first, the width is approximated by a piece-wise constant function in the ζ1 and ζ2 domains. In this case, the elasticity equation (8) can be discretizedin obtainp i = Ciwi + pitip + ∆σi, (13)5¯ij¯ij1wij1 wij wij+1⇣ 12⇣Figure 4: Discretization of the width w and normalized particle concentration φ¯.where p i and wi denote respectively pressure and width vectors (defined at all grid points and time ti), pitipis the pressure correction at the tip (defined only at the tip points and time ti),(Ci)jk = −E′l1∆ζ4pi((zj − zk)2 − 14 l21∆ζ2)−1, zk < 0,(Ci)jk = −E′(l1 + l2)∆ζ8pi((zj − 12 l2∆ζ)(zj + 12 l1∆ζ))−1, zk = 0,(Ci)jk = −E′l2∆ζ4pi((zj − zk)2 − 14 l22∆ζ2)−1, zk > 0is the elasticity matrix (depends on time through l1 and l2), while ∆σi is the term that comes from thepresence of stress barriers (again, this depends on time through l1 and l2). Here ∆ζ denotes the mesh sizeassociated with the discretized coordinates ζ1 and ζ2 (the same element size is used for both ζ1 and ζ2),while zj refers to the location of the jth element of wi. The tip pressure is added as an unknown since theaccuracy of pressure at the tip is low due to the singular nature of the kernel in (8). By using the backwardEuler scheme to approximate the time derivative, equations (12a) and (13) can be written in a discretizedform aswi −wi−1ti − ti−1−Biwi +A(wi, φ¯i)Ciwi +A(wi, φ¯i)(p itip + ∆σi)= Si−1/2, (14)where Bi is the matrix that accounts for the “moving mesh terms” coming from the time derivatives in (11),A(wi, φ¯i) approximates the flux divergence term, while Si−1/2 accounts for the source and leak-off terms.Central differences are used to calculate matrix Bi, while A(wi, φ¯i) is calculated using1l1,2(∂qsz∂ζ)ij= JAj+1/2p ij+1 − p ij(l1,2∆ζ)2− JAj−1/2p ij − p ij−1(l1,2∆ζ)2= Ajm p im, (15)whereJAj±1/2 = −wij±1/2µ′((wij±1/2)2Qs(φ¯ij±1/2) + a2φ¯ij±1/2D),and wij±1/2 = 12 (wij±1 +wij), while φ¯ij±1/2 are defined at the midpoints, see Fig. 4. Note that the pressureis defined at the same points as the width, while coefficients EAj±1/2 share the mesh with φ¯ij±1/2. Thediscretized equation (14) approximates the corresponding differential equation inside the domain, and thusdoes not capture the boundary conditions. To find the discretized equations for the boundary nodes, oneneeds to integrate equation (6a) over the tip elements and use the boundary conditions (9). This providestwo equations, which can be written in the general form asp itip = p itip(V1, V2,wi, φ¯i). (16)If the correction for the tip were not necessary, then (16) would allow us to find the unknown tip velocitiesV1 and V2 and close the system of equations (14). However, since the pressure at the tip elements cannotbe computed accurately, one has to impose two additional conditions. One possibility is to assume that thewidth in the tip elements should follow the appropriate asymptotic solution, as used in [11], and use thecorresponding asymptotic formulas for the tip velocity. In this case, two additional equations areV1 = V¯1(wi, C ′,K ′), V2 = V¯2(wi, C ′,K ′), (17)6where functions V¯1 and V¯2 depend on the regime of propagation of the hydraulic fracture. Note that theproppant cannot occupy the near-tip elements due to the “blocking” functions, introduced in (3), for thisreason, “classical” asymptotic solutions can be used. For instance, in the viscous regime, i.e. in situationwhen C ′ = K ′ = 0, one hasV1 =(wi2)3β3m(l1∆ζ)2, V2 =(wiNζ−1)3β3m(l2∆ζ)2(18)where βm = 21/3 ·35/6, while wi2 and wiNζ−1 are the values of the width for the second and the penultimatenodes. Note that wi1 = wiNζ = 0 due to the boundary conditions in (9). More information about asymptoticsolutions for the KGD fracture can be found in [19].Finally, at each time step, equations (14) and (16) are solved iteratively using the appropriate expressionfor the tip velocities (17). Then, the fracture footprint is updated usinglij = li−1j + (ti − ti−1)Vj , j = 1, 2.Proppant transport. As indicated earlier in this section, first, the fracture propagation is determined fora given time step, and then (12b) is solved numerically to update the particle distribution over the fracturelength. Equation (12b) has the form of a conservation law, which is both nonlinear and heterogeneous. Todeal with such a problem, a finite volume method with cell-centered flux functions [20] is used. In particular,the moving mesh term is integrated by parts, and (12b) is discretized as follows(wφ)i − (wφ)i−1ti − ti−1+ V1,2l1,2(wφ)i−1 +DE±1/2q˜pzi−1= Si−1/2p , (19)where Si−1/2p is the source term, matrix D is a finite difference operator, which employes a central differencescheme to approximate the flux derivative. i.e. (DA)j−1/2 = Ai−1j −Ai−1j−1, E±1/2 is a shift operator, whichshifts the flux function evaluation either to the left or to the right by a half grid cell size, i.e. E±1/2Aj =Aj±1/2. Note that (wφ)i are defined at the grid points that correspond to φ¯ij±1/2 (see Fig. 4), so that(wφ)ij±1/2 = 12 (wij +wij±1)φ¯ij±1/2. The most difficult challenge in using (19) is to determine how to use theshift operator appropriately to have a stable scheme.An in-depth analysis of the nonlinear conservation laws (and the associated numerical schemes) entailsfinding the characteristics, which depend on the derivative of the nonlinear flux, i.e. ∂q˜pz/∂(wφ¯) (see e.g. [20,21]). Unfortunately, the proppant flux q˜pz depends on w, which, in turn, is functionally dependent upon φ¯through (6a), so that exact evaluation of this derivative is not trivial. For this reason, it is assumed that thevariation of w with respect φ¯ is small (∂w/∂φ¯ ≈ 0), in which case the differentiation of (6a) with respect toφ¯ yieldsqsz = const. ⇒ vsz =qszw = const., (20)where vsz is an average velocity of the slurry, and “const.” means constant with respect to φ¯. The implicationof equation (20) is that the velocity of the slurry may be taken as constant during the evaluation of thederivative of the flux, i.e.∂q˜pz∂(wφ¯) = B(wa)w vsz∂∂(wφ¯)[Qˆp(wφ¯w ,wa)]− B(wa)a2wµ′ (ρp−ρf)g ∂∂(wφ¯)[Gˆp(wφ¯w ,wa)]− V1,2l1,2ζ1,2. (21)One may also interpret the assumption ∂w/∂φ¯ ≈ 0 and consequently (20) in a different way. In the numericalscheme, first the width profile is updated, and then the equation for particle concentration is solved. Oncethe width is calculated, it cannot be changed while solving the proppant transport equation (within the sametime step), hence w = const., which consequently implies (20). It is interesting to observe that accordingto Fig. 2, derivatives of the functions Qˆp and Gˆp appearing in (21) can be large in magnitude for φ¯ ≈ 1.7In particular, it can be shown that their absolute values behave like w2/a2  1. Note, however, thatQˆs = O(a2/w2) and consequently vsz=O(a2/w2) for such particle concentrations, so that the product of thederivative of Qˆp and vsz is O(1) in terms of the small parameter a/w. Since a2 multiplies the gravitationalsettling term in (21), the fact that the derivative of Gˆp is O(w2/a2) similarly does not make the whole termlarge in magnitude. This shows that even in the limit of very small particles, the flux derivative (21) isbounded. The width of the fracture can also be small, but, clearly, the differentiation in (21) cancels w andleads to some final value of the flux derivative for small w.Since (12b) has the form of a conservation law, equation (21) allows us to calculate the velocity of thenonlinear wave, which is then used to find the sign of the “wind” and utilize it in the numerical scheme.One of the best options is to use the Godunov scheme, which, however, requires the solution of the Riemannproblem [21]. Unfortunately, finding the solution of the Riemann problem may be challenging and requiressignificant amount of the computation time (since the proppant flux depends on functions that are computednumerically), for this reason an approximate Riemann solver is used. To assist with the construction of thenumerical scheme, it is noted that q˜pz = q˜pz (wφ¯, w), and that the shock velocities between the elements canbe defined asV shj =q˜pz(wij+1/2φ¯ij+1/2, wij+1/2)− q˜pz(wij−1/2φ¯ij−1/2, wij−1/2)wij+1/2φ¯ij+1/2 − wij−1/2φ¯ij−1/2, (22)where wij±1/2 = 12 (wij + wij±1). Fig. 5 shows the algorithm for determining the proppant flux based on thesign of the “wind” (21) and the shock velocity (22). There are three cases, where: (a) the “wind”, calculatedwij(a)wij(b)wijV shj > 0V shj < 0V shj = 0(c) ⇢¯ij1/2 ¯ij+1/2 ¯ij1/2 ¯ij+1/2¯ij1/2 ¯ij+1/2⇣ @q˜px@(w¯)⌘ij1/2> 0⇣ @q˜px@(w¯)⌘ij+1/2> 0⇣ @q˜px@(w¯)⌘ij+1/2< 0⇣ @q˜px@(w¯)⌘ij1/2< 0⇣ @q˜px@(w¯)⌘ij1/2⇣ @q˜px@(w¯)⌘ij+1/2< 0q˜pxij = q˜px(wij1/2¯ij1/2, wij1/2)q˜pxij = q˜px(wij1/2¯ij1/2, wij1/2)q˜pxij = 12hq˜px(wij1/2¯ij1/2, wij1/2)+q˜px(wij+1/2¯ij+1/2, wij+1/2)iq˜pxij = q˜px(wij+1/2¯ij+1/2, wij+1/2)q˜pxij = q˜px(wij+1/2¯ij+1/2, wij+1/2)⇣ ⇣⇣zzzz zzzzzzzzzzzzzFigure 5: Schematics of the algorithm for approximating the proppant flux at the point that corresponds towij : (a) if the “wind” at both neighbouring points is positive, then use the “left” value, (b) if the “wind”at both neighbouring points is negative, then use the “right” value, and (c) if the direction of the “wind” isdifferent for neighbouring points, then use the sign of V shj to determine value of the flux.according to (21), is positive for both neighbouring points, (b) the “wind” is negative for both neighbouringpoints, and (c) the direction of the “wind” is different for the neighbouring points. In the latter case, theshock velocity in (22) is used to determine the value of the proppant flux, see Fig. 5.8Since equation (19) represents an explicit scheme, stability poses a restriction on the magnitude of thetime step ti − ti−1. In other words, the time step has to be reduced to make sure that the Courant-Friedrichs-Lewy (CFL) condition [21] is satisfied. On the other hand, since the numerical scheme for thecrack opening (14) is implicit, there is no restriction on the time step for solving (14). To allow for arbitrarylarge time steps ti−ti−1 for the whole numerical algorithm, the time step ti−ti−1 is subdivided into smalltime steps when solving (19), each of which satisfy the CFL condition. This decomposition enables us to uselarge time steps for the stiff part of the problem (6a) and (8), which would require a time step restriction∆t ∼ O(∆z3), while the less stiff equation (6b) can be treated explicitly as it only involves a CFL condition∆t ∼ O(∆z).It should also be noted that the equations are first scaled and then solved numerically. Once the solutionis obtained, then the dimensions of the parameters are recovered. More information about scaling for theKGD fracture can be found in [19].3.3 Numerical examplesThis section is devoted to numerical examples that highlight the effects produced by the proppant. Due tothe considerable number of parameters that influence the result, it is instructive to specify a reference set ofparameters, and then change just some of them if needed. All numerical simulations in this section start attstart = 1 s, assume a symmetric crack of length l1 = l2 = 1 m, and take the opening being elliptic with themaximum width wmax = 5×10−4 m. The fracture is then propagated until tpr = 1000 s using pure fluid, andthereafter the proppant is introduced, so for t > tpr, the mixture of the proppant and the fluid is used. Thesimulations run until tend, which is specified uniquely for each calculation. The input volume concentrationof particles is taken φ¯0 = 0.2, but note that φ¯0 is the scaled concentration, so that the true concentrationis φmφ¯0. Other parameters used for the calculations are E′ = 25×109 Pa for the plane strain modulus,µ′ = 1.2 Pa·s for the viscosity (times 12), Q0 = 2×10−4 m2/s for the inlet flux, C ′ = 5×10−5 m/s1/2 forthe leak-off coefficient, K1c = 106 Pa·m1/2 for the fracture toughness, a = 4×10−4 m for the particle radius,∆ρ = 1300 kg/m3 for the difference in mass densities between the proppant and the fluid and g = 9.8 m/s2for the gravitational acceleration. A stress barrier is assumed to be symmetric, located lσ = 10 m from theinlet, and to have a magnitude ∆σ = 2.5× 106 Pa.The problem with no leak-off is considered first. Fig. 6 plots pictures of the fracture width and proppantconcentration at different time instants, as well as pressure and length histories. First of all, without theleak-off, the proppant does not accumulate rapidly in the tip regions, while, at the same time it can reachthe bottom tip due to the gravitational settling. Even though the proppant changes the viscosity notably,its uneven distribution does not influence the symmetry much. This is because high viscosity has the mostinfluence in the regions with high pressure gradients, i.e. near fracture tips. So, there is almost no influenceof proppant before it reaches the tip region and starts to accumulate there. This hypothesis is supportedby the length histories shown in Fig. 6. Indeed, noticeable asymmetry is induced only at t ≈ 4500 s, whichcorresponds to the time of developing of the bottom proppant plug. Note that the kink on the pressurehistory at t ≈ 110 s corresponds to reaching the stress barrier. There is no visible pressure change due tothe injection of the proppant as well as due to the developing of the bottom plug.To investigate the tip screen-out effect, Fig. 7 shows the results of a simulation for the reference parameterset, which includes the leak-off. There are three features in the pressure diagram. First, there is a kink thatcorresponds to the time when the fracture reaches the stress barrier. Then, there is a characteristic pressureincrease due to the developing of the bottom proppant plug at t ≈ 1800 s. After that, when the proppantreaches the top fracture tip, it causes an additional pressure rise at t ≈ 2000 s. These observations are alsosupported by the length histories which clearly indicate the slower fracture growth due to the stress barriers,initiation of the asymmetry at t ≈ 1800 s, as well as the slower propagation of the top fracture tip afterthe formation of the proppant plug there. Also note that the fracture width is noticeably affected by theformation of the proppant plugs.9−100 −50 0 50 10000.0050.010.015z [m]w[m]−100 −50 0 50 10000. [m]φ¯0 1000 2000 3000 4000 50001234t [s]p[MPa]0 1000 2000 3000 4000 5000050100150t [s]l[m]  l1l2Stress barrierStress barrierBottom plugFigure 6: The fracture width and proppant distribution (top pictures) for the reference parameter set and noleak-off at different time instants tev =500 s (no proppant at this time), 1100 s, 3000 s, and 5000 s. Bottompictures show the histories of pressure at the inlet and the lengths l1 (distance from the inlet to the bottomtip) and l2 (distance from the inlet to the top tip).10−40 −20 0 20 4000.0050.010.015z [m]w[m]−40 −20 0 20 4000. [m]φ¯0 1000 2000 30001234t [s]p[MPa]0 1000 2000 3000010203040t [s]l[m]  l1l2Stress barrierBottom plugTop plugStress barrierBottom plugTop plugFigure 7: The fracture width and proppant distribution (top pictures) for the reference parameter set and noleak-off at different time instants tev =500 s (no proppant at this time), 1100 s, 1800 s, and 3500 s. Bottompictures show the histories of pressure at the inlet and the lengths l1 (distance from the inlet to the bottomtip) and l2 (distance from the inlet to the top tip).11xzl(t)h(x, t) HFracture footprintStress barriers zyw(x, z)gFigure 8: Schematics of the P3D fracture.4 Numerical solution for a P3D fracture4.1 Problem formulationTo highlight the versatility of the proppant transport model, a simple multidimensional case is considered,namely the pseudo-3D (P3D) model for hydraulic fracture propagation with symmetrical stress barriers [22].Fig. 8 shows the schematics of the fracture. The fracture is propagating between two symmetric stressbarriers, where an additional stress ∆σ further compresses the fracture in the regions |z| > H/2. Thefracture tip is assumed to be vertical, the horizontal length of the fracture is denoted by l(t), while theheight of the fracture is h(x, t), see Fig. 8. Other assumptions of the model are: i) the fluid pressure isuniform over the height of the fracture, i.e. does not depend on z, p=p(x), which implies the symmetry ofthe fracture, ii) a plain strain elasticity condition exists in any vertical (y, z) plane, and iii) leak-off occursonly in the reservoir layer (|z|<H/2) and follows the Carter’s leak-off model [18].To facilitate the development of the appropriate proppant transport model, it is useful to formulate the2D equations for both fracture width and the particle concentration, so that∂w∂t +∂qsx∂x +∂qsz∂z +C ′H(H−2z)H(2z+H)√t− t0(x)= Q(z)H δ(x),∂wφ¯∂t +∂qpx∂x +∂qpz∂z =φ¯0Q(z)H δ(x), (23)where the fluxes are given in (2), Q(z) is a source density that is distributed over the vertical coordinatez, and H is a Heaviside step function. The boundary conditions for (23) require the vanishing of all fluxesalong the fracture boundary, as well as the appropriate asymptotic behaviour of the width near the verticaltips, see [11] for details. Following [22], since the pressure is assumed to be uniform over the height, theelasticity equations can be solved to obtainw(x, z) = 2E′ (p(x)−∆σ)χ+4∆σpiE′{χ arcsin(Hh)− zln∣∣∣Hχ+ 2zψHχ− 2zψ∣∣∣+ H2 ln∣∣∣χ+ ψχ− ψ∣∣∣}, (24)where χ =√h2 − 4z2, ψ =√h2 −H2, while E′ as in (8) is the plane-strain Young’s modulus. Again, asfollows from [22], the application of the toughness propagation criterion leads top = ∆σ[1 +√2piHKIc∆σ√Hh −2pi arcsin(Hh)], (25)12where KIc is mode I fracture toughness. Formulas (24) and (25) apply in the regions where h > H. Whenh = H, an elliptic fracture width profile is used instead of (24), i.e. w = 2(E′)−1χp(x). To formulate theP3D model, one also needs to introduce average width, flux, and total inlet flux as follows:w¯ = 1H∫ h/2−h/2w dz, q¯sx =1H∫ h/2−h/2qsx dz, Q0 =1H∫ h/2−h/2Q(z) dz. (26)Note that φ¯ is an average concentration over the thickness (i.e. in the out-of-plane y direction) of the fracture,while w¯ and q¯sx are respectively the width and flux, averaged over the height (i.e. z direction) of the fracture.With (26) in mind, equation (23a) can be integrated over z to obtain∂w¯∂t +∂q¯sx∂x +C ′√t− t0(x)= Q0H δ(x), (27)whereq¯sx = −∂p∂x1H∫ h/2−h/2[w3µ′ Qs(φ¯) + a2wµ′ D(φ¯)]dz. (28)Relation (24) can also be integrated to obtainw¯ = HE′(√ pi2HKIc( hH)3/2+ ∆σ√h2H2−1), h > H, (29)Given the fact that p can be expressed as a function of w¯ from (25) and (29), equation (27) can be solved ifthe variation of φ¯ versus z is provided (so that the integral in (28) can be calculated). Once solved for w¯,(29) can be used to find h, which finally allows us to obtain w from (24) and (25). In other words, knowingthe average width w¯ is enough to “restore” the fracture width profile w. This property allows us to obtainthe 2D fracture footprint by solving the one-dimensional problem governed by (27). Unfortunately, sucha useful “restoring” property does not apply for the proppant. Indeed, having only the average proppantconcentration is not sufficient to “restore” the concentration profile, as there are many different (physicallyadmissible) particle distributions that can have the same mean value. For this reason, a 2D proppanttransport model has to be considered.To assist the solution of the proppant transport equation, the vertical (or z) component of the slurry fluxhas to be computed first. Formally, due to the assumptions of the model, there is no pressure gradient inthe vertical direction, which implies no flux in the vertical direction. This, however, should be interpretedin a sense that the vertical flux is small compared to the horizontal flux, but not necessary zero. To find thevertical component of the slurry flux, one can integrate (23a) to determineqsz =∫ z−h/2[Q(z)H δ(x)−C ′H(H−2z)H(2z+H)√t− t0(x)− ∂w∂t −∂qsx∂x]dz. (30)Finally, the system of equations that describes the P3D problem with proppant is∂w¯∂t +∂q¯sx∂x +C ′√t− t0(x)= Q0H δ(x),∂wφ¯∂t +∂qpx∂x +∂qpz∂z =φ¯0Q(z)H δ(x), (31)where q¯sx is given in (28), the relations between w, w¯, p, and h are given by (24), (25) and (29), while theproppant fluxes areqpx = −B(wa)w3µ′∂p∂xQp(φ¯),qpz = B(wa)Qˆp(φ¯, wa)qsz −B(wa)a2wµ′ (ρp−ρf)gGˆp(φ¯, wa), (32)13where g is the gravitational acceleration constant. The boundary conditions for (31) areq¯sx|x=l = 0, w¯|x=l = 0. (33)Note that the boundary conditions at the top and bottom sides of the fracture are accounted for in (30).Also, the blocking function B restricts the presence of the particles near the fracture tip, so that the zero-proppant-flux boundary condition is always satisfied automatically.4.2 Numerical algorithmTo facilitate the numerical calculations, first, the problem parameters are scaled aslˆ = ll∗, hˆ = hh∗, tˆ = tt∗, wˆ = ww∗, pˆ = pp∗, qˆsx =qsxq∗, qˆsz =l∗qszh∗q∗,Kˆ = K1cK∗, Cˆ = C′C∗, aˆ = aa∗, gˆ = gg∗, (34)where all “hat” quantities are dimensionless, while the scaling factors are computed asl∗ =H4∆σ4Q0µ′E′3, h∗ = H, t∗ =∆σ5H6Q20µ′E′4, w∗ =H∆σE′ , p∗ = ∆σ, q∗ =Q0H ,K∗ = ∆σ√2Hpi , C∗ =Q0E′µ′1/2H2∆σ3/2 , a∗ =H∆σE′ , g∗ =Q0µ′E′3∆ρH4∆σ3 . (35)The biggest advantage of this scaling lies in the fact that it highlights the number of independent parametersthat govern the problem. For this problem, there are four of such quantities: Kˆ, Cˆ, aˆ and gˆ. Note that thescaling (34) is done slightly differently from [22].To aid the solution of the moving boundary problem, a moving mesh in both x and z directions isintroduced. The following scaled coordinates are introduced asξ = xlˆ, ζ = zhˆ, (36)where lˆ = lˆ(tˆ) and hˆ = hˆ(ξ, tˆ). In this case, the derivatives transform as∂(·)∂tˆ∣∣∣x,z= ∂(·)∂tˆ∣∣∣ξ,ζ− V ξlˆ∂(·)∂ξ∣∣∣t,ζ−(∂hˆ∂tˆ− V ξlˆ∂hˆ∂ξ) ζhˆ∂(·)∂ζ∣∣∣t,ξ,∂(·)∂x∣∣∣t,z= 1lˆ∂(·)∂ξ∣∣∣t,ζ− ∂hˆ∂ξζlˆhˆ∂(·)∂ζ∣∣∣t,ξ, (37)∂(·)∂z∣∣∣t,x= 1hˆ∂(·)∂ζ∣∣∣t,ξ,where V = dlˆ/dtˆ is the velocity of the crack tip. By substituting (37) into (31), and simplifying the result,one may write∂ ˆ¯w∂tˆ− V ξlˆ∂ ˆ¯w∂ξ +1lˆ∂ ˆ¯qsx∂ξ +Cˆ√tˆ− t0(lˆξ)= 1lˆδ(ξ),∂hˆlˆwˆφ¯∂tˆ+ ∂q˜px∂ξ +∂q˜pz∂ζ = φ¯0hˆQˆ(ζ)δ(ξ), (38)whereq˜px = hˆ(qˆpx − V ξwˆφ¯), q˜pz = lˆqˆpz −∂hˆ∂tˆlˆζwˆφ¯− ∂hˆ∂ξ ζq˜pxhˆ, (39)14q˜xq˜zqˆxqˆzV ⇠✓@hˆ@ tˆ⇣hˆ⇣⇣⇠lˆ⇠Figure 9: The element in the physical domain (left) and the computational domain (right).are the fluxes that account for the moving mesh terms. Note that similar changes of variables can be appliedto (23a), in which case the integral (30) for the calculation of the flux transforms toqˆsz =(∂hˆ∂tˆwˆ + ∂hˆ∂ξq˜sxhˆlˆ)ζ + 1lˆ∫ ζ−1/2[hˆQˆ(ζ)δ(ξ)− hˆlˆ CˆH(1−2hˆζ)H(2hˆζ+1)√tˆ− t0(lˆξ)− ∂hˆlˆwˆ∂tˆ− ∂q˜sx∂ξ]dζ, (40)whereq˜sx = hˆ(qˆsx − V ξwˆ).Note that the fluxes in (39) can also be derived from physical considerations. Fig. 9 shows the elements in thephysical (on the left) and computational (on the right) domains. The fluxes in the picture are generic and canbe applied to either the slurry or proppant. The element in the computational domain is rectangular and doesnot move, while the corresponding element in the physical domain is distorted and moves horizontally withthe velocity V ξ and vertically with velocity ∂hˆ/∂tˆ ζ. The angle θ can be found from tan(θ) = −∂h/∂ξ ζ/lˆ.Since the fluxes through the sides of the element should be preserved, one can write(qˆpx − V ξwˆφ¯)hˆ∆ζ = q˜px∆ζ,(qˆpz −∂hˆ∂tˆζwφ¯+ (qˆpx − V ξwˆφ¯) tan(θ))lˆ∆ξ = q˜pz∆ξ,which allows us to recover (39).To close the system of equations, it is noted thatqˆsx = −1lˆ∂pˆ∂ξ[wˆ3Qs(φ¯) + aˆ2wˆD(φ¯)],ˆ¯qsx = −hˆlˆ∂pˆ∂ξ∫ 1/2−1/2[wˆ3Qs(φ¯) + aˆ2wˆD(φ¯)]dζ,qˆpx = B( wˆaˆ)Qˆp(φ¯, wa)qˆsx, (41)qˆpz = B( wˆaˆ)Qˆp(φ¯, wa)qˆsz −B( wˆaˆ)aˆ2wˆgˆGˆp(φ¯),where the pressure derivative is∂pˆ∂ξ =dpˆdhˆdhˆd ˆ¯wd ˆ¯wdξ = Y (hˆ)d ˆ¯ωdξ , (42)whereY (hˆ) = dpˆdhˆdhˆd ˆ¯w =4√hˆ− 2Kˆ√hˆ2−1pihˆ2(3Kˆ√hˆ2−1 + 2√hˆ) . (43)15Finally, to “restore” the fracture opening based on the average width, one can rewrite (24), (25) and (29) aswˆ(hˆ) = 4pi[Kˆχˆ√hˆ− hˆζln∣∣∣ χˆ+ 2ζψˆχˆ− 2ζψˆ∣∣∣+ 12ln∣∣∣ hˆχˆ+ ψˆhˆχˆ− ψˆ∣∣∣],ˆ¯w(hˆ) = Kˆhˆ3/2 +√hˆ2−1, hˆ > 1, (44)where χˆ =√1− 4ζ2 and ψˆ =√hˆ2 − 1. Relations (44) allow us to calculate a function wˆ( ˆ¯w), i.e. “restore”the fracture width. When hˆ = 1, the latter equations combine to yieldωˆ = 4pi√1−4ζ2 ˆ¯w, hˆ = 1.The numerical algorithm for the P3D geometry is somewhat similar to that for the KGD fracture, inthat it is divided into two main parts: i) calculating the fracture propagation and ii) updating the proppantconcentration inside the fracture. Since the equation for fracture propagation is very similar to that for KGDfracture, a similar algorithm is used to update the fracture footprint. The average width ˆ¯w is approximatedby a piece-wise constant function, the time derivative is approximated by backward differences, and theintegral in (41) is approximated using midpoint rule, in which case (38a) is reduced to a system of algebraicequations that is solved iteratively. The big difference, however, lies in the absence of a pressure singularitynear the right tip, in which case the velocity of the crack tip is calculated based on the zero flux condition.Equation (38b) is solved numerically using a finite volume method and a generalization of the 1-dimensionalalgorithm shown in Fig. 5. The analog of the condition (20) isˆ¯qsx = const.,i.e. the average flux stays constant during the differentiation with respect to φ¯. Another difference comesfrom the fact, that a line source has to be used, and, moreover, that the line cannot go all the way to thefracture boundary since the proppant cannot be there due to the blocking function. The distribution ofthe intensity of the source is taken proportional to the cube of the fracture width, and contained insidew > (2N+ 12 )a (N =3 is used for all calculations). In this case there is no proppant in the prohibited areasand the source is concentrated near the centre of the fracture height. Other details about the numericalscheme are analogous to the 1D model and omitted for brevity.4.3 Numerical examplesThis section covers several numerical examples that highlight the effects associated with the presence ofproppant. First, it should be noted that the numerical code (without proppant) was tested against thesolutions in [22], and the results showed good agreement. There are two main effects associated with thepresence of the proppant, that are considered in the examples: i) gravitational settling and ii) tip screen-out.All figures in this section that plot the fracture footprint and show the proppant concentration in colour.Note that the maximum value for φ¯ is 1, since the concentration is scaled by φm. For all figures, simulationsstart at tstart = 1 s, with l = 1 m, and an elliptic width profile with a maximum opening wmax = 10−3 m.At time tpr the proppant is introduced, and at tend the simulation stops. Note that the effect of the initialcondition decays with time, so, as long as the total initial volume of the fracture is sufficiently small, there ispractically no effect associated with the initial condition. The reference problem parameters are H = 25 mfor the width of the reservoir layer, µ′ = 1.2 Pa·s for the shear viscosity (times 12), E′ = 25 × 109 Pa forthe plane strain modulus, Q0 = 2×10−2 m3/s for the injection rate, ∆σ = 2.5×106 Pa for the magnitude ofthe stress barriers, ∆ρ = 1300 kg/m3 for the difference between proppant and carrying fluid mass densities,K1c = 106 Pa·m1/2 for the fracture toughness, C ′ = 5×10−5 m/s1/2 for the leak-off coefficient, a = 4×10−4 mfor the particle radius, g = 9.8 m/s2 for the gravity constant, φ¯0 = 0.2 for the proppant concentration atthe inlet, and tpr = 1000 s for start time of proppant injection. At the same time, different values of tend areused. For the all figures in this section, the parameters are assumed to be taken from the above referenceset, except those, which are specified directly.16It is important to recognize the presence of the time scale associated with the gravitational settling. Theasymptotic behaviour of the function Gˆp(φ¯) (see [12]) can be used to estimate the settling velocity, whilethe settling time can be calculated based on the vertical distance the proppant needs to travel, which is 1/2in the scaled formulation. Combining these assumptions and the last equation in (41), the settling time canbe estimated asts =12vˆsettl= 12 83 aˆ2gˆt∗ =3∆σ4H416∆ρa2gQ0E′3,where vˆsettl is the dimensionless settling velocity calculated for small particle concentrations. This settlingtime needs to be compared to the proppant injection duration. Since the proppant is first injected at tpr,the duration of the injection is tend − tpr. In this case, it might be useful to introduce the ratio between twotime scales which determines whether the gravitational settling is significant or notGs =16∆ρa2gQ0E′3(tend−tpr)3∆σ4H4 . (45)If the parameter Gs  1, then the settling occurs before the end of the fracturing job, while if Gs  1, thenthe gravity does not affect the proppant distribution much. It is interesting to note that if tend−tpr =const.,the viscosity does not enter (45), so that changing the viscosity of the carrier fluid alone cannot be usedto alter the settling pattern. This counterintuitive phenomenon can be understood in the following way: ahigher viscosity leads to a slower vertical settling velocity, however, at the same time, the horizontal velocitybecomes smaller too. Since both settling and horizontal velocities are proportional to the inverse of theviscosity, the direction of the velocity vector does not change, and so the proppant pattern is unaffected.However, if the design fracture length is regarded as fixed, then the total treatment time tend becomes afunction of viscosity, and then Gs will no longer be independent of µ′.In addition, since high powers of E′, ∆σ and H appear in the scaling parameters (34), and in particularin (45), the fracture footprint and proppant distribution become very sensitive to the values of E′, ∆σ andH. In other words, if any of those quantities is measured inaccurately, the prediction of the model can beunreliable.To illustrate the importance of the parameter Gs, Fig. 10 plots the distribution of φ¯ for different valuesof Gs. As can be seen from Fig. 10, Gs = O(1) leads to some skewness of the proppant distribution, whichis nearly identical (as it should be) for the top pictures. The bottom left picture corresponds to Gs  1,and so the effect of the gravity is significant. The bottom right picture corresponds to Gs  1, and so theeffect of the gravitational settling is minimal. One can also see the effect of the “blocking” functions, as theproppant cannot sink all the way to the bottom of the fracture and the maximum concentration in the plugdoes not increase beyond unity (the maximum allowed concentration).Crack tip screen-out is another very important consequence of the presence of proppant. It is importantto recognize that particles can reach the crack tip even without leak-off. First, the proppant flows faster bya factor 1.2 for small concentrations, which can be concluded from the asymptotic behaviour of Qˆp, see [12].This happens because of the nonuniform distribution of particles over the width of the channel. A similarthing happens in the vertical direction, the proppant tends to flow in the areas where the fracture is wide,while the fluid flows everywhere inside the fracture. So, the proppant gets concentrated near the centre,which again implies that, on average, it flows faster than the slurry. To illustrate this phenomenon, Fig. 11plots the fracture footprints and the proppant distributions at different time instants t = 800 s, t = 1200 s,t = 3000 s and t = 4500 s for C ′ = 0 (other parameters are taken from the reference set). The variationof the pressure at the inlet, the length of the fracture, and the height at the inlet versus time are alsoshown. As can be seen from the pictures, the length of the fracture is below 100 m right before the proppantinjection starts, and by a length of 300 m, the proppant is already in the tip region. The proppant travelsapproximately 250 m, while the fracture grows by 200 m during the same time period, which shows thatthe proppant is faster by approximately 25%. There is no plug formation (in the x direction), however, atthese times and the proppant is distributed almost uniformly inside the fracture. At the same time, thereis a plug in the vertical z direction due to gravitational settling. The variation of the pressure, the length,and the height is smooth, although there is a small kink at t = 1000, which corresponds to the beginning ofproppant input.17x [m]z[m]  0 50 100 150 200 250 300−60−40−20020406000. [m]z[m]  0 100 200 300 400 500−60−40−20020406000. [m]z[m]  0 50 100 150 200−150−100−5005010015000. [m]z[m]  0 50 100 150 200 250 300−60−40−20020406000. 10: Fracture footprint and the proppant concentration φ¯ indicated by colour calculated for thereference parameters and tend = 4000 s, except C ′ = 0 (top left), C ′ = 0 and µ′ = 0.24 Pa·s (top right),C ′ = 0 and ∆σ = 1.5 × 106 Pa (bottom left) and C ′ = 0 and H = 35 m (bottom right). The gravitationalsettling parameters are Gs = 0.67 for both top pictures (notable settling), Gs = 5.15 for the bottom leftpicture (significant settling), and Gs = 0.17 for the bottom right picture (almost no settling).18x [m]z[m]  0 50 100 150 200 250 300−60−40−20020406000. [m]z[m]  0 50 100 150 200 250 300−60−40−20020406000. [m]z[m]  0 50 100 150 200 250 300−60−40−20020406000. [m]z[m]  0 50 100 150 200 250 300−60−40−20020406000. 1000 2000 3000 4000 50000.511.522.5t [s]p[MPa]0 1000 2000 3000 4000 50000100200300400t [s]l[m]0 1000 2000 3000 4000 500020406080100120t [s]h[m]Figure 11: Fracture footprint with the proppant concentration φ¯ indicated by colour calculated for thereference parameters and C ′ = 0 at t = 800 s (top left), t = 1200 s (top right), t = 3000 s (centre left) andt = 4500 s (centre right). The case with t = 4500 s corresponds to Gs = 0.89. Bottom pictures show thepressure at the inlet, the length of the fracture, and the height of the fracture at the inlet versus time.19x [m]z[m]  0 20 40 60 80 100 120 140−40−200204000. [m]z[m]  0 20 40 60 80 100 120 140−40−200204000. [m]z[m]  0 20 40 60 80 100 120 140−40−200204000. [m]z[m]  0 20 40 60 80 100 120 140−40−200204000. 500 1000 1500 2000 2500 30000.511.522.5t [s]p[MPa]0 500 1000 1500 2000 2500 3000050100150t [s]l[m]0 500 1000 1500 2000 2500 300020406080100t [s]h[m]Figure 12: Fracture footprint with the proppant concentration φ¯ indicated by colour calculated for thereference parameters at t = 800 s (top left), t = 1200 s (top right), t = 1700 s (centre left) and t = 3000 s(centre right). The case with t = 3000 s corresponds to Gs = 0.45. Bottom pictures show the pressure atthe inlet, the length of the fracture, and the height of the fracture at the inlet versus time.20To promote the accumulation of proppant in the tip region, leak-off needs to be introduced. Fig. 12plots the fracture footprint and the proppant distribution at different time instants t = 800 s, t = 1200 s,t = 1700 s and t = 3000 s (all parameters are taken from the reference set). The variation of the pressureat the inlet, the length of the fracture, and the height at the inlet versus time are also shown. As can beseen, the leak-off significantly retards the fracture propagation, and even at t = 800, the fracture is 50%shorter compared to that in Fig. 11. Once the proppant is introduced, it reaches the crack tip much faster,accumulates there, and significantly slows further fracture propagation. After the plug is formed, only asmall amount of fluid can penetrate and so the fracture starts to grow noticeably in the vertical direction.5 SummaryThis paper applies a model for proppant transport, which is based on an empirical constitutive law for themixture of a viscous fluid with spherical particles, to two hydraulic fracturing problems, namely with theKGD and the P3D geometries. In the adopted formulation, the slurry flux has two terms, one Poiseuille-law-type term with an effective viscosity (which goes to infinity as the concentration reaches a critical value),and a Darcy-law-type term, where the average velocity is proportional to the particle size squared and thepressure gradient. The flux of the particles also has two terms, one proportional to the slurry flux, andanother related to gravitational forces. The first term describes the advective motion, while the secondterm describes gravitational settling. The proppant transport model is first applied to the KGD fracturegeometry with stress barriers. The numerical simulations show that the developed model is able to capturethe initiation and further growth of a proppant plug in the crack tip region, which leads to tip screen-out.The gravitational settling introduces asymmetry, leads to faster screen-out at one side of the crack, and maystop the propagation there for some time. The proppant transport model is then implemented with the P3Dfracture geometry. Despite the fact that the P3D model reduces to the solution of a 1D equation, the proppanttransport cannot be treated in a similar fashion, and requires the numerical solution of a 2D problem. Twomain objectives in the analysis of the numerical solutions include the gravitational settling and tip screen-out. A dimensionless parameter, which controls the magnitude of the particle settling, is introduced. Oneinteresting fact is that this parameter does not directly depend on the viscosity of the fluid. It is furthershown that the particles can reach the tip of the fracture even without leak-off. This occurs due to the factthat the proppant is concentrated near the centre of the channel, and thus, on average, gets transportedfaster than the carrying fluid. When leak-off is introduced, the proppant reaches the crack tip region notablyfaster and accumulates there forming a plug. Once the plug is developed, only a small amount of fluid canpenetrate the plug (due to the Darcy-law-type term), which switches fracture growth predominantly to thevertical direction. The main drawbacks of the model include: its inability to capture asymmetry causedby gravitational settling; and the rigid plug property, in which the proppant can sustain some stress oncethe fracture tends to close. Both issues cannot be implemented since one of the key assumptions of theP3D model – uniform pressure over the height of the fracture (and the resulting solution of the elasticityequation), would not be satisfied. This could be overcome only by adding a proppant transport module intoa fully planar 3D hydraulic fracture propagation model, which is a challenging problem for future research.AcknowledgementsThe authors gratefully acknowledge the support of the British Columbia Oil and Gas Commission and theNSERC discovery grants program.References[1] U. Frank and N. Barkley. Remediation of low permeability subsurface formations by fracturing en-hancements of soil vapor extraction. J. Hazard. Mater., 40:191–201, 2005.[2] A.S. Abou-Sayed, D.E. Andrews, and I.M. Buhidma. Evaluation of oily waste injection below the21permafrost in prudhoe bay field. In Proceedings of the California Regional Meetings, Bakersfield, CA,Society of Petroleum Engineers, Richardson, TX, pages 129–142, 1989.[3] R.G. Jeffrey and K.W. Mills. Hydraulic fracturing applied to inducing longwall coal mine goaf falls. InPacific Rocks 2000, Balkema, Rotterdam, pages 423–430, 2000.[4] M.J. Economides and K.G. Nolte, editors. Reservoir Stimulation. John Wiley & Sons, Chichester, UK,3rd edition, 2000.[5] S.A. Khristianovic and Y.P. Zheltov. Formation of vertical fractures by means of highly viscous fluids.In Proc. 4th World Petroleum Congress, volume 2, pages 579 – 586, 1955.[6] D.I. Garagash and E. Detournay. Near tip processes of a fluid-driven fracture. ASME J. Appl. Mech.,67:183–192, 2000.[7] J.I. Adachi and E. Detournay. Self-similar solution of a plane-strain fracture driven by a power-lawfluid. Int. J. Numer. Anal. Meth. Geomech., 26:579–604, 2002.[8] E. Detournay and D. Garagash. The tip region of a fluid-driven fracture in a permeable elastic solid.J. Fluid. Mech., 494:1–32, 2003.[9] E. Detournay. Propagation regimes of fluid-driven fractures in impermeable rocks. Int. J. Geomech.,4:1–11, 2004.[10] J. Adachi, E. Siebrits, A. Peirce, and J. Desroches. Computer simulation of hydraulic fractures. Int. J.Rock Mech. Min. Sci., 44:739–757, 2007.[11] A. Peirce and E. Detournay. An implicit level set method for modeling hydraulically driven fractures.Comput. Methods Appl. Mech. Engrg., 197:2858–2885, 2008.[12] E.V. Dontsov and A.P. Peirce. Slurry flow, gravitational settling, and a proppant transport model forhydraulic fractures. http://hdl.handle.net/2429/50000, 2014.[13] F. Boyer, E. Guazzelli, and O. Pouliquen. Unifying suspension and granular rheology. Phys. Rev. Lett.,107:188301, 2011.[14] S.A. Boronin and A.A. Osiptsov. Two-continua model of suspension flow in a hydraulic fracture. DokladyPhysics, 55:199–202, 2010.[15] E. Chekhonin and K. Levonyan. Hydraulic fracture propagation in highly permeable formations, withapplications to tip screenout. Int. J. Rock Mech. Min., 50:19 – 28, 2012.[16] D. Eskin and M.J. Miller. A model of non-newtonian slurry flow in a fracture. Powder Technology,182:313–322, 2008.[17] J.R. Lister. Buoyancy-driven fluid fracture: the effects of material toughness and of low-viscosityprecursors. J. Fluid. Mech., 210:263–280, 1990.[18] E.D. Carter. Optimum fluid characteristics for fracture extension. In Howard GC, Fast CR, editors.Drilling and production practices, pages 261–270, 1957.[19] J.I. Adachi. Fluid-Driven Fracture in Permeable Rock. PhD thesis, University of Minnesota, 2001.[20] D. Bale, R. J. LeVeque, S. Mitran, and J. A. Rossmanith. A wave-propagation method for conservationlaws and balance laws with spatially varying flux functions. SIAM J. Sci. Comput., 24:955–978, 2002.[21] R.J. LeVeque. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, 2002.[22] J. I. Adachi, E. Detournay, and A. P. Peirce. An analysis of classical pseudo-3D model for hydraulicfracture with equilibrium height growth across stress barriers. Int. J. of Rock Mech. and Min. Sci.,47:625–639, 2010.22


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