Numerical modelling of proppant transport in hydraulic fractures

The distribution of proppant injected in hydraulic fractures significantly affects the fracture conductivity and well performance. The proppant transport in thin fracturing fluid used during hydraulic fracturing in the unconventional reservoirs is considerably different from fracturing fluids in the conventional reservoir due to the very low viscosity and quick deposition of the proppants. This paper presents the development of a threedimensional Computational Fluid Dynamics (CFD) modelling technique for the prediction of proppant-fluid multiphase flow in hydraulic fractures. The proposed model also simulates the fluid leak-off behaviour from the fracture wall. The Euler-Granular and CFD-Discrete Element Method (CFD-DEM) multiphase modelling approach has been applied, and the equations defining the fluid-proppant and inter-proppant interaction have been solved using the finite volume technique. The proppant transport in hydraulic fractures has been studied comprehensively, and the computational modelling results of proppant distribution and other flow properties are in good agreement with the published experimental study. The parametric study is performed to investigate the effect of variation in proppant size, fluid viscosity and fracture width on the proppant transport. Smaller proppants can be injected early, followed by larger proppants to maintain high propping efficiency. This study has enhanced the understanding of the complex flow phenomenon between proppant and fracturing fluid and can play a vital role in hydraulic fracturing design.


Introduction
In the last decade, the advancements in horizontal drilling and multistage hydraulic fracturing have resulted in considerable progress in the hydrocarbon production [Warpinski, Mayerhofer, Vincent et al. (2009)]. Both these techniques are closely related to geomechanics, i.e., in order to create a multi transverse hydraulic fracture, wells are drilled horizontally in the direction of minimum horizontal in-situ stress. One conventional method of generating multiple fractures is the plug and perf method, in which multiple fractures are created in stages, and each stage consists of a cluster of perforations ]. The highly pressurised fluid is injected at sufficiently high rates to initiate and propagate the fractures. It is followed by the injection of proppant laden fluid, to ensure that the fracture remains open against the geomechanical stresses when the fracturing fluid pressure is reduced [Economides and Nolte (2000); Gaurav, Dao and Mohanty (2010)]. The proppant distribution inside the fracture aids to maintain a conductive path for the flow of hydrocarbon from surrounding tight rock [Sharma, Sircar and Gupta (2018)]. In the hydraulic fracturing of unconventional reservoirs, slick water is mainly used as the fracturing fluid. Unlike conventional fracturing fluids, due to the very low viscosity of slick water and negligible chemical additives, the tendency of suspending the proppant particles dramatically decreases [Sahai, Miskimins and Olson (2014)]. This results in early deposition of the proppants compare with the fracturing fluids in conventional reservoirs [Alotaibi and Miskimins (2015)]. The hydraulic fracturing in an unconventional reservoir is considered successful when the long multiple hydraulic fractures are created with uniform proppant distribution resulting in the flow of hydrocarbon fluids economically [Gu and Mohanty (2014)]. Proppant transport plays a vital role in hydraulic fracturing. After injection, proppant particles are dispersed in the slurry making it a multiphase flow problem [Kostenuk and Browne (2010)]. The proppant follows the fluid path, momentum exchange occurs between the two phases and the proppant settles down away from the wellbore due to the gravity forces. This results in two primary physical mechanisms, namely fluid-proppant interaction and proppantproppant or proppant-wall interaction [Daneshy (2011)]. Proppant distribution inside the fractures is a complex phenomenon, and there are many factors affecting the proppant transport process, including properties of proppant and fracturing fluid, fracture propagation, and fluid leak-off [Gadde, Liu, Norman et al. (2004)]. Many experiments have been reported in the literature to investigate the proppant transport and distribution in laboratory-scale fracture slots [Boronin and Osiptsov (2010); Kostenuk and Browne (2010); Sahai, Miskimins and Olson (2014)]. Some of the earliest experiment was performed by Kern et al. [Kern, Perkins and Wyant (1959)] and Baree et al. [Barree and Conway (1994)], which mainly focussed on evaluating the effect of proppant and fracturing fluid properties on proppant distribution. Later, Liu et al. [Liu and Sharma (2005)], and Woodworth et al. [Woodworth and Miskimins (2007)] extended the proppant transport study to the detailed analysis of fracture geometry and large-scale laboratory experiment. More recently, Sahai et al. [Sahai, Miskimins and Olson (2014)], Alotaibi et al. [Alotaibi and Miskimins (2015)], Tong et al. [Tong, Singh and Mohanty (2017)] studied proppant distribution in complex fracture geometries, fracture network, slickwater and foam fracturing fluids. The experimental studies reported in the literature, although provides a good understanding of the variation of proppant properties on proppant deposition behaviour, but due to the laboratory limitations and assumptions, the results can only be qualitatively used. Hence, in addition to the experimental work, numerical simulation methods can aid in a detailed understanding of multiphase flow physics. Two fundamental methods to numerically solve the proppant transport equations in fracturing fluid flow can be categorised into Eulerian-Lagrangian method and the Eulerian-Eulerian method [Kong, McAndrew and Cisternas (2016); Tsai, Fonseca, Lake et al. (2012)]. The advantages and limitations of both these models are discussed as follows.
In the Eulerian-Lagrangian method, the equations for the primary phase or the continuous phase is solved using mass and momentum conservation equations. The equation for dispersed phase or proppants is solved using Newton's second law of motion by tracking their motion ]. Thus tracking the motion of individual proppant in fracturing fluid slurry with inter-particle interaction forces result in this technique providing accurate proppant distribution results with the cost of being computationally very demanding [Zhang, Gutierrez and Li (2017)]. This limitation impedes the Eulerian-Lagrangian method for use in industrial-scale fractures. The Eulerian-Lagrangian technique can be broadly grouped into two categories based on proppant volume fraction in the slurry, causing inter-proppant and proppant-fluid interaction. The two categories are the Discrete Particle Method (DPM) and Computational Fluid Dynamics-Discrete Element Method (CFD-DEM). The DPM model is mainly used when the proppant volume fraction in the slurry is low (approximately<10%) so that the inter-proppant interaction forces are not dominant and can be neglected ]. Conversely, the CFD-DEM model uses a soft-sphere approach to capture the detailed proppant-proppant and proppant-wall interaction, which makes them optimal for higher volume fraction proppant slurry flow. As mentioned earlier, the detailed proppant physics involved in the CFD-DEM model demands considerable computational time and makes it unattractive for using it in industrial-scale fractures [Deng, Li, Ma et al. (2014); Patankar and Joseph (2001); Snider (2001); Wu and Sharma (2016)]. The alternative methods for solving proppant transport equations widely used in the literature is the Eulerian-Eulerian method or Eulerian-Granular method. The motion of both the phases, i.e., continuous phase and dispersed solid phase in solved by mass and momentum conservation equations [Clifton and Wang (1988); Liu (2006); Kong, McAndrew and Cisternas (2016); Li, Zhang and Lu (2018); Roostaei, Nouri, Fattahpour et al. (2018)]. The interaction between the phases is modelled using the Kinetic Theory of Granular Flow (KTGF) [Clifton and Wang (1988)]. The Eulerian-Granular method offers computationally faster results, and some researchers have demonstrated that it provides reasonable results, however the particle-wall interaction is not very well captured that immensely impacts the proppant flow dynamics [Clifton and Wang (1988); Liu (2006); Kong, McAndrew and Cisternas (2016); Li, Zhang and Lu (2018); Roostaei, Nouri, Fattahpour et al. (2018)]. For numerical modelling of proppant distribution in hydraulic fractures much work has been done [Patankar and Joseph (2001); Gadde, Liu, Norman et al. (2004); Liu (2006); Wu and Sharma (2016); Kong, McAndrew and Cisternas (2016); Zhang, Gutierrez and Li (2017); Roostaei, Nouri, Fattahpour et al. (2018)]. However, the existing models ignore the fluid leak-off from the fracture wall in the proppant transport equations for simplifying the model. The amount of proppant suspended in the slurry is highly governed by the carrying fluid and the rate of fluid leak-off from the hydraulic fracture to the surrounding porous rock. As soon as the fracturing fluid slurry is injected in the hydraulic fracture, it leaks into the surrounding rock depending upon the reservoir characteristics and leaves the remaining proppants in the slurry that have an inclination to deposit and form a proppant bed. Ignoring the fluid leak-off effects cannot provide the true estimate of the proppant distribution. To the best of our knowledge, the existing literature studies mentioned earlier have mainly solved the proppant dynamics in planar fracture geometry with an assumption of no fluid leak-off from the fracture wall. In this paper, the current work aims to develop a threedimensional (3D) Computational Fluid Dynamics (CFD) model to include the fluid leakoff rate from the fracture walls and model the proppant micro-mechanics to numerically investigate the proppant transport in hydraulic fractures using the advanced modelling techniques namely Eulerian-Granular method and Eulerian-Lagrangian method. Firstly, a proppant transport model was developed and compared with the experimental results. It was followed by a base case that models the fluid leak-off from the fracture walls and proppant transport in a planar fracture. Subsequently, a parametric study was performed to critically examine the role of proppant properties (proppant size), fluid properties (fluid viscosity) and geomechanical properties (fracture width) on proppant transport. In addition, the Eulerian-Granular model has been compared with the more accurate Discrete Element Model for proppant flow.

Methodology
In the present study, two different numerical modelling techniques are used to study proppant transport and distribution in hydraulic fractures, namely Eulerian-Granular model and Computational Fluid Dynamics-Discrete Element Model (CFD-DEM), described in the following sections. The key objective in the present study is to provide a detailed understanding of the proppant transport considering the effect of fluid leak-off from the fracture wall in a planar fracture in the unconventional reservoir. Some of the assumptions underlying the current model are as follows: First, the model is a small scale because of the large simulation time in the CFD-DEM model. Second, the slurry is Newtonian fluid. Third, the fracture geometry is constant and assumed as a cuboid; no dynamic fracture propagation is considered in this study.

Two-fluid model or eulerian-granular model
The Two-Fluid model or Eulerian-Granular model is a multiphase flow model in which both phases are defined as a continuous phase. This means the flow governing equations (continuity and momentum equations) are solved separately for each phase. The primary phase is fluid, and the second phase is defined as granular phase (solid phase). The particleparticle collision or inter-particle interaction is explicitly modelled using a collision model, kinetic theory of granular flows (KTGF) and frictional models [Kong, McAndrew and Cisternas (2016); Reuge, Cadoret, Coufort-Saudejaud et al. (2008)]. The particle-fluid interaction is defined by interphase exchange coefficients and is modelled using the empirical models [Burns, Frank, Hamill et al. (2004); Kong, McAndrew and Cisternas (2016); Reuge, Cadoret, Coufort-Saudejaud et al. (2008)]. The governing momentum equation for the granular phase includes additional terms to define the properties for granular flow such as solid pressure and solid stress tensor terms from the application of the kinetic theory of granular flows (KTGF) [Jenkins and Savage (1983); Savage and Jeffrey (1981)]. A key parameter in the KTGF models for solids phase stress is a parameter known as granular temperature. The granular temperature provides a measure of the kinetic energy associated with solid particles velocity fluctuations. The granular temperature is a function of the fluctuating velocity of the particles and is obtained using an algebraic equation described in detail in the next section. One of the underlying assumptions in this method is that the position of each phase is a function of phase volume fraction (interpenetrating continua).

Flow governing equations
The governing equations that refer to the equation of conservation of mass and momentum for granular-liquid coupled flow are described below. The equations are based on an assumption of isothermal and incompressible condition for the fracturing fluid. The detailed derivation of these equations can be found in Banerjee et al. [Banerjee and Chan (1980)], Versteeg et al. [Versteeg and Malalasekera (2007)] and Jakobsen [Jakobsen (2014)].
The mass conservation equation is given by: where α represents volume fraction, ρ refers to the density, v refers to velocity, refers to mass source term and subscript i refers to phase (liquid or solid) For the fracturing fluid, the conservation of momentum equation is given by: where g refers to acceleration due to gravity, M ls ������⃗ = M sl ������⃗ refers to the interfacial momentum exchange between fluid and granular phase, subscript and refers to liquid and granular phase respectively, refers to the momentum source term and τ l � is the fluid phase stressstrain tensor given by: where λ l and µ l refer to the bulk viscosity and dynamic viscosity of continuous phase (fracturing fluid) respectively. For the granular phase, the conservation of momentum equation is given by: where p s refers to granular phase pressure and τ s � refers to the stress-strain tensor for granular phase. The interfacial momentum exchange term and the stress terms are described in detail below.

Interfacial momentum exchange
Burns et al. [Burns, Frank, Hamill et al. (2004)] described that due to the fracturing fluid velocity variation, the motion of granular phase is changed due to the transfer of momentum between the two phases. It is given by Eq. (6) which is a combination of the drag force F �⃗ drag , lift force F l ���⃗ and virtual mass force F vm �������⃗ .
The lift force is because of the fluid velocity gradient. The virtual mass force refers to the 304 FDMP, vol.16, no.2, pp.297-337, 2020 force required to accelerate the fluid surrounding the particle [Sankaranarayanan, Shan, Kevrekidis et al. (2002)]. Issa et al. [Issa and Oliveira (1993)] and Ekambara et al. [Ekambara, Sanders, Nandakumar et al. (2009)] concluded in their study that these forces are negligible when the ratio of granular density to fluid density is greater than 1. Thus, only the effect of the drag force is accounted for in this study.

Drag force modelling
The drag force is described by the Eq. (7). Numerous drag force models are available for multiphase flow modelling that differs in the definition of inter-phase momentum exchange coefficient, K ls or K sl .
v l ���⃗ − v s ���⃗ is the relative velocity between the phases. Gidaspow [Gidaspow (1994)] proposed a drag force model which provides the flexibility to use it for a wider application range based on the proppant volume fraction. Gidaspow drag model is used in the present study as described by Eq. (8): where d s represents the granular phase diameter and C D refers to the drag coefficient and calculated by Eq. (9).
where Re s refers to the Reynolds number of the granular phase and calculated by:

Stresses model for the proppant phase
Savage et al. [Savage and Jeffrey (1981)] described that the solid stress for the granular phase, τ s � (in Eq. (5)) is based on the kinetic theory of granular flow (KTGF) models as expressed in Eq. (11) where λ s and µ s refer to the bulk viscosity and dynamic viscosity of the granular phase respectively and I ̿ is the unit tensor.

Granular temperature
As mentioned earlier, the granular temperature provides a measure of the kinetic energy associated with solid particles velocity fluctuations. The granular temperature is a function of the fluctuating velocity of the particles as described in Eq. (12). It can be obtained using Eq. (13) or granular energy transport equation.
where Θ , , and Φ refers to the granular temperature, velocity fluctuation of granular phase, and granular energy transfer respectively. Θ , Θ and α denotes granular energy dissipation rate, diffusion coefficient and granular phase volume fraction respectively. Alternatively, the granular temperature can be obtained using an algebraic expression in Eq. (14) proposed by Wachem et al. [Wachem, Schouten, Bleek et al. (2001)] that assumes a steady-state solution of the granular energy neglecting the convection and diffusion terms.
2.2.5 Granular phase pressure model Granular phase pressure, P s , is a function of normal force due to particles motion and can be calculated using a correlation from Lun et al. [Lun, Savage, Jeffrey et al. (1984)] given by Eq. (15).
where refers to the restitution coefficient from particles collision. The restitution coefficient varies from 0 to 1 for a perfectly elastic collision to a perfectly inelastic collision. The value of the restitution coefficient assumed in the current study is 0.9. The detailed explanation of the granular phase pressure model can be found in Suri et al. [Suri, Islam and Hossain (2019)].
where refers to the friction angle and refers to the friction pressure. The detailed explanation of granular shear viscosity model can be found in Suri et al. [Suri, Islam and Hossain (2019)].

Computational fluid dynamics-discrete element model (CFD-DEM)
CFD-DEM is based on Eulerian-Lagrangian method as explained earlier. Unlike other Eulerian-Lagrangian methods, for instance, Discrete Phase Model (DPM) which is applicable only for the low volume fraction of particles (<10%), the Discrete Element Method (DEM) can be used when the higher volume fraction of particles is present (>10%). Thus, CFD-DEM can accurately model the multiphase flow where the inter-particle interaction is imperative, such as proppant flow in the fracturing fluid. Cundall et al. [Cundall and Strack (1979)] proposed the DEM method, and it was later coupled with CFD by other researchers to study fluid-solid flow modelling ]. In this approach, the primary phase is solved using a conventional Eulerian method meaning continuity and momentum equations are solved using CFD. However, the solid phase is solved using DEM by tracking every dispersed particle, thus it is a computationally expensive technique. Particles are tracked by calculating and tracking the mass, velocity, and forces acting on a particle using Newton's second law of motion. This is referred to as tracking in the Lagrangian frame in the DEM method [Zhang, Gutierrez and Li (2017)]. Finally, the drag forces and interphase momentum exchange terms are used to model the interaction, energy dissipation and coupling of both the phases, i.e., continuous and discrete phases. The method is based on the following two assumptions: 1) Particles are considered rigid. However, particle-particle deformation is present. A simple force-displacement law is used to take into account particle overlapping and particle-particle deformation.
2) Newton's second law of motion is used to predict each particle motion. In order to account for accurate particle micro-mechanics and particle collision, it is assumed that after the collision, the two particles deform and defined by the overlap displacement of the particles. This approach is called the soft-sphere approach that outlined an accurate contact model and is explained later in the following section.

The governing equations for the particles
The distribution of discrete phase particle motion is calculated by integrating the force balance on the particle, which is written in a Lagrangian reference frame. Using Newton's second law of motion, the governing equations of the particle motion can be defined as follows: The above equations can be re-written in the following form as The velocity and spatial location of discrete particles are calculated using Eqs. (20)-(21) respectively. The term F �⃗ other refers to other forces such as forces of the collision, cohesion, electrostatic forces, lift forces, magnetic forces and virtual mass forces. The collision model is described in the next section and the forces of cohesion, electrostatic forces, lift forces, magnetic forces and virtual mass forces are not considered in the present study as explained earlier. The variable τ r is the droplet or particle relaxation time is given by: is the drag force per unit particle mass, v l ���⃗ and v p ����⃗ is the fluid and particle velocity respectively, µ is the fluid viscosity, ρ and ρ p are the fluid and particle density respectively, d p is the particle diameter, and Re is the Reynolds number, defined as: The Navier-Stokes equations (mass and momentum conservation equations) of the continuous phase are described below: where S m and refers to the mass and momentum source term respectively, and it takes into account the particle motion, F �⃗ represents external body force term, τ � represents stress tensor, g is the acceleration due to gravity and ρ is the density. In the CFD-DEM method to ensure the numerical stability and converged solution, usually the time step for discrete phase DEM modelling is smaller than continuous phase CFD modelling. This is done to capture the particle micro-mechanics correctly. In the present study, the time step used for the CFD continuous phase flow simulations is 1.0e-3 seconds, and for the discrete phase, DEM simulation is 1.0e-6 seconds, which is three order of magnitude lower. The DEM time step was selected based on the numerical stability and convergence of the solution using the Eqs. (31)-(32) and Eq. (35) described in the next section. Cundall et al. [Cundall and Strack (1979)] proposed the "soft sphere" approach in order to model the collision forces of granular phases in the DEM method. F �⃗ other term in Eq. (20) and Eq. (22) accounts for these forces. The granular collision forces are calculated by the deformation, resulting from the overlap between pairs of spheres (Fig. 1). The springdashpot collision model is used in the present study for modelling inter-particle collision. The force exerted on the proppant due to the collision by another proppant particle is given by Eq. (27) 1 ���⃗ = ( + ( ⃗ 12 . ⃗ 12 )) ⃗ 12 (27) By Newton's third law of motion, the force on the second proppant particle is given by 2 ����⃗ = − 1 ���⃗ (28) where ⃗ 12 , is the unit vector, γ is the damping coefficient, represents the overlap, represents the spring constant, ⃗ 12 is the relative velocity [Issa and Oliveira (1993)] and is given by

Collision model
= � 2 + 2 (32) is the coefficient of restitution for particles collision, which can vary from 0 to 1 corresponding to from perfectly inelastic to a perfectly elastic collision. The inelastic collision with =0.9 is used in this study. is the spring constant and can be estimated by the Eq. (34).
The spring constant value of 1000 is used in the present study. The particle timestep can be estimated from Eq. (35) to get an accurately resolved collision. Thus, the particle time step used was 1.0e-6.

Numerical modelling parameters
Proppant transport and distribution were investigated in a hydraulic fracture using the CFD technique in ANSYS FLUENT. The geometry or computational domain used in the current study is, as shown in Fig. 2, with dimensions 1.5 m (length)×0.5 m (height)×0.05 m (width). In order to obtain a mesh independent solution, a mesh sensitivity analysis was done with by varying the mesh sizing parameter 0.0025 m, 0.005 m, 0.0075 m, and 0.01 m. The results of the mesh sensitivity analysis were compared against the proppant volume fraction and proppant axial velocity vs fracture height Fig. 3. The results from the mesh sensitivity study, suggest that the mesh size of 0.005 m (300×100×10 elements) was reasonably able to provide the mesh independent, numerically converged and computationally efficient solution. Next, appropriate boundary conditions, cell zone conditions, and simulation properties were defined. A velocity inlet boundary condition is used at the inlet where fluid and proppants are injected at 0.5 m/s. The particle size distribution is assumed to be of uniform diameter 1 mm. All the walls shown in Fig. 2 are assumed as no-slip stationary walls. In order to mimic the fluid leak-off into the surrounding porous rock, the fluid leakage effect is modelled through the fracture sidewalls with the help of a user-defined function (UDF). The momentum and mass source terms are defined and included in the governing transport equations through UDF been written in C++, which is interpreted by the CFD solver, ANSYS FLUENT 18.1. In order to obtain the fluid leakage rate, an explicit CFD study was carried out to calculate the water leaking off rate along the fracture sidewall. The underlying equations describing the source terms and UDF used to model the fluid leakoff is explained in detail in Suri et al. [Suri, Islam and Hossain (2019)]. The fluid leak-off profile along the fracture length to a surrounding porous rock with porosity 5% and permeability 1 mD used in the current study is shown in Fig. 4.

Figure 4:
Fluid leak-off rate at fracture wall along the fracture length The pressure-based solver with gravitational effects was used to solve the governing proppant transport equations as described earlier. In order to model the turbulence, the Shear Stress Transport (SST) k-ω model [Menter (1993)] was used due to its applicability in adverse pressure gradients and separating flow [Versteeg and Malalasekera (2007)]. Lastly, the transient phenomenon was used to investigate the variation of proppant transport and distribution with time. Two different models were studied for proppant transport in fracture-Eulerian-Granular model ( Thus, in the current study, the proppant volume fraction of 20% is used to model the inter-proppant interaction. The viscosity of the granular phase is calculated from the Gidaspow et al. [Gidaspow, Bezburuah and Ding (1991)] correlation. The primary role of granular viscosity is used to consider the frictional losses. The frictional viscosity refers to the shear viscosity based on the viscous-plastic flow and is calculated using the Johnson et al. [Johnson and Jackson (1987)] correlation. The packing limit defines the maximum volume fraction of the Percetange of injected water mass leak-off rate Non-dimensional distance along fracture length granular phase. For the uniform proppant size, this value is equal to 0.63. Friction packing limit refers to a threshold volume fraction at which the frictional regime becomes dominant, and friction packing limit of 0.6 is used.
In the Eulerian-Granular method, the drag force used to model the interaction between the two phases is based on Gidaspow drag law [Gidaspow (1994)]. The collision between the proppant particles is modelled using the restitution coefficient, as explained earlier. On the other hand, for CFD-DEM modelling, the time step used for continuous phase CFD modelling was 1.0e-3 s, and for discrete phase was 1.0e-6 s, which is three orders of magnitude higher. Thus, this approach is computationally expensive. The time step was selected based on the numerical stability and convergence of the solution using the Eqs. (33)-(34) and Eq. (37) described earlier. To accurately model the inter-particle collision, DEM collision with the spring-dashpot model was used in the normal and tangential direction as explained earlier. The reflect DPM boundary condition used at walls so that the particles will reflect after the collision with the wall. Finally, the Phase-coupled SIMPLE algorithm is used as a solution method for pressurevelocity coupling [Versteeg and Malalasekera (2007)]. The discretisation of momentum, volume fraction, and turbulent kinetic energy was solved by the second-order upwind scheme.

Comparison with the experimental results
The present simulation model was compared against the experimental study of Tong et al. [Tong and Mohanty (2016)]. The simulation was performed with the geometry similar to the experimental setup. All the modelling parameters are presented in Tab. 2, which are similar to experimental parameters. Eulerian-Granular multi-phase flow model was used. Fracturing fluid (water, in this case) along with the proppant is injected at the inlet. Fig. 5 shows a comparison of experimental and simulation results at three different time periods after the start of injection. The results between numerical modelling study and experimental study were compared by calculating the non-dimensional proppant equilibrium height and non-dimensional proppant bed length at the centre and plotted in Fig. 6. The results from Figs. 5 and 6 suggests a reasonable match between the experimental study and the current model. The average error calculated is 5.8% and 6.2% for nondimensional proppant equilibrium height and non-dimensional proppant bed length respectively, which suggests an overall good match and the simulation model can be used to perform further analysis of proppant distribution in the slickwater fracturing fluid.

Results from fluid leak-off modelling
Fluid leak-off is one of the critical phenomena that govern the proppant suspension in the slurry. As the fracturing fluid slurry is injected in the fracture, the fracturing fluid leaks off from the fracture wall to the surrounding porous rock at a rate depending upon the reservoir characteristics. The remaining proppants in the slurry have a tendency to deposit and form proppant bed at the fracture bottom. The higher leak-off rate can result in a greater flow of thin fracturing fluid to the surrounding reservoir rock, leaving behind the proppant in the remaining slurry and consequently early deposition of the proppants. The fluid leak-off depends on the reservoir characteristics (porosity and permeability). The simplest model to take into account the fluid leak-off is defined by Carter [Carter (1957)], which describes the leak-off rate depends on a mathematical constant and elapsed time. Leak off effects plays a vital role in shale reservoirs where due to the use of thin fracturing fluid, the ability to suspend the proppants is considerably low. Furthermore, greater fluid leak-off from the Experimental study Current Numerical Study fracture wall will increase the rate of proppant bed formation and early fracture tip screen out. The fracture tip screen out will then inhibit any further proppant transport into the fracture, and the unpropped section of the fracture will close down, resulting in loss of fracture conductivity. The fracturing fluid leak-off from the fracture wall is ignored in the existing numerical proppant transport modelling studies, as explained earlier, resulting in inaccurate flow and transport properties of proppants. The effect of fluid leak-off from the fracture wall in the proppant distribution was investigated by comparing it with a simulation case with no fluid leak-off effect. The results were compared based on the variation of proppant volume fraction and proppant horizontal velocity, as shown in Figs. 7 and 8. The proppant volume fraction and proppant horizontal velocity were calculated at two different longitudinal locations from the inlet x=0.25 m and x=1.2 m and results were plotted against the fracture height at t=2 s (Figs. 7 and 8). The results show that as the fluid leaks off the fracture wall, leaving the proppants in the remaining slurry, it increases the tendency for the proppants to settle at fracture bottom and forms a bed. Thus, comparatively greater proppant bed height is noticed in the fluid leakoff case against the no leak-off the case. Furthermore, comparing the proppant horizontal velocity suggests that the lower value of the horizontal velocity of proppants is noticed in the fluid leak-off case against the no fluid leak-off case. This can be explained by the greater tendency of the proppants to settle down at the fracture bottom and thus to have lower horizontal transport velocity. The comparison study of fluid leak-off effect with no fluid leak-off effect suggests that the proppant bed height will be under predicted by 10-50% if the leak-off effects are ignored in the proppant transport model and can significantly impact the hydraulic fracture design.

Effect of proppant size
The proppant size was varied, keeping all the other parameters constant, and simulation run was performed. The three cases of variation in proppant diameter studied are diameter=0.3 mm, 0.5 mm and 1 mm. Figs. 9 and 11 are the contour plots of proppant volume fraction and proppant horizontal velocity respectively at fracture mid-plane for different time step and all the three cases of variation in proppant sizes. It shows the difference in proppant distribution inside the fracture with time. It can be interpreted from Fig. 9 that greater particle deposition is noticed for proppants with greater size, or in other words, the greater size proppants tends to settle more quickly. This is due to as the proppant size increases, it increases vertical settling velocity, thus as the particles get larger, the tendency for deposition increases. On the other hand, the smaller size proppants have a lower settling velocity in the vertical direction and occupy greater volume in the suspending region. Next, the proppant volume fraction and proppant horizontal velocity were calculated at two different locations in the longitudinal direction from the inlet at 0.25 m and 0.8 m. The results were plotted with the fracture height to investigate the advancement of proppant volume fraction with time ( Figs. 10 and 12). The results from Fig. 10 show that the proppant volume fraction is identical for all the cases at the beginning, but later with time, the smaller size proppant particles (0.3 mm and 0.5 mm) are more suspended and fill a larger volume of the fracture, while the bigger size particles show greater deposition. The results from Fig. 12 show that the proppant horizontal velocity profiles initially at t=0.5 s is similar for all the three cases. However, later with the progression of time, close to the wellbore (x=0.25 m from inlet), the greater size proppants are moving with higher velocity. However, away from the wellbore (x=0.8 m from inlet), the greater size proppants appear to move slowly compared to the smaller size proppants. This velocity lag of greater size proppants away from the wellbore can be attributed to the higher drag forces experienced in comparison to the smaller size proppants. Consequently, the smaller size proppants are transported to a more considerable distance with the fluid flow. The reverse flow contributing to the negative velocity in the velocity profiles is due to the proppants after colliding with the fracture wall, moves in the reverse direction. Thus, the greater reverse flow velocity of smaller size proppants observed in Fig. 12 explains that a higher number of smaller size proppants colliding the fracture wall with higher velocity and rebounding back in the suspension layer. The parametric study of the proppant distribution to particle diameter suggests that it can play a significant role in optimising fracture conductivity. One effective approach, for low viscosity fluid-like slick water, could be injecting the smaller size proppant particles first in the slurry displaced by larger size proppant particles for distributing the higher volume of the fracture with proppants. To understand the effect of multisize proppant injection, smaller proppants with diameter 0.3 mm were injected for 0-1 s followed by 0.5 mm proppant diameter for 1-2 s and 1 mm proppant diameter for 2-3 s respectively . It can be observed that as the smaller proppant has a greater tendency for suspension and are injected prior to the larger proppant, thus multisize proppant injection results in improved distribution and can lead to more uniform fracture conductivity. Tsai et al. [Tsai, Fonseca, Lake et al. (2012)] also reported similar observations in their study with different proppant size.  , vol.16, no.2, pp.297-337, 2020 Time=2.5 s legend

Effect of fluid viscosity
In the next study, the fluid viscosity was varied, keeping all the other parameters constant, and simulation run was performed. The three cases of variation in fluid viscosity studied are µ=1 cP, 10 cP and 100 cP. Figs. 15 and 17 are the contour plots of proppant volume fraction and proppant horizontal velocity respectively at fracture mid-plane for different time step showing all the three cases of variation in fluid viscosity. It can be interpreted from the contour plots that particle deposition is much dependent on the fluid rheological property. With the increase in the viscosity of the fracturing fluid, the percentage of the proppant deposition decreases substantially. This is due to the increasing viscous forces that provide more flow resistance and drag force on the proppant particles. As a consequence, higher viscosity fracturing fluid prevents proppant particles from depositing. This observation is evident in Fig. 15 that for the simulation case of 100 cP fracturing fluid viscosity, the proppants deposition is considerably low. Next, the proppant volume fraction and proppant horizontal velocity were plotted with the fracture height and the advancement of proppant volume fraction with time at the twodifferent vertical planes was analysed (Figs. 16 and 18). The results from Fig. 16 show that the proppant volume fraction is identical for all the cases at the beginning, but later with time, low viscosity fracturing fluids result in greater proppant deposition compared with high viscosity fracturing fluids. The results from Fig. 18 show that the proppant velocity profile is significantly dependent on the fracturing fluid viscosity. With the progress of time, closer to the wellbore (0.25 m from inlet), the proppant flow in high viscous fracturing fluid lags behind low viscous fracturing fluid. This can be attributed to the higher viscous resistance force provided by the high viscosity fracturing fluid, which promotes the suspension ability of the proppants and retards the proppant deposition. This parametric study results in an important conclusion that the proppant transport, distribution and settling is substantially dependent on the fracturing fluid viscosity. Highly viscous fracturing fluids impede the proppant deposition.  , vol.16, no.2, pp.297-337, 2020 Figure 15 , vol.16, no.2, pp.297-337, 2020 @

Effect of fracture width
In the next study, the fracture width was varied, keeping all the other parameters constant, and simulation run was performed. The three cases of variation in fracture width studied are width=0.01 m, 0.05 m and 0.005 m. Figs. 19 and 21 are the contour plots of proppant volume fraction and proppant horizontal velocity respectively at fracture mid-plane for different time step and shows all the three cases of variation in fracture width. Two important observations can be made. Firstly, the difference in proppant distribution is less pronounced initially at time=0.5 s, and the effect of the fracture width is visible only at later times as the flow progresses. And secondly, the smaller width leads to greater wall resistance to the flow. As a result, the proppant in the lower fracture width case tends to deposit quickly, leading to a greater height of the dune formation. On the other hand, for the greater fracture width results in lower wall resistance, leading to the proppant particles travelling farthest and covers maximum volume. Next, the proppant volume fraction and proppant horizontal velocity were plotted with the fracture height and the time evolution of proppant volume fraction, and proppant horizontal velocity at the two-different vertical cross-sections was analysed respectively (Figs. 20 and 22). An important observation that can be noticed in Fig. 20 is that for the cross-section of 0.25 m from the inlet and lower fracture width case (0.005 m), higher proppant deposition and suspension characteristics is observed. This is because the lower fracture width due to greater wall resistance tends to form higher particle dune. On the contrary, for the crosssection of 0.8 m from the inlet and lower fracture width case (0.005 m), lower proppant deposition characteristics are noticed. This is due to unlike higher fracture width case, the proppant in lower fracture width do not tend to spread to the higher volume of the fracture, thus resulting in lower concentration away from the wellbore. The results from Fig. 22 show that the fracture width plays a significant role in the proppant velocity profile. The proppant velocity for the lower fracture width is considerably lower compared with greater fracture width. At t=1.5 s, away from the wellbore (0.8 m from inlet), the proppant horizontal velocity is almost three times higher for fracture width=0.05 m compared with fracture width of 0.01 m. This can be attributed to the lower fracture width results in greater wall resistance and consequently, higher proppant deposition. 326 FDMP, vol.16, no.2, pp.297-337, 2020

Comparison of eulerian-granular method with discrete element model
To understand the difference in proppant distribution between Eulerian-Granular method and CFD-DEM method, a separate study was carried out using an inlet velocity of 0.1 m/s and proppant volume fraction of 0.15. All the other simulation parameters were the same as described earlier in Tab. 1. Fig. 23 shows the comparison of Eulerian-Granular and DEM methods. The DEM Model shows the accurate proppant location, as it captures the complete particle micromechanics and tracks the individual particle and is computationally very expensive. On the other hand, Eulerian-Granular method provides proppant volume fraction, which can act as a substitute for the proppant position. One of the most significant advantages of using Eulerian-Granular method in proppant transport is that it is computationally economical compared with the DEM model. Fig. 23 compares both the approaches and shows that the particle distribution rate at the suspension layer and fracture bottom in Eulerian-Granular method is significantly different from the DEM method. This can be explained by the different ways in which particleparticle and particle-wall interaction is captured in both these approaches. DEM method models the particle motion explicitly with a detailed inter-proppant and proppant-wall interaction and tends to capture the physical phenomenon close to reality. On the other hand, the Eulerian-Granular method is based on KTGF and considers the granular particles as continuous media. Thus, it describes more fluid-like behaviour for the proppants and results in higher particle distribution rate at fracture bottom. Proppants modelled using DEM has a greater tendency to collide and suspend in the slurry, resulting in transporting proppant to a longer length, whereas, in the Eulerian-Granular method, the proppant tends to settle quickly and form relatively greater proppant bed.
Figs. 24 and 25 shows the quantitative comparison of Eulerian-Granular and CFD-DEM methods. The proppant volume fraction and axial proppant velocity were plotted with the fracture height, and the time evolution of proppant volume fraction and axial proppant velocity at the vertical cross-section of x=0.25 m and x=0.8 m from inlet was analysed. It can be observed from Fig. 24 that the proppant suspension layer for the DEM method is considerably greater at both the cross-section suggesting that the proppant transport using DEM methods tends to suspend greater proppants and transport proppants to a long distance away from the wellbore, compared with the Eulerian-Granular method. This can be explained by explicit treatment of frictional velocity, viscosity and inter-particle interaction in CFD-DEM method, which provides an accurate prediction of proppant distribution inside the fracture. Whereas, as explained earlier, the Eulerian-Granular method considers the granular particles as a continuous phase with a high viscosity. Thus, it describes more fluid-like behaviour for the proppants. On close observation from the Fig. 25, which is the plot of proppant velocity vs fracture height, it shows that at the fracture bottom, the greater proppant velocity observed in Eulerian-Granular method relates to the greater proppant distribution noticed in the contour plot (Fig. 23). Further, the proppant velocity for the DEM method is relatively lower compared with the Eulerian-Granular method. This can again be explained by the difference in proppant physics between the two methods. The Eulerian-Granular method uses KTGF and considers the granular particles as continuum media with a high viscosity. Thus, it describes more fluid-like behaviour for the proppants resulting in relatively greater proppant velocity at fracture bottom. In general, from the comparison of both the techniques, it can be interpreted that Eulerian-Granular method provides a reasonable approximation to the proppant particle physics inside the fracture. Considering the considerable simulation time required for the DEM method, and applicability for upscaling the model to field-scale hydraulic fractures, this comparison study suggests that Eulerian-Granular method can be used for practical problems of petroleum engineering interests for proppant distribution and settling.  , vol.16, no.2, pp.297-337, 2020 Figure 25: Comparison of proppant velocity for Eulerian-Granular vs CFD-DEM model

Conclusion
Proppant transport study in hydraulic fractures was conducted using the advanced numerical flow modelling methods, namely, Eulerian-Granular method and Discrete Element Method in commercial computational fluid dynamics (CFD) software, ANSYS FLUENT. A user-defined function was defined in order to mimic and model the fluid leakoff rate in the porous reservoir through the hydraulic fracture. It was established by adding the momentum and mass source terms in the flow governing equations.
The results were validated with the reported experimental study and show good agreement. The parametric study was performed to understand the proppant settling and transport mechanism by the variation in proppant properties (proppant diameter), fluid property (fluid viscosity) and geomechanical property (fracture width). The results show that proppant distribution is significantly affected by these properties. Small diameter proppant tends to remain suspended in the slurry, and larger diameter proppant tends to settle down quickly. Secondly, highly viscous fluids prevent the proppants from depositing due to the significant increase in the drag forces and proppants with lower fracture width tends to form deposition dune quickly. Finally, the comparison of the Eulerian-Granular method was made with the DEM method. The Eulerian-Granular method provides an approximate match with the DEM; however, the particle distribution rate in the Eulerian-Granular method is relatively higher than the DEM method. This was explained by the different ways in which particle-particle interaction is captured, and particle physics is handled in both these approaches. DEM provides a more accurate particle physics, but the computational time required is significantly higher. Considering the considerable simulation time required for the DEM method, and applicability for upscaling the model to field-scale hydraulic fractures, the current study suggests that Eulerian-Granular method can be used for practical problems of petroleum engineering interests for proppant distribution and settling. The current study has enhanced the understanding of complex proppant transport phenomenon in hydraulic fractures with fluid leak-off by capturing the proppant-fracturing fluid interaction and interparticle physics accurately using the advanced computational methods.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.