Coupled Hydraulic Fracture and Proppant Transport Simulation

Abstract

This paper focuses on the study of proppant transport mechanisms in fractures during frac-packing operation. A multi-module, numerical proppant, reservoir and geomechanics simulator has been developed, which improves the current numerical modeling techniques for proppant transport. The modules are linked together and tailored to capture the processes and mechanisms that are significant in frac-pack operations. The proposed approach takes advantage of a robust and sophisticated numerical smeared fracture simulator and incorporates an in-house proppant transport module to calculate propped fracture dimensions and concentration distribution. In the development of software capability, the propped fracture geometry and proppant concentration, which are the output of the proppant module, are imported to the hydraulic fracture simulator through mobility modification. Complex issues of proppant transport in fractures that are addressed in the literature and captured by the current model are: hindered settling velocity (terminal velocity of proppant in the injection fluid), the effect of fracture walls, proppant concentration and inertia on settling (due to extra drag forces applied on particles, compared to single-particle motion in Stokes regime in unbounded medium), possible propped fracture porosity and also mobility change due to the presence of proppant, and fracture closure or extension during proppant injection. A sensitivity analysis is conducted using realistic parameters to provide guidelines that allow more accurate predictions of the proppant concentration and fluid flow. The main objective of this study is to link a numerical hydraulic fracture model to a proppant transport model to study the fracturing response and proppant distribution and to investigate the effect of proppant injection on fracture propagation and fracture dimensions.

1. Introduction

In the numerical simulation of hydraulic fracturing in a reservoir, the actual size of a fracture can be modeled if it extends over several grid-blocks [1,2]. Since a fracture has high permeability with a very small width, compared to the model dimensions, having one grid refinement for the fracture and another grid refinement for the reservoir [3,4] requires a large degree of refinement for the entire model which poses serious numerical stability problems and requires long run-time of the computer. To overcome these limitations, maximize numerical stability and reduce run time, [4,5] proposed the use of one common grid system for both the fracture and the reservoir, where the fracture is modelled by increasing the permeability of the grids that contain the fracture. However, the proppant transport simulation requires explicit modelling of the fracture. Proppant transport simulation is mostly carried out by explicit modeling of the fracture. However, exceedingly small grids that are used to discretize the fracture cannot be used everywhere in the model. Therefore, the idea of smeared fracture modelling combined with proppant transport appears promising. In the numerical scheme proposed in this paper, the approach described in [4,5] is used to treat the fracture in the reservoir module, and the proppant module simulates an explicit fracture with realistic dimensions. In other words, the benefits of the previous two approaches have been combined in the current modular simulation. A third module is considered that calculates stresses and displacements during fracture formation and propagation to include the geomechanical effects.

Modelling proppant transport has been around since the 1970s. The modellings can be categorized based on the treatment of the fracture or based on the mathematical formulation that describes the transport phenomenon. Considering the transport formulations, the modeling approaches are either simplified proppant transport and settling, mixture model or kinetic theory for granular material. Since the approach in this paper falls within the mixture model, the literature review is forced on the work pertinent to the mixture model.

Typically, in mixture-type modelling, slurry flow inside the fracture is described by the lubrication theory. An essential element of this theory is the averaging of the variables such as concentration and velocities across the fracture’s width. Consequently, it implies that the concentration of proppant is constant across the fracture aperture in any cross-section, or its variation is insignificant, which can be neglected. This simplification resolves a significant modeling complication: since the fracture width is typically so much smaller than the other dimensions, discretizing an actual size fracture across its width for a field-scale model requires a large number of elements to keep the element aspect ratio reasonably small. Among other assumptions of this modeling, which are removed in later simulations, the no-slip condition between the proppants and the fracturing fluid is an important assumption. In addition, diffusion of proppant is assumed to be negligible, and thus the front of the proppant concentration profile remains sharp [6]. Moreover, fluid flow is assumed to be incompressible, and the proppant particles are small compared to the fracture aperture.

Reference [7,8,9] adopted Settari et al.’s approach and integrated commercial reservoir, fracture, and proppant simulators. Reference [6] assumed incompressible fracturing fluid and proppant particles, and simulated proppant transport in pseudo-3D hydraulic fractures. He further assumed that the only mechanism that accounts for the slippage between proppant and carrying fluid is settling.

Reference [10,11] published other important work in this area. These two publications serve as a validation to their commercial proppant simulator GOHFER. Based on the theory on GOHFER, the Stokes equation has been modified to account for the effect of proppant concentration with one equation. However, the wall effects and the inertial effects of the proppant have not been considered. The same approach has been used to account for the differences in horizontal velocities between the proppant and the flowing fluid. In References [10,11], the fracture simulator uses an integral equation for the width of the fracture and fluid pressure (the pressure is obtained by integrating the displacement for the point load over the surface of the fracture). Simulations have been carried out using GOHFER [12,13] under these limitations.

The simulations performed at the University of Texas at Austin by [14,15,16,17] are among other well-established works in proppant simulation, and they have provided successive improvements in their models. The first few models assumed analytical PKN fracture geometry until [14] incorporated a purely numerical hydraulic fracture simulator into their model. However, the model was based on elasticity neglecting plastic deformations in the formation.

Reference [18] combined the formulation of steady viscous flow with the flow of spherical proppants in Khristianovich–Zheltov–Geertsma–De Klerk (KGD) and pseudo-3D fractures. They were able to capture the formation of proppant bed and its growth. It is concluded that even without the leak-off, proppants may still advance to the fracture’s tip. However, their model is not capable of modeling fracture closure and expected asymmetry caused by gravitational settling.

Reference [19] integrated ITASCA’s PFC3D software, based on the Discrete Element Method (DEM), with a numerical fluid flow model to simulate proppant embedment in shale hydraulic fracturing. They discussed the fracture closure and shale hydration effects on proppant embedment. Reference [20] applied the DEM for hydraulic fracturing of two-layer formations in which the primary layer is under high stress while the upper layer is under lower stress. They proposed guidelines on the optimal proppant size and the fracture’s pinching as it propagates towards the upper layer.

Many of the published numerical models for simulating proppant transport use an analytical hydraulic fracture model. In most of these works, the fracture width is calculated from PKN, KGD, P3D or PL3D models [6,15,16,17,21,22]. In the most recent works, an adaptive re-meshing technique is used to couple a full 3D fracture model with a proppant transport model. However, the fracture model is fully elastic [14]. These models neglect plastic deformations of the medium and assume plasticity has minor effects on fracturing. Besides, conventional simulators do not include the deformation resulting from the interaction between stress and fluid flow response in a porous medium.

This paper presents the linkage between the three modules of fluid flow simulator, geomechanics module and proppant transport simulator. The main algorithms implemented in the numerical code, such as adaptive re-meshing of the propagating fracture, moving boundary conditions, time-stepping scheme and mobility averaging in the simulator, are discussed in this paper.

2. Hydraulic Fracture Module

The hydraulic fracture numerical simulator consists of two modules: fluid flow and geomechanics. The fluid flow simulator is the host or master module in the iteratively coupled model, which means it is run at the beginning of each time-step, and it triggers (calls) the geomechanics module to calculate stresses and displacements. It is well known that the orientation of a fracture is determined by the in-situ stress field: the hydraulic fracture will propagate perpendicular to the minimum principal stress. In all the simulations, it is assumed that the minimum principal stress is horizontal. Therefore, the fracture plane is vertical, normal to the direction of the minimum in-situ stress.

There are many different fracture initiation and propagation criteria in the literature. In this work, geomechanical analysis based on effective stresses in the rock is used to determine the fracture initiation and propagation. Assuming compressive stresses are negative, fracture initiation and propagation in this work are assumed to occur when the tensile stress exceeds the rock’s tensile strength:

$${\sigma }_{tip}+p\ge T$$

(1)

where σ_tip is the stress at the tip of the fracture in the minimum principal stress direction, p is pore pressure inside the fracture and T is the rock tensile strength. This equation means that fracture will initiate through a grid-block if the minimum principal effective stress at that grid-block exceeds the tensile strength of the rock. This criterion is checked at every time-step in the geomechanical module to calculate the size of the fracture.

The total length of the fracture is determined by the sum of the length of all the elements that satisfy Equation (1). The height of the fracture is assumed to be equal to the pay thickness. The width, on the other hand, can be calculated from the geomechanical analysis. If the initiation criterion is met in any grid block, then the fracture’s width at the corresponding location is determined from the nodal displacement of that grid block in the direction normal to the minimum principal stress:

$$w(x,y)=2u_y(x,y)$$

(2)

where w is fracture width and u_y is displacement perpendicular to the minimum principal stress.

In the smeared fracture simulator, only one common grid system is considered for both the reservoir and the fracture. If the fracture is modeled based on its actual dimensions, it will require many time-steps in the simulation and requires large computational resources. Alternatively, the permeability and porosity of the fracture are smeared within a grid-block. Based on the fluid flow cubic law, the average permeability of the fractured grid-blocks can be written as:

$$k_{avg}=\frac{k_m(\Delta y-w)+k_{f}w}{\Delta y}$$

(3)

where k_m is the matrix permeability, and Δy is the grid size. During fluid injection into the formation, the porosity of the rock changes. It is desirable to express this change in porosity in terms of the volumetric strains so that the direct outputs of the geomechanics module can be used to obtain the new porosity value in each time-step. The porosity change can be calculated from:

$${\varphi }_{new}=\frac{{\epsilon }_v+{\varphi }_0}{1+{\epsilon }_v}+\frac{L_{f_i}w}{\Delta x\Delta y}$$

(4)

where ε_ν is the volumetric strain, ϕ₀ is the initial porosity, and L_fi is the fracture length in grid-block i. The first term in this equation comes from the deformation of the rock, and the second term accounts for fracture porosity. By treating porosity in this way, fracture closing or opening, which in turn results in a fracture volume change, can be captured as a change in fracture porosity. Further details about the formulation of the hydraulic fracture module can be found in [23,24,25,26].

3. Proppant Transport Modules

Several algorithms have been implemented in the model to capture different phenomena occurring during the injection. To better understand the coupling of the modules, these algorithms are explained in this section before presenting the last step of our development: linkage of the modules of the numerical tool.

3.1. Mathematical Formulation

According to the law of conservation of mass, proppant mass entering the system minus proppant mass leaving the system plus proppant source or sink in an increment of time will result in a mass change in the system. Mass balance equation of the proppant can be expressed as:

$${\rho }_p\frac{\partial [cw]}{\partial t}+\frac{\partial [{\rho }_p{u}_{p}cw]}{\partial x}+\frac{\partial [{\rho }_p{v}_{p}cw]}{\partial y}-c{q}_{inj}=0$$

(5)

where ρ_p is the proppant density, w is the fracture width, u_p and ν_p are horizontal and vertical proppant velocities, respectively, c is proppant concentration and q_inj is the injection flow rate. This form of mass conservation is categorized as pure advection, assuming no diffusion. The 5th order WENO scheme and 4th order Runge-Kutta method are used to solve this equation. The ν_p term contains the settling velocity of the proppant. This parameter is calculated from Stokes law with corrections for the effect of inertia, fracture walls and concentration.

On the other hand, by adding the mass balance equation of the proppant and the injection fluid, the mass balance equation of slurry becomes:

$$\frac{\partial w}{\partial t}+\frac{\partial }{\partial x}(u_{sl}w)+\frac{\partial }{\partial y}(v_{sl}w)-q_{inj}=0$$

(6)

where $u_{sl}$ and $v_{sl}$ are slurry horizontal and vertical velocities defined as:

$$u_{sl}=c{u}_p+(1-c)u_f$$

(7)

$$v_{sl}=c{v}_p+(1-c)v_f$$

(8)

where $u_f$ and $v_f$ are fluid horizontal and vertical velocities, respectively. In hydraulic fracturing applications, the cubic law describes the motion of a fluid as a relationship between flow velocity (or momentum) and pressure [27]:

$$u_{sl}=-\frac{w^2}{12{\mu }_{sl}}\frac{\partial p}{\partial x}$$

(9)

similarly, for the y direction:

$$v_{sl}=-\frac{w^2}{12{\mu }_{sl}}\frac{\partial (p-{\rho }_{sl}gy)}{\partial y}$$

(10)

where $u_{sl}$ and ${\rho }_{sl}$ are the equivalent viscosity and density of the slurry, respectively. By applying the cubic law to the slurry mass balance equation, one obtains:

$$\frac{\partial }{\partial x}(\frac{w^3}{12{\mu }_{sl}}\frac{\partial p}{\partial x})+\frac{\partial }{\partial y}(\frac{w^3}{12{\mu }_{sl}}\frac{\partial (p-{\rho }_{sl}gy)}{\partial y})-q_{inj}=\frac{\partial w}{\partial t}$$

(11)

The conventional central finite difference scheme is used to solve this elliptic PDE and obtain the fluid pressures.

The slurry viscosity is affected by the proppant concentration. We assume that the increase in viscosity can be described by the following expression [10]:

$$\mu =\frac{{\mu }_0}{{(1-\frac{c}{c^{*}})}^{1.82}}$$

(12)

where c* is the saturation concentration, and μ₀ is the initial viscosity of the clean fluid. When the proppant concentration reaches the saturation concentration, a pack of proppant is formed inside the fracture, which behaves like a porous medium. Experimentally determined values of saturation concentration vary between 0.52 (loose packing) and 0.65 (dense packing) (for example, see [28,29]. A value of 0.6 for saturation concentration is assumed, which corresponds to a dense packing structure.

3.2. Proppant Entry Requirement

When the injection of the proppant starts, the dimensions of the fracture are known. As explained in Section 3.1, the hydraulic fracture simulator (linked reservoir and geomechanics modules) provides the initial fracture length and width based on the effective stresses in the rock and nodal displacements, respectively. However, the proppant cannot enter the fracture if the width is not 3 to 4 times greater than the largest proppant diameter. This phenomenon is well known in the literature and is called “fracture entry requirement” [27].

In all of the simulations, it was assumed that a proppant particle could enter the fracture as long as its diameter is equal to the fracture width. Hence, the length of the fracture that is transferred to the proppant simulator is shortened to the point that its width is equal to proppant diameters. It is this new shortened length that is discretized and used in the proppant simulator module. Figure 1 shows this process schematically.

3.3. Mesh Design and Adaptive Re-meshing of Different Modules

Fluid flow and geomechanics modules are incorporated in the hydraulic fracture simulator. Although the two modules can have different meshing strategies, the same mesh is used for both to prevent the need for mappings between the two modules. However, proppant simulator grids are designed to be much finer than the hydraulic fracture (or geomechanics) simulator grids to model the fracture explicitly. Therefore, the fracture grid that may contain proppant should be merged into the bigger reservoir grid. The term “sub-grids” is used to refer to the proppant simulator grids since they are much smaller than the reservoir and geomechanics grids. Figure 2 shows the schematic of the mesh design of the simulators. In all the simulations of this paper, based on the input parameters which determine the numerical stability of the model, each reservoir simulator grid is divided into 10 sub-grids for proppant transport simulation.

The width of the fracture, as mentioned before, is calculated from the geomechanics module. Although the proppant simulator assumes uniform properties along the fracture’s width and does not solve mass balance along the width, a gradual width reduction towards the tip of the fracture is considered in calculating fracture permeability, according to the cubic law. Fracture width in each sub-grid is calculated by linear interpolation between the width coming from the geomechanics module.

As the fracture propagates, its length and width will change while the height is assumed to be constant. At the start of each time-step, the mesh for the proppant simulator is adaptively changed to include newly fractured grids. However, the dimension of the grids is kept the same during the entire simulation time. This approach provides two advantages. First, the stable time-step for the proppant simulator (based on Courant-Friedrichs-Lewy condition) will change only due to changes in velocity and not changes in fracture length. Second, there is no need to interpolate the result of the previous time-step from the old mesh to the new time-step with the new mesh. Figure 3 shows how the adaptive re-meshing technique has been implemented in the model.

The pore pressures should also be linearly interpolated between larger reservoir grids and proppant sub-grids in addition to interpolating width between larger geomechanics grids and proppant sub-grids. Figure 4 is a schematic of the interpolation process.

It should be noted that although the length and height of the sub-grids are kept constant, the width might change during each time-step. Consequently, the volume of each sub-grid is not constant during the simulation. Obviously, proppant concentration is inversely proportional to the volume of the sub-grids:

$$c=\frac{V_p}{V_{subgrids}}$$

(13)

where V_p denotes proppant volume. As can be seen from the above equation, any changes in the sub-grids’ volume will change the concentration. Therefore, before each time-step, the concentration of the previous time-step should be adjusted to account for the change of the sub-grids volume. The old volume of each sub-grid is assumed to be:

$$V_{subgrid}^{Old}=dx.dy.w^{Old}$$

(14)

while the new volume of the sub-grid is:

$$V_{subgrid}^{New}=dx.dy.w^{New}$$

(15)

Hence, the adjusted proppant concentration would be:

$$c^{New}=\frac{V_p}{V_{subgrid}^{New}}=c^{old}\times \frac{V_{subgrid}^{Old}}{V_{subgrid}^{New}}=c^{old}\times \frac{w^{Old}}{w^{New}}$$

(16)

3.4. Averaging Mobility

As the concentration of proppant increases, the viscosity of the slurry is also increased. On the other hand, the fracture’s width and, consequently, the permeability, along its length, will be different. Therefore, the ratio of permeability to viscosity or mobility of the slurry will be a variable during the simulation.

The mobility that is needed in the reservoir simulator should be obtained using an averaging technique. Before running the simulators, each reservoir grid block is ensured to contain a certain number of proppant sub-grids (for example, in every and each reservoir grid-blocks, there exist 10 sub-grids). After running the proppant simulator, the proppant concentration and the dimensions of the sub-grid are known for each sub-grid. If sub-grids’ permeability and viscosity are assumed to be k₁, k₂, …, k_n, and µ₁, µ₂, …, µ_n, respectively, with n being the total number of sub-grids in a reservoir grid-block, considering Figure 5, the total pressure drop along all sub-grids will be equal to:

$$\Delta p_t={\displaystyle \sum _{i=1}^n\Delta p_i}$$

(17)

substituting Darcy’s equation for the pressure drop in Equation (17):

$$\frac{q_t\Delta X}{dy{w}_{avg}}{(\frac{\mu }{k})}_{fr}=\frac{q_{1}dx}{dy}{(\frac{\mu }{wk})}_1+\frac{q_{2}dx}{dy}{(\frac{\mu }{wk})}_2+\dots +\frac{q_{n}dx}{dy}{(\frac{\mu }{wk})}_n$$

(18)

where $\Delta X$ is the length of the coarse reservoir grids:

$$\Delta X={\displaystyle \sum _{i=1}^n\Delta x_i}$$

(19)

and $w_{avg}$ is the average width of the first and last sub-grids:

$$w_{avg}=\frac{w_i+w_{i+1}}{2}$$

(20)

However, the flow rate passing through each sub-grid would be equal:

$$q_t=q_1=\dots =q_n$$

(21)

therefore:

$$\frac{\Delta X}{w_{avg}}{(\frac{\mu }{k})}_{fr}=dx[{(\frac{\mu }{wk})}_1+{(\frac{\mu }{wk})}_2+\dots {(\frac{\mu }{wk})}_n]$$

(22)

or:

$${(\frac{\mu }{k})}_{fr}=\frac{w_{avg}dx}{\Delta x}[{(\frac{\mu }{wk})}_1+{(\frac{\mu }{wk})}_2+\dots {(\frac{\mu }{wk})}_n]$$

(23)

Since the permeability of each sub-grid is related to the width according to the cubic law, one finally obtains:

$${(\frac{k}{\mu })}_{fr}=\frac{\Delta X}{w_{avg}dx}[\frac{1}{{\displaystyle \sum _{i=1}^n{(\frac{\mu }{w^3})}_i}}]$$

(24)

Now that the average mobility of the sub-grids is found, a parallel averaging should be used between matrix and fracture mobilities. According to Figure 6, the total flow rate passing through a reservoir grid is the sum of flow rates through matrix and fracture:

$$q_t=q_m+q_f$$

(25)

Using Darcy’s law in Equation (25):

$${(\frac{k}{\mu })}_{avg}\frac{\Delta Y\Delta Z\Delta p_t}{\Delta X}={(\frac{k}{\mu })}_{fr}\frac{w_{avg}\Delta Z\Delta p_{fr}}{\Delta X}+{(\frac{k}{\mu })}_m\frac{(\Delta Y-w_{avg})\Delta Z\Delta p_m}{\Delta X}]$$

(26)

But the pressure drops would be the same across the fracture and matrix:

$$\Delta p_t=\Delta p_f=\Delta p_m$$

(27)

Finally, the average mobility of the grids containing proppant will be:

$${(\frac{k}{\mu })}_{avg}=\frac{1}{\Delta Y}[{(\frac{k}{\mu })}_{fr}{w}_{avg}+{(\frac{k}{\mu })}_m(\Delta Y-w_{avg})]$$

(28)

If the fracture partially penetrates a grid, the procedure of finding the total mobility is the same, with another series averaging between the portion of the grid that is fractured and the portion that is not fractured. Injection of proppant may cause a significant increase in the viscosity of the slurry, which in turn decreases the slurry mobility to very low values. Our formulation accounts for the change of width and proppant concentration at the sub-grid scale.

3.5. Fracture Closure

When the mobility of the injection fluid is decreased due to the presence of proppant particles, a large pressure drop might occur inside the fracture. Consequently, the length and width of the fracture will change. As shown in Figure 7a, two outcomes should be expected to occur:

(1)

If no proppant exists in some portions of the fracture, it may close completely, and the width is equal to zero. Therefore, the fracture length will be shortened.

(2)

On the other hand, portions of the fracture that contain proppant will not close completely with zero width since the proppant prevent the fracture from closing. However, the width of the fracture will be reduced to the “propped fracture width.” At this state, the fracture width cannot decrease any further.

The first case can be easily implemented in the numerical code since an adaptive re-meshing technique, along with a moving boundary, has been used in the code. However, the amount of fracture closure on proppant and resulting width reduction in the second scenario should be calculated.

The reduction of fracture width will modify proppant concentration, while the volume of the proppant inside the sub-grids will be constant. Before the reduction of fracture width occurs, the volume of proppant in the sub-grids will be:

$$V_{p1}=c_{1}dxdy{w}_1$$

(29)

It is assumed that after fracture closes on proppant, the concentration in the sub-grid will reach the saturation concentration. Hence, the proppant volume can be obtained from:

$$V_{p2}=c_{sat}dxdy{w}_{propped}$$

(30)

where $w_{propped}$ is fracture propped width. Since the volume of proppant will not change during the closure,

$$V_{p1}=V_{p2}$$

(31)

or:

$$w_{propped}=\frac{c_1{w}_1}{c_{sat}}$$

(32)

This calculation process is shown schematically in Figure 7b. The value of the propped width is calculated for all the sub-grids at the beginning of each time-step. It determines the maximum reduction of y-displacement in the geomechanics module. When the total y-displacements in two adjacent geomechanics nodes reach the critical y-displacement (half of the propped width), the nodes are fixed in the y-direction, assuring that no more y-displacement will occur. There is, however, one issue with this approach that needs to be resolved. Due to the injection of slurry, the fracture may reopen, the nodes previously fixed should be freed. This can be done by looking at the pore pressures inside the grids. For every grid-cell, the pore pressure and corresponding fracture width are saved in a table in the numerical code. The propped width will also correspond to certain pore pressure. This pore pressure is called “critical pore pressure.” As soon as the pore pressure of the grid reaches the critical pore pressure, two nodes of the grid will be freed.

There are other methods of simulating fracture closure in which the stiffness matrix is changed to ensure nodal y-displacement is a constant number ($\raisebox{1ex}{$w_{propped}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$). The approach is called the penalty approach [30,31]. The penalty approach could not be implemented in current work since the commercial software FLAC2D used in this study does not allow modification of the global stiffness matrix.

3.6. Saturation Concentration and Proppant Pack Permeability

There is a maximum concentration of proppant that can fill up a fracture. This maximum concentration is an empirical parameter and depends on the density of packing. In the current simulations, it was assumed that the maximum proppant packing concentration is 0.6, which corresponds to a dense pack.

Before the proppant concentration reaches its maximum value, the fracture’s permeability is calculated from the cubic law. However, this law does not apply to a stationary pack of particles. If the concentration in any element reaches the maximum value, it is assumed that it behaves like a porous medium with a permeability calculated from the Kozeny-Carman equation [32]:

$$k=\frac{{(1-c)}^3{d}_p^2}{180c^2}$$

(33)

The coefficient of 1/180 was suggested by [33] for uniform spherical particles. The permeability of the fractured elements does not suddenly jump to the value calculated by Equation (33), as this will cause severe numerical instability. Instead, the permeability is reached to this value over several iterations. The procedure is explained by [23] in detail.

4. Verification of Finite Conductivity Fracture

The study discussed in this paper forms the basis of a multi-module numerical hydraulic fracture and proppant transport simulator. The numerical hydraulic fracture simulator development, its verification and validation and mathematical challenges in solving first-order hyperbolic proppant transport partial differential equations have already been published.

Modeling a field frac-and-pack problem is very complicated due to the lack of data and variability of material properties; hence, it is not plausible to validate the model using actual field data. Therefore, the verification and validation of the model were carried out on each module of the simulator separately. [23] described the hydraulic fracture module of this work. The first validation of the hydraulic fracture module was later presented by [24] against hydraulic fracturing laboratory experiments on cohesionless sands. Additionally, the hydraulic fracture module was verified against [34] analytical solution for pressures around a fracture [26]. Reference [25] applied the simulator to actual field hydraulic fracturing in oil sand during cold water injection. They discussed the advantages of the model over other simulators in weak/unconsolidated formations.

The first step in developing the proppant simulator module is obtaining the mathematical solution for the transport equation. The solution schemes are described by [35]. The verification of the proppant transport module is performed by simulating the proppant injection into fixed slots, and the results are compared with commercial software COMSOL [36].

When proppant is injected inside the fracture, the mobility (defined as k/μ) of the slurry is reduced due to an increase in viscosity. The effect of the decrease in mobility on pressure variations along the fracture is similar to a reduction in permeability. In other words, permeability decrease, and viscosity increase both reduce fluid mobility and consequently, the conductivity of the fracture will decrease.

To verify the average mobility calculations explained in Section 3.4, the numerical model results are compared with analytical solutions by [37], and the numerical work of [38], for static fractures. The [37] model contains a vertical fracture with finite conductivity at the center of an isotropic, horizontal, infinite, circular drainage area. They assumed that the cylindrical reservoir slab was bounded by impermeable strata. As shown in Figure 8, a finite conductivity fracture intersects the wellbore, and a slightly compressible fluid is injected inside the wellbore. Reference [38] simulated an explicit fracture in a finite square reservoir with the same abovementioned assumptions. In both studies, the results are presented in terms of dimensionless pressure and dimensionless time defined in Equations (34) and (35):

$$t_{Df}=\frac{kt}{\mu c\varphi x_f^2}$$

(34)

$$p_{wD}=\frac{2\pi kh}{q\mu B}(p_{wf}-p_{init})$$

(35)

It should be noted that the results of these two cases from the literature deviate in the long term since [37] assumed an infinite acting reservoir. Figure 9 shows a comparison between the analytical work of [37] and [38] numerical work for a full penetrating fracture with 0.2π conductivity.

A 2D model with input parameters of

Table 1. Input Parameters for Finite Conductivity Fracture Simulation.

Parameter	Value
Porosity	0.25
Permeability	0.9869 × 10−13 m2 (0.1 Da)
Fluid Viscosity	1 cp (0.001 Pa.s)
Fluid Compressibility	1 × 10−10 Pa−1
Injection Flow Rate	8 × 10−5 m3/s
Reservoir Dimensions (Drainage Area)	300 m by 300 m (90,000 m2)
Initial Reservoir Pore Pressure	5 MPa
Fracture Conductivity	0.2π
Fracture Half Length	150 m
Reservoir Grid Size	Variable
Fracture Sub-grid Size	0.1 of Reservoir Grid Size
Fracture Width	0.001 m
Slurry Viscosity	107.5 cp (0.1075 Pa.s)

was built to verify the averaging scheme explained in Section 3.4.

It is assumed that the half-length of the fracture is 150 m, and its width is constant at 0.001 m, with a proppant concentration of 0.5. Assuming a maximum packing concentration of 0.66, this gives a slurry viscosity of 107.5 cp, according to Equation (12). Each reservoir grid is divided into 10 sub-grids in all the simulations. According to Equation (28), the average mobility of the slurry becomes 9.3 × 10⁻⁹ m²/Pa·s. This is equivalent to a fracture conductivity of 0.2π.

Figure 10 shows the comparison between the result of this study and [37] analytical solution and [38] numerical results. A good match is obtained except at the beginning of time. According to the mesh sensitivity shown in Figure 10, accurate modeling of early transient flow requires a very finer grid, comparable to wellbore radius, near the wellbore to capture wellbore boundaries.

5. Numerical Algorithm of Coupling the Three Modules

In general, the coupling between the simulators falls into three main categories. It can be classified as fully decoupled, fully coupled, or partially decoupled. In the first two coupling approaches, an explicit fracture is simulated in the model. In the third approach, however, some effects of the fracture are included.

In general, in a fully decoupled approach, the fracture equations are solved independently of the reservoir response. In other words, the fracture geometry is obtained by numerical or analytical methods, and its effect is represented in the reservoir model using several methods, ranging from an overall leak-off to permeability modification in the reservoir block containing the fracture to increasing the wellbore radius [39,40]. Although this approach is computationally efficient, it is nevertheless a simplistic approximation of the reservoir. The main disadvantage of this type of coupling is its limited range of applicability in field cases.

In fully coupled models, as the name implies, the fracture and reservoir equations are simultaneously solved at all times. The propagation of fracture is simulated in a fixed reservoir grid. Although the reservoir fluid flow is not simplified in this approach, the implementation of a numerical solution is cumbersome and lacks flexibility. Simultaneous formulation of flow and geomechanics results in large matrices. Therefore, the model will be computationally time-consuming. Moreover, in fully coupled models, an iterative procedure is needed in the same fashion as iterative coupling, and this approach does not reduce the complexity of the problem.

The third approach falls in between the other two approaches in which the fracture propagation and reservoir fluid flow are solved independently. However, the results from each module are transferred to the other simulator to improve the outcome. In this approach, which is called partially coupled or partially decoupled, new fracture grids can be easily generated. In principle, it can be conveniently attached to any type of reservoir or fracture simulator. This greatly increases the flexibility and range of application of the method. In modular coupling, information is transferred between different modules and iterations are used to obtain a converged solution. Such an approach can even use a highly advanced commercial reservoir and geomechanics simulators.

The intent here is not to discuss fully coupled models since the focus is using a multi-module simulator, where the modules are linked together through partially coupled principals. This approach seems to be effective when considering the choices of the available geomechanics codes outside petroleum engineering. An in-house fluid flow simulator has been developed using MATLAB and used as the coordinating and linking module. The flow simulator is linked to ITASCA’s commercial stress simulator FLAC to develop the numerical hydraulic fracture simulator. After the start of proppant injection, a third module, which is the proppant transport simulator, is added to the modules.

Models can be built based on solving the fluid flow and stress equations in different modules. The model here consists of three separate simulators: fluid flow, geomechanics and proppant. In modular simulators, different strategies are applied to link the modules. In “one-way coupling,” pore pressures, which ae the output of reservoir simulator, are transferred to the geomechanics module, but no information is transferred back to the reservoir simulator [41]. In other words, the geomechanics module does not attempt to improve the fluid flow solution. This old-style type of coupling can lead to large errors when porosity strongly depends on the flow. In the “loose coupling,” reservoir and geomechanics modules are run sequentially, and the solution from each module is transferred to the other one [41]. However, this transfer of information happens once in each time step, and no iteration is performed. Such coupling cannot represent complex constitutive plasticity models of the formation. Loose coupling is also called “explicit coupling” or “sequential coupling.” The type of coupling that is used here is “iterative coupling.” In this type of coupling, iteration is carried out in each time step, between the fluid flow and geomechanics modules until specific convergence criteria are met. In each iteration, the previous guess of the permeability and porosity is used to solve the flow equation. The corresponding change of pore pressure is used to calculate new deformations and stresses, which in turn provide a new update of permeability. Iterative coupling, when converged, gives an equivalent solution to a fully coupled model, while it is much more flexible and less computationally demanding.

In all the current simulations, three modules, namely, fluid flow, geomechanics and proppant simulators, have been used. The coupled model is a new stress-dependent hydraulic fracture simulator, in which, an in-house finite-difference fluid flow simulator is coupled to the finite difference code FLAC2D. Following the definition of [42], such a treatment is called a partially coupled simulation. The coupling between the modules relies on the stress state computations inside the reservoir, induced by injection and pore pressure increase. In the next step, the calculated strains, and displacements in the geomechanical module are used to modify the reservoir permeability for use in fluid flow simulation. The partial coupling has the advantage of being both flexible and computationally cost-effective compared to a fully coupled scheme. In this section, the coupling of the stress-dependent reservoir simulator, geomechanical simulator, and proppant transport simulator is described in detail.

Any existing sophisticated modeling tool for each component can be used in the linkage that has been implemented. The linkage has been termed a partially decoupled or partially coupled approach since there is a separate module for stress and fluid flow that solves the equations separately in each time increment.

The critical aspect in the partially decoupled modeling is the modelling of fracture in the reservoir simulator. The numerical tool consists of three computational modules. The geomechanics module gives the stresses and fracture geometry based on the smeared approach [43]. The reservoir module solves the fluid flow equations in the reservoir and includes fracture by assigning a high permeability to the grids that contain the fracture. Finally, the proppant module simulates the proppant transport until all the proppant becomes immobile when the concentration reaches the saturation concentration, and the fracture is filled up with proppant. Figure 11 shows a schematic of the coupling process.

The proppant simulator solves the hyperbolic partial differential equations described in Section 3.1. The grids for the proppant simulator are generated dynamically based on fracture growth in the geomechanics module. If a fixed number of grid-blocks is assigned to the entire fracture in all the time-steps, the sizes of the cells become bigger as the simulation progresses, and this reduces the accuracy significantly. Therefore, a constant size is assigned to the elements for all time-steps. In this way, as the computations proceed, the number of grid-blocks becomes larger, and no stability problems occur. Besides, since the position of the grid-blocks is the same for all time-steps, no interpolation is required. If the length of the fracture in the current time-step is L_fⁿ and the assumed length of the proppant cells is Δx_p, then the number of grid-blocks that are added to the simulation model will be:

$$N_{NewElements}=\frac{L_f^n-L_f^{n-1}}{dx}$$

(36)

Since the structured mesh is used to discretize the moving boundary problem, there is no need to check the smoothness of the mesh or regularity of the elements. The incremental addition of new elements effectively prevents the element size from becoming too large.

As shown in Figure 11, there are two iteration loops in the coupled algorithms. The first iteration loop is between fluid flow and geomechanics simulators and provides the fracture dimension. The second iteration loop is between the proppant transport simulator and the hydraulic fracture package (fluid flow and geomechanics) and calculates the proppant concentration distribution and mobility.

The proppant simulator’s primary result is the distribution of proppant concentration in the fracture, which can be translated into the independent grid of reservoir simulator by reducing the fractured grid-block mobility.

Because permeability of the grid-block, which contains the fracture, has a strong dependency on the displacements in the geomechanics module, the proppant calculations must be done in small time-steps after proppant injection starts. The stable time-step in the explicit reservoir simulator is larger than the same in the proppant simulator. However, both of these stable times are very small. Therefore, the hydraulic fracture module is run for several time-steps at a constant concentration, and the result is transferred to the proppant module. Next, the proppant transport simulator is run for several time-steps at constant fracture dimensions until the total time reaches that of a reservoir simulator. Figure 12 shows the numerical algorithm of the linked modules.

Reference [26] discussed different coupling algorithms applied to multi-modular simulations. The implemented system has several advantages over a fully coupled model:

(1)

Since fracture and reservoir use different grids, better resolution of the fracturing process can be obtained. Any mesh sensitivity analysis with the reservoir model can be performed without the need to regenerate the fracture description.

(2)

The proppant simulator can be coupled with any fracture simulator, either smeared or explicit, numerical, or analytical. This greatly increases the versatility of this approach and its range of applications.

(3)

Due to the smeared approach, the computational efficiency is greatly improved, and the memory requirements, especially for the fracturing simulation, is significantly reduced.

(4)

Simulation of the leak-off during the fracture growth occurs naturally in the smeared approach, simplifying the numerical modeling scheme in this respect.

(5)

The proppant simulator requires tiny time-steps for numerical stability considerations. Decoupling it from the reservoir and geomechanics allows small time-steps to be assigned to the proppant simulator, while bigger time-steps can be assigned to the reservoir simulator. The proppant module is run for several time-steps until it catches up to the other simulators in time. Due to this severe time-step limitation, decoupling is inevitable.

It should be noted that the hydraulic fracture and proppant transport simulator have been verified or validated in the authors’ previous publications [24,26,35,36,43].

6. Numerical Simulations

Several numerical examples are carried out in this section to investigate the effects of proppant injection on hydraulic fracture shape and propagation. Since many parameters can influence the results, it is more informative to define a base case with reference parameters, and sensitivity analyses are carried out by changing one parameter at a time. In this simulation, the fracture is allowed to propagate for 2000 s (33.3 min) of pad injection. At this time, a proppant laden slurry is introduced, and injection continues until the whole fracture is packed with proppant. The boundary conditions and input parameters of the reservoir and geomechanics simulators are shown in Figure 13a, and

Table 2. Input Parameters for Dynamic Fracture Model.

Parameter	Value
Model Dimensions	100 m by 50 m
Injecting Fluid Viscosity	0.001 Pa.s.
Formation Initial Porosity	0.2
Formation Initial Horizontal Permeability	0.1 Da
Formation Initial Vertical Permeability	0.1 Da
Initial Reservoir Pressure	1.6 MPa
Young’s Modulus	1.785 GPa
Poisson’s Ratio	0.3
Vertical Principal Stress	4 MPa
Maximum Horizontal Principal Stress	6.7 MPa
Minimum Horizontal Principal Stress	3 MPa
Cohesion	1.185 MPa
Friction Angle	20
Dilation Angle	22
Injection Time	2000 s

. The additional proppant related parameters and boundary conditions of proppant simulator are described here in Figure 13b, and Table 2.

For the case of dynamic fracture propagation, a 50 m by 50 m model is built. Table 2 presents the input parameters of the model. It is assumed that the formation rock follows the Mohr-Coulomb plasticity model, and water is used as the injection fluid.

It is assumed that two barriers exist at the top and bottom of the reservoir, which block fluid flow outside the reservoir. It should also be noted that the proppant transport simulator needs two sets of boundary conditions since slurry and proppant mass balances are solved separately in this module. The boundary condition specified at the right-hand side of the model also moves as the location of the fracture tip changes due to fracture extension or closure. Proppants and fluid properties in the simulation are described in

Table 3. Input Parameters for Proppant Module.

Parameter	Value
Model Dimensions	Varies According to Fracture Propagation
Injecting Fluid Initial Viscosity	0.001 Pa.s.
Proppant Diameter	0.0006 m
Proppant Density	2600 kg/m3
Saturation Concentration	0.6
Slurry Viscosity Exponent	1.82
Slurry Injection Start Time	2000 s

.

Figure 14 shows the fracture length, fracture width and proppant front advance during the injection time. The fracture length reduces as the proppant is introduced at 2000 s, which means fracture has partially closed due to the increase of slurry’s viscosity. The partial closure occurs since no proppant is present near the fracture’s tip to prevent its closure. However, as injection continues, the length of the fracture starts to increase again. After some time, this propagation stops again because the mobility of injecting fluid is decreased significantly. If the closure algorithm explained earlier is not implemented, the length of the fracture will become zero. However, the algorithm keeps track of the proppant front advancement and does not let full closure to occur. When the proppant front reaches the end of the fracture, no more fracture propagation is observed, and the fracture length is stabilized.

The situation for the width of the fracture is very different. Increasing viscosity of the injection fluid causes an increase in the pressure inside the fracture. Consequently, it is expected that the width at the injection point increases at a higher rate. This fact is clearly shown in Figure 15 after proppant injection starts. Even when the proppant front reaches the end of the fracture, the width, unlike length, increases since the concentration and consequently slurry viscosity are still increasing.

Figure 15 shows the changes in fracture width at different locations along the fracture with time. The first observation is the change in the rate of width variation after the start of proppant injection. In sub-grids closer to the injection point, the rate of width increase is enhanced while in sub-grids further away, this rate is reduced. As the sub-grids become closer to the fracture tip, the width might become zero, which means total fracture closure in this area. However, as shown in Figure 15, some width development occurs after closure, which indicates fracture reopening during proppant injection.

Figure 16 and Figure 17 show proppant concentration contours at different times during slurry injection. Again, the fracture’s length changes until the proppant front reaches the fracture’s end when further length variation stops. A bed of proppant with maximum packing concentration is formed at the bottom of the fracture. The size of this bed determines the efficiency of the frac-packing operation. If a large bed is formed, it can block the transport of proppant to the near tip area of the fracture and result in premature screen out.

At this point, it is necessary to check the global proppant mass balance and use it to verify the accuracy of the numerical tool. The mass of proppant that is being injected into the fracture can be easily obtained as:

$$m_{injected}=q_{inj}{c}_{inj}{\rho }_{p}t$$

(37)

where $q_{inj}$ is the slurry volumetric flow rate, $c_{inj}$ is proppant injection concentration, ${\rho }_p$ is proppant density and t is injection time. On the other hand, the amount of proppant inside the fracture can be obtained using:

$$m_{inside}={\displaystyle \sum _{i=1}^N{c}_i{V}_i{\rho }_p}$$

(38)

where c is the amount of proppant concentration in each sub-grid, and V is sub-grid volume. The difference between these two parameters serves as an error check on mass balance. This error always exists in numerical simulations. Figure 18 shows the mass of proppant injected and inside the fracture and their difference, which is called mass balance error expressed in percentage. The error is negligible compared to the amount of injected mass of proppant.

Sensitivity Analysis

A sensitivity analysis is conducted to investigate the effect of proppant density, proppant diameter and injection fluid viscosity. In all the simulations, the initial fracture geometry before proppant injection is assumed to be the same as in the base case model, i.e., the proppant is introduced after 2000 s of pad injection.

Table 4. Input Parameters for Sensitivity Analysis Simulations.

Simulation Number	Proppant Density (kg/m3)	Proppant Diameter (m)	Viscosity of Injection Fluid (Pa.s)
1-1	1100	0.0006	0.01
1-2	1800
1-3	2600
1-4	3900
2-1	2600	0.0001	0.01
2-2		0.0006
3-1	2600	0.0006	0.001
3-2			0.01
3-3			0.4
3-4			1

shows the parameters in each simulation.

Figure 19 shows the proppant concentration distribution at 2194.3 s of injection for different proppant densities ranging from 1100 kg/m³ to 3900 kg/m³ (the maximum density that can be found in ceramic proppants with ultra-high strength). The terminal settling velocity of particles in an infinite medium depends on the drag force that is exerted on the particle. Extensive experiments have been conducted to measure the amount of this drag force as a function of Reynolds number [44,45,46,47]. The experimental data are usually plotted on a log-log graph to present the so-called “standard drag curve” (Figure 19).

Obtaining the terminal velocity from Figure 19 requires trial and error since the particle velocity term is repeated in the drag coefficient and Reynolds number definitions. To remove the time-consuming trial and error procedure in the simulator, we exploited the definition of another dimensionless number which is called Galileo Number [48] defined as:

$$Ga=\frac{3}{4}{C}_D{\mathrm{Re}}^2=\frac{d_p^{3}g{\rho }_f({\rho }_p-{\rho }_f)}{{\mu }^2}$$

(39)

where C_D is the drag coefficient, Re is Reynolds number, d_p is proppant diameter, ρ_f and ρ_p are fluid and proppant densities respectively and μ is viscosity. As Equation (39) shows, Galileo number is independent of terminal velocity and removes the necessity of performing trial and error. Galileo number can be obtained from the physical properties of carrying fluid and proppants. The experimental data of drag coefficient in Figure 19 are converted to Galileo number and displayed in Figure 20. This figure is explicitly incorporated into the simulator to account for the effect of inertia on settling velocity.

In the next step, the settling velocity obtained from Figure 20 is adjusted for the effect of concentration and fracture walls by:

$$v_{Adj}=v_{set}f(c)f(w)$$

(40)

where $v_{set}$ is obtained from Figure 20 based on Reynolds number, $v_{Adj}$ is settling velocity adjusted for the effect of concentration and fracture walls, c is the proppant concentration, w is fracture wall, $f(c)$ is the voidage function and $f(w)$ is the wall factor defined as:

$$f(w)=\frac{v_w}{v_{set}}$$

(41)

$$f(c)=\frac{F_{Dt}}{F_{Ds}}$$

(42)

where $v_w$ is settling considering the effect of fracture walls and $F_{Dt}$ and $F_{Ds}$ are terminal drag forces in multi-particle and single-particle systems, respectively.

In this study, the following relationship is used to calculate $f(w)$ [49].

$$f(w)=R_0+{\displaystyle \sum _{i=1}^{19}{R}_i{\lambda }^i}$$

(43)

where $R_0$ and $R_i$ are coefficients provided in

Table 5. Coefficients for Equation (43).

Coefficient	Triangular Cylinder	Square Cylinder	Parallel Plates
R0	1.0000000	1.0000000	1.0000000
R1	−0.1524694 × 101	−0.1923777 × 101	−0.4027060 × 100
R2	−0.9356945 × 101	0.1649393 × 101	−0.8435362 × 101
R3	0.6788950 × 102	−0.1153624 × 102	0.3487996 × 102
R4	−0.1634936 × 103	0.2682020 × 102	−0.2359584 × 102
R5	0.6563649 × 102	0.1367386 × 102	−0.1193919 × 103
R6	0.1929998 × 103	−0.5060226 × 102	0.1362242 × 103
R7	0.4729873 × 102	−0.1042480 × 103	0.1601959 × 103
R8	−0.3751033 × 103	0.1170802 × 103	−0.4106427 × 101
R9	0.2752887 × 103	0.2395431 × 103	−0.3171554 × 103
R10	−0.1190656 × 104	−0.1757552 × 103	−0.1989548 × 103
R11	0.6542166 × 103	0.1097079 × 102	0.4608181 × 102
R12	0.2075038 × 104	−0.2409061 × 103	0.3581750 × 103
R13	0.3268518 × 103	−0.9802373 × 102	0.2128604 × 103
R14	−0.2014247 × 104	0.1344775 × 103	0.2338137 × 103
R15	−0.2869939 × 104	0.2037087 × 103	−0.8912624 × 102
R16	0.1553641 × 104	0.9298401 × 102	−0.6472198 × 102
R17	0.1483049 × 104	−0.1153909 × 103	−0.2528621 × 103
R18	0.2078562 × 104	0.2791606 × 103	0.2000882 × 103
R19	−0.2211987 × 104	−0.3434466 × 103	0.2897953 × 103

and $\lambda$ is the ratio of particle diameter to fracture’s width.

Additionally, $f(c)$ or the voidage function is obtained from [50] expression:

$$f(c)={(1-c)}^{-\beta }$$

(44)

where exponent βchanges with the Reynolds number:

$$\beta =\{3.7-0.65e^{\raisebox{1ex}{${(1.5-{\log }_{10}\mathrm{Re})}^2$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\}$$

(45)

Any other voidage function or wall factor can be implemented in the simulator. The details of various proposed $f(c)$ and $f(w)$ in the literature is explained in detail by [51].

On the other hand, the convection velocity can be expressed as:

$$v_{convection}=\frac{w^2}{12\mu }\frac{\partial ({\rho }_{sl}gy)}{\partial y}$$

(46)

where ${\rho }_{sl}$ is the slurry density defined as:

$${\rho }_{sl}=\frac{100}{[\frac{c_m}{{\rho }_p}+\frac{1-c_m}{{\rho }_f}]}$$

(47)

and $c_m$ is the proppant mass concentration.

According to Equation (46), increasing proppant density increases the settling velocity directly. In addition, proppant with higher densities create a higher density slurry, and convection velocity may increase. As a result, a bigger layer of proppant bed is expected to form at the bottom of the fracture. The effect of density is more pronounced here due to lower viscosity that is assigned to carrying fluid.

In Figure 21, the blue area displays regions with zero proppant concentration or no sand; the red area displays the regions where proppant concentration has reached its maximum (saturation concentration) set in the model. A value of 0.6 for saturation concentration is assumed, which corresponds to a dense packing structure. Technically the red region displays the sandbank or proppant bed.

The suspending proppant is shown by the greenish/yellowish region in Figure 21 in which sand concentration is below saturation concentration. The sand flow can be seen in Figure 21 in front of the sandbank in which the bank extends. Sand particles roll or saltate and create a second bed with lower height attached to the main proppant bed.

Although very different proppant distribution is obtained, the fracture geometry in all simulations remains the same as shown in Figure 14 and Figure 15. The reason is that the concentration in most part of the fracture is nearly constant at 0.3, except for the proppant bed that develops at the bottom of the fracture. This gives a very similar value of viscosity and, consequently, mobility, which is the main fracture geometry controlling parameter. The same fracture propagation behavior was also observed in simulations with varying proppant diameters (Figure 22).

According to Equation (46), the proppant diameter does not affect the convection velocity. This is particularly true when the fracture aperture is several orders of magnitude larger than the proppant diameter. However, settling velocity is proportional to the proppant diameter. Therefore, increasing proppant diameter will create a bigger bed of proppant compared to the same order increase in density. Although having larger proppant will result in a proppant packed fracture with higher permeability (Equation (33)), large proppant may not go into a significant portion of fracture due to minimum width for entry requirement.

Based on the simulations’ results, the carrying fluid’s viscosity has the most effect on the fracture geometry and proppant distribution. Figure 23a,b show how the length and width of the fracture changes when the carrying fluid’s viscosity is equal to 400 cp and 1000 cp. In both cases, shortly after proppant injection starts, fracture shortens and widens due to an increase in mobility. The final propped fracture length is larger, and its width is narrower when a lower-viscosity fluid is used.

Figure 24 shows the proppant distribution of two viscosity cases at 2312.2 and 2411.2 s. Viscosity significantly changes both settling and convection. A better proppant placement is obtained by having a higher viscosity fluid. However, at the same time, a shorter fracture is to be expected.

In another simulation, the slurry viscosity is decreased to 1cp. A perfect example of tip screen-out is obtained in this simulation. It is apparent from Figure 25 that low viscosity results in longer and narrower fractures.

It is interesting to note that the proppant front never reaches the fracture tip. The reason is more evident when examining the proppant concentration distribution in Figure 26. It can be concluded from Figure 26 that a large bed of settled proppant is created at the bottom that covers a huge portion of the fracture. This is expected to occur as the carrying capacity of the fluid is reduced when viscosity is reduced. However, the size of the bed is so large that it blocks nearly the fracture’s full height before the proppant reaches the tip. In Figure 25, the proppant front is depicted in the middle section of the model, and it can be seen in Figure 26 that it never reaches the tip of the fracture during the simulated injection time.

7. Conclusions

In this paper, a hydraulic fracture and proppant transport simulator that includes the reservoir and geomechanical effects have been formulated. The main objective is to link the numerical hydraulic fracture model to a proppant transport model to study the fracturing response and proppant distribution to investigate the effect of proppant injection on fracture propagation and dimensions.

The model considers variable fracture permeability, poroelastic effects, fracture closure on proppant, and adaptive mesh design to reflect the physics of proppant transport more realistically compared to the existing numerical simulators. The fluid flow, mechanical and proppant equations are iteratively coupled in a multi-module numerical tool with adaptive time-stepping and re-meshing. The significance of the current coupling development is its flexibility in integrating any commercial fracture, fluid flow or geomechanics simulators efficiently. A framework is developed that allows the simulator to seamlessly capture proppant settling, fracture closure on proppant, fracture width variation, proppant accumulation into a packed bed between fracture walls and tip screen-out.

A series of sensitivity analyses is carried out to confirm that the simulation results are plausible and investigate the effect of different controlling parameters on proppant distribution inside hydraulic fractures. The sensitivity analysis shows the impact of proppant density, proppant diameter and fluid viscosity on hydraulic fracture behavior. Viscosity is by far the most important parameter on proppant transport, which can significantly affect fracture geometry and proppant placement. This sensitivity analysis is useful in testing the robustness of the developed numerical tool, and it provides a better insight into the parameters affecting the treatment operation. Moreover, the model’s errors are identified through sensitivity analysis when an unexpected relationship between inputs and outputs is found.

The application of the newly developed model provides valuable information for field frac-packing operation and optimization. The effect of different design parameters in proppant transport operation could be assessed before the actual operation. Currently, there is significant uncertainty in the effect of physical proppant and fluid parameters on the final proppant distribution. Therefore, the proposed numerical tool will increase the understanding of the relationship between fluid and proppant properties and the final distribution of these particles, which in turn determines the conductivity of the propped fracture, leading to a reduction of uncertainty and more realistic production forecast especially for reservoirs under improved or enhanced oil recovery scheme as found in heavy oil and oil sands projects.