Multiphase flow models for hydraulic fracturing technology

The technology of hydraulic fracturing of a hydrocarbon-bearing formation is based on pumping a fluid with particles into a well to create fractures in porous medium. After the end of pumping, the fractures filled with closely packed proppant particles create highly conductive channels for hydrocarbon flow from far-field reservoir to the well to surface. The design of the hydraulic fracturing treatment is carried out with a simulator. Those simulators are based on mathematical models, which need to be accurate and close to physical reality. The entire process of fracture placement and flowback/cleanup can be conventionally split into the following four stages: (i) quasi-steady state effectively single-phase suspension flow down the wellbore, (ii) particle transport in an open vertical fracture, (iii) displacement of fracturing fluid by hydrocarbons from the closed fracture filled with a random close pack of proppant particles, and, finally, (iv) highly transient gas-liquid flow in a well during cleanup. The stage (i) is relatively well described by the existing hydralics models, while the models for the other three stages of the process need revisiting and considerable improvement, which was the focus of the author’s research presented in this review paper. For stage (ii), we consider the derivation of a multi-fluid model for suspension flow in a narrow vertical hydraulic fracture at moderate Re on the scale of fracture height and length and also the migration of particles across the flow on the scale of fracture width. At the stage of fracture cleanaup (iii), a novel multi-continua model for suspension filtration is developed. To provide closure relationships for permeability of proppant packings to be used in this model, a 3D direct numerical simulation of single phase flow is carried out using the lattice-Boltzmann method. For wellbore cleanup (iv), we present a combined 1D model for highly-transient gas-liquid flow based on the combination of multi-fluid and drift-flux approaches. The derivation of the drift-flux model from conservation olaws is criticall revisited in order to define the list of underlying assumptions and to mark the applicability margins of the model. All these fundamental problems share the same technological application (hydraulic fracturing) and the same method of research, namely, the multi-fluid approach to multiphase flow modeling and the consistent use of asymptotic methods. Multi-fluid models are then discussed in comparison with semi-empirical (often postulated) models widely used in the industry.


Introduction
The technology of hydraulic fracturing of a hydrocarbon bearing underground formation is based on injecting a fluid laden with rigid particles under a high pressure (up to several hundred bar) into the well to create fractures in the porous medium, which are filled with particles. After the end of pumping, fractures closed on packed granular material provide high-conductivity channels to transport hydrocarbons from reservoir to the well and all the way up to the surface. The well may be vertical (when a single bi-wing fracture is formed) or near-horizontal with several perforation clusters providing reservoir contact (the so-called multi-stage fracturing in low-permeability formations). The latter case gives rise to several transversal fractures. With respect to different stages of the hydraulic fracturing technology, we consider four classes of multiphase flows that can be modelled within the multi-continua (or multi-fluid) approach [2]. In a more detail, we distinguish the following classes: (i) the flow of suspension of fluid with particles in a circular pipe at high Reynolds numbers during pumping, (ii) the flow of suspension in a narrow vertical hydraulic fracture at moderate Re during pumping [3,4], (iii) suspension filtration through a packed of proppant particles in a closed fracture during cleanup [5], and (iv) multiphase gas-liquid flow with admixture of rigid particles in a circular pipe during well startup, cleanup and testing in a wide range of the Reynolds numbers [6]. We discuss the advantages and limitations of the multi-fluid approach based on the simulation examples from each of the four classes of multiphase flows, in comparison with simplified semi-empirical approaches, e.g. the drift-flux model for well flows, the effective-fluid model for suspension transport in fractures, and the deep-bed filtration model. The talk ends up with recommendations for future research on the topic. A detailed review of the state of the art in numerical modeling of hydraulic fracturing can be found in [7], and the most recent review with a specific focus on fluid mechanics of hydraulic fracturing is [1].

Two-fluid model for suspension transport in a fracture
In this section, we present in a lumped form the two-fluid model for suspension flow and sedimentation in a hydraulic fracture. This model was derived in the multi-continua approach [2] using asymptotic methods in the lubrication approximation. The key assumption is that the width-to-length ratio of the fracture is a small parameter: ε = w/L ≪ 1. The cross-flow particle concentration profile is assumed uniform. In the width-averaged variables, the equations are as follows [3]: Here, Cartesian coordinate system Oxy is introduced in the cell plane, so that y-axis (with the basis vector e 2 ) is vertical and origin O is located in the bottom left corner of the computational domain; C is the particle volume fraction; w(x, y, t) is the width of the fracture (in hydraulic fracturing simulators, the width is available from solving the geomechanics problem of fracture growth [7]); V f , V p , V s are the width-averaged velocities of the fluid and the particles, and the settling velocity of particles relative to the fluid; v l is the velocity of fluid leak-off through the porous walls; differential operator '∇' acts in the (x, y) plane as we applied the averaging procedure along the cell width. The flow scales are as follows: L is the cell length, U is the scale of the injection velocity, d is the cell width scale, ρ f 0 is the fracturing fluid density, µ 0 and τ 0 are the fracturing fluid plastic viscosity and yield stress, respectively; g is the gravity acceleration; Bu is the Buoyancy number. Boundary and initial conditions for the hyperbolic equation for concentration (1) are given by: Boundary conditions for the pressure equation (2) are as follows: The particle settling velocity is given by an empirical formula, which is the generalization of the well-known Richardson-Zaki expression, taking into account the fact that particles slow down and stop completely when reaching the packed bed on the bottom. The existing effectivefluid models of suspension flows [8] contain an assumption that the volume-averaged suspension velocity is governed by the Poiseuille law, while in the present two-fluid model it is shown based on the derivation from the conservation laws that the Poiseuille law governs the mass-averaged velocity of the carrier fluid (3). Also, earlier models contained an assumption in the algebraic expression for the particle velocity (3) that the particles settle relative to the volume-averaged velocity of the suspension [8], and not relative to the fluid, as is the case on the two-fluid approach [3]. As a result, in contrast to the existing models the two-fluid model includes an additional term −∇(wCV s ) in the right-hand side of the pressure equation (2), which takes into account the two-speed effects.

Cross-flow inertial migration of particles in a fracture
Particle migration in fractures is essentially a multi-scale problem, spanning from the lift force on a single particle settling in a slot [9], through the inertial igration in a suspension flow in the horizotnal section of a fracture [10,11] to the effects of cross-flow migration on the global transport and sedimentation in the entire fracture [12].
In thi seciton, we will consider in detail the inertial migration of particles in a dilute suspension flow through the entry region of a plane channel, which is important in application to modeling of proppant migration in the near-wellbore zone of the fracture. Within the two-fluid approach, an asymptotic one-way coupling model of the dilute suspension flow in the entry region of a channel is constructed in [10]. The carrier phase is a viscous incompressible Newtonian fluid, and the dispersed phase consists of identical non-colloidal rigid spheres. In the interphase momentum exchange, we take into account the drag force, the virtual mass force, the Archimedes force, and the inertial lift force with a correction factor due to the wall effect and an arbitrary particle slip velocity. The channel Reynolds number is high and the particle-tofluid density ratio is of order unity or significantly larger unity. The solution is constructed using the matched asymptotic expansions method. The problem of finding the far-downstream cross-channel profile of particle number concentration is reduced to solving the equations of the two-phase boundary layer developing on the channel walls. The full Lagrangian approach is used to study the evolution of the cross-flow particle concentration profile.  The system of equations of the two-phase boundary layer are as follows [10]: , Here the latter terms in the right-hand side of the momentum conservation equations for the particulate phase correspond to the virtual mass force and the Archimedes force. Additionally, φ is the Blasius function from the well known solution of the boundary layer problem [13]. External flow is uniform, hence the longitudinal pressure gradient is zero. Boundary conditions take the form:  (1) and small pore channels (2).
The inertial migration in the entry region of a plane channel (a circular pipe) results in particle accumulation on two symmetric planes (an annulus) distanced from the walls, with a non-uniform concentration profile between the planes (inside the annulus) and particle-free layers near the walls. When the particle-to-fluid density ratio is of order unity, an additional local maximum of the particle concentration on inner planes (an inner annulus) is revealed. The inclusion of the corrected lift force makes it possible to resolve the non-integrable singularity in the concentration profile on the wall, which persisted in all previously published solutions for the dilute-suspension flow in a boundary layer [14]. The numerical results are compared with the tubular pinch effect observed in experiments, and a qualitative analogy is found.
In this case the trajectories shown in Fig. 3 are substantially different from the case of a dusty gas. The reason is that in the case of a suspension the particle-to-fluid density ratio is of order unity ξ ∼ 1, and the terms due to the Archimedes force and the virtual mass force should be retained in the momentum conservation laws. These terms result in the formation of an additional local maximum in the cross-flow concentration profile (see Fig. 4,(a)).

Suspension filtration in proppant packings in a closed fracture
Suspension flow in the propped fracture is described within the three-continua approach: suspended partices (solid particles carried by the flow of a fluid within the porous space), trapped particles (solid particles, which are sedimented in pores), and carrier fluid (viscous incompressible Newtonian fluid). Suspended particles are characterised by the phase density ρ mob p and the mass-averaged velocity U mob p ; the phase of trapped particles is characterized by the density ρ sed p ; the carrier phase is characterised by the mass-averaged velocity U f and density ρ f . Particles have constant substance density ρ 0 p , and the fluid -a constant substance density ρ 0 f . Multi-continua modeling of suspension flows is applicable if the following hierarchy of scales exists: particle diamter is significantly smaller than the diameter of pore channels, while being significantly larger than the mean free path of the fluid molecules [2]. Large pore channels are the pore space excluding trapped particles and the void space between those. Small pore channels are pore space between trapped particles (Fig. 5, b).
Density of the particle phases are related to the porosity and the particle concentration by the following relations [15]: where C is the particle volume fraction in the pore space available for the fluid to flow (in large pore channels), σ -the volume concentration of trapped particles in the entire volume of the porous medium, C max is the maximum packing concentration (random close packing), ϕ t -the porosity of the medium formed by the trapped particles and the matrix (large and small channels in total), ϕ c -the porosity of the medium formed by the large pore channels; ϕ 0initial porosity of the medium. The mass conservation equation for particles and fluid is written as [19]: Here j = 0 in a plane flow and j = 1 in a radial flow, q s is the rate of trapping and mobilization of particles in the matrix, determined by the formulas [16]: Here λ is the colmatation coefficient (intensity of particle trapping) [20], α -the coefficient of mobilization of trapped particle, U s is the filtration velocity in large and small pore channels, U crit is the critical suspension velocity, above which trapped particles are mobilized.
In large pore channels with the permeability k(σ), the suspension is flowing with the viscosity µ(C). In small pore channels with the permeability k s (σ), only the clean fluid is flowing with the dynamic viscosity µ 0 . The momentum conservation equations (Darcy laws) for the carrier phase, suspended particles and the suspension as a whole are obtained by volume-averaging in large and small pore channels: Here U is the suspension filtration velocity in large and small pores, U filtr f is the fluid filtration velocity in large and small pore channels, U p = U mob p Cϕ c is the volume-averaged velocity of suspended particles in large pore channels, p -pressure.
The work [17] suggests the following relation between the permeability in large pore channels and the volume fraction of trapped particles: Permeability in small pore channels: Here k s0 is the permeability of small pore channels during full packing of porous space (σ = ϕ 0 C max ), determined from the Kozeny-Carman formula [21]; d is the diameter of particles in suspension.
Using the definitions introduced above, the conservation equations can be reformulated as: Classical deep-bed filtration models of suspension flow do not take into account the fact that the particles trapped in pores form a secondary porous medium (a random close packing of sedimented particles) with the permeability smaller than the permeability of the matrix. Filtration through the close packing of trapped particles is not considered in the classical models.
Below for comparison we present the classical model of deep-bed filtration in large pore channels without considering the flow in small pore channels in the packed bed of trapped particles (ϕ c = ϕ t = ϕ 0 − σ, U = U s , k s = 0) [17,22,15,24,18]: Experiments on colmatation of core sample [23] have shown that the model (13) can be improved by introducing the following expression:λ = λ 0 (1 + βσ). Here, λ 0 is the original colmatation coefficient. Thus, the classical model (13) contains two tuning parameters λ 0 β.
Two versions of boundary conditions are considered: -(i) pressure is specified at the entry and at the outlet of the core sample (a): r = r 0 : p = p 0 , C = C 0 ; r = L : p = 0. -(ii) the filtration velocity is specified at the inlet, and pressure -at the outlet (b): r = r 0 : U s = U 0 , C = C 0 ; r = L : p = 0.
Initial conditions: t = t 0 : ϕ = ϕ 0 , σ = 0. Permeability of mixed proppant pakcings from LBM simulations. Closure relations for permeability of the random close packing of proppant in a hydraulic fracture can be found from direct numerical simulation supported by proper conductivity experiments (see, for example, LBM simulations in mixed packing of proppant particles of various shape in [25]). Granular materials are important in many applications in petroleum and water supply industries. For example, during hydraulic fracturing, proppant particles are injected into a fracture to avoid its complete closure. In gravel packing, particulate material is injected into a well and placed in the borehole between a screen and the formation to prevent sand production from the reservoir. Granular mixtures of particles are also used in various other branches of industry, particularly in chemical and civil engineering. Because of their industrial importance, the properties of such materials have been the subject of extensive studies. The proppant packs are generated using the sequential deposition algorithm (Coelho et al., 1997). The pack is created by deposing particles into a box with a rectangular cross-section wL, where w is the fracture width (Fig. 6) and H ≫ L. The length of the pack L is chosen equal to 6.4l max , where l max is the maximum semi-axis of the particle. The binary samples are periodical in the longitudinal direction, while in the transverse direction they are either periodical or have flat walls. For the pack located between impermeable walls, the ratio of the fracture width to the proppant size w/l max is in the range from 1.6 to 6.4. Using the lattice Boltzmann method developed in (Keehm, 2003) for slow single-phase flows in porous media, we conducted a series of numerical simulations to obtain the permeability tensors of different proppant packs. In our simulations, we consider the three-dimensional steady-state flow of a viscous incompressible Newtonian fluid through a granular packing. A series of numerical simulations is then carried out to obtain the permeability tensors for various random packs. As an example, permeability tensor K/100D for the pack of ellipsoids (ϕ = 0.65): 41.09 -0.54 -0.07 -0.53 42.34 0.41 -0.05 0.41 43. 52 We considered separately the case of a spatially periodical pack and the case of a pack confined between two impermeable walls. The results of numerical simulations demonstrate that the pack of cylinders provides the largest porosity and the highest permeability, while the pack of spheres has the lowest porosity and the smallest permeability. The porosity and permeability of the pack of ellipsoids are intermediate. All permeability tensors obtained in simulations are fairly isotropic. All numerical results are well approximated by the non-dimensional permeabilityporosity relation, where permeability is scaled by the equivalent radius squared: A comparison with lab data is shown in Fig. 7. More details on the topic can be found in the paper [25].

Gas-liquid flows in a well during cleanup
In this section we will discuss the first-principles derivation of simplified drift-flux equations for gas-liquid flows in a well [26]. We consider a transient isothermal flow of a gas-liquid mixture in a long circular pipe with a variable inclination angle to the horizon. The flow is assumed to be axisymmetric and non-swirling. The liquid is a continuous carrier phase. The gas is a dispersed phase present in the form of monodisperse spherical bubbles suspended in the carrier fluid. The gas is compressible, and the liquid is incompressible. The cross-flow migration of the bubbles and their merging are not considered; however, we take into account the resulting nonuniform cross-flow profile of the gas volume fraction, which is formed as a result of the migration of the bubbles. The pressure difference in the bubbles and in the surrounding liquid due to surface tension is neglected. The bubble size is assumed to be much smaller than the spatial scales of variation of the fluid velocity, and the Reynolds numbers of the flow around individual bubbles are small. The two-phase flow is considered within the approach of interpenetrating and interacting continua [2]. The problem is described by conservation laws in differential form, written for gas and liquid phases. The mass transfer between the phases is absent.
The balance laws of mass and momentum in the differential form for gas and liquid take the form [2]: Here, the indexes i, j = g, l (i ̸ = j) denote gas and liquid, respectively, α i , ρ i , and v i are the volume fractions, densities and velocities of the phases, p i and τ i are the pressures and viscous stress tensors in each phase, and g is the gravity force acceleration. The momentum exchange between the phases is described by the term ±n b F ij , where F gl = F is the force exerted on a single bubble by the fluid, F lg = −F lg , and n b is the number concentration of the bubbles. For simplicity, the following calculations are performed for a vertical pipe, although the results may be generalized to the case of an inclined pipe, excluding flows in near-horizontal pipes.
It is assumed that the chaotic motion of the bubbles can be neglected, and the deviation of the bubble velocity from the mass-averaged velocity of the dispersed phase v g is small, hence the stress tensor in the dispersed phase can be neglected [3]. The presence of the dispersed phase affects the stress tensor of the carrier phase. On the other hand, the bubbles of a compressible gas travel with a velocity different from that of the fluid, and the bubble volume fraction varies. Accordingly, the condition ∇v l = 0 is no longer true. In this sense, the averaged liquid phase is compressible, in contrast to the liquid as a material. Therefore, the stress tensor of the carrier phase is written as for a viscous compressible fluid, with the coefficients of shear viscosity µ and bulk viscosity ζ dependent on the gas volume fraction: Here, e l is the strain rate tensor, and I is the unit tensor. To determine the dependences µ (α g ) and ζ (α g ) is a separate problem, which is usually solved for neutrally-buoyant particles, neglecting the phase slip. In what follows, we assume that µ (0) = µ 0 , where µ 0 is the viscosity of the pure fluid, and ζ (α g ) ∇v l → 0 as α g → 0.
The bubble radius is bounded by a limiting value R c above which the bubble loses the spherical shape, becomes unstable, and is fragmented into smaller bubbles: Here, σ is the gas-liquid surface tension. Under the above assumptions, the total force exerted on a bubble by the fluid is a superposition of the Stokes F s , the Archimedes F A , the added mass F am , and the Basset-Boussinesq F BB forces. Accordingly, the forces on a bubble can be written in the form: For a transient flow in a well, the equations of the drift-flux model, known in the literature and implemented in the commercial reservoir simulator ECLIPSE (Schlumberger), take the form [27,28]: Here, A is the pipe cross-section, v m = α g v g + α l v l is the volume-averaged velocity of the mixture, ρ m = α g ρ g + α l ρ l is the mixture density, f = f(α g , v m , p) is the friction coefficient, d is the pipe diameter, and θ is the angle between the pipe axis and the vertical.
In the applications, a quasi-steady-state variant of the drift-flux model is widely used, in which the time-derivative in the momentum equation for the mixture is neglected and the total pressure difference is equal to the sum of the terms responsible for the gravity force, friction, and acceleration [27]. In the literature, this model is referred to as the no-pressure wave model, since it does not take into account the disturbance propagation with the transport velocity. In contrast to this quasi-steady-state formulation, we retain the time derivative of velocity in the momentum equation to take into account highly transient effects.
In the literature, formula (18) [2] (where C 0 = C 0 (α g , v m , p) is the profile parameter which takes into account a non-uniform cross-flow profile of the bubble volume fraction and the carrierphase velocity, and v d = v d (α g , v m , p) is the drift velocity) is called the drift-flux model relation.
There is also another known formulation of the drift-flux model [29], with the mixture momentum equation written as ∂ ∂t (α g ρ g v g + α l ρ l v l ) + ∂ ∂z (α g ρ g v 2 g + α l ρ l v 2 l + p) = Q l + Q g , where Q i are the source terms for each phase, and the drift-flux relation written as The asymptotic equations are derived in the long-channel approximation ε ≪ 1, similar to the boundary-layer approximation, the narrow-channel approximation for fracturing flows [3], and the lubrication approximation for thin-film flows [30].
Our analysis demonstrates that the drift-flux model [27,28] in the form (16)- (18) for the present flow configuration strictly follows from the conservation laws in the limited number of cases (i)-(iii) and represents essentially the effective-fluid model. The closure relations published in open literature are obtained from a calibration against a large body of experimental data [27,28] for the governing parameters satisfying at least one of conditions (i)-(iii). At the same time, our analysis indicates that the drift-flux model [29] in the form (16), (19), and (20) is more general because it follows from the balance laws without any additional assumptions, besides the requirement of noninertial slip εSt ≪ 1.

Conclusions
We present a thorough review of a family of closely-related multi-scale models of multiphase flow constructed to describe all stages of the technology of hydraulic fracturing. These models were derived by the author from conservation laws within the multi-fluid approach using perturbation methods. The models cover suspension injection into the well to create primary fractures, suspension transport and sedimentation in a fracture, cross-flow inertial migration of particles in a horizontal section of a hydraulic fracture, filtration of suspensions of non-colloidal particles in a random close packing of proppant particles in the closed fracture, and gas-liquid flows in a well during cleanup and startup. In particular, our review covers: 1. A novel two-fluid model of suspension transport in a hydraulic fracture. It is shown that the model is different from the existing effective-fluid model of suspension flow by additional terms due to two-speed effects. These terms are important at high buoyancy numbers Bu, which corresponds to the case of high-rate slick-water fracturing. In the case of conventional fracturing with high-viscosity cross-linked gels, the models match. This model is then generalized to the case of a Bingham fluid flow and validated against four different sets of experiments on slumping, Saffman-Taylor fingering in Newtonian fluids, fingering and channeling in yield-stress fluids, and