8/10/2019 SPE 152541_12_

1/20

8/10/2019 SPE 152541_12_

2/20

2 SPE 152541

In recent years, several types of hydraulic fracturing simulators dedicated to shale reservoirs were developed (Meyer et al.

2011, Xu et al., 2010), in an attempt to account for the complex interaction between hydraulic fractures with a discrete

natural fracture network. In the work presented in this paper, the hydraulic fracturing simulations were conducted with the

UFM (Unconventional Fracture Model) which explicitly considers discrete and complex sets of natural fractures, stress

shadow effect and proppant transport (Weng et al. 2011, Kresse et al., 2010). The challenge for production forecast is that the

resulting network is often far too complex to be exported into a conventional reservoir simulator. It is only very recently that

complex fracturing simulators and reservoir simulators for shale have been coupled (Cipolla et al., 2011), but this is still a

technology under development. Today, stimulation engineers do not have a simple tool to forecast the impact a hydraulicfracturing treatment will have on the production of a shale gas reservoir.

This paper primarily presents a workflow to estimate the potential production from a simulated fracturing treatment in a

shale reservoirs. The motivation of this work was to build a production model, inside the fracturing simulator, that could

estimate the production with minimal additional input from the user. Therefore, this new tool is aimed primarily at

stimulation engineers seeking immediate feedbacks on their treatment design.

The workflow consists of the extrapolation of the results from the UFM simulation which are processed with a novel

production model. The model considers explicitly the flow inside the DFN (Discrete Fracture Networks) formed by the

hydraulic fracture network, and the flow from the matrix to the fractures as a source term given by an analytical solution.

This solution is based on the solution of the pressure diffusion equation inside the matrix which assumes a constant pressure

inside the fracture, an invalid assumption under realistic conditions where the bottom hole pressure may take years to diffuse

inside a low conductivity fracture network. This issue is solved with an original method that extends the validity of the semi-

analytical model to the full range of fracture conductivities for application in real cases. This is done by calculating a time

delay based on the local mass balance of each matrix block through an inversion algorithm. This model has been validated

against a commercial reservoir simulator to illustrate its capability to provide accurate results for a large range of fracture

conductivity. The paper presents the details of the model and discusses the benefits and limitations of this approach compared

with standard reservoir simulators. In the last part of the paper, we illustrate the type of insight this new tool provides

regarding the design of a treatment (proppant selection and pumping rate).

Modeling production from a hydraulic fracture networkSimulating the production from a complex fracture network is a difficult problem that has been approached previously in

different manners. The differences of properties between the fracture and the rest of the reservoir (matrix) lead to the

development of dual-porosity models (blatt and Zeltov, 1960) which consider two coarse grids connected to each other, one

for the fracture network and another for the matrix. Examples of this method applied to hydraulic fracturing can be found in

Zang et al. (2009) and Cipolla et al. (2010). This method needs some averaging of properties (especially for the fracture

network) and simplifications in order to model the mass transfer term between the two media. It is useful for naturally

fractured reservoirs on a large scale but not appropriate for the -wellbore effect of the fracture network created by hydraulicfracturing. The opposite approach is to have one medium that contains both the fracture and the reservoir, but in this case the

numerical grid required has to be extremely refined near the fractures, leading to extreme computational cost. Also, this

method often needs flexibility on the gridding (unstructured mesh generation) that requires complex algorithms. An example

of this technique is given in Cipolla et al. (2011) that coupled the UFM simulator with an unstructured reservoir simulator to

estimate production. A middle ground is the use of dual porosity equations on a DFN (Basquet et al., 2003). With this

approach the difficulty to simulate the flow from the matrix to the fractures remains. One of the main issues, in particular for

compressive reservoir fluids like gas, is to consider the production history from each matrix block into the DFN. One way is

to grid the matrix block at the cost of having additional unknowns in the system of equations to solve (Gong et al, 2010), or

use an analytical solution (Xu et al., 2011) for the matrix blocks. Analytical solutions can be derived from several techniques,

like for example the use of a Greens function like in the point source methods (Lin and Zhu, 2011), or obtained by a Laplace

transform of the continuity equation (Cinco-Ley et al, 1982, Xu et al., 2011), like in this paper. If this solution considers a

transient fracture pressure, a complex expression is obtained that requires numerical integration in time. But if we consider a

constant fracture pressure, then we can obtain an expression of the flow rate between the matrix and the fracture that is linearwith respect to pressure. This simple and convenient solution is valid only for very conductive fractures where the pressure

variations inside the DFN are negligible (at constant wellbore pressure). The method presented in this paper extends the

validity of this analytical solution to the full range of fracture conductivity for general application of the model in hydraulic

fracturing in a naturally fractured reservoirs.

Description of the workflowThis workflow has been developed to minimize the user intervention between the simulation of hydraulic fracturing and the

visualization of the production simulation. It comprises five steps:

1. Simulation of the hydraulic fracturing of the well. This will generate the data describing the resulting DFN.2. Automated extrapolation of the UFM data to create a DFN in a format that can be used by the production model

(Fig. 1). This format considers a unique averaged value for each property at each fracture branch. The fracture

8/10/2019 SPE 152541_12_

3/20

SPE 152541 3

branches are defined as the plane that connects two points of fracture intersection, or a fracture intersection and a

fracture tip. These properties are :

Spatial coordinates at the extremity of the branch

The averaged conductivity

The averaged height

The averaged reservoir pressure at the branch location

The averaged reservoir permeability at the branch location

Figure 1: Example of the creation of the simulated hydraulic fracture network from the UFM into an equivalent DFN network as usedin the present model (red dots are the connection points, blue l ines are the branches)

The description of the DFN by nodes (intersections) and branch is necessary because the present model calculates

the pressure at the nodes and uses the branch to both connect the nodes and calculate the production from the

adjacent matrix blocks.

3. Automated evaluation of the depth of drainage is carried out in the matrix block in front of each fracture branch(Fig. 2). This has to be done in such way that the total and actual volume of a given matrix block (not in contact

with the reservoir boundaries) may be drained.

Figure 2: Example of the calculation of the matrix depth to be depleted on each side of all branches (dashed red lines are the depthof the matrix block)

This calculation is similar to the method of Xu et al. (2011) to calculate equivalent block length for orthogonal

fracture network. The Fig 3illustration explains the definition of the equivalent block length.

Hydraulic Fracturing Simulation

Fluid

Slurry

Bank

Fracture branch

Fracture tip or

intersection

Automated

Data Export

Production network

Depth of matrix

block to be depleted

Fracture branchFracture tip or

intersection

8/10/2019 SPE 152541_12_

4/20

4 SPE 152541

Figure 3: Calculation of the equivalent block depth

To give an example, for a square matrix block surrounded by four fractures of equal length, we assume that each

quarter of the block can be depleted by the fracture branch it is in contact with. Because we consider a linear flow

from the matrix to the fracture (to be detailed later), we assume that this quarter of the matrix block has the length of

the fracture branch. Therefore, the depth of this quarter of the block has to be equal to one-fourth of the block

length for the total volume to be depleted to be the same. For more complicated block shapes, the technique is morecomplex.

4. This is the only step that cannot be automated. The user is required to define production parameters, such as:

Bottom hole Pressure (BHP)

Reservoir fluid viscosity at reservoir conditions

Reservoir fluid compressibility at reservoir conditions

Duration over which production is to be simulated

5. Simulation of the production and visualization of the results. The user can not only visualize the production declineand the cumulated production, but also the dynamics of the pressure field in both the fracture network and the matrix

blocks (Fig. 4).

Figure 4: Example of visualization of the production data in time (140 days). Top left: Production rate. Bottom left: Cumulatedproduction. Right: Reservoir pressure (colored squares) and pressure inside the fracture network (blue lines).

Description of the method to calculate the productionThe description of the method comes in four parts. The first one describes the main equations and their analytical solutions.

Fracture branch 1

Volume of the matrix

block depleted byfracture branch 1

Equivalent block depthto be depleted

Linear Flow

ApproximationL

L

L/4

0.0+00

5.0+0

1.0+0

1.5+0

2.0+0

2.5+0

3.0+0

3.5+0

4.0+0

0 200 400 00 00 1000

100000

100000

300000

500000

00000

00000

1100000

1300000

1500000

0 200 400 00 00 1000

8/10/2019 SPE 152541_12_

5/20

SPE 152541 5

The details of the calculation are provided in Appendix 1. The second part explains the limitation of the analytical part of the

model to high-conductivity fractures and provides an example with a detailed description of what happens in a single branch

for both high- and low-conductivity cases. From this example, the third part explains the innovative numerical algorithm that

addresses this issue and extends the validity of the analytical model.

The assumptions of the model are:

- Single phase flow- Constant bottom hole pressure

- Linear flow from the matrix to the fracture- No desorption mechanism

Governing equations and analytical solutions

The continuity equation for compressible fluid in porous media is applied to both the matrix and the fractures. Inside the

fracture network, the continuity equation can be re-written as follows:

, = , ( 1)Qmfis the flow rate from the matrix to the reservoir, Qfis the flow rate inside the fracture, is the fluid density and xfis the

axis along the fracture. We consider that the fracture permeability (conductivity divided by the width) is so large that the

transient term of the continuity equation may be neglected over the time scale considered for production simulation (from

days to years). We assume Darcy flow inside the fracture network.

, = 2 ( 2)Pf is the pressure inside the fracture, C the conductivity, T the temperature. The function m is the pseudo pressure (Al-

Hussainy, 1966).

= 2 ( 3)Inside the matrix, the continuity equation for a compressible fluid takes the following form.

= ( 4)Pmis the pressure inside the matrix, kmis the matrix permeability, ctis the fluid compressibility, is the viscosity, Zis the

volume factor and mis the porosity of the matrix. For simplification, Eq. 4can be re-written as follow.

= 1 ( 5)ais defined in Eq. 6.

=

( 6)

To calculate Qmf, we solve Eq. 5, withxmthe coordinate along an axis orthogonal to the fracture.

8/10/2019 SPE 152541_12_

6/20

6 SPE 152541

Figure 5: Illustration of the coordinates in the fracture and matrix block.

The solution of Eq. 5is found by using a Laplace transform (Jeannot, 2009; Bello, 2009), as explained in Appendix 1

, = _ 2 1 , 1

,

,

( 7)

This solution being linear in pressure, Eq. 2 is easy to integrate and solve.

= _ + , , + , , ( 8)

= 1 , 1

,

,

( 9)

Knowing the pressure distribution inside the network, the production rate can be calculated with Darcys law.

The advantage of this calculation is that it is simple to implement and carries no constraint on the time step. But it has a

major underlying limitation: the flow from the matrix is based on the assumption that the pressure inside the fracture stays

constant. If the production is done at a constant BHP and the conductivity is extremely high, this assumption is correct. But in

reality, only a fraction of the fracture branches of the network are likely to have high conductivities (Cipolla et al., 2009). In

the following we illustrate this issue and explain the proposed solution.

Analysis of the validity of the analytical solution

To study the validity of the analytical model for different values of fracture conductivity, we analyze the evolution of the

pressure and production at a single fracture branch of a complex network, consisting of two sets of equally spaced parallel

fractures as shown in Fig. 6:

xf

xm0

L

-L

wellbore

Branch 1Matrix block 1

Qmf,1 : Flow rate frommatrix block1to

fracture branch 1

DFN (Discrete Fracture Network)

8/10/2019 SPE 152541_12_

7/20

SPE 152541 7

Figure 6: Description of the single branch in the DFN that will be analyzed.

For a very high conductivity in the fracture network (2,500 mD.ft) in a reservoir of only 0.001 mD, the BHP is almost

instantaneously diffused inside the network and from there the variation of pressure inside the DFN can be neglected,

compared to the pressure differential between the initial reservoir pressure and the BHP (Fig. 7).

Figure 7: Pressure inside the DFN at two different time of production for high-conductivity DFN

This is illustrated in Fig. 8where the pressure inside the selected fracture branch 1 (see Fig. 6) can be considered constant

during three years of production.

Figure 8: Pressure and initial time during three years of production in case of high conductivities

The consequence of this almost constant pressure in the DFN is shown in Fig. 9where the cumulated volume produced (from

the matrix block 1 into fracture branch 1) converges towards the maximum volume recoverable as defined by mass balance

(or initial volume of gas in place). Because we are considering compressible fluids, in this example, the measure of volume is

done at constant pressure condition (BHP).

Figure 9: Cumulated production from branch 1 versus time in case of high-conductivities

This convergence means that the analytical solution verifies mass balance in this case of a high-conductivity DFN. The same

analysis for a low-conductivity DFN (12.5 mD.ft) has a different conclusion. In Fig. 10, it is shown that the pressure in the

Day 1 Day 21

Initial reservo ir pressure

0

200

400

00

00

1000

1200

1400

0

0.2

0.4

0.

0.

1

0 500 1000 1500

( 0)

0

4000

000

12000

1000

0 500 1000 1500

1

Maximum recoverable volume

8/10/2019 SPE 152541_12_

8/20

8 SPE 152541

DFN varies a lot compared to the pressure range of the problem ([BHP, initial reservoir pressure]).

Figure 10: Pressure inside the DFN at two different times of production for low-conductivity DFN

This pressure variation is seen on the pressure recorded in the fracture branch 1 versus time (Fig. 11).

Figure 11: Pressure and initial time during 3 years of production in case of low conductivities

This variation of the pressure in the DFN means that the assumption of constant pressure boundary condition in the analyticalsolution is not valid. The consequence is that the calculated flow rate from the matrix is underestimated and the mass balance

is incorrect as illustrated in Fig. 12. .

Figure 12: Cumulated production from branch 1 versus time in case of low conductivities

The low diffusivity in the fracture network results in a delay in the production of the block depending on how far (or how

connected) it is from the wellbore. This observation is the starting point for the method to extend the validity of the analytical

solution to low conductivity fractures.

Extension of the validity of the analytical solution

The idea to extend the validity of this analytical solution is to modify the initial time t0in such a way that the volume

produced so far from the matrix block is equal to the volume that would have been produced according to the analytical

Day 1 Day 110

Initial reservoir pressure

0

200

400

00

00

1000

1200

1400

0

0.2

0.4

0.

0.

1

0 500 1000 1500

( 0)

Variation of the boundary condition

0

4000

000

12000

1000

0 500 1000 1500

1

Error on volume

Maximum recoverable volume

8/10/2019 SPE 152541_12_

9/20

SPE 152541 9

solution under current pressure conditions inside the DFN. By doing this search at every time step and on each side of each

fracture branch, the analytical solution is forced to satisfy mass balance. The search for the t0begins by defining the objective

function Fto minimize.

,, = , , , , ,,

= 0 ( 10)Mtot is the volume produced at time t from the matrix block on the side k of the fracture branch. It is compared to theintegration of the flow rate from the matrix over the length of the fracture branch and from the initial time t0,kto t. The search

fort0,ksuch that Fis equals to zero, we use the iterative algorithm of Newton-Raphson as described in Eq. 11..

,= , ,, , ,, , ,, = ,, , ( 11)

The derivative of the function F0,kis calculated by a numerical gradient. If t0,kmeets its time boundaries the optimization uses

the bisection method. This optimization algorithm is very efficient because the solution from the previous time step is used as

the initial guess for the next iteration. From a numerical point of view, the calculation of the approximated volume requires

integration in time, which is the most CPU intensive part of the simulation. Fortunately, the optimization algorithm is applied

for each side of each branch with minimal dependencies between the variables, making this part of the algorithm a perfect

candidate for parallel computing.To illustrate the mechanism behind this approach, we use the same analysis as before, on a single fracture branch of a

DFN with low conductivity (12.5 mD.ft).

Figure 13: Pressure and initial time during 3 years of production in the case of low conductivities

Fig.13 shows the computed pressure inside the fracture and the initial time t0,k updated with the proposed method. The

increase of t0,kwith time is necessary to sustain the flow rate from the matrix and the cumulated production (see Fig. 14).

Figure 14: Cumulated production from branch 1 versus time in case of low conductivities

We can see in this figure that the method significantly reduces the error on the mass balance, because the cumulated

production converges close to the maximum recoverable volume. This example illustrates the ability of the method to extend

the validity of the analytical solution. Fig. 15 shows the distribution of pressure and the initial time as calculated by the

0

200

400

00

00

1000

1200

1400

0

0.2

0.4

0.

0.

1

0 500 1000 1500

( 0)

0

4000

000

12000

1000

0 500 1000 1500

1

Error on volumeMaximum recoverable volume

8/10/2019 SPE 152541_12_

10/20

10 SPE 152541

algorithm on the entire fracture network at different time steps.

Figure 15: Pressure and initial time (or delay) in the reservoir at different times of production. The pressure column shows thepressure inside the reservoir blocks (colored marks) and the pressure inside the fracture network (blue lines). The initial time

column shows the initial time for each block calculated by the algorithm.

Fig. 15shows that the greater the distance from the perforations (center of the grid), the smaller the initial time. The greater

the distance is from the perforations, the longer it takes for the BHP to diffuse up to that location.

Validation of the method.To illustrate the capabilities of the production model referred to as UPM (Unconventional Production Model) in the

following, simulations have been compared with those from a commercial reservoir simulator. The permeability of the

reservoir is 0.001 mD with a porosity of 8%, the initial reservoir pressure is 4000 psi and the bottom-hole pressure is 1000

psi. The DFN is made up of 13 identical fractures in each orthogonal direction with a vertical well in the middle. In thisexample the volume factorZ and the gas viscosity were constant and equal to 1 and 0.02 cP, respectively.

10 days

40 days

120 days

Pressure Initial timeTime

8/10/2019 SPE 152541_12_

11/20

SPE 152541 11

Figure 16: Reservoir and DFN used in a comparison between simulations done with a commercial reservoir simulator and thepresented method.

The conductivities of the DFN vary from 0.025 mD.m to 2500 mD.m. Fig. 18-19compares the two simulators.

Figure 18: Comparison of the production rate between commercial reservoir simulator and the UPM.

Eclipse UPM

1.+02

1.+03

1.+04

1.+05

1.+0

1.+0

1 10 100 1000 10000

0,025 .

0,25 .

2,5 .

25 .

250 .

2500 .

0,025 .

0,25 .

2,5 .

25 .

250 .

2500 .

8/10/2019 SPE 152541_12_

12/20

12 SPE 152541

Figure 17: Comparison of the cumulated production between Eclipse and the UPM.

These comparisons show reasonably good agreement between the two simulators, in particular in the case of low

conductivity where the algorithm to update the initial time plays a major role. To illustrate the importance of the algorithm

correcting the initial time, Fig. 18compares the simulation results for the case of fracture conductivity equal to 25 mD.m.

Figure 18: Comparison of the rate (left) and cumulated production (right) between commercial reservoir simulator, the UPM and theUPM without delay.

In Fig. 18, UPM without delay means that the UPM simulator was using only the analytical part of the model, with a

constant initial time equal to 0. When the fracture conductivity is increased, the difference between the Eclipse simulation

and the UPM simulation without delay is reduced.

Similar validations as the ones presented in Fig. 16Error! Reference source not found.and Fig. 17have been done with

lower reservoir pressure and lower drawdown, as well as pressure dependant volume factor zand viscosity. In all cases the

results had the same or better accuracy than in Fig. 16Error! Reference source not found.and Fig. 17.

Sensitivity analysisThis part illustrates how the workflow presented in this paper can be used to investigate the relationship between the

production and the reservoir properties or the treatment design. We present the application of the production model to a

realistic case to do a sensitivity analysis on the stress anisotropy of the reservoir and the proppant size. The pumping schedule

and the reservoir properties used in this case are described in Table 1and Table 2.

1.+02

2.+0

4.+0

.+0

.+0

1.+0

1.+0

1 10 100 1000 10000

0,025 .

0,25 .

2,5 .

25 .

250 .

2500 .

0,025 .

0,25 .

2,5 .

25 .

250 .

2500 .

0

100000

200000

300000

400000

500000

00000

00000

0 500 1000 1500

( )

0.0+00

2.0+0

4.0+0

.0+0

.0+0

1.0+0

1.2+0

0 500 1000 1500

( )

8/10/2019 SPE 152541_12_

13/20

SPE 152541 13

Table 1: Pumping schedule

Stage Pumping rate (bbl/min) Volume (gal) Proppant concentration (ppa)

Slick water 80 100,000 0

Slick water + proppant (100

mesh (diameter = 6.4e-3 in))80 300,000 1

Slick water + proppant

(40/70 (diameter = 1e-2 in)) 80 600,000 1

Table 2: Reservoir properties

Gross height 400 ft

Reservoir pressure 3393 psi

porosity 5%

permeability 0.001 md

Youngs modulus 4 Mpsi

Bottomhole pressure 1000 psi

Minimum horizontal stress 4538 psi

Production and Hydraulic Fracture Surface AreaIn the first set of simulations, we adjust the maximum far-field horizontal stress so that UFM simulations produce hydraulic

fracture networks with different complexities. When a hydraulic fracture is stopped by a natural fracture, it requires more

pressure in the hydraulic fracture tip to open the natural fracture and propagate inside it. Consequently, it becomes more

difficult to create a complex hydraulic fracture network. In Fig. 19, the stress anisotropy varies between 1% and 6%, and the

same treatment produces very different networks.

Figure 19: Simulations of the same hydraulic fracture treatment with different horizontal stress anisotropy (minimum stressconstant). Colors illustrate the proppant distribution in lbm/ft

2.

As expected, low-stress anisotropy produces a more complex network. The UPM is applied to these four networks and gives

the production forecast shown in Fig. 20.

1% 3%

5% %

8/10/2019 SPE 152541_12_

14/20

14 SPE 152541

Figure 20: Flow rate and cumulated production at surface condition for reservoirs with different stress anisotropies (Figure 19).

Fig. 20shows that the productivity of the well decreases as the maximum horizontal stress increases. The reason is that in the

four cases of our example, the propped and total surface area of the hydraulic fracture also decreases, despite the fact that the

averaged propped conductivity increases with increased anisotropy. This means that in this case, the surface area of the

hydraulic fracture has more impact on production than the actual fracture conductivity. This observation is probably true only

if the propped fracture conductivities are high enough (at least several mD.ft). In these simulations, the unpropped fracture

conductivity is 0.001 mD.ft. Further discussions on the relation between productivity and fracture network conductivities can

be found in Cipolla et al (2009).

Table 3: Results from simulations in Fig. 19.

Stress

anisotropy

Total hydraulic fracture

surface area (106ft2)

Total propped hydraulic

fracture surface area (106ft2)

Average propped

conductivity (mD.ft)

1% 11.45 1.34 39.75

3% 8.39 1.33 58.02

5% 5.42 0.47 98.95

6% 5.08 0.61 120.63

Production and proppant sizeIn the following sensitivity analysis, we modify the proppant diameter to investigate how it affects reservoir production. The

reservoir is the same as in the example above with a stress anisotropy of 1%. In this example, the total injected mass of fluid

and proppant are the same as in the previous example. The treatment consists of a pad of 100000 gallons, followed by

900,000 gallons of slurry with a proppant concentration of 1 ppa. The proppant has a specific gravity of 2,694, but its

diameter is increased from 4.10-3

in to 12.10-3

in.

0.+00

1.+0

2.+0

3.+0

4.+0

0 1000 2000 3000 4000

1%

3%

5%

%

0.+00

1.+0

2.+0

3.+0

0 1000 2000 3000 4000

1%

3%

5%

%

8/10/2019 SPE 152541_12_

15/20

SPE 152541 15

Figure 21: Simulations of the same hydraulic fracture treatment in the same reservoir but with different proppant diameters..

Fig. 21shows that, as one might expect, the proppant is spread much further into the network with the smaller size proppant

than the larger ones. This leads to a larger propped hydraulic fracture surface area for the smaller proppant. Fig. 21shows the

production forecast from the hydraulic fracture network in Fig. 20.

Figure 22: Flow rate and cumulated production at surface condition for the injection of different proppant diameter (Fig. 21).

4.103 .103

.103 10.103

12.103

0.0+00

5.0+0

1.0+0

1.5+0

2.0+0

2.5+0

3.0+0

3.5+0

0 1000 2000 3000 4000

43

3

3

103

123

0.0+00

5.0+05

1.0+0

1.5+0

2.0+0

2.5+0

0 1000 2000 3000 4000

43

3

3

103

123

8/10/2019 SPE 152541_12_

16/20

16 SPE 152541

Fig. 22shows that the production for this reservoir increases with the proppant size. One reason could be that large proppant

produces higher fracture conductivity, which compensates for the smaller propped surface area. This can be observed in

Table 4which compares the propped surface area with the corresponding average conductivity for each simulation.

Table 4: Results from simulations in Fig. 20.

(103

)

(10

2)

(.)

4 4.42 11.5

3.4 14.0

2.5 1.51

10 1.3 2.0

12 1.43 31.35

An interesting observation from the comparison of Fig. 19and Fig. 20, is that the treatment that combines 100 mesh and

40/70 proppant in Fig. 19(stress anisotropy of 1%) has a better cumulated production than any of the similar treatments with

a single type of proppant.

Discussion on the advantages and limitations of the UPMThese two examples of sensitivity analysis demonstrate how a simple tool like the UPM can help the stimulation engineer in

optimizing the fracture treatment design. It can help define an appropriate mix of proppant type, volume of fluid, or injectionrate. Nevertheless, the UPM is a simplified model that relies on several assumptions:

- The flow from the matrix to the fracture is based on a predefined block volume. Each one of these blocks is assumedto be depleted by a unique associated fracture branch. Therefore the model simplifies the matrix flow. In particular it

ignores the flow exchange between the blocks. The consequence is that once a block becomes partially depleted, the

hydraulic fracture can act as a flow barrier. Therefore, the model forecast is generally more accurate for early

production stages.

- For now, the model only considers single phase flow. In reality, the well should produce a significant part of thestimulation fluid in early times (two-phase flow), and some reservoirs produce a mix of gas and oil (three-phase

flow).

- The model does not yet consider desorption as it occurs in shale reservoir.- The fracture network is 2D (invariant in one of the three spatial dimensions)

SummaryThe paper presents a simple workflow that combines a hydraulic fracturing simulator for unconventional reservoirs with a

new semi-analytical production model. This new tool offers an opportunity to obtain a quick estimate of the performance of a

fracturing job design. The production model (UPM) is explained and a set of validation cases against a reservoir simulator

Eclipse is presented. The paper also presents two examples of sensitivity analysis using this simulation workflow, to

investigate the influence of fracture network geometry and proppant size on production.

AcknowledgementThe authors would like to thank Schlumberger for permission to publish this paper.

Nomenclaturect= compressibility (Pa

-1)

C=conductivity (m2.m)

Fo,k= objective function (m

3

)H= fracture height (m)

km= matrix permeability (m2)

L = maximum length of drainage (m)

m= real gas pseudo pressure (Pa/s)

M = molar mass (kg/mol)

Pm= matrix pressure (m2)

Pm0= initial matrix pressure (m2)

Pf= initial fracture pressure (m2)

Qtot= total flow rate from a matrix block into a fracture (m3/s)

Qmf= local flow rate from a matrix block into a fracture (m3/s)

Qf= flow rate inside the matrix (m3/s)

R= universal gas constant (J/mol/K)

8/10/2019 SPE 152541_12_

17/20

SPE 152541 17

t= time (s) to,k= initial production time (s)

T = Temperature (K)

xf= coordinate in the fracture (m)

xm= coordinate in the matrix (m)

Z= volume factor

= viscosity (Pa.s)

m= porosity

SI Metric Conversion Factorsft x 3.048 e-01 = mpsi x 6.894 757 e+00 = kPa

Mscf x 2.8317 e-03 = m3

ReferencesAl-Hussainy, R., Ramey, H. J., Crawford, P. B. The Flow of Real Gases Through Porous Media, Journal of Petroleum Technology, page

624-636, May 1966.

Barenblatt, G.I., Zheltov, I.P., and Kochina, I.N., "Basic Concepts in the Theory of Homogeneous Liquids in Fissured Rocks", J. AppliedMath. Mech. (U. S. S. R.), (1960).

Basquet, R., Jeannin, L., Lange, A., Bourbiaux, B., (2003). "Gas Flow Simulation in Discrete Fracture Network Models". Paper SPE 79708presented at the SPE Reservoir Simulation Symposium, 3-5 February 2003, Houston, Texas

Bello, R.O., Rate Transient Analysis in Shale Gas Reservoirs with Transient Linear Behavior, PhD Thesis, 2009.Cinco-Ley, H., Samaniego, F.. Pressure Transient Analysis for Naturally Fractured Reservoirs, SPE paper 11026 presented at the Annual

Fall Technical Conference and Exhibition held in New Orleans, LA, Sept 26, 1982.Cipolla, C. L., Lolon, E. P., Mayerhofer, M. J., The Effect of Proppant Distribution and Un-Propped Fracture conductivity on Well

Performance in Unconventional Gas Reservoirs. SPE paper 119368 presented at the SPE Hydraulic Fracturing Conference held inThe Woodlands, Texas, USA, 19-21, January 2009.

Cipolla, C. L., Lolon, E. P., Mayerhofer, M. J., Reservoir Modeling and Production Evaluation in shale-Gas Reservoirs, SPE paper13185 presented at the International Petroleum Technology Conference held in Doha, Qatar, December 7th 2009.

Cipolla, C.L., Fitzpatrick, T., Williams, M.J., and Ganguly U.K. Seismic-to-Simulation for Unconventional Reservoir Development, SPEpaper 146876 presented at the SPE Reservoir Characterization and Simulation Conference and Exhibition held in Abu Dhabi, UAE, 9-11 October 2011.

Cipolla, C.L., Williams, M.J., Weng, X., Mack M., and Maxwell, S. Hydraulic Fracture Monitoring to Reservoir Simulation: Maximizing

Value. SPE paper 133877 presented at the Annual Technical conference and Exhibition held in Florence, Italy, 19-22 Septembre2010.

Gong, B., Guan, Q., Douglas, C., Yuan, S., Detailed Modeling of the Complex Fracture Network of Shale Gas Reservoirs. SPE paper142705 presented at the SPE Midle East Unconventional Gas Conference and Exhibition held in Muscat, Oman, 31 January 2011.

Kresse, O., Cohen, C., Weng, X., Wu, R., Gu, H. Numerical Modeling of Hydraulic Fracturing in Naturally Fractured Formations. Paperpresented at the 45th US Rock Mechanics / Geomechanics Symposium held in San Francisco, CA, June 26-29, 2011.

Lin, J., Zhu, D. Modeling Well Performance for Fractured Horizontal Gas Wells. SPE paper 130794 presented at the International Oiland Gas Conference and Exhibition in China, 8-10 June 2010, Beijing, China.

Li, J., Zhan, H., Huang, G. Application of the Linearized Governing Equation of Gas Flow in Porous Media. Transport in PorousMedia, vol. 87, no3, pp. 815-834, 2011.

Jeannot, Yves. Thansfert Thermique, Textbook, Ecole des Mines de Nancy, 2009. http://www.thermique55.com/principal/thermique.pdfMeyer, B.R. and Bazan, L.W. 2011. A Discrete Fracture Network Model for Hydraulically-Induced Fractures: Theory, Parametric and

Case Studies. Paper SPE 140514 presented at the SPE Hydraulic Fracturing Conference and Exhibition, Woodlands, Texas, January

24-26.Xu, W., Thiercelin, M., Ganguly, U., Weng, X., Gu, H., Onda, H., Sun, J., and Le Calvez, J. 2010. Wiremesh: A Novel Shale Fracture

Simulator. Paper SPE 132218 presented at the CPS/SPE International Oil & Gas Conference and Exhibition in China held in Beijing,China, 810 June

Xu, W., Li, J., Du, M.C.. "Quick Estimate of Initial Production from Stimulated Reservoirs with Complex Hydraulic Fracture Network".Paper SPE 146753 presented at the SPE Annual Technical Conference and Exhibition held in Denver, Colorado, USA, 30 October2011.

Weng, X., Kresse, O., Cohen, C., Wu, R., Gu, H. Modeling of Hydraulic Fracture Network Propagation in Naturally Fractured

Formation. SPE paper 140253 presented at the SPE Hydraulic Fracturing Conference and Exhibition, Woodlands, Texas, 2011.Wu, Y.S., Pruess, K. and Persoff, P. Gas flow in Porous Media with Klinkenberg Effects. Transport in Porous Media, 32: 117-137, 1998.Zhang, X., Du, C., Deimbacher, F., Crick, M. and Harikesavanallur, A. Sensitivity Studies of Horizontal Wells with Hydraulic Fractures

in Shale Gas Reservoirs. SPE paper presented at the International Petroleum Technology Conference held in Doha, Qatar, 7-9December 2009.

Appendix 1 Equation, implementation and algorithm for the presented method.

8/10/2019 SPE 152541_12_

18/20

18 SPE 152541

Pressure profile in the matrix: solution for constant fracture pressure

The first assumption of our model is that the gas behavior can be described by the following real gas model:

= inside the matrix ( 12) = inside the fracture network ( 13)The basis equation for the linear gas flow inside the matrix block is

= ( 14)With= 1 1 ( 15)The following definition of a pseudo pressure simplifies the resolution of previous equation

= 2

( 16)

Equation ( 14)then becomes

= 1 ( 17)With = ( 18) = _= _ ( 19)Boundary conditions:

= , = = , = _ ( 20) 0, = 0 ( 21) , = 0 ( 22)Equations ( 17)to ( 22)are solved using the Laplace transform. Details of the calculation can be found in Jeannot (2009).

The solution for the pseudo pressure inside the matrix is

_ = _ 1 2 + 1 2

+ 2 + 1 + 2

( 23)

Flow rate from the matrix to fracture with constant fracture pressure

The flow rate from the matrix to the fracture is given by Darcys law

, = , ( 24)Lkcorrespond to the maximum length of drainage on the side kof the fracture.

8/10/2019 SPE 152541_12_

19/20

SPE 152541 19

= _ 1 ( 25)

Which gives

, = _ 2 1 , 1

, ,

( 26)

The diffusion coefficient ais constant in this equation, thus it has to be calculated with averaged values of the viscosity and

compressibility which are pressure dependent.

= ( 27)Comparison with reservoir simulation showed that the following expression give a good approximation, even with a large

pressure drawdown of the same order than the reservoir pressure, as shown in the validation section.

= _ + 2 ( 28) = _ + 2 ( 29)A discussion on the average value of the compressibility to consider in this analytical solution can be found in Li et al. (2011)

and Wu et al. (1998).

Flow inside the fracture branch between intersections iandjThe mass balance inside the fracture network is described by the following equation

, = , ( 30)With Qmfthe flow rate (m

3/s) from the matrix to the fracture and Qfthe flux (m

2/s) from the fracture. Under the assumption of

a real gas this equation becomes:

2 = , ( 31)With the following boundary conditions:

= 0 _ = = 0 _ = , _ = , = _ = = _ = , _ = , ( 32)WithLfthe length of the fracture between two nodes. If we write

= _ ( 33)And introduce equation (26)into equation (31), we obtain

= 0 ( 34)

8/10/2019 SPE 152541_12_

20/20

20 SPE 152541

with

= 1 , 1

,

,

( 35)

The solution to equation ( 34)is

= _ + , , + , , ( 36)The flow rate at intersection ifrom branch i,jis

,, = ,2 ( 37)

When introducing equation ( 8)in equation ( 37)

,, = ,2 ,+ 2 , ( 38)

Mass balance at intersection between fractures

,,_ = 0 ( 39)

WithN_ijthe number of branches reaching this intersection.

,,_ = 2 , ,, +, 2 , ,, , _ = 0

( 40)

This can be rearranged as follow

, , +, , ,_ ,

, ,_

= , , +, 2

,

,

_

( 41)