Model development of proppant transport through hydraulic fracture network and parametric study

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

Contents lists available at ScienceDirect

Journal of Petroleum Science and Engineering journal homepage: www.elsevier.com/locate/petrol

Model development of proppant transport through hydraulic fracture network and parametric study ⁎

Oliver Changa, , Robert Dilmoreb, John Yilin Wanga a The Pennsylvania State University, Petroleum and Natural Gas Engineering, Department of Energy and Mineral Engineering and EMS Energy Institute, 202 Hosler Building, University Park, PA 16802, USA b US Department of Energy, National Energy Technology Laboratory, 626 Cochrans Mill Road, P.O. Box 10940, Pittsburgh, PA 15236-0940, USA

A R T I C L E I N F O

A BS T RAC T

Keywords: Proppant transport Hydraulic fracture network Shale gas reservoirs

A model capable of simulating proppant transport through hydraulic fracture network is developed and summarized in this paper. The proppant transport model (PTM) is able to capture multiple proppant transport patterns, including suspension, saltation and creeping. These patterns are first identified, and then quantified to establish proppant transport equations. The governing equations are programed into a three-dimensional, finite-difference model to simulate the proppant transport process. The PTM is coupled to a previously developed hydraulic fracture network propagation model, which updates essential input parameters such as fracture geometry, velocity distribution and pressure profile for each step. In every step, the proppant transport model extracts values for these parameters and solves the mass transport equations for all three patterns. Finally, the PTM generates proppant concentration, fracture conductivity and distribution throughout the created fracture network and predicts, at the end of the designed treatment the propped stimulated reservoir volume (PSRV) – a critical indicator of long-term stimulation effectiveness for hydraulically fractured oil/gas reservoirs. Parametric Studies of several important treatment, operational, reservoir, and geomechanical parameters are done in this paper to illustrate the impact of each factor on the PSRV.

1. Introduction Hydraulic fracturing is one of the most effective and efficient methods to enhance oil and gas recovery from unconventional reservoirs. Pre-existing natural fractures are commonly present in the unconventional reservoirs. Due to this distinct character, hydraulic fracturing treatment generates a hydraulic fracture network instead of a bi-wing planar fracture in the formation. Proppant transport behavior in the hydraulic fracture network system is more complicated and is not able to be simulated using conventional methods. The complexity of the fracture network causes the evaluation of the effectiveness of the treatment very difficult. Traditionally, engineers use the concept of stimulated reservoir volume (SRV) to evaluate the effectiveness of the hydraulic fracturing treatment. The quantification of the SRV is often estimated by microseismic mapping. However, the micro-seismic is not capable of indicting the actual propped-open portion of the fracture network, which actually contributes to oil and gas recovery. Ahn et al. (2014) firstly stated the concept of the propped-open stimulated reservoir volume (PSRV) that represents the actual propped-open portion of the created fracture network. The PSRV is the actual effectively stimulated

⁎

volume, which is contributing to oil and gas recovery. Weng et al. (2011) developed a model able to simulate the interaction between the hydraulic fracture and pre-existing fracture network. The proppant transport function in this model is not described and simplified. Mack et al. (2014) conducted proppant transport experiments in a set of intersecting transparent narrow slots to mimic the scenario of hydraulic fracture network. The experiment described some observation and phenomenon of proppant transport behavior in the narrow slot network. Recently, more researches investigating particle transport in hydraulic fractures are conducted. Blyton et al. (2015) proposed a CFDDEM type method to model particle transport in hydraulic fractures. The integrated model is able to capture detailed effects of particle – particle, particle – fracture interactions. The research suggests that the propped open length by proppants could be shorter than expected due to low proppant velocity because of concentration effect. Shiozawa et al. (2016) compared simulation results from Pseudo-3D and fully 3D models. The result indicates that the P3D results are consistent with the fully 3D results and with better efficiency. The research also investigates the effect of tip screen out and fracture closure. The results suggest that tip screen out is not significant in shale reservoirs, but the settling of proppants are phenomenal and is detrimental to proppant

Corresponding author.

http://dx.doi.org/10.1016/j.petrol.2016.12.003 Received 25 March 2016; Received in revised form 17 October 2016; Accepted 3 December 2016 0920-4105/ © 2016 Elsevier B.V. All rights reserved.

Please cite this article as: Chang, O., Journal of Petroleum Science and Engineering (2016), http://dx.doi.org/10.1016/j.petrol.2016.12.003

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

O. Chang et al.

Fig. 1. Hydraulic fracture network geometry at the end of 60-min treatment (Ahn et al., 2014).

Propagation Model, and (2) Proppant Transport Model. The Hydraulic Fracture Propagation Model is coupled to the Proppant Transport Model to capture the proppant placement during the injection period. The Hydraulic Fracture Propagation Model developed by Ahn et al. (2014) is a finite-difference based model with capability to simulate fracture propagation in a naturally-fractured shale gas reservoir. This model captures the interaction between hydraulic and natural fractures, fluid flow in matrix and geo-mechanics of the rock. Complex fracture propagation can be generated as a function of different input parameters. The HFPM assumes two sets of uniformly spaced, orthogonal pre-existing natural fractures in the reservoir. The procedure starts with the injection of fluid into the perforation with corresponding pressure increase. Using an iterative method, fluid flow equation and geo-mechanical equation are solved at the same time. Once the pressure exceeds the initiation criterion of fracture propagation, the grid block then is converted into a fracture block. If the treatment continues, fracture propagates through the pre-existing fracture; if the fracture encounters an intersection, crossing and lateral initiation criteria are calculated to govern the propagation. Depending on different input parameters, the fracture either propagates across the intersection without opening secondary fractures, initiates opening of the secondary fracture, or both. Thus, for a designed treatment schedule and other scenario parameter input assumptions, the HFPM generates a hydraulic fracture network in the reservoir. Fig. 1 illustrates results from a case in which the HFPM treated at 17 bbl/min rate, with fluid viscosity of 2 cp, a total volume of 42,840 gallon is injected for half of the simulated cluster. The treatment information is provided in the appendix section. The proppant transport models developed in this article captures fluid flow velocity, fracture geometry and other essential input parameters throughout the fracture network domain that are output from the HFPM as in each time step, and uses that as input in simulation of particle transport dynamics.

placement. Han et al. (2016) investigated proppant transport in complex fracture geometry also using CFD method. The authors established several simplified scenarios for intersecting fractures such as T-junction and crossing geometries to mimic a part of the natural fracture network. The results indicate that a great portion of proppants remain in the major fracture instead of entering the natural fractures. The research also indicates that lighter proppants, higher fluid viscosity, and higher pumping rate increase the amount of proppants entering into the natural fractures and vertical propped coverage. Most of the works focused on proppant transport in a bi-wing planar fracture, which is not suitable for naturally-fractured unconventional reservoirs. The previous experiments or numerical studies are done under fixed fracture geometry, which neglects the effect of changing fracture geometry and fluid flow velocity. Some of the previous researches use CFD-DEM method to model particle movements in the fractures. The computational burden makes it really difficult to simulate proppant transport in a larger scale and efficient application. The physics governing proppant transport in fracture network is not identified and quantified either. Also, there is no published model with multiple governing physics that simulates proppant transport in a hydraulic fracture network dynamically. From our previous research, Chang et al. (2016) summarized and also developed the governing physics of proppant transport through hydraulic fracture network. The author provided quantification of each governing physics such as advection, turbulent dispersion, settling, bedload transport, and the uneven fraction proppant mass flux at fracture intersection. All the major governing physics were identified and quantified. In this paper, we developed a three-dimensional, finitedifference, and multi-mechanistic proppant transport model to quantify the proppant distribution and concentration in the hydraulic fracture network in naturally fractured unconventional reservoirs. This is the first model that includes comprehensive governing physics of proppant transport in naturally fractured system. With the model, one is able to quantify the actual propped-open stimulated reservoir volume (PSRV) and evaluate the effectiveness of the treatment.

2.2. Mechanisms of proppant transport and the proppant transport model (PTM) 2. Model development Proppant transport in slick water is a complex process that is controlled by multiple physical phenomena and proceeds occurs under different flow regimes simultaneously, as described in an article by Chang et al. (2016, In preparation). Four regions can be identified that are associated with proppant transport in the fracture. The first zone near the tip of the propagating fracture is the pad solution that

2.1. Hydraulic Fracture Propagation Model (HFPM) In order to model proppant distribution in complex hydraulic fracture network, multiple physics needed to be captured. Two models are coupled to capture the physics – (1) Hydraulic Fracture 2

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

O. Chang et al.

contained no proppant. The second zone is the proppant slurry suspension. The third zone is the moving bed layer underlaying the suspension in which creep, saltation and repation occur. The last zone is the non-moving bed where the impact of flow is negligible – located at the bottom of the fracture. The mass transport equations governing the particle behavior in each of these zones are developed in a differential form as discussed in the following sections. The Proppant Transport Model will then iteratively solve the proppant concentration within a given time step in the fracture network domain. Once done with a time step the model saves the proppant distribution and the Hydraulic Fracture Propagation Model will propagate the fracture network again for another given time step. The process will continue until the scheduled pressure pumping treatment duration is achieved.

vs′=vs′fw ∅s 5.5=

where qCs is the negative source term in terms of proppant concentration change per volume, qCp is the proppant concentration rate, vs is settling velocity,vs′ is the correlated settling velocity, ∅s is the ratio of liquid volume to total slurry volume, fw is wall factor, accounting for the impedance of particles settling as a result of fractional force from the fracture wall, ρp is density of proppant, ρl is density of fluid, and Vb is the fracture block volume.

2.3.2. Governing equation of saltation and creep The concentration of proppant in the settled bed also changes with time and location in the fracture network. The bed receives settling proppant from the suspension and also from bedload transport such as creeping and saltation. The change in proppant concentration of the settled bed through time can be determined in a differential form.

2.3.1. Governing equation of suspension flow Turbulent eddies and advective fluid flow both contribute to proppant mass transport during hydraulic fracturing treatment. The governing advection-turbulent dispersion equation of slurry can be expressed as:

∂ 2C

∂Cpbed ∂t

∂ 2C

DC ∂C ∂C ∂C +u +v +l =Kx 2 +Ky 2 +Kz 2 Dt ∂x ∂y ∂z ∂x ∂y ∂z

(1)

q A

qCpsalt =mbsalt / Vb=

(2)

2(1−ν 2 ) h f E

( pf −σn )

where

(10)

Vb

τ* = τ* −1, crit

(s − 1) gdp 3 wf ρp

(11)

* is dimensionless τ * is dimensionless shear stress, τcrit

critical shear stress, s =

ρp ρf

, and σ=2. 45

* τcrit (ρs / ρf )0.4

Proppant concentration creep flux can be determined according to Wang and Zheng (2004).

(3)

qCpcreep =mbcreep / Vb

where σn is normal stress acting perpendicular to fracture element, and w is the fracture width, E is Young's moduli, ν is Poisson's ratio, and pf is the fluid pressure in the fracture. The fracture propagation model is capable of capturing horizontal fluid flow velocity; vertical fluid flow is neglected. Hence, the equation is modified as:

DC ∂C ∂C ∂ 2C ∂ 2C ∂ 2C +u +v =Kx 2 +Ky 2 +Kz 2 Dt ∂x ∂y ∂x ∂y ∂z

(9)

qb wf

⎞ ⎛ 1 qb = 0.635r τ * ⎜1 − ln(1 + σr ) ⎟ ⎠ ⎝ σr

where q is the local flow rate in the fracture, A is the fracture cross sectional area, vf is fluid flow velocity. The elasticity equation describes the pressure/width relationship. It is also assumed that plane-strain condition in the z-direction and consequently, the displacement in z-direction are zero.

w (x, t )=

=qCpsaltin−qCpsaltout +qCpcreepin−qCpcreepout +qCs

where Cp bed is proppant concentration of the settled bed, qCp salt is proppant concentration flux per unit width donated from saltation. qCp creep is proppant concentration flux per unit width resulted from creeping. Proppant concentration flux due to saltation can be expressed as the following term according to Yalin (1963).

where C is proppant concentration, K is turbulence diffusion coefficient factor, u is fluid flow in x direction, v is fluid flow in y direction and l is fluid flow in z direction. The fluid flow velocity is calculated from the pressure and fracture geometric data generated in each time step from the Hydraulic Fracturing Propagation Model. The fluid flow equation within the fracture is expressed as:

vf =

(8)

18μ

qCp=proppant concentration rate

2.3. Governing equations of PTM

∂ 2C

fw gdp2 ∅s 5.5 (ρp − ρl )

0

mbcreep=

∫−δ

c

⎡ ρf 2 ⎛V ⎞ ⎢ ρp vcreep ⎜ s ⎟ dy=⎢ Vs 2 ⎝ Vt ⎠ ⎣⎢ 2ρp Vt gεm

⎤ εm − εd ⎥ 4 v* cdg ⎥⎥ ⎦

(12)

where εm is static friction coefficient, εd is dynamic friction coefficient, c is experimentally determined constant ≅ 1.5, δc is the thickness of creep layer, v* is shear velocity, and Vs is the volume fraction of the layer.

(4)

Vt

For a given time step ∆t , proppant settles vertically toward the bottom of the fracture. The settled amount of proppant is treated as a negative source term that draws particles out of the cell. The equation is modified as:

DC ∂C ∂C ∂ 2C ∂ 2C ∂ 2C +u +v +qCs−qCp=Kx 2 +Ky 2 +Kz 2 Dt ∂x ∂y ∂x ∂y ∂z

qCs=

vs=

vs′∆t C ∆zVb

2.3.3. Finite-difference approximation method The purpose of the work is to predict proppant distribution throughout the hydraulic fracture network over time. A finite–difference method is used for proppant transport modeling. A unit cell of a fracture block is illustrated in Fig. 2. A mass transfer equation can be established according to the state of flow, location of the block, and other input parameters representing conditions in each unit cell. The governing equation of the system can be written in terms of proppant concentration. For this model, finite-difference scheme is used for the approximation. For suspension transport, the numericalized form is expressed as:

(5)

(6)

gdp2 (ρp − ρl ) 18μl

(7) 3

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

O. Chang et al.

of major dependent variables at the interface of the grid blocks using harmonic average as interpolation function. The transmissibility term for turbulent dispersion is expressed as: In x-direction:

⎡ KA Txturb i + 1 , j, k =⎢⎢ x x 2 ∆x ⎣

⎤ ⎥ ⎥ i + 1 , j, k ⎦ 2

⎤ ⎡ 2Kx i, j, k Kx i +1, j, k Ax i, j, k Ax i +1, j, k ⎥ =⎢ ⎣ Kx i, j, k Ax i, j, k ∆x i +1, j, k + Kx i +1, j, k Ax i +1, j, k ∆x i, j, k ⎦

(17)

In y-direction:

⎡ ⎢ Ky Ay Tyturb i, j + 1 , k =⎢ 2 ∆y ⎣

⎡ ⎤ 2Ky i, j, k Ky i, j +1, k Ay i, j, k Ay i, j +1, k ⎥ =⎢ ⎢⎣ Ky i, j, k Ay i, j, k ∆y i, j +1, k + Ky i, j +1, k Ay i, j +1, k ∆y i, j, k ⎥⎦

Fig. 2. A discretized fracture block with particle mass flux.

Cin, j+1 ,k

−

Cin, j, k

∆t

Vi, j, k +vfx

Cin+1, j, k

−

Cin−1, j, k

2∆x

Vi +vfy

Vi +qCsin, j+1 V −qCp n +1 Vi, j, k =Kx Ax , k i, j , k i, j , k

+Kz Az

∆t

i, j , k

Cin, j+1 +1, k

2Cin, j+1 ,k

−

+

Kx Ax

Cin, j+1 −1, k

⎡ KA Tzturb i, j, k + 1 =⎢⎢ z z 2 ∆z ⎣

∆x n +1 − 2Cin, j+1 Cin, j+1 +1, k , k + Ci, j −1, k n +1 n +1 Cin, j+1 , k +1 − 2Ci, j, k + Ci, j, k −1

∆y

+1 n +1 n +1 Cin+1, j, k − 2Ci, j, k + Ci −1, j, k

+Kz Az

∆x n +1 − 2Cin, j+1 , k + Ci, j, k −1 (14)

∆z

Since the advection term is directional, an upstream weighting scheme is applied to obtain accurate approximation of the entering and exiting proppant concentration flux. For the following expressions, we assume the fracture fluid is flowing toward the X (j+1) and Y (i+1). If the fluid flow direction is opposite, the code recognizes the fluid direction and adjusts the setting of upstream weighting.

⎤ ⎡ 2vfx i, j, k vfx i +1, j, k Ax i, j, k Ax i +1, j, k ⎤ ⎥ =⎢ ⎢⎣ vfx i, j, k Ax i, j, k + vfx i +1, j, k Ax i +1, j, k ⎥⎦

(20)

⎤ ⎡ 2vfy i, j, k vfy i, j +1, k Ay i, j, k Ay i, j +1, k ⎤ ⎥ =⎢ ⎢⎣ vfy i, j, k Ay i, j +1, k + vfy i, j, k Ay i, j +1, k ⎥⎦

(21)

⎡ Txadv i + 1 , j, k =⎢vfx Ax ⎣ 2

i + 1 , j, k ⎥ ⎦ 2

⎡ Tyadv i, j + 1 , k =⎢vfy Ay ⎣ 2

i, j + 1 , k ⎥ 2 ⎦

Once the transmissibility term is determined, the residual equation can be re-written in the following form. n +1 n +1 n +1 +1 R (Ci, j, k )=Txturb i + 1 [Cin+1, j, k − Ci, j, k ]+ Txturb i − 1 [Ci −1, j, k − Ci, j, k ]+ Tyturb j + 1 2

2

2

n +1 n +1 n +1 [Cin, j+1 +1, k − Ci, j, k ]+ Tyturb j − 1 [Ci, j −1, k − Ci, j, k ]+ Tz turb k + 1 2

n Cin, j+1 , k − Ci, j, k

∆t +Ky Ay

= K Ax

V i, j , k

Cin, j+1 +1, k

−

x 2Cin, j+1 ,k

+1 n +1 n +1 Cin+1, j, k − 2Ci, j, k + Ci −1, j, k

n +1 [Cin, j+1 , k +1− Ci, j, k ]+ Tz turb k − 1 2

∆x

+

Cin, j+1 −1, k

∆y

+Kz Az

2

∆z

+Kz Az

∆x

n +1 n +1 Cin, j+1 , k +1 − 2Ci, j, k + Ci, j, k −1

∆z

−

∆t

V i, j , k

2

(22)

Once these finite-difference equations are solved and the new proppant concentration of the suspension calculated for a time step, the program starts to calculate the mass flux donated by saltation and creep. The approximation of saltation and creeping can be presented as follows:

n +1 n +1 Cin, j+1 +1, k − 2Ci, j, k + Ci, j −1, k

∆y

Cpbed in, j+1 − Cpbed in, j, k ,k

1 − v Ax [Cin, j, k −Cin−1, j, k ] 2 fx

∆t

Vi, j, k =qCpsalt n +1 Vi, j, k −qCpsalt n +1 Vi, j, k +qCpcreep n +1 Vi, j, k i −1, j, k

−qCpcreep n +1 Vi, j, k i, j , k

1 − v A [Cin, j, k −Cin, j −1, k ]−qCsin, j+1 V +qCp n +1 Vi, j, k , k i, j , k i, j , k 2 fy y n Cin, j+1 , k − Ci, j, k

∆t

i, j , k

The residual equation can be written as:

+Ky Ay

n Cin, j+1 , k − Ci, j, k

Vi, j, k +qCp n +1 Vi, j, k

(15)

n +1 n +1 +1 Cin+1, j, k − 2Ci, j, k + Ci −1, j, k

2

n +1 [Cin, j+1 , k −1− Ci, j, k ]−

n +1 n +1 n +1 n +1 +1 Vi, j, k −Txadv i − 1 [Cin−1, j, k − Ci, j, k ]− Tyadv j − 1 [Ci, j −1, k − Ci, j, k ]− qCs i, j, k

n +1 n +1 Cin, j+1 , k +1 − 2Ci, j, k + Ci, j, k −1

1 1 − v Ax [Cin, j, k −Cin−1, j, k ]− v A [Cin, j, k −Cin, j −1, k ]−qCsin, j+1 V +qCp n +1 Vi, j, k , k i, j , k i, j , k 2 fx 2 fy y

R (Cin, j+1 , k )=K x Ax

(19)

Similar to turbulent dispersion, the transmissibility term for advection can be also expressed in a harmonic form in the following equations.

(13)

∆z

⎤ ⎥ ⎥ i, j , k + 1 ⎦ 2

⎤ ⎡ 2Kz i, j, k Kz k +1 Az i, j, k Az i, j, k +1 ⎥ =⎢ ⎢⎣ Kz i, j, k Az i, j, k ∆z i, j, k +1 + Kz i, j, k +1 Az i, j, k +1 ∆z i, j, k ⎥⎦

∆y

Cin, j+1 , k +1

(18)

In z-direction:

2∆y

1 1 Vi, j, k + v Ax [Cin+1, j, k −Cin−1, j, k ]+ v A [Cin, j, k −Cin, j −1, k ] 2 fx 2 fy y

+qCsin, j+1 V −qCp n +1 Vi, j, k = , k i, j , k +Ky Ay

−

Cin, j −1, k

n +1 n +1 +1 Cin+1, j, k − 2Ci, j, k + Ci −1, j, k

+Ky Ay

n Cin, j+1 , k − Ci, j, k

Cin, j, k

⎤ ⎥ ⎥ 1 i, j + , k ⎦ 2

i, j , k

i −1, j, k

(23)

Bedload transport is an effect observed in the settled bed domain. Therefore, the bedload mass transport is separated from the slurry suspension domain. The bedload transport term is coupled explicitly to the model on the settled bed region. Once the velocity distribution is obtained, the mass due to bedload transport can be calculated. The proppant concentration of the settled bed domain is then re-distributed depending on the bedload transport flux at each grid.

(16)

Transmissibility term is often used to more accurately capture mass transfer between discretized blocks. The transmissibility terms are defined as the magnitude of mass transfer corresponding to the change 4

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

O. Chang et al.

Fig. 3. Top view when slurry encounters an intersection in the fracture network. S1, S2 and S3 represent the portion of slurry entering the branches (Chang et al., 2016).

2.3.4. Slurry fractional flow at intersections of fracture network When proppant bearing slurry flows through an intersection of fracture network, depending on the properties of proppant and fluid, and the fluid flow velocity field, proppant mass fluxes entering the branches are not proportional to the fluid fluxes. For most of the cases, proppants prefer to enter the branch in the same direction instead of turning into intersected fracture branches. Therefore, a mathematical model is developed to predict the unequal fractional flow of proppant entering each encountered branch (Chang et al., 2016). An illustrative top view when slurry reaches an intersection is shown in Fig. 3. The model is coupled to the PTM to obtain accurate mass fluxes that enter fracture network branches during the treatment. The proportion entering each branch can be expressed as follows. The parameter S1 represents the characteristic length of proppant mass flux entering the major fracture branch 1. The parameter S2 represents the characteristic length of the proppant mass flux entering the secondary fracture branch 2 perpendicular to the major fracture. In this demonstration case, S2 is equal to S3 due to symmetry.

⎡ mlog (bvlx t + m ) ⎤ tsy S1=S3=⎢vlx t − ⎥ |0 andS2=wf 2−S1−S3 ⎣ ⎦ b

The factors are then applied to the advection transmissibility terms to obtain the proportional allocation of the upstream mass flux to each branch.

F1=

S1 , S1 + S2 + S3

F2=

S2 S3 , and F3= S1 + S2 + S3 S1 + S2 + S3

(25)

Once the proppant mass flux at each junction is calculated, the ratio of the mass flux versus total flux can be obtained. The ratio is treated as a multiplier to the transmissibility term of the advection term of the slurry domain that corrects the apparent total flux to the actual proppant mass flux. 2.3.5. Slurry front tracking scheme using time of flight The propped stimulated reservoir volume (PSRV) is defined as the total reservoir volume of matrix block connected with propped-open fractures; PSRV is therefore highly sensitive to the accurate prediction of the slurry front. Generally, grid block sizes can be reduced to more accurately capture the profile of modeled parameters of interest. However, the proppant transport model is coupled to a fracture propagation model and a reservoir model that is already computationally expensive; resorting to grid size reduction to improve accuracy would quickly become computationally intractable. To avoid such excessive calculation, a “time of flight” method was developed and applied to accurately predict slurry front without reducing grid block size. This time-of-flight scheme is often used to track front in streamline simulation of gas or polymer flooding. First the streamline is constructed by using traditional Cartesian blocks or other analytical methods and assumed to be unchanging through a relatively longer period. The velocity field along the streamline is then constructed according to the pressure field solved by finite-difference type of flow model in the reservoir. To capture the front of injected fluid, BuckleyLeverett flow model along a streamline is introduce to obtain the saturation front (Buckley and Leverett, 1942). The front is called a Buckley-Leverett shock front and the saturation profile following the

(24)

wf 1

where tsy= v , CD is drag coefficient of a round particle ranging from 2 0.44 to 1.5, m is the mass of a proppant particle, vlx is the fluid flow velocity entering secondary branch in x direction, d is the diameter of a proppant particle, ρf is the density of fracture fluid, wf 2 , is the width of the major fracture, wf 1 is the width of the secondary fracture, v2 is the fluid flow velocity in the major fracture, v1 is the fluid flow velocity in π the secondary fracture, and b = CD 8 d 2ρf . The fraction that enters each branch can be obtained as a ratio of the proportion that enters a branch to the total flux. F1 represents the proportion of proppant mass flux entering the major fracture. F2 and F3 represent the proportions of proppant mass flux entering the secondary fracture branches. The PTM reads the velocity field and fracture geometry and calculates the fractional factors at each intersection. 5

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

O. Chang et al.

blocks. To smoothen the discrete velocity profile, a piece-wise linear approximation of the velocity field in each direction, for each block, is used.

vx =vi −1, j +mx (x−x 0 ); mx =

vy=vi, j −1+m y ( y−y0); m y=

(vi, j − vi −1, j ) ∆x

(26)

(vi, j − vi, j −1) ∆y

(27)

After the velocity field is smoothed, the travel distance of a slurry front along the direction within the grid block can be expressed in the following equation for a given period of time.

Fig. 4. Interpolation of discretized velocity profile.

X=

1 ln[vxexp (mx ∆t −vx0 )]+X0 mx

(28)

Y=

1 ln[vyexp (m y ∆t −vy0 )]+Y0 my

(29)

By implementing this time of flight scheme, the location of a front can be predicted for each time step throughout the fracture network domain. The front indicates the maximum travel distance the slurry can reach. The concentration profile behind the front can be solved numerically using turbulence advection-dispersion, and mass balance equations stated in the previous section. Different from traditional streamline simulation using time of flight, the proppant transport model updates the discretized velocity field for every time step from the HFPM instead of a fixed velocity field in a longer global time step. Therefore, both a dynamic and accurate front and concentration profile

front follows a Buckley-Leverett wave form. Similarly, to capture the slurry front and proppant concentration behind the front, the time of flight scheme is implemented in the proppant transport model. First, in the fracture network domain, the discretized velocity and pressure profile are solved and obtained from the HFPM. To make the velocity field solution smoother (i.e., to capture a more accurate velocity field), an interpolation technique is applied to the fracture domain. Fig. 4 is a bird's-view of a discretized velocity field in an intersection of fracture network domain where the velocities are solved on the edge of the grid

Fig. 5. Flow chart of the proppant transport model.

6

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

O. Chang et al.

criteria are reached.

along the fracture network can be obtained by combining an analytical front tracking technique and a numerical mass transport solution, without need for extensive grid refinement that requires excessive computational time. Fig. 5 provides the flow chart of the proppant transport model. The model first requires input of fluid flow properties, fracture geometry, pressure distribution, initial proppant concentration and distribution and other essential parameters from the Hydraulic Fracture Propagation Model. Fluid flow velocity, wall shear stress, diffusion factors and other parameters related to particle movement are calculated from the input parameters in the Proppant Transport Model. User-defined values can be selected to describe a scenario of interest. The Jacobian matrix of residual equations for each cell is formulated and calculated using initial input parameters. The matrix is iteratively solved until the summation of difference of each cell is smaller than the assigned tolerance. If it converges for a time step, mass balance check will be calculated to ensure proppant mass injected is the same as proppant mass remaining in the fracture network. Proppant concentration and distribution is updated after convergence. The updated profiles are sent back to the HFPM for the calculation of the next time step. The entire process continues until the designed total treatment time is reached.

J (Ci, j, k )[∆Ci, j, k ]=−F (Ci, j, k ) ∂F ∂CC ∂F ∂CW

∂F ∂CC ∂F ∂CW

⋱

⋱

⋱

⋱ ⋱

∂F ∂CN

⋱

⋱ ⋱ ⋱

⋱

∂F ∂CA

⋱

⋱

⋱

⋱

F1,2,1

⋱

⋮ ∆C2,1,1

⋱ ⋱

⋱

⋮ ∆C1,1,2

⋮ F1,1,2

⋱

⋱

⋮

⋮

⋱ ⋱ ⋱

⋱

=

⋮ F2,1,1

2.3.7. Matrix storage schemes The configuration of the Jacobian matrix for both Proppant Transport Model and fracture closure function contains a number of zeros which consumes CPU memory and therefore affects the speed of calculation for cases where the domain has a very large number of grid blocks. For the Proppant Transport Model the grid block size is divided into smaller blocks in order to capture the front of the slurry. In addition, the complex multi-phase fracture closure function will both create extensive populates sparse matrix as a coefficient matrix. In order to overcome potential memory limitations, that the size of the matrix increases logarithmically and interrupts the modeling, a matrix size reducing technique is introduced to the original Jacobian matrix. The sparse matrix storage scheme is implemented for the original Jacobian matrix system. Several different storage techniques can be used to reduce the matrix size. The sparse matrix operation techniques are summarized in the book “Matrix Computations3rded” by Golub and Van Loan (1994). The Compressed Sparse Row (CSR) is a matrix storage technique that transforms a matrix with equal row and column into three storage vectors. The first vector records indices defining where each row starts. The second vectors stores the coordinate indices of the column location of the corresponding value. The third vector captures the values of non- zero components of the original Jacobian matrix. The compression significantly reduces the size of the matrix by eliminating zeros and creates pointers to locate each non- zero element. To implement the compression of matrix with more flexibility, instead of using CSR in which the compression records the starting location of each row, the coordinate list (COO) compression method can be used to achieve the same goal. The COO method transforms the Jacobian matrix into three columns: the first vector records the column location indices, the second vector stores the row indices, and the third column records non-zero elements of the matrix. After compression, the matrix is transformed into three vectors with the following expression.

1 1

2 1

4 1

1 2

∂F ∂CC

∂F ∂CE

∂F ∂CS

∂F ∂CW

…… …… ……

This three vector configuration is more compatible with the PTM algorithm because the same coordinating indices of row and column are always used, so there is no need to re-coordinate the transformed matrix to the solver system.

⋱ ⋱

∆C1,2,1

⋱

⋱

Column index Roll index Matrix Coefficient

⋱ ⋱

⋱

k +1=C n +1k +∆C n +1k +1until Cin, j+1 ∑ −F (Ci,j,k ) < ε ,k i, j , k i, j , k

⋱ ⋱ ⋱

⋱ ⋱

⋱

⋱

F1,1,1

⋱

⋱ ⋱

∂F ∂CB

⋱

⋱ ⋱

⋱

⎤ ⎥ ⋮ ⎥⎥= ∂F (Cm, n, l ) ⎥ ∂Cm, n, l ⎦ ∂F ∂CE

⋱

∂F ∂CA

0

∂F ∂CS

⋱

∆C1,1,1

∂F ∂CB

∂F ∂CS

⋱

∂F ∂CN

2.3.6. Solution methodology and matrix operation In the proppant transport system, models are coupled together to model the dynamic proppant distribution within the fracture network domain. Each model has its own mathematical characteristic in terms of residual equation. The numerical solution methodologies for the Proppant Transport Model and the Fracture closure function are introduced in the following sections. According to the mathematical characteristics of the two models, the system of equations has nonlinearity. The Proppant Transport Model has moderate non-linearity due to the turbulent diffusion coefficient terms with proppant concentration as a function of particle suspension concentration. The fracture closure function exhibits strong non-linearity due to the transmissibility term. In the transmissibility term, volumetric factor, relative permeability, width of the fracture and others parameters are functions of dependent variables of pressure and saturation. To address this nonlinearity, the system of equations is written into generalized NewtonRaphson form to iteratively approach the solution of the system of nonlinear equations. The Jacobian matrix of the Proppant Transport Model is constructed using generalized Newton-Raphson procedure. For a three dimensional, single proppant component model, a hepta-diagonal Jacobian matrix is constructed such that each diagonal line is the partial differential of the residual equations with respect to the concentration change in a direction.

⎡ ∂F (Ci, j, k ) ⎢ ∂Ci, j, k ⋯ J (Ci, j, k )=⎢ ⋮ ⋱ ⎢ ⎢ 0 ⋯ ⎣

∂F ∂CE

3. Model validation For validation, the model developed herein is simplified into a onedimensional case and benchmarked with the analytical solution provided in Thongmoon et al. (2008). Results of that comparison are

To solve the change of concentration through time and location, the residual vector and vector of concentration variation are lined up and iteratively approach the solution for each row until the convergence 7

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

O. Chang et al.

Table1 Input settings of the base case.

Fig. 6. Model validation with analytical equation in 1D.

shown in Fig. 6. The results generated from the PTM are compared to the analytical solution. The error between the simulated results and the analytical solution is defined as: (concentration simulated–concentration analytical)/(concentration analytical)×100%. The average error between the PTM results with 30 grid blocks is 10.1% compared to the analytical solution. The average error between the PTM results with 120 grid blocks is 0.15% compared to the analytical solution. In the frontal area, the absolute difference is actually very small but the relevance is larger, which is caused by the tolerance level. This causes the average error to be larger. Since the absolute error is very small, it has minimal impact on mass balance accuracy. Finer gridding size is able to reduce the difference between the finite-difference model and the analytical solution to negligible level. The analytical solution of the differential equation can be expressed as:

⎛x − v t⎞ ⎛ x + v t ⎞⎤ ⎛ xv ⎞ 1⎡ x x ⎟⎟ +exp ⎜ x ⎟ erfc ⎜⎜ ⎟⎟ ⎥ C (x, y)= ⎢erfc ⎜⎜ ⎝ Kx ⎠ 2 ⎢⎣ ⎝ 4Kx t ⎠ ⎝ 4Kx t ⎠ ⎥⎦

Parameters

Unit

Base case

Proppant size Fluid viscosity Injection rate Horizontal stress difference Natural fracture spacing initial reservoir pressure reservoir permeability porosity sig (max-min) Young's modulus Poisson's ratio Fracture height Slurry concentration

mesh cp bbl/min psi ft psi md % psi psi fractional ft lb/gal

40/70 5 17 100 20×10 4000 0.0001 6 100 3.00E+06 0.2 100 1

reference case, as the average normalized proppant packing of the settled bed. The height of the settled proppant bed is also very close to the reference case. The PTM model does not capture the fluid flow velocity passing through the perforation. However, the PTM captures the turbulent dispersion near the perforation to model the entrance effect. Due to the above reason, the proppant concentration near the perforation is a little different. However, a short distance away from the perforation, where the entrance effect of perforation is mitigated, the PTM generates accurate proppant concentration and distribution prediction compared to the reference model. 4. Results and discussion 4.1. Results and analysis of the base case In order to compare the effectiveness by tuning the input engineering and natural parameters, a base case is set up as a reference. Table 1 shows the base case settings. In this set of simulation studies, a vertically confined fracture is generated by the HFPM that the fracture propagates only in the pay zone without height growth. In the following figures of results, the 3-D results generated from the model are converted into 2-D results on the X-Y plane due to the graphing algorism in the program. The results are averaged vertically (in zdirection) and presented on the figures. Fig. 8 and Fig. 9 show proppant concentration and fracture conductivity of the base case. In the simulation, pad fluid is injected for the first 20 min and followed by slurry for 60 min. The total treatment period for pressure pumping is 80 min. Proppant concentration is ranging from 0 to 0.12 lbm/ft2 . 73% of the created fracture

(30)

where Kx =0. 01, vx =0. 1, t=5, and dx=0. 1. The PTM is further validated with the results provided by Shiozawa and McClure (2016). With the same input provided from the paper, the PTM is able to generate a very identical result with the published reference (Fig. 7). The normalized proppant concentration generated from the PTM in the suspension is around 0.1, which is very close to the reference. The normalized packing of the settled proppant bank from the reference case is not a single value. In order to calculate the height of the settled proppant bank, we use 0.88, by taking the spatial average of the normalized proppant packing of the settled bed of the

Fig. 7. Model validation with the result provided by Shiozawa and McClure (2016).

8

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

O. Chang et al.

Fig. 8. Proppant concentration of the base case.

Fig. 10. Fracture conductivity of case1 (20/40 proppant).

Fig. 9. Fracture conductivity of the base case.

Fig. 11. Fracture conductivity of case1 (20/40 proppant).

network volume is propped-open. The PSRV for this case is 5,620,000 ft 3. The average fracture conductivity is 6.2 md-ft. The fracture conductivity is calculated by the following equation.

Cf =k f whf

4.3. Effect of proppant size In the case1, 20/40 mesh proppant is modeled to compare with the one modeled in the base case using smaller 40/70 mesh proppant and the results are shown in Fig. 10 and Fig. 11. In the figures proppant concentration reaches up to 0.27 lbm/ft2 and fracture conductivity to 324 md-ft. In this case, settling of the proppant is much significant than the base case using smaller sized proppant. The average size of 20/40 proppant is 600 µm that causes the settling terminal velocity to be around 7 times larger than the 40/70 (220 µm) case. Substantial proppant settlement to the settled bank increases the overall proppant concentration in a unit cell. A huge portion of proppant concentration is contributed by the accumulation of proppant than in the suspension for this case. The average proppant concentration is higher than the base case by 30 times. Higher proppant concentration, due to substantial proppant settlement, corresponds to larger healed width and higher proppant pack permeability is achieved because of largersized proppant is selected. The above two reasons are accounted for the results. However, in the other hand, the size of the PSRV decreases from 5,620,000 ft 3 to 321,500 ft 3 and only 7% of the total created fracture volume is propped-open. Severe setting is also accounted for the outcome for this case. The proppant particles are not capable of being transported and carried into the further part of the fracture network before settling to the bottom. Therefore the propped-open area is restricted remaining a large portion of created fracture volume vacant. Average fracture conductivity and proppant concentration is increased by several times by using larger proppant. However, the PSRV is severely reduced. A possible implication can be brought that by using larger proppant the less total amount of hydrocarbon can be extracted from the reservoir but the production rate could be higher due to higher fracture conductivity.

(31)

where Cf is fracture conductivity (md-ft), k f is proppant pack permeability (md), and whf is healed fracture width (ft). The PSRV is calculated as the formation volume contacted by the propped fracture. 4.2. Effect of engineering factors Several engineering factors such as proppant size and fluid viscosity impact the effectiveness of hydraulic treatment. In this set of numerical study, different proppant size and fracture fluid viscosity is modeled to evaluate its impact on PSRV. Table 2 is the input parameters for the case study. The yellow shaded sections are the controlling parameter. Table 2 Input settings of parametric study investigation engineering factors. Parameters

Unit

Base case

Case1

Case2

Proppant size Fluid viscosity Injection rate Horizontal stress difference Natural fracture spacing initial reservoir pressure reservoir permeability porosity sig (max-min) Young's modulus Poisson's ratio Fracture height Slurry concentration

mesh cp bbl/min psi

40/70 5 17 100

20/40 5 17 100

40/70 25 17 100

ft psi md % psi psi fractional ft lb/gal

20×10 4000 0.0001 6 100 3.00E+06 0.2 100 1

20×10 4000 0.0001 6 100 3.00E+06 0.2 100 1

20×10 4000 0.0001 6 100 3.00E+06 0.2 100 1

9

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

O. Chang et al.

Table 3 Input settings of parametric study investigation natural factors.

Fig. 12. Proppant concentration of case2 (25cp).

Parameters

Unit

Base case

Case3

Case4

Proppant size Fluid viscosity Injection rate Horizontal stress difference Natural fracture spacing initial reservoir pressure reservoir permeability porosity sig (max-min) Young's modulus Poisson's ratio Fracture height Slurry concentration

mesh cp bbl/min psi

40/70 5 17 100

20/40 5 17 100

40/70 25 17 100

ft psi md % psi psi fractional ft lb/gal

20×10 4000 0.0001 6 100 3.00E+06 0.2 100 1

200×100 4000 0.0001 6 100 3.00E+06 0.2 100 1

20×10 4000 0.0001 6 300 3.00E+06 0.2 100 1

4.5. Effect of natural factors In this section, impact of natural factors such as difference in in-situ stress and natural fracture spacing are modeled and evaluated. The highlighted parts listed in Table 3 are the controlling parameters tuned to investigate the impact due to the nature of the reservoirs. In the following paragraphs, quantitative evaluation is made to evaluate the effectiveness of a hydraulic fracturing applied to different reservoirs processing varying properties. 4.6. Effect of difference in horizontal stress Fig. 13. Fracture conductivity of case2 (25cp).

Difference in horizontal stress is an important factor affecting fracture geometry and fluid flow field in the fracture network. Results of the case are shown in Fig. 14 and Fig. 15. In case 3, difference in horizontal stress is increased to 300 psi. Proppant concentration ranges from 0 to 0.07 lbm/ft2 . Fracture conductivity reaches up to 14 md-ft with average conductivity of 4 md-ft. 51% of the created fracture volume is contacted by proppant. The volume of the PSRV is 2,342,000 ft 3. Observation can be made from the results that the fracture network geometry is significantly different from the base case with 100 psi in horizontal stress difference. Direction of fracture propagation is elongated toward the maximum stress direction with longer major fractures and less lateral secondary fractures opened. Confined by the maximum in-situ stress, the absolute size of PSRV is decreased. In addition, the percentage of the propped-open volume is also decreased from 73% to 51%. The reason is because the created fracture network has an elongated geometry with 1.5 times longer in the maximum in-situ stress direction. However, by using the same proppant and fracture fluid viscosity, proppant particles has a limited reachable length. Therefore, the slurry is not capable of contact the additional fracture volume created in the maximum in-situ stress direction outside the reachable length, leaving a huge portion the

4.4. Effect of fracture fluid viscosity The effect of fracture fluid viscosity is also examined by the parametric study. In this case, a fracture fluid with higher viscosity (25cp) is modeled to investigate its impact and compared to the base casing with lower viscosity 5(cp). The results are shown in Fig. 12 and Fig. 13. Proppant concentration ranges from 0 to 0.06 lbm/ft2 and fracture conductivity reaches to 12 md-ft. The average fracture conductivity is reduced to 1.5 md-ft. The PSRV increase from 73–89% of the total created fracture volume. However, the exact volume of the PSRV is decreased to 1,541,500 ft 3. Several observations can be made based on the results. Settling velocity is decrease by the fluid viscosity. The percentage of terminal velocity reduction has an exact inverse relation with the proportion of the increase of fluid viscosity. In other words, in this case, the fluid viscosity is brought up 5 times higher, and therefore, the terminal velocity is then reduced to 1/5 of its original value. With the terminal velocity is reduced, fracture fluid is capable of carrying proppant further into the fracture network. Therefore we can observe that 89% of the fracture volume is contacted by proppant. Proppant concentration and average fracture conductivity is reduced accordingly because proppant is widely distributed in the suspension instead of restricted by a limited length. In the contrast, although the percentage contact by proppant is increased, the exact volume of propped-open volume is less than the base case. The reason is because fracture fluid with higher viscosity generates wider but shorter fracture network. Also more secondary lateral fractures are opened because of higher viscosity. Therefore the size of the PSRV is reduced but the intensity is enhanced. The results suggest that higher fluid viscosity mitigates the settling and achieves better efficiency. However, the total size of the PSRV is reduced but the intensity increases in the other way. Larger-sized treatment could be considered to stimulate larger reservoir volume since high fracture network intensity and favorable proppant transport ability can be achieved.

Fig. 14. Proppant concentration of case3 (300 psi).

10

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

O. Chang et al.

Fig. 15. Fracture conductivity of case3 (300 psi).

Fig. 17. Fracture conductivity of case4 (200×100 ft).

fracture network volume in propped. Moreover, very few lateral secondary fractures are opened due to high in-situ stress which causes low fracture network density. Due to this distinct fracture propagation behavior, to enhance the effectiveness of the treatment in terms of PSRV, smaller proppant may be encouraged to increase the maximum travel length of slurry and the ability to enter narrower secondary fractures.

size of PSRV is favorable due to larger contacted area which may imply preferable recovery.

4.8. Effect of bedload transport

Shale reservoirs possess various characteristics and natural fracture spacing is a critical factor affecting the PSRV. Fig. 16 and Fig. 17 show the proppant concentration and fracture conductivity of case 4. Proppant concentration reaches its highest at 0.1 lbm/ft2 . Fracture conductivity ranges from 0 to 21 md-ft with an average of 7.7 md-ft. 21% of the created fracture volume is propped-open. The exact size of the PSRV is 14,003,600 ft 3. Fracture network spacing has a great impact on the PSRV. Fracture geometry has a significant difference from the base case. The coverage of the fracture network increases almost an order in both x and y direction because fracture fluid creates fracture network through pre-existing fractures with huge spacing. The percentage of propped-open area decreases drastically from 73% to 21% due to limited slurry travel length. Although the extension of the fracture network is large, the maximum reachable length is controlled by the terminal settling velocity. Unlike case3, most of the unpropped area is in the maximum stress direction. In this case, not only the maximum in-situ stress direction but also the minimum in-situ stress direction is not contacted by the slurry leaving the two sides unpropped. In contrast, although the propped-open percentage is low, the exact PSRV volume increases due to large fracture spacing and higher fluid flow velocity that enhance the travel length of slurry. The density of the fracture network is also low due to the spacing. Smaller-sized proppant may be suggested in this case to maximize spatial transportability of proppant to maximize the effectiveness. The

Bedload transport, saltation and creeping, contribute to mass transport by removing settled particles from the bank toward downstream direction. Mass transport due to bedload transport is modeled in this section to quantify its impact on final proppant concentration and fracture conductivity distribution in the fracture network. Fig. 18 and Fig. 19 show the result when mass transport contributed by bedload transport is taken into consideration. Proppant concentration reaches 0.06 lbm/ft2 . Fracture conductivity achieves 13 md-ft with an average of 2.8 md-ft. 81% of the generated fracture volume is propped open. The actual size of the PSRV is 6,201,500 ft 3. Fig. 20 is the comparison of the base case and case 5. The red line indicates the furthest location where the proppant is able to reach. 795 ft is reached in the base case and 896 ft is reached in case 5. In this case, bed load transport contributes to 100 ft incremental propped length in the major fractures by remobilizing the settled proppant. The ratio of the propped area is increased to 81%. The average fracture conductivity is decreased due to the removal of settled proppant. In addition, proppant concentration gradient of the proppant in the fracture network is less deep compared to the base case. Bedload transport is capable of flattening the uneven distribution of proppant in the fracture network and enhances homogeneity in terms of fracture conductivity. Especially at the near wellbore region, bedload transport is effective in removing excessive settled proppant which could potentially cause a screen out event. To summarize the section, bed load transport is capable to increasing the size of PSRV and the effectiveness of the treatment. The process also decreases the risk of screen out at the near wellbore region by removing overbuilding proppant bank.

Fig. 16. Proppant concentration of case4 (200×100 ft).

Fig. 18. Proppant concentration of case5 (with bedload transport).

4.7. Effect of difference in natural fracture spacing

11

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

O. Chang et al.

Fig. 21. Proppant concentration of case6 (fractional slurry flow at intersections).

Fig. 19. Fracture conductivity of case5 (with bedload transport).

Fig. 22. Fracture conductivity of case6 (fractional slurry flow at intersections).

the concentrated distribution of proppant around the center fracture. the proportion of the propped-open volume is decreased because of the same mechanism that restricts the lateral transport ability of slurry. Uneven fractional slurry flow at intersections has a great impact on final proppant distribution. Proppant is placed in a centralized manner in the middle leaving the rest of the fracture network either stimulated by diluted slurry or unpropped. The behavior may decrease the final recoverable hydrocarbon volume due to reduced PSRV size and unfavorable fracture conductivity in rest of the fracture network. The risk of encountering screen out may be increased due to limited exit of proppant where proppant bridging and excessive bank building are favorable to occur under higher proppant concentration scenario. Fig. 20. Comparison of fracture conductivity between the (a) base case and (b) case5.

4.10. Effect of shut-in time

4.9. Effect of uneven fractional flow of slurry at fracture network intersections

After pressure pumping, the well is often shut-in for a period to transfer into production phase. Proppant concentration and distribution in the fracture network is changing depending on the shut-in schedule. In this section, the proppant transport model is linked to the fracture closure function in the HFPM to capture proppant distribution and fracture closure in the same time. Fig. 23 and Fig. 24 show the results after the well is shutin for 10 h where proppant in suspension has been closed on by fracture walls. In order to quantify the effect of shut-in, the results only show the fracture conductivity contributed by proppant in the suspension. Proppant concentration reaches 0.011 lbm/ft2 and fracture conductivity reaches 2.3 md-ft at the perforation. The average conductivity is 8 md-ft. 24% of the created volume is propped-open. The size of PSRV is 2,038,500 ft 3. After 10 h of shut-in, proppant concentration in the suspension and the size of PSRV decrease drastically due to settling of proppant. Almost two third of overall fracture conductivity is lost due to settling of proppant. Proppant in suspension in major fractures, except for a little area neat the perforation, is completely lost before being trapped by fracture walls. Fig. 25 is a zoomedin view of fracture conductivity in the fracture network near the perforation. We can clearly observe that proppant concentration in the fractures

According to Chang et al. (2016) uneven proportion of proppant enters different fracture branches depending on fracture network geometry, fluid flow field, and particle properties. The PTM is coupled to a fractional flow module to calculate the proportion of slurry entering the fracture braches in downstream directions. The effect of fractional slurry flow is modeled in the section to investigate its impact on PSRV. Fig. 21 and Fig. 22 present the results when fractional slurry flow is considered. Proppant concentration reaches 0.27 lbm/ft2 at the highest. Fracture conductivity achieves 54 md-ft at the perforation with an average value of 9.2 md-ft. 51% of created fracture volume is propped-open. The actual size of PSRV is 3,874,000 ft2 . Proppant distribution in fracture network has a huge difference compared to the base case. Proppant distribution is concentrated in the major fracture in the middle with high proppant concentration where very few proppant turns into secondary fractures. Besides the major fracture in the middle, most of the propped-open fractures are only stimulated by very little amount of proppant. The average fracture conductivity is increased due to 12

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

O. Chang et al.

viscosity could be used to mitigate proppant settlement as well as the shut-in time could also be decreased by immediate flowback to expedite the closure time to maintain a desirable PSRV connectivity. 5. Conclusion In this article, a three-dimensional, finite-difference proppant transport model (PTM) in hydraulic fracture network is developed. The model includes multiple particle transport physics such as turbulent dispersion, advection, saltation and creeping to more accurately simulate proppant transport process. Several numerical and mathematical techniques are also applied to improve the accuracy and efficiency of the model. The PTM is able to generate proppant concentration, distribution and final propped stimulated reservoir volume (PSRV) in a hydraulic fracture network when treatment schedule, engineering and geomechanical properties are given. With its wide applicability, the effectiveness of a hydraulic fracturing treatment in a naturally fractured reservoir can be evaluated by utilizing the model for a variety of scenarios. This model can serve as a tool to understand factors contributing to stimulation effectiveness and optimize fracture treatment design in naturally-fractured reservoirs. The impacts of several important factors are also evaluated in the article. The effect of each factor is summarized as follows.

Fig. 23. Proppant concentration of case7 (effect of shut-in).

1. Proppant size has first order impact on the PRSV. Larger proppant should only be used at the end of the treatment, otherwise the PSRV decreases significantly. 2. Higher fluid viscosity reduces proppant settling, the PSRV is increased due to longer travel length. 3. Large difference in in-situ stress results in an elongated fracture network, leaving a great portion unpropped. The efficiency of the PSRV decreases. 4. Uneven fractional flow of proppants at intersections reduces proppant concentration in the minimum in-situ stress direction significantly. 5. Shut-in time is critical to the PSRV. Potential huge loss of fracture conductivity in the major fractures can happen if permeability of the formation is low.

Fig. 24. Fracture conductivity of case7 (effect of shut-in).

Disclaimer This project was funded in part by the Department of Energy, National Energy Technology Laboratory, an agency of the United States Government, through a support contract with URS Energy & Construction, Inc. Neither the United States Government nor any agency thereof, nor any of their employees, nor URS Energy & Construction, Inc., nor any of their employees, makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

Fig. 25. Zoomed-in figure of fracture conductivity of case7 (effect of shut-in).

aligning in the maximum horizontal stress direction is almost completely lost. However, proppant concentration in the secondary fractures in the minimum horizontal stress direction is maintained to a certain value without severe loss in conductivity. The secondary fractures are subjected to maximum in-situ stress that the width is narrower. Hence, more proppant in suspension is trapped due to earlier closure of the secondary fractures. Despite that the secondary fractures maintain favorable fracture conductivity, the pathways of hydrocarbon production are impeded due to the absence of proppant in the major fractures resulted from late closure. Recovery could be undesirable because of poor connectivity of propped-open fracture network even if a favorable density of stimulated fractures is achieved. Smaller-sized proppant or higher fluid

Acknowledgment This work was performed as part of a collaborative research effort with the National Energy Technology Laboratory's Office of Research and Development, with support from the RES contract DE-FE0004000.

13

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

O. Chang et al.

Appendix A See Table A1. Table A1 Treatment schedule for the example case in Fig. 1. Parameters

Unit

Base case

Fluid viscosity Injection rate Horizontal stress difference Natural fracture spacing initial reservoir pressure reservoir permeability porosity sig (max-min) Young's modulus Poisson's ratio Pay height

cp bbl/min/per half cluster psi ft psi md % psi psi fractional ft

2 17 100 20×10 4000 0.0001 6 100 3.00E+06 0.2 100

is Presented at SPE, Low Perm Symposium, Denver, Colorado, May, 5–6. Mack, Mark, Sun, Juan, Khadilar, Chandra. 2014. Quantifying Proppant Transport in Thin Fluids-Theory and Experiments. Paper SPE 39428 presented at the SPE Hydraulic Fracturing Technology Conference and Exhibition, Woodland, Texas, Febuary, 4 - 6 Shiozawa, S., McClure, M., 2016. Comparison of Pseudo-3D and Fully-3D Simulations of Proppant Transport in Hydraulic Fractures, Including Gravitational Settling, Formation of Proppant Banks, Tip-Screen Out, and Fracture Closure, Paper SPE179132-MS is present at SPE Hydraulic FracturingTechnology Conference, The Woodland, Texas, USA, February, 9–11. Thongmoon, M., Tangmanee, S., Mckbbin, R., 2008. A comparison of splines interpolations with standard finite difference methods for one-dimensional Advection-diffusion equation. Int. J. Mod. Phys. 19 (8), 1291–1304. Wang, Z., Zheng, X., 2004. Theoretical prediction of creep flux in aeolian sand transport. J. Powder Technol. 139, 123–128. Weng, X., Kresse, O., Cohen, C., Wu, R., 2011. Modeling of Hydraulic-Fracture-Network Propagation in a Naturally Fractured Formation, SPE production and operations. 26 (4), 368–380. Yalin, M.S., 1963. An Expression for Bed-Load Transportation. Proc. ASCE 89, 221–250.

References Ahn, C.A., Chang, O., Dilmore, R., Wang, J.Y., 2014. A Hydraulic Fracture Network Propagation Model in Shale Gas Reservoirs: Parametric Studies to enhance the Effectiveness of Stimulation. Paper SPE-2014-1922580-MS is presented URTC conference, Denver, Colorado, August, 25–27. Blyton, C.P.G., Gala, D.P., Sharma, M.M., 2015. A Comprehensive Study of Proppant Transport in a Hydraulic Fracture. Paper SPE-174973-MS is Presented at SPE, ATCE, Houston, Texas, September, 28–30. Buckley, S., E., Leverett, M., C., 1942. Mechanism of fluid displacements in sands. Trans. AIME 146, 107–116. Chang, O., Dilmore, R., Wang, J.Y., 2016. Physics of proppant transport in hydraulic fracture network. ASME J. Energy Resour. Technol. ( Rev.). Golub, G., H., Van Loan, C., F., 1994. Matrix Computation 3rd ed. The Johns Hopkins University Press, Baltimore, USA. Han, J., Yuan, P., Huang, X., Zhang, H., Sookprasong, A., Li, C., Dai, Y., 2016. Numerical Study of Proppant Transport in Complex Fracture Geometry. Paper SPE-180243-MS

14