doi.org/10.26434/chemrxiv.8206214.

v1

Production Prediction of Hydraulically Fractured Reservoirs Based on

Material Balances
Raul Velasco, Palash Panja, Milind Deo

Submitted date: 30/05/2019 • Posted date: 31/05/2019

Licence: CC BY-NC-ND 4.0
Citation information: Velasco, Raul; Panja, Palash; Deo, Milind (2019): Production Prediction of Hydraulically
Fractured Reservoirs Based on Material Balances. ChemRxiv. Preprint.

Transient flow is dominant during most of the productive life of unconventional wells in ultra-low permeabilities
resources such as shales. As a result, traditional reservoir performance analysis such as the conventional
material balance have been rendered inapplicable. A new novel semi-analytical production predictive tool
based on application of material balance on a transient linear flow system is developed.

File list (1)

Production Prediction.pdf (1.37 MiB) view on ChemRxiv download file

 Production Prediction of Hydraulically Fractured Reservoirs Based on Material Balances
Raul Velascoa, Palash Panja a,b, *, Milind Deo b
a
Energy & Geoscience Institute, b Department of Chemical Engineering,
University of Utah, Salt Lake City, Utah, the United States

* Corresponding author
ppanja@egi.utah.edu

Abstract
Transient flow is dominant during most of the productive life of unconventional wells in ultra-low
permeabilities resources such as shales. As a result, traditional reservoir performance analysis such
as the conventional material balance have been rendered inapplicable. A new novel semi-analytical
production predictive tool based on application of material balance on a transient linear flow
system is developed. This method considers two important regions during transient production of
oil reservoirs: the saturated region where gas evolves and flows with oil, and the undersaturated
region where only oil flows. These regions evolve in extent as pressures diffuse from the fractures
into the reservoir while production takes place. First, the proposed semi-analytical method and
concept are introduced and validated using numerical simulations. After verification is done, the
method is applied to three field cases including black and volatile oil examples where production
prediction is highlighted. We observed that application of this method on tight reservoirs around
the U.S. shows very good forecasting capabilities. Input requirements for this multiphase method
are PVT data, relative permeability ratios, and average pressure estimations. Unlike conventional
predictive tools, this method is based on material balance application to several evolving flow
zones during production and can accommodate for variable bottom-hole pressures. We present a
method simple enough for quick application on a spreadsheet, but comprehensive enough to
capture important multiphase reservoir physics not present in most traditional production analysis.
There are several applications of this tool, ranging from production performance evaluation to
prediction of hydrocarbon production.

Introduction
One of the main objectives regarding reservoir engineering is to analyze production performance
such as fluid rates, cumulative production, pressure declines, and to determine hydrocarbon
estimated ultimate recovery (EUR). Numerous methods have been developed and modified since
the early studies of reservoir behavior (Millikan,1926, Coleman et al.,1930). These studies can be
divided into few broad categories such as decline curve and type curve analysis, semi-analytical
methods, surrogate models, material balance methods and numerical simulation.
The decline curve analysis method has been used widely due to its quick predictive capabilities
based on available production data. DCA has a long history after Arps (Arps,1945) developed the
first mathematical model where a curve is fitted on a plot of flow rate vs time. The fitted curve
then is used to forecast fluid rates as well as cumulative production and eventually used to
determine EUR. Many new and modified versions of Arp’s method (Mead,1956, Mannon,1965,
Fetkovich,1980) were developed consequently. Most decline curve analyses are suitable for single
phase flow in conventional reservoir with high permeabilities. However, Arp’s original model has
 2

been unsuccessful for analysis of production from unconventional and tight permeability reservoirs
(Duong,2010). Arp's empirical formula is constructed on the assumption that the reservoir has a
boundary dominated flow. Unconventional reservoirs have predominantly linear flow with long
transient flow regimes due to the very low matrix permeability and presence of induced fractures.
These resulted in value of decline exponent (b) greater than one which overestimated production.
As alternatives, a series of techniques that embody the advantages of the DCA have been
developed in recent years to address unconventional reservoirs. A few examples include the Power
Law Exponential (PLE) (Ilk et al.,2008), Stretched-Exponential model(SE) (Valko,2009), Logistic
Growth Analysis (LGA) (Clark et al.,2011)and Duong Method (Duong,2010). Rate transient
analysis (Palacio and Blasingame,1993) using type curves was successfully applied in analysis of
gas production data from tight formations. However, multiphase flow of oil and gas cannot be
captured using this type curves. Concept of transient productivity index (PI) (Medeiros et al.,2010)
was first introduced to analysis production data from fractured and unfractured horizontal wells.
The better estimations are obtained with longer data which includes multiple flow regime. One of
the disadvantages of these techniques, however, is the sensitivity to variation in production data
due to changes in flowing bottom-hole pressures. Calculations of EUR also vary with the selection
of fitted curves, a problem that only be alleviated by fitting more production data. Various methods
such as Arp's, PLE, SE, RTA and transient PI method were compared in estimating oil-in-place in
fractured tight formation (Kabir et al.,2011). Dynamic material balance was developed in
estimating original-gas-in-place when production rate is variable (Mattar and Anderson,2005).
This method is only applicable during boundary-dominated flow in gas reservoir (Ismadi et
al.,2011). It fails to predict original-hydrocarbon-in-place for multi-phase (oil and gas) flow at
transient state conditions. Therefore, the dynamic material balance method is not an efficient tool
for unconventional reservoirs.
Analytical methods based on continuum fluid mechanics have been developed for several
geometries and production conditions. Researchers (Levine and Prats,1961, Miller,1962,
Fetkovich,1980, El-Banbi and Wattenbarger,1998, Wattenbarger et al.,1998, Clarkson and
Beierle,2010) have developed expressions derived from the conservation of mass and Darcy’s law
to tie the physics of fluid flow and surface production. Analysis using these expressions is made
universal by use of dimensionless groups that combine production, completion, operational,
petrophysical, and fluid parameters. The equations provided by governing partial differential
equations give us great insight into the physics of fluid. In spite of the single-phase assumption in
these equations, conventional and unconventional production analysis has been successful as long
as the well has a stable production life. This approach, however, is not applicable for complex and
heterogeneous reservoirs.
State-of-the-art numerical simulations is the to-go approach for highly complex multi-phase
systems. In many cases, simulation has become the ultimate reservoir performance tool. The
drawback, however, is the extensive data required for a simulation model to be meaningful.
Required input data for unconventional wells include seismic data interpretation, petrophysical
properties, fluid data, relative permeabilities from cores, etc. Any uncertainty in the input data can
lead to misinterpretation of the simulation results. As a consequence, reservoir simulations tend to
be expensive and time consuming for extensive projects.
Surrogate based models (Amorim and Schiozer,2012, Dahaghi et al.,2012) such as polynomial
response surface (Parikh,2003, Li and Friedmann,2005, Carreras et al.,2006, Panja and Deo,2016)
and artificial intelligence method such as least square support vector machine (E. El-
Sebakhy,2007, Fatai Adesina Anifowose,2011, Panja et al.,2016) and artificial neuron network are
 3

gaining interest due to high computational efficiency. In surrogate models, important parameters
(Panja et al.,2016) are selected based on their hierarchy to affect the production. The main
disadvantage of this approach is that the models are developed based on simplified reservoir
geometries.
The material balance method is a very useful tool to estimate initially hydrocarbon in place,
forecast reservoir performance and to calculate individual contributions from different drive
mechanisms. The general material balance equations which are based on hydrocarbon volume
balance are initially presented by Schilthuis (Schilthuis,1936). Material balance comprises three
variables namely average reservoir pressure, recovery factor and cumulative GOR. Other variables
such as compressibilities and fluid PVT properties are function of average pressure.
For prediction of production performance, many smart approaches (Tarner,1944,
Muskat,1945, Tracy,1955) were tried in order to solve the material balance equation progressively
by small pressure steps in solution-gas-drive reservoir. These methods are applicable to well
depleted saturated reservoirs.
Another use of material balance equation is to calculate initial hydrocarbon in place by using
field production data. For this approach, the equation is rearranged in the form of straight line. The
straight line form was first proposed by Havlena and Odeh (Havlena and Odeh,1963) by re-
arranging the material balance equations. The linear form is convenient for various calculations
such as initial oil in place, gas cap size, water influx and other drive mechanisms.
The material balance method is advantageous for requiring only production data and fluid PVT
properties. The existence of kerogen in shale plays that might alter the PVT properties(Pathak et
al.,2015) are not considered in this study. Fluid properties are determined at average reservoir
pressures i.e., a uniform pressure throughout the entire reservoir. Application of material balance
using average reservoir pressure is suitable for high permeability conventional reservoirs. The
assumption of average reservoir pressure leads to erroneous predictions if the variation of pressure
is significant across reservoir. This is the case for transient flow in ultra-low permeability (10 nD
to 1000 nD) reservoirs such as shales.
Unconventional reservoirs such as shales, tight formation etc. behave differently from
conventional reservoirs (Orangi et al. , Panja et al.,2016, Panja and Deo,2016). Long transient
state prevails in ultra-low permeable reservoirs where the total average pressure does not go below
bubble point for long time. Conventional material balance fails to capture unconventional
production behavior such as produced GOR, liquid flowrates etc. Modification of conventional
material balance is essential to address this problem. In this paper, we approach transient flow by
dividing the reservoir into meaningful regions; a similar concept was introduced before for
condensates (Fevang and Whitson,1996). The difference between conventional material balance
and the proposed method is depicted in Figure 1.
 4

Figure 1 – Pressure distribution for a transient flow system

The detailed description of different regions is provided in the next section.

A simple spreadsheet for rapid and simple calculations was developed for ultra-low
permeability, hydraulically fractured reservoirs during transient and boundary dominated linear
flow. The results from the model are compared with fine grid reservoir simulations using IMEX
(Computer Modeling Group, Calgary, Canada).
In the following sections, we introduce the method’s core concept, provide simulation
verification, and present application to three field cases. Since this method is based on material
balances, the same assumptions apply.

Methodology
We begin by defining the system and determining the important components of multi-phase flow
behavior in porous media. The system may be of any regular geometry, but for the purpose of
application to shale plays we define a horizontal well with a number of vertical and equally spaced
fractures draining from a constant thickness pay zone as shown by Figure 2.
For practical purposes, most analysis techniques assume that the fractures along the horizontal
well behave identically. By taking advantage of symmetry, the entire length of the horizontal well
is condensed into one representative fracture draining from one of its sides as shown by the red
volume in Figure 2.
 5

Figure 2 – Idealized multi-stage horizontal well model

Reservoir Segmentation
Depending on the fluids originally present in the reservoir, initial reservoir pressure and flowing
bottom-hole pressures, a number of different multi-phase configurations evolve as fluids drain into
the fracture. As a working example, we consider an initially undersaturated reservoir producing at
a bottom-hole pressure lower than the bubble point pressure. As a consequence, transient flow
develops two important regions: the undersaturated region and the saturated region.
Both regions are a consequence of the pressure behavior inside the system as described in
Figure 3 (a). It should be noted, however, that even though gas evolves in the saturated portion of
the system, not all of the gas flows. Gas critical saturations determine the point at which gas starts
flowing with the oil as described in Figure 3 b).

(a)
 6

(b)
Figure 3 – System segmentation for an initially undersaturated reservoir during transient
flow: (a) Typical pressure profile describing evolution of gas near the fracture (b) Typical
gas saturation profile describing immobile and mobile gas portions inside the saturated
region
Figures 3(a) and (b) show a key parameter: the bubble point position,𝑥𝑥𝑏𝑏 , defined as the distance
from the fracture at which the pressure is at the bubble point. During transient flow, the bubble
point position (xb) travels away from the fracture, effectively increasing the size of the saturated
region, and decreasing the undersaturated portion of the reservoir.

Undersaturated Region (𝑷𝑷𝒃𝒃 ≤ 𝑷𝑷 � 𝑼𝑼 ≤ 𝑷𝑷𝒊𝒊 )

This region is characterized by the flow of oil and water without gas saturation. The inner boundary
of this zone is located at bubble point position, 𝑥𝑥𝑏𝑏 , and outer boundary is equal to half fracture
spacing, 𝐿𝐿𝑓𝑓 . By defining a pressure gradient for a slightly compressible fluid, the average
undersaturated pressure, 𝑃𝑃�𝑢𝑢 , is calculated as shown below.
𝐿𝐿𝑓𝑓
1
���
𝑃𝑃𝑈𝑈 = � 𝑝𝑝 𝑑𝑑𝑑𝑑 (1)
𝐿𝐿𝑓𝑓 − 𝑥𝑥𝑏𝑏
𝑋𝑋𝑏𝑏

where 0 ≤ 𝑥𝑥𝑏𝑏 ≤ 𝐿𝐿𝑓𝑓 and 𝑥𝑥𝑏𝑏 is a function of time.

As the bubble point position progresses, the undersaturated region shrinks while providing the
saturated region a larger hydrocarbon production potential. We define the hydrocarbon
contribution from the undersaturated region to the saturated portion as the bubble point position
travels away from the fracture. To calculate this contribution, we must first determine the “excess”
volume of fluid in the undersaturated region as oil expands and the region’s volume shrinks. The
excess volume in the undersaturated region due to oil expansion and region shrinkage is evaluated
as:
 7

2𝑥𝑥𝑓𝑓 ℎ𝐵𝐵𝑜𝑜 (𝑃𝑃�𝑆𝑆 ) 𝑆𝑆𝑜𝑜 (𝑃𝑃�𝑈𝑈′ )𝜙𝜙(𝑃𝑃�𝑈𝑈′ ) ′

𝑆𝑆𝑜𝑜 (𝑃𝑃�𝑈𝑈 )𝜙𝜙(𝑃𝑃�𝑈𝑈 )
𝐸𝐸𝑢𝑢 = � �𝐿𝐿𝑓𝑓 − 𝑥𝑥𝑏𝑏 � − �𝐿𝐿𝑓𝑓 − 𝑥𝑥𝑏𝑏 �� (2)
5.615 𝐵𝐵𝑜𝑜 (𝑃𝑃�𝑈𝑈′ ) 𝐵𝐵𝑜𝑜 (𝑃𝑃�𝑈𝑈 )

where the apostrophe denotes values evaluated at the previous step and porosity (if not
constant) can be determined with: 𝜙𝜙(𝑝𝑝) = 𝜙𝜙𝑖𝑖 �1 + 𝑐𝑐𝑓𝑓 �𝑝𝑝 − 𝑝𝑝𝑟𝑟𝑟𝑟𝑟𝑟 �� for a linear elasticity
compaction model or a pore volume lookup tables for a non-linear elasticity model. The detailed
derivation of equation 2 is provided in Appendix B.
Due to material balance, the undersaturated excess volume, Eu, can no longer be contained in
the undersaturated region and accounts for the contribution from the undersaturated region to the
saturated region. The oil and gas contribution from the undersaturated to the saturated region are
determined as:

𝐶𝐶𝑢𝑢→𝑠𝑠,𝑜𝑜 = 𝐸𝐸𝑢𝑢 (3.a)

𝐵𝐵𝑔𝑔 (𝑃𝑃�𝑆𝑆 ) (3.b)

𝐶𝐶𝑢𝑢→𝑠𝑠,𝑔𝑔 = 𝐸𝐸𝑢𝑢 [𝑅𝑅𝑠𝑠𝑠𝑠 − 𝑅𝑅𝑠𝑠 (𝑃𝑃�𝑆𝑆 )]
𝐵𝐵𝑜𝑜 (𝑃𝑃�𝑆𝑆 )

Once 𝑥𝑥𝑏𝑏 = 𝐿𝐿𝑓𝑓 , the undersaturated region disappears terminating the oil contribution to the
saturated region as shown by equation (2).

Saturated Region (𝑷𝑷𝒘𝒘𝒘𝒘 ≤ 𝑷𝑷 � 𝑺𝑺 ≤ 𝑷𝑷𝒃𝒃 )

In this region, gas evolution occurs as pressures develop below the bubble-point. The pressure at
the fracture boundary is the flowing bottom-hole pressure and the pressure at the undersaturated-
saturated interface (𝑥𝑥𝑏𝑏 ) is the bubble point pressure. The average pressure in this region is
calculated as:
𝑥𝑥𝑏𝑏
1
𝑃𝑃�𝑆𝑆 = � 𝑝𝑝 𝑑𝑑𝑑𝑑 (4)
𝑥𝑥𝑏𝑏
0

In this case, as the bubble point position progresses, the saturated region grows by tapping into
the undersaturated oil potential. The oil contribution from the undersaturated region has been
determined earlier. As the average pressure in this region changes, so do the oil and gas volumes
at reservoir conditions. The oil and gas volumes at reservoir conditions are calculated as follows:

𝐵𝐵𝑜𝑜 (𝑃𝑃�𝑆𝑆 )
𝑉𝑉𝑠𝑠,𝑜𝑜 = 𝑉𝑉𝑠𝑠′,𝑜𝑜 (5.a)
𝐵𝐵𝑜𝑜 (𝑃𝑃�𝑆𝑆′ )

𝐵𝐵𝑔𝑔 (𝑃𝑃�𝑆𝑆 ) 𝐵𝐵𝑔𝑔 (𝑃𝑃�𝑆𝑆 ) (5.b)

𝑉𝑉𝑠𝑠,𝑔𝑔 = 𝑉𝑉𝑠𝑠′,𝑔𝑔+ 𝑉𝑉
𝑠𝑠
′
,𝑜𝑜 [𝑅𝑅𝑠𝑠
�𝑆𝑆′ ) − 𝑅𝑅𝑠𝑠 (𝑃𝑃�𝑆𝑆 )]
(𝑃𝑃
𝐵𝐵𝑔𝑔 (𝑃𝑃�𝑆𝑆′ ) 𝐵𝐵𝑜𝑜 (𝑃𝑃�𝑆𝑆′ )
Once new saturated volumes are determined, the saturated region total volumes are calculated
at the new step as:
 8

𝑉𝑉𝑠𝑠,𝑜𝑜 = 𝑉𝑉𝑠𝑠,𝑜𝑜 + 𝐶𝐶𝑢𝑢−𝑠𝑠,𝑜𝑜 (6.a)

𝑉𝑉𝑠𝑠,𝑔𝑔 = 𝑉𝑉𝑠𝑠,𝑔𝑔 + 𝐶𝐶𝑢𝑢−𝑠𝑠,𝑔𝑔 (6.b)

Similarly to the undersaturated excess volume, the saturated excess volume is the difference
between the saturated volume expansion and the size of the region.

2𝑥𝑥𝑓𝑓 ℎ𝑥𝑥𝑏𝑏 𝜙𝜙(𝑃𝑃�𝑆𝑆 )

𝐸𝐸𝑆𝑆 = 𝑉𝑉𝑠𝑠,𝑜𝑜 + 𝑉𝑉𝑠𝑠,𝑔𝑔 − (7)
5.615

Similar to the undersaturated region procedure, the excess fluid in the saturated region is
contribution from the region to the surface, i.e. production. Because the excess fluid is calculated
as a combination of oil and gas barrels at reservoir conditions, we need a function that determines
the portion of the excess fluid that flows as oil and the portion that flows as gas. Fractional flow
has been used widely for CO2 and waterflooding analyses (Willhite,1986). We introduce the
fractional flow concept for this case as:

1
𝑓𝑓𝑔𝑔 =
𝑘𝑘𝑟𝑟𝑟𝑟 (𝑠𝑠𝑔𝑔 )𝜇𝜇𝑔𝑔 (𝑃𝑃�𝑆𝑆 ) (8.a)
1+
𝑘𝑘𝑟𝑟𝑟𝑟 (𝑠𝑠𝑔𝑔 )𝜇𝜇𝑜𝑜 (𝑃𝑃�𝑆𝑆 )

𝑓𝑓𝑜𝑜 = 1 − 𝑓𝑓𝑔𝑔 (8.b)

𝑉𝑉𝑠𝑠,𝑔𝑔
where gas saturation is simply calculated as: 𝑠𝑠𝑔𝑔 = 𝑉𝑉
𝑠𝑠,𝑜𝑜 +𝑉𝑉𝑠𝑠,𝑔𝑔
Note that besides oil and gas, water flow can be easily implemented by addition of 𝑓𝑓𝑤𝑤 .
Surface production is readily calculated as shown below:

𝐸𝐸𝑆𝑆
𝑄𝑄𝑜𝑜 = 𝐶𝐶𝑠𝑠→𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠,𝑜𝑜 = 𝑓𝑓𝑜𝑜 (9.a)
𝐵𝐵𝑜𝑜 (𝑃𝑃�𝑆𝑆 )

𝐸𝐸𝑆𝑆 (9.b)
𝑄𝑄𝑔𝑔 = 𝐶𝐶𝑠𝑠→𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠,𝑔𝑔 = 𝑓𝑓𝑔𝑔 + 𝑄𝑄𝑜𝑜 𝑅𝑅𝑠𝑠 (𝑃𝑃�𝑆𝑆 )
𝐵𝐵𝑔𝑔 (𝑃𝑃�𝑆𝑆 )

Equations (1) to (9) close the system of equations and solve for cumulative oil and gas
production by advancing 𝑥𝑥𝑏𝑏 → 𝐿𝐿𝑓𝑓 . When the reservoir becomes fully saturated, production
continues by declining the average pressure in the saturated region. The total reservoir average
pressure can be expressed by combination of equations (1) and (4):
𝑋𝑋𝑏𝑏 Lf
1 𝑥𝑥𝑏𝑏 𝑃𝑃�𝑆𝑆 + (𝐿𝐿𝑓𝑓 − 𝑥𝑥𝑏𝑏 ) ���
𝑃𝑃𝑈𝑈
���
𝑃𝑃𝑅𝑅 = �� 𝑝𝑝 𝑑𝑑𝑑𝑑 + � 𝑝𝑝 𝑑𝑑𝑑𝑑 � = (10)
𝐿𝐿𝑓𝑓 𝐿𝐿𝑓𝑓
0 Xb

Saturated and undersaturated average pressures depend on the particular case study. Estimates
can be found in the Appendix.
 9

Region Integration
Once all regions and PVT properties have been defined within a system, all contributive volumes
from all regions have to be accounted for up until the producing region. In cases where we
encounter high gas oil ratios at high bubble point pressures, gas buildup is substantial, especially
once the reservoir is fully saturated. In such cases, it is encouraged that the saturated region be
split into a number of sequential saturated portions. The volume excess and volume contribution
equations presented earlier still apply simply by introducing a new travelling variable similar to
𝑥𝑥𝑏𝑏 . In an ideal situation, we wish to introduce variable 𝑥𝑥𝑠𝑠𝑠𝑠𝑠𝑠 defined as the distance where gas
saturation reaches its critical point. However, 𝑥𝑥𝑠𝑠𝑠𝑠𝑠𝑠 introduces more unknowns to the system and
has been found to be of little effect in our experience. Similarly, Whitson and Fevang (Fevang and
Whitson,1996) concluded that critical condensate saturations had no effect on condensate gas
deliverability.
For practical purposes, we choose to divide the saturation region into parts that grow equally
as 𝑥𝑥𝑏𝑏 grows. An example of a saturated region partition is shown Figure 4.

Figure 4 – Initially undersaturated reservoir segmentation with two saturated portions

evaluated at a 𝚫𝚫𝒙𝒙𝒃𝒃

If the saturated region is divided into equal portions, each portion’s volume changes and
contributions to other portions are calculated as:

′
𝐵𝐵𝑜𝑜 (𝑃𝑃�𝑖𝑖 )
𝑉𝑉𝑖𝑖,𝑜𝑜 = 𝑉𝑉𝑖𝑖,𝑜𝑜 (11.a)
𝐵𝐵𝑜𝑜 (𝑃𝑃�𝑖𝑖′ )

𝐵𝐵𝑔𝑔 (𝑃𝑃�𝑖𝑖 ) 𝐵𝐵𝑔𝑔 (𝑃𝑃�𝑖𝑖 ) (11.b)

′
𝑉𝑉𝑖𝑖,𝑔𝑔 = 𝑉𝑉𝑖𝑖,𝑔𝑔 + 𝑉𝑉 ′
𝑠𝑠,𝑜𝑜 [𝑅𝑅𝑠𝑠 (𝑃𝑃�𝑖𝑖′ ) − 𝑅𝑅𝑠𝑠 (𝑃𝑃�𝑖𝑖 )]
𝐵𝐵𝑔𝑔 (𝑃𝑃�𝑖𝑖′ ) 𝐵𝐵𝑜𝑜 (𝑃𝑃�𝑖𝑖′ )
 10

𝑥𝑥 (12)
2𝑥𝑥𝑓𝑓 ℎ � 𝑛𝑛𝑏𝑏 � 𝜙𝜙(𝑃𝑃�𝑖𝑖 )
𝐸𝐸𝑖𝑖 = 𝑉𝑉𝑖𝑖,𝑜𝑜 + 𝑉𝑉𝑖𝑖,𝑔𝑔 −
5.615

𝐵𝐵𝑜𝑜 (𝑃𝑃�𝑖𝑖+1 ) (13.a)

𝐶𝐶𝑖𝑖→(𝑖𝑖+1),𝑜𝑜 = 𝑓𝑓𝑖𝑖,𝑜𝑜 𝐸𝐸𝑖𝑖
𝐵𝐵𝑜𝑜 (𝑃𝑃�𝑖𝑖 )

𝐵𝐵𝑔𝑔 (𝑃𝑃�𝑖𝑖+1 ) 𝐵𝐵𝑔𝑔 (𝑃𝑃�𝑖𝑖+1 ) (13.b)

𝐶𝐶𝑖𝑖→(𝑖𝑖+1),𝑔𝑔 = 𝑓𝑓𝑖𝑖,𝑔𝑔 𝐸𝐸𝑖𝑖 + 𝐶𝐶𝑖𝑖→(𝑖𝑖+1),𝑜𝑜 [𝑅𝑅𝑠𝑠 (𝑃𝑃�𝑖𝑖 ) − 𝑅𝑅𝑠𝑠 (𝑃𝑃�𝑖𝑖+1 )]
𝐵𝐵𝑔𝑔 (𝑃𝑃�𝑖𝑖 ) 𝐵𝐵𝑜𝑜 (𝑃𝑃�𝑖𝑖+1 )

where 𝑖𝑖 → (𝑖𝑖 + 1) denotes fluid contribution from portion 𝑖𝑖 to portion (𝑖𝑖 + 1)

Equations (11) through (13) are adaptations of the equations presented earlier. They are used
for each saturated portion up until the last portion where the excess fluid contribution becomes the
production as shown in equation (9). Because this is a direct method, a simple spreadsheet can be
created by implementing all equations into the system and allowing 𝑥𝑥𝑏𝑏 to travel horizontally. In
summary, spreadsheet steps for an initially undersaturated, linear flow case are shown in Figure 5.
Note that water flow can easily be implemented to the methodology by a similar fashion to oil
expansion and contribution to neighboring regions. Relative permeabilities would then need to be
integrated using Stone II model.

Figure 5 – Method flowchart where “n” is the number of partitions in the saturated region
 11

Validation: Comparison with numerical simulation

The proposed segmentation approach is compared to numerical simulation as a mean of
validation for a black oil and a volatile oil model. We consider a simple, vertical and infinite
conductivity fracture model that penetrates the entire thickness of the formation as shown in Figure
2. A numerical model was built and ran using a finite difference simulator with a logarithmic cell
distribution (Panja et al.,2013) around the fracture to capture steep pressure declines as shown in
Figure 6.

Figure 6 –Pressure distribution of a hydraulic fracture draining from the matrix at an

arbitrary time

Along the horizontal well, a significant difference in reservoir pressure is observed. Pressure
along the fracture decline to the bottom-hole pressure of 1000 psi while pressure just a few feet
away is at the initial conditions as shown in Figure 6. The main consideration in conventional
material balance is the calculation of fluid properties at average reservoir pressure which is
determined at the moment when the production data is collected. It is also assumed that the flow
is boundary dominated i.e., reservoir has been depleted enough to produce from entire reservoir.
To measure the reservoir pressure, well is shut in for hours or days for equilibration of pressure
throughout the reservoir. On the other hand, in case of low permeability reservoir, initial pressure
front reaches only few feet far from the fracture even after few years. In some cases, initial pressure
front even never reaches boundary in the entire life of the well. Therefore, average pressure which
represents the entire reservoir is not applicable for partially depleted reservoir such as tight
formations. Under these circumstances, conventional material balance calculation cannot be
applied.
 12

Model key parameters for both cases are shown in Table 1 and Figure 7.

Table 1 – Reservoir simulation parameters

Reservoir Parameter Black Volatile Oil
Oil
Thickness, ft 200 200
Hydraulic fracture half-length, ft 500 500
Fracture spacing, ft 200 200
Initial pressure, psi 5500 6500
API gravity 36.5 45
Bubble point pressure, psi 1615 5015
Initial solution gas-oil ratio, scf/stb 745 2000
Initial oil formation volume factor, bbl/stb 1.08 1.80
Oil formation volume factor at the bubble point, bbl/stb 1.20 1.87
Initial oil viscosity, cp 0.47 0.19
Bottom-hole pressure, psi 1000 1000
Matrix porosity 0.08 0.08
Matrix permeability, md 0.0005 0.0005
Critical gas saturation, fraction 0.01 0.005
krg at connate liquid, fraction 0.87 0.74
kro at connate gas, fraction 0.95 0.74
ng 1.8 1.86
nog 2 4

Figure 7 – Gas formation volume factor and viscosity for the black oil and volatile oil cases
The entire study is focused on the two phase flow of oil and gas. Therefore, the initial water
saturation is kept equal to one and water saturation is zero.
Similar to numerical simulation gridding, we must make a choice of Δxb values. Fortunately,
we can afford small Δxb steps since the calculations are direct and no iterative solvers are
necessary. Choosing a bubble point distance step size of 0.5 feet and a saturated reservoir pressure
step size of 5 psi, we built a multiphase cumulative production spreadsheet with PVT functions
 13

and relative permeability ratios based on the steps shown in Figure 5. Note that bottom-hole
pressures decline from initial pressures to their constant values in table 1.

Black Oil Case

Figure 8 shows a good agreement in terms of oil and gas cumulative production for the black
oil case considering that the method does not require of permeability as an input.

(a) (b)
Figure 8 – (a) Cumulative oil and gas production and (b) cumulative gas oil ratio
simulation results compared to the proposed method for a black oil case

If the number of hydraulic fractures is uncertain and their extent unknown, these values may
be assumed to be one. Essentially, simplifying reservoir dimensions to unity leads to calculation
of recovery factors instead of cumulative production. The number and extent of hydraulic fractures
can be estimated by iterating initial guesses until the calculated production matches historical
production. This procedure can be used to estimate the original oil in place in the stimulated
reservoir volume.
The proposed method predicts a maximum recovery factor of 29.5% if the reservoir average
pressure declines close to the bottom-hole pressure (assuming no other production enhancement
techniques are applied). This value matches the simulation recovery factor under similar
conditions. The average pressures are also compared in Figure 9.
 14

Figure 9 – Comparison of average reservoir pressure obtained from flow simulation and
proposed model in black oil case
The average pressure from the proposed model matches very well with the average pressure
from commercial simulator. The match is almost perfect at higher pressure side due to bigger
undersaturated region with less gas evolution.
We found that the reservoir becomes fully saturated at an average pressure of almost 1500 psi.
Note that even though the saturation pressure is 1615 psi, the entire reservoir becomes fully
saturated only when xb reaches the boundary of the reservoir. When this happens, there is no point
in the reservoir with a pressure value above the bubble point, hence the average reservoir pressure
must be less than the bubble point pressure. Volumetric integration of the pressure profile results
in an average pressure of 1468 psi which is just below the value predicted by the proposed method.

Figure 8 (a) shows a very good agreement in terms of cumulative oil and gas production
proving that accurate results can be achieved with only two saturated regions. The reservoir
becomes fully saturated at an average pressure of around 4070 psi compared to 4130 psi obtained
from simulation. The excess volumes of fluid of 2 segments for saturated section (Es1 and Es2)
and 1 segment for undersaturated section (Eu) are calculated using equation 12 and 2 respectively
as shown in Figure 10
 15

Figure 10 – The changes in excess volume in two saturated segments and the single
undersaturated segment in black oil case
Boundary (xb) of saturated and undersaturated sections moves away from fracture with time.
Consequently, the average pressure of each segment decreases and this causes the reduction of
excess volumes for both section. Excess volume in undersaturated section i.e., Eu is always greater
than excess volume in saturated sections (Es) because of higher average pressure in undersaturated
section. At higher pressure the formation volume factor as well as the solution gas oil ratio is
higher.

Accuracy can be enhanced by partitioning the saturated region many times to capture steep
saturation changes, especially for cases with high gas oil ratios. This practice, however, turns this
methodology into a framework akin to a finite-difference simulator. In an effort to keep the
workflow as simple as possible, we partitioned the saturated region into two equal portions.
The pressure and gas saturation profiles for the simulation model are shown in Figure 11 to
draw comparison with the conceptual Figure 3 (b).
 16

Figure 11 – Gas saturation and pressure profiles from simulation

Volatile Oil Case

Two distinct production trends are shown in Figure 8(a) and Figure 12 (a).

(a) (b)
Figure 12 – (a) Cumulative oil and gas production and (b) cumulative gas oil ratio simulation
results compared to the proposed method for a volatile oil case
Figure 8 (b) and Figure 12 (b) shows a cumulative GOR plot that features an initial plateau
higher than the initial solution gas-oil ratio. This behavior is seen in many wells in the Bakken and
other tight oil linear flow reservoirs and discussed in the literature (Jones, 2016).

In the black oil case, when the reservoir becomes fully saturated (xb=Lf) both oil and gas
production see a steep production slope versus pressure where cumulative oil production increases
by about 70% in a span of 500 psi. In the volatile oil case, production after the reservoir becomes
fully saturated increases by over 100% in the span of 3000 psi. This difference in behavior is a
direct result of different fluid expansion characteristics (black oil versus volatile oil) from the PVT
data. The absolute maximum recovery with no production enhancement is calculated to be 22%
 17

which matches the output from simulation. A comparison of average pressure from proposed
model and commercial simulator is given in Figure 13.

Figure 13 – Comparison of average reservoir pressure obtained from flow simulation and
proposed model in volatile oil case
Similar to black oil case, the average pressure calculated using the proposed model matches very
well with results from commercial simulator in case of volatile oil.
The excess volumes of 2 saturated segments (Es1 and Es2) and 1 undersaturated segment (Eu)
for volatile oil case are shown in Figure 14

Figure 14 – The changes in excess volume in two saturated segments and the single
undersaturated segment in volatile oil case
 18

The decrease in excess volume is less in undersaturated segment in volatile oil case because of
the fact that the higher amount of dissolved gas is evolved which sustains the pressure of the
section.

Field Applications
We begin with the analysis of an Eagle Ford shale horizontal well that produced for over a year
(Well A). PVT data as well as other parameters are shown in table 2. Bottomhole pressure data as
well as production as shown in Figure 15.

(a) (b)
Figure 15-Data for Well A (a) bottom-hole pressure and (b) fluid production.

Table 2 – Well A parameters

Parameters Values
Initial Pressure, psi 5500
Bubble Point Pressure, psi 2800
Initial solution GOR, scf/stb 500
Initial Boi, rb/stb 1.4
Number of hydraulic fractures 60
Estimated fracture half-length, ~400
ft
Formation thickness, ft 180-200
Average porosity, fraction 0.05

Similarly, to the simulation examples, a spreadsheet was set up with a bubble point distance
step size of 0.5 feet and a saturated reservoir pressure step size of 5 psi. Relative permeability
ratios were obtained from in table 1 where critical gas saturation is 0.05 for this case. Production
data as a function of average reservoir pressure is not available for this well. To overcome this
obstacle, we plotted cumulative oil versus cumulative gas and compared that to the proposed
method’s results as shown in Figure 16 (a). Cumulative oil and cumulative gas are calculated from
proposed model using the PVT data, initial conditions, operational parameters, and completion
 19

parameters for the field. Predicted cumulative fluid production versus average pressure are shown
in Figure 16 (b).

(a) (b)
Figure 16-Well A (a) predicted cumulative production using different fracture half-lengths
and (b) fluid cumulative production versus average reservoir pressure.
As shown in Figure 16 (a), we plotted cumulative fluid production for different hydraulic
fracture sizes (half-lengths of 200, 400, and 500 ft). Assuming that the reported fracture half-
lengths of 400 ft are representative of the system, we can calculate OOIP volumetrically for the
stimulated portion of the reservoir. As shown in Figure 16(a), the predicted production for a
fracture half-length of 200 feet does not match the historical production. Based on this the
minimum average half-length fracture is 200 ft, and the SRV OOIP is at least 2.4 MMSTB. In this
fashion, uncertain hydraulic fractures can be determined by plotting cumulative data versus the
results from this method. Only future production will decidedly corroborate the reported 400 ft
half-lengths and the actual OOIP. Using the methodology, we are also able to predict a 10.3%
maximum oil recovery regardless of fracture dimensions. This is due to high bottom-hole pressures
and can be improved by the addition of gas lift.
We now apply the methodology to two prolific shale/tight formation plays, namely Bakken in
North Dakota /Montana and Eagle Ford in Texas where information is scarce. Field parameters
are given in Table 3.

Table 3 – Key field parameters

Parameters Bakken Eagle Ford
(Well B) (Well C)
Initial Pressure, psi 4600 7000
Flowing Bottom-Hole Pressure, psi ~970 ~1000
Bubble Point Pressure, psi 2800 3644
Initial solution GOR, SCF/STB 700 1577
Initial Boi, rb/STB 1.44 1.70

As shown in Figure 17, production data from both wells match the predicted data using
proposed the model. This shows good applicability of the method to multi-phase flow in field
cases. Calculation of original oil in place (OOIP) is the direct implication of this results if the
geometry of the reservoir is unknown. If field data is available along with initial reservoir
conditions, a unit-size system can be created using the proposed method. The unit-size results
 20

would then be multiplied by a guess of OOIP until a match like the ones shown in Figure 16 (a)
and 17 is achieved.

(a) (b)
Figure 17-Comparison of field data with the proposed method for a (a) Bakken well and an
(b) Eagle Ford well.

Conclusions
We have introduced a quick and direct method based on material balances with the aim to analyze
multiphase reservoir performance during transient and boundary dominated flow. The method
makes use of a traveling variable that divides the reservoir into regions depending on the
developing phases. This technique was verified by comparing oil and gas production using several
numerical simulation models, two of which were presented here for black and volatile oils.
Field applications have been explored and performance analysis suggested even when data is
scarce. One clear outcome from the use of this method is the estimation of OOIP which entails
knowledge of fracture surface area for horizontal wells. OOIP is based upon the stimulated
reservoir volume (SRV) defined by the dimensions of hydraulic fractures. The method relies on
the availability of PVT data, relative permeability ratios, and assumptions regarding average
pressures. Application of this method is not limited to geometry and technically not limited to type
of fluids. Condensates, for instance, may be coupled into the method by introduction of Rv data.
The proposed methodology does not require absolute permeabilities as input. The price we pay is
that the current framework does not account for time.

Nomenclature
Symbol Description Units
‘ Apostrophe denoting parameter is evaluated at the previous step -
2
α Diffusivity coefficient ft /sec
φ Porosity -
γ Constant in equation A-1 1/ft
�
𝑃𝑃 Average pressure psi
�
𝑃𝑃𝑆𝑆1 Average pressure of first saturation part psi
 21

𝑃𝑃�𝑆𝑆1 Average pressure of second saturation part psi

𝑃𝑃�𝑆𝑆2 Average pressure of second saturation part psi
���
𝑃𝑃𝑅𝑅 Average reservoir pressure psi
𝑃𝑃�𝑆𝑆 Average pressure in saturation region psi
𝑃𝑃�𝑆𝑆′ Average pressure in saturation region at previous step psi
���
𝑃𝑃𝑈𝑈 Average pressure of the undersaturated region psi
𝑃𝑃�𝑈𝑈′ Average pressure of the undersaturated region at previous step psi
Volume of gas contributed by ith portion to (i+1)th portion in
𝐶𝐶𝑖𝑖→(𝑖𝑖+1),𝑔𝑔 rb
saturated region
Volume of oil contributed by ith portion to (i+1)th portion in
𝐶𝐶𝑖𝑖→(𝑖𝑖+1),𝑜𝑜 rb
saturated region
Volume of gas from nth portion or last portion of saturated region to
𝐶𝐶𝑛𝑛→𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠,𝑔𝑔 scf
surface
Volume of oil from nth portion or last portion of saturated region to
𝐶𝐶𝑛𝑛→𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠,𝑜𝑜 stb
surface
𝐶𝐶𝑠𝑠→𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠,𝑔𝑔 Volume of oil from saturated region to surface scf
𝐶𝐶𝑠𝑠→𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠,𝑜𝑜 Volume of oil from saturated region to surface stb
Volume of gas contributed by undersaturated region to saturated
𝐶𝐶𝑢𝑢→𝑠𝑠,𝑔𝑔 rb
region
Volume of oil contributed by undersaturated region to saturated
𝐶𝐶𝑢𝑢→𝑠𝑠,𝑜𝑜 rb
region
µ Viscosity cp
µg Gas viscosity cp
µo Oil viscosity cp
Bg Gas formation volume factor rb/scf
Bo Oil formation volume factor rb/stb
c1 Constant in equation A-4 psi/ft
cf Formation compressibility 1/psi
ct Total compressibility 1/psi
ES Excess volume in the saturated region rb
Eu Excess volume in the undersaturated region rb
f(xb) Pressure gradient at the bubble point pressure psi/ft
fg Fractional gas flow -
fo Fractional oil flow -
fw Fractional water flow -
h Formation thickness ft
k Absolute permeability mD
krg Gas relative permeability -
kro Oil relative permeability -
 22

Limits of integration if the saturated region is divided into smaller

l1, l2 ft
portions
Lf Half fracture spacing ft
P Pressure psi
Pb Bubble point pressure psi
Pi Initial reservoir pressure psi
Pref Reference Pressure psi
Pwf Flowing Bottom Hole pressure psi
Qg Volume of gas produced at surface scf
Qo Volume of oil produced at surface stb
Rs Solution Gas to Oil Ratio scf/stb
Rsi Initial solution gas to oil ratio scf/stb
Sg Gas saturation -
Sgc Critical gas saturation -
So Oil saturation -
Vi,g Volume of gas in ith portion of saturation region rb
′
𝑉𝑉𝑖𝑖,𝑔𝑔 Volume of gas in ith portion of saturation region at previous step rb
Vi,o Volume of oil in ith portion of saturation region rb
′
𝑉𝑉𝑖𝑖,𝑜𝑜 Volume of oil in ith portion of saturation region at previous step rb
Vs,g Volume of gas in saturation region rb
′
𝑉𝑉𝑠𝑠,𝑔𝑔 Volume of gas in saturation region at previous step rb
Vs,o Volume of oil in saturation region rb
𝑉𝑉𝑠𝑠′,𝑜𝑜 Volume of oil in saturation region at previous step rb
xb Bubble point distance from fracture ft
𝑥𝑥𝑏𝑏′ Bubble point distance from fracture at previous step ft
xf Fracture half length ft
xs2 Distance of boundary of first saturation portion ft
xsgc Distance where gas saturation reaches its critical point ft

References
Amorim, Tiago C A De and Denis Jose Schiozer (2012). Risk Analysis Speed-Up With Surrogate
Models. SPE Latin America and Caribbean Petroleum Engineering Conference. Mexico City,
Mexico, Society of Petroleum Engineers. SPE153477.
Arps, J. J. (1945). "Analysis of Decline Curves." Transactions of the AIME SPE-945228-G.
Carreras, Patricia Elva, Scott Edward Turner and Gwendolyn Tharp Wilkinson (2006). Tahiti:
Development Strategy Assessment Using Design of Experiments and Response Surface Methods.
SPE Western Regional/AAPG Pacific Section/GSA Cordilleran Section Joint Meeting.
Anchorage, Alaska, USA, Society of Petroleum Engineers.
 23

Clark, Aaron James, Larry Wayne Lake and Tadeusz Wiktor Patzek (2011). Production
Forecasting with Logistic Growth Models. SPE Annual Technical Conference and Exhibition.
Denver, Colorado, USA, Society of Petroleum Engineers. SPE-144790-MS.
Clarkson, Christopher R. and Joshua Beierle (2010). Integration of Microseismic and Other Post-
Fracture Surveillance with Production Analysis: A Tight Gas Study. SPE Unconventional Gas
Conference. Pittsburgh, Pennsylvania, USA, Society of Petroleum Engineers. SPE-131786-MS.
Coleman, Stewart, H. D. Wilde, Jr. and Thomas W. Moore (1930). "Quantitative Effect of Gas-oil
Ratios on Decline of Average Rock Pressure."
Dahaghi, Amirmasoud Kalantari, Soodabeh Esmaili and Shahab D. Mohaghegh (2012). Fast Track
Analysis of Shale Numerical Models. SPE Canadian Unconventional Resources Conference.
Calgary, Alberta, Canada, Society of Petroleum Engineers. SPE 162699.
Duong, Anh N. (2010). An Unconventional Rate Decline Approach for Tight and Fracture-
Dominated Gas Wells. Canadian Unconventional Resources and International Petroleum
Conference. Calgary, Alberta, Canada, Society of Petroleum Engineers. SPE-137748-MS.
E. El-Sebakhy, T. Sheltami, S. Al-Bokhitan, Y. Shaaban, I. Raharja, Y. Khaeruzzaman (2007).
Support Vector Machines Framework for Predicting the PVT Properties of Crude-Oil Systems,
Kingdom of Baharin, 15th SPE Middle East Oil & Gas Show and Conference.
El-Banbi, Ahmed H. and Robert A. Wattenbarger (1998). Analysis of Linear Flow in Gas Well
Production. SPE Gas Technology Symposium. Calgary, Alberta, Canada, Society of Petroleum
Engineers. SPE-39972-MS.
Fatai Adesina Anifowose, AbdlAzeem Oyafemi Ewenla, Safiriyu Ijiyemi (2011). Prediction of Oil
and Gas Reservoir Properties using Support Vector Machines, Bangkok, Thailand, International
Petroleum Technology Conference,.
Fetkovich, M. J. (1980). "Decline Curve Analysis Using Type Curves." Journal of Petroleum
Technology 32(06).
Fevang, Øivind and C.H. Whitson (1996). "Modeling Gas-Condensate Well Deliverability." SPE
Reservoir Engineering 11(4): 221-230.
Havlena, D. and A.S. Odeh (1963). The Material Balance as an Equation of a Straight Line.
Ilk, Dilhan, Jay Alan Rushing, Albert Duane Perego and Thomas Alwin Blasingame (2008).
Exponential vs. Hyperbolic Decline in Tight Gas Sands: Understanding the Origin and
Implications for Reserve Estimates Using Arps' Decline Curves. SPE Annual Technical
Conference and Exhibition. Denver, Colorado, USA, Society of Petroleum Engineers.
Ismadi, Danu, C. Shah Kabir and A. Rashid Hasan (2011). The Use of Combined Static and
Dynamic Material-Balance Methods in Gas Reservoirs. SPE Asia Pacific Oil and Gas Conference
and Exhibition. Jakarta, Indonesia, Society of Petroleum Engineers: 16.
Kabir, Shah, Faisal Rasdi and B. Igboalisi (2011). "Analyzing Production Data From Tight Oil
Wells." Journal of Canadian Petroleum Technology 50(05): 48-58.
Levine, J. S. and M. Prats (1961). "The Calculated Performance of Solution-Gas-Drive
Reservoirs." Society of Petroleum Engineers Journal 1(3).
Li, B. and F. Friedmann (2005). Novel Multiple Resolutions Design of Experiment/Response
Surface Methodology for Uncertainty Analysis of Reservoir Simulation Forecasts. SPE Reservoir
Simulation Symposium. The Woodlands, Texas, 2005,. Society of Petroleum Engineers Inc. SPE
92853.
Mannon, Robert W. (1965). Oil Production Forecasting By Decline Curve Analysis, Society of
Petroleum Engineers.
 24

Mattar, L. and D. Anderson (2005). Dynamic Material Balance(Oil or Gas-In-Place Without Shut-
Ins). Canadian International Petroleum Conference. Calgary, Alberta, Petroleum Society of
Canada: 10.
Mead, Homer N. (1956). Modifications to Decline Curve Analysis, Society of Petroleum
Engineers.
Medeiros, Flavio, Basak Kurtoglu, Erdal Ozkan and Hossein Kazemi (2010). "Analysis of
Production Data From Hydraulically Fractured Horizontal Wells in Shale Reservoirs." SPE
Reservoir Evaluation & Engineering 13(03): 559-568.
Miller, G. Frank (1962). "Theory of Unsteady state Influx of Water in Linear Reservoirs." Journal
of the Institute of Petroleum 48.
Millikan, C. V. (1926). "Gas-oil Ratio as Related to the Decline of Oil Production, with Notes on
the Effect of Controlled Pressure." Transactions of the AIME SPE-926147-G.
Muskat, Morris (1945). "The Production Histories of Oil Producing Gas‐Drive Reservoirs."
Journal of Applied Physics 16(3): 147-159.
Orangi, Abdollah, Narayana Rao Nagarajan, Mehdi Matt Honarpour and Jacob J. Rosenzweig
Unconventional Shale Oil and Gas-Condensate Reservoir Production, Impact of Rock, Fluid, and
Hydraulic Fractures, Society of Petroleum Engineers.
Palacio, J. C. and T. A. Blasingame (1993). Decline-Curve Analysis With Type Curves - Analysis
of Gas Well Production Data. Low Permeability Reservoirs Symposium. Denver, Colorado,
Society of Petroleum Engineers: 1.
Panja, Palash, Tyler Conner and Milind Deo (2013). "Grid sensitivity studies in hydraulically
fractured low permeability reservoirs." Journal of Petroleum Science and Engineering 112: 78-87.
Panja, Palash, Tyler Conner and Milind Deo (2016). "Factors Controlling Production in
Hydraulically Fractured Low Permeability Oil Reservoirs." International Journal of Oil, Gas and
Coal Technology (In Press).
Panja, Palash and Milind Deo (2016). "Factors That Control Condensate Production From Shales:
Surrogate Reservoir Models and Uncertainty Analysis." SPE Reservoir Engineering and
Evaluation SPE-179720-PA.
Panja, Palash and Milind Deo (2016). "Unusual Behavior of Produced Gas Oil Ratio in Low
Permeability Fractured Reservoirs." Journal ofPetroleumScienceandEngineering 144: 76-83.
Panja, Palash, Manas Pathak, Raul Velasco and Milind Deo (2016). Least Square Support Vector
Machine: An Emerging Tool for Data Analysis. SPE Low Perm Symposium. Denver, Colorado,
USA, Society of Petroleum Engineers. SPE-180202-MS.
Parikh, Harshal (2003). Reservoir Characterization Using Experimental Design And Response
Surface Methodology Master of Science, Texas A&M University.
Pathak, Manas, Milind Deo, Palash Panja and Levey Raymond (2015). The Effect of Kerogen-
Hydrocarbons Interaction on the PVT Properties in Liquid Rich Shale Plays. SPE/CSUR
Unconventional Resources Conference. Calgary, Alberta, Canada. SPE-175905-MS.
Schilthuis, Ralph J. (1936). "Active Oil and Reservoir Energy."
Tarner, J (1944). "How different size gas caps and pressure maintenance programs affect amount
of recoverable oil." Oil Weekly 144.
Tracy, G. W. (1955). "Simplified Form of the Material Balance Equation." Transaction of AIME
204.
Valko, Peter P. (2009). Assigning value to stimulation in the Barnett Shale: a simultaneous analysis
of 7000 plus production hystories and well completion records. SPE Hydraulic Fracturing
Technology Conference. The Woodlands, Texas, Society of Petroleum Engineers: 19.
 25

Wattenbarger, Robert A., Ahmed H. El-Banbi, Mauricio E. Villegas and J. Bryan Maggard (1998).
Production Analysis of Linear Flow Into Fractured Tight Gas Wells. SPE Rocky Mountain
Regional/Low-Permeability Reservoirs Symposium, Society of Petroleum Engineers.
Willhite, G.P. (1986). Waterflooding. Dallas, texas, USA, Society of Petroleum Engineers.

Appendix: A
Average Pressure Estimations

Average pressure approximations rely heavily upon the particular case study and should be derived
individually. In the present work, we used modified versions of known pressure equations for
single phase linear flow as presented by Miller (Miller,1962) .These equations were initially
derived from the study of heat transfer with different boundary conditions and accommodated for
linear flow in infinite and bounded systems. The exact solution to the single phase, incompressible
flow is calculated as:
∞
𝑝𝑝(𝑥𝑥, 𝑡𝑡) − 𝑝𝑝𝑤𝑤𝑤𝑤 4 𝑒𝑒𝑒𝑒𝑒𝑒(−𝛾𝛾 2 𝛼𝛼 𝑡𝑡) 𝑠𝑠𝑠𝑠𝑠𝑠(𝛾𝛾𝛾𝛾)
= �� � (A-1)
𝑝𝑝𝑖𝑖 − 𝑝𝑝𝑤𝑤𝑤𝑤 𝜋𝜋 2𝑛𝑛 − 1
𝑛𝑛=1

(2𝑛𝑛−1) 𝜋𝜋 𝑘𝑘
where, 𝛾𝛾 = , and 𝛼𝛼 = 𝜙𝜙𝜙𝜙𝑐𝑐
2 𝐿𝐿𝑓𝑓 𝑡𝑡
The exact pressure profile for the undersaturated portion of a transient system can be
determined with by modifying 𝑥𝑥 and 𝐿𝐿𝑓𝑓 and iteratively solving for the bubble point pressure
distance. The average pressure can be determined by analytical integration as shown in equation
(1). Undersaturated average pressure calculations using this equation were done with a simple
Newtonian solver, however, the process takes away from the proposed direct method and requires
substantial additional information such as time. We found that an approximation to this equation
must comply with the following conditions:

𝑃𝑃�𝑈𝑈 (𝑥𝑥𝑏𝑏 = 0) = 𝑃𝑃𝑖𝑖

𝑃𝑃�𝑈𝑈 �𝑥𝑥𝑏𝑏 = 𝐿𝐿𝑓𝑓 � = 𝑃𝑃𝑏𝑏

An approximation based on these conditions can be used for a number of cases as:

𝑥𝑥𝑏𝑏 𝑚𝑚
𝑃𝑃�𝑈𝑈 (𝑥𝑥𝑏𝑏 ) = 𝑝𝑝𝑤𝑤𝑤𝑤 + (𝑝𝑝𝑖𝑖 − 𝑝𝑝𝑤𝑤𝑤𝑤 ) �1 − � (A-2)
𝐿𝐿𝑓𝑓 𝑚𝑚

where 𝑚𝑚 = 1 for linear pressure decline cases.

The linear form of equation (A-2) is a decent approximation to the analytical solution for
several cases. This a preferred way to estimate undersaturated average pressures since it requires
to iterative methods.
The saturated portion of the reservoir can be estimated using a variety of ways. In our analysis,
we start at the diffusivity equation for a linear system:
 26

𝜕𝜕𝜕𝜕 𝜕𝜕 2 𝑝𝑝
= 𝛼𝛼 2 (A-3)
𝜕𝜕𝜕𝜕 𝜕𝜕𝑥𝑥
𝜕𝜕𝜕𝜕
For a saturated portion of a linear system, we estimate the pressure profile by assuming that 𝜕𝜕𝜕𝜕
is constant for a particular value of 𝑥𝑥𝑏𝑏 . This condition has been used before for pseudo-steady state
1 𝜕𝜕𝜕𝜕
flow. Letting 𝛼𝛼 𝜕𝜕𝜕𝜕 = 𝑎𝑎, rearranging terms and integration gives:

𝜕𝜕𝜕𝜕
= 𝑎𝑎𝑎𝑎 + 𝑐𝑐1 (A-4)
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
At this point, pseudo-steady state dictates that �
𝜕𝜕𝜕𝜕 𝑥𝑥=𝑥𝑥𝑏𝑏
= 0, however, we deviate from this
condition because the undersaturated region contributes to the flow. Hence, the pressure gradient
𝜕𝜕𝜕𝜕
at 𝑥𝑥𝑏𝑏 is not zero, but a function of 𝑥𝑥𝑏𝑏 . Letting 𝜕𝜕𝜕𝜕 � = 𝑓𝑓(𝑥𝑥𝑏𝑏 ) and solving for 𝑐𝑐1:
𝑥𝑥=𝑥𝑥𝑏𝑏

𝜕𝜕𝜕𝜕
= 𝑎𝑎𝑎𝑎 + 𝑓𝑓(𝑥𝑥𝑏𝑏 ) − 𝑎𝑎𝑥𝑥𝑏𝑏 (A-5)
𝜕𝜕𝜕𝜕

Integration and application of bottom-hole boundary condition yields:

𝑥𝑥 2
𝑝𝑝(𝑥𝑥) = 𝑝𝑝𝑤𝑤𝑤𝑤 + 𝑥𝑥𝑥𝑥(𝑥𝑥𝑏𝑏 ) + 𝑎𝑎 �𝑥𝑥𝑥𝑥𝑏𝑏 − � (A-6)
2

where 𝑥𝑥𝑏𝑏 is a constant at each step.

There is one more boundary condition that must be met: 𝑝𝑝(𝑥𝑥𝑏𝑏 ) = 𝑝𝑝𝑏𝑏 , this can be achieved by
solving for 𝑎𝑎 or by scaling the pressure function. Both methods yield the same solution.
Keeping in mind that the saturated region may be divided in n parts, application of scaling and
equation (1) gives us an average pressure expression for a single phase incompressible fluid. The
saturated, region, however contains important amounts of gas that grow during transient flow and
affect production substantially once the reservoir is fully saturated. By application of square-
pressures analogy:

2 (𝑙𝑙13 − 𝑙𝑙23 )
2
�𝑝𝑝𝑏𝑏 − 𝑝𝑝𝑤𝑤𝑤𝑤 2 � (𝑙𝑙1 − 𝑙𝑙22 )𝑓𝑓(𝑥𝑥𝑏𝑏 ) + �(𝑙𝑙12 − 𝑙𝑙22 )𝑥𝑥𝑏𝑏 − � (A-7)
� 3
𝑃𝑃�𝑆𝑆 = 𝑃𝑃𝑤𝑤𝑤𝑤 2 +
(𝑙𝑙1 − 𝑙𝑙2 )𝑥𝑥𝑏𝑏 2𝑓𝑓(𝑥𝑥𝑏𝑏 ) + 𝑥𝑥𝑏𝑏

where, 𝑙𝑙1 and 𝑙𝑙2 are the limits of integration if the saturated region is divided into smaller
portions. If the saturated region consists of only one cell, then 𝑙𝑙1 = 𝑥𝑥𝑏𝑏 and 𝑙𝑙2 = 0.
Recalling that 𝑓𝑓(𝑥𝑥𝑏𝑏 ) is the pressure gradient at the bubble point pressure as the saturated region
grows, we invoke two conditions:
𝜕𝜕𝜕𝜕 𝑝𝑝𝑏𝑏 − 𝑝𝑝𝑤𝑤𝑤𝑤
� =
𝜕𝜕𝜕𝜕 𝑥𝑥𝑏𝑏 →0 𝑥𝑥𝑏𝑏
 27

𝜕𝜕𝜕𝜕
� =0
𝜕𝜕𝜕𝜕 𝑥𝑥𝑏𝑏 =𝐿𝐿𝑓𝑓

According to these conditions, 𝑓𝑓(𝑥𝑥𝑏𝑏 ) can be approximated as:

𝑝𝑝𝑏𝑏 − 𝑝𝑝𝑤𝑤𝑤𝑤 �𝑝𝑝𝑏𝑏 − 𝑝𝑝𝑤𝑤𝑤𝑤 �

𝑓𝑓(𝑥𝑥𝑏𝑏 ) = − 𝑥𝑥𝑏𝑏 (A-8)
𝑥𝑥𝑏𝑏 𝐿𝐿2𝑓𝑓

Equation (A-2) was used for the majority of cases, whereas a few required the iterative method
for the estimation of the undersaturated average pressure. However, Equation (A-7) worked very
well for all cases ranging from black oils to volatile oils.

Appendix: B
Calculation of Excess Volume

Oil occupying the reservoir volume extended 'x' distance from fracture face is calculated as

2𝑥𝑥𝑓𝑓 ℎ 𝑥𝑥 𝑆𝑆𝑜𝑜 (𝑃𝑃�)𝜙𝜙(𝑃𝑃�)

𝑂𝑂𝑂𝑂𝑂𝑂 𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣 𝑎𝑎𝑎𝑎 𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟 𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐 = (𝑟𝑟𝑟𝑟) (B-1)
5.615

𝑃𝑃� denotes the average reservoir pressure in the section where the calculation is performed.
Conversion of the volume at reservoir conditions to stock tank conditions using formation volume
factor is shown in equation B-2

2𝑥𝑥𝑓𝑓 ℎ 𝑥𝑥 𝑆𝑆𝑜𝑜 (𝑃𝑃�)𝜙𝜙(𝑃𝑃�)

𝑂𝑂𝑂𝑂𝑂𝑂 𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣 𝑎𝑎𝑎𝑎 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑘𝑘 𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡 𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐 = (𝑠𝑠𝑠𝑠𝑠𝑠) (B-2)
5.615 𝐵𝐵𝑂𝑂 (𝑃𝑃�)

A simplified version of Figure 4 is shown in Figure B.1 to understand the volume change in
saturated portion due to pressure depletion below bubble point in the undersaturated portion.
 28

Figure B.1- Calculation of excess volume in the saturated portion

The length of saturated segment is increased from x'b to xb after a small time period. The
average pressures of the saturated and undersaturated segments at current time step are P �S and P
�U
� ′
respectively. The average pressure of undersaturated portion at previous time step was 𝑃𝑃𝑈𝑈 . From
the figure B.1, it is evident that the volume of oil in undersaturated zone at previous time step was

2𝑥𝑥𝑓𝑓 ℎ �𝐿𝐿𝑓𝑓 − 𝑥𝑥𝑏𝑏′ �𝑆𝑆𝑜𝑜 (𝑃𝑃�𝑈𝑈′ )𝜙𝜙(𝑃𝑃�𝑈𝑈′ )

Oil volume in undersaturated segment at previous time step = (B.3)
5.615𝐵𝐵𝑜𝑜 (𝑃𝑃�𝑈𝑈′ )

Similarly the volume of oil in undersaturated zone at current time step is calculated as shown
in Equation B.4
� 𝑈𝑈 �𝜙𝜙�𝑃𝑃
2𝑥𝑥𝑓𝑓 ℎ �𝐿𝐿𝑓𝑓 − 𝑥𝑥𝑏𝑏 �𝑆𝑆𝑜𝑜 �𝑃𝑃 � 𝑈𝑈 �
Oil volume in undersaturated segment at current time step =
� 𝑈𝑈 �
(B.4)
5.615 𝐵𝐵𝑜𝑜 �𝑃𝑃

The change in oil volume in these two consecutive time steps is calculated as
 29

Change in Oil volume in undersaturated zone

� 𝑈𝑈 �𝜙𝜙�𝑃𝑃
2𝑥𝑥𝑓𝑓 ℎ �𝐿𝐿𝑓𝑓 − 𝑥𝑥𝑏𝑏 �𝑆𝑆𝑜𝑜 �𝑃𝑃 � 𝑈𝑈 � 2𝑥𝑥𝑓𝑓 ℎ �𝐿𝐿𝑓𝑓 − 𝑥𝑥𝑏𝑏 �𝑆𝑆𝑜𝑜 �𝑃𝑃
� 𝑈𝑈 �𝜙𝜙�𝑃𝑃
� 𝑈𝑈 �
= −
5.615 𝐵𝐵𝑜𝑜 �𝑃𝑃 � 𝑈𝑈 � 5.615 𝐵𝐵𝑜𝑜 �𝑃𝑃 � 𝑈𝑈 �

′ ′
� 𝑈𝑈 � 𝜙𝜙 �𝑃𝑃
2𝑥𝑥𝑓𝑓 ℎ ⎡𝑆𝑆𝑜𝑜 �𝑃𝑃 � 𝑈𝑈 � � 𝑈𝑈 �𝜙𝜙�𝑃𝑃
𝑆𝑆𝑜𝑜 �𝑃𝑃 � 𝑈𝑈 � ⎤
⎢ ′
= � 𝐿𝐿𝑓𝑓 − 𝑥𝑥𝑏𝑏 � − � 𝑈𝑈 �
� 𝐿𝐿𝑓𝑓 − 𝑥𝑥𝑏𝑏 ⎥
�
5.615 ⎢ 𝐵𝐵 �𝑃𝑃 �′ � 𝐵𝐵𝑜𝑜 �𝑃𝑃 ⎥
⎦ (B.5)
𝑜𝑜 𝑈𝑈
⎣

The section (as marked by light green in Figure B.1) where this change occurred has now
become the part of saturated region and the reservoir condition has also been changed. The change
in oil volume calculated in equation B.5 is at stock tank conditions. This volume which is referred
as excess volume (Eu) is converted back to reservoir conditions which is basically the conditions
of saturated regions now. Therefore, the formation volume factor of oil at average pressure at
saturated region at current time step is used in the conversion as shown in equation B.6

2𝑥𝑥𝑓𝑓 ℎ𝐵𝐵𝑜𝑜 (𝑃𝑃�𝑆𝑆 ) 𝑆𝑆𝑜𝑜 (𝑃𝑃�𝑈𝑈′ )𝜙𝜙(𝑃𝑃�𝑈𝑈′ ) ′

𝑆𝑆𝑜𝑜 (𝑃𝑃�𝑈𝑈 )𝜙𝜙(𝑃𝑃�𝑈𝑈 )
𝐸𝐸𝑢𝑢 = � �𝐿𝐿𝑓𝑓 − 𝑥𝑥𝑏𝑏 � − �𝐿𝐿𝑓𝑓 − 𝑥𝑥𝑏𝑏 �� (B.6)
5.615 𝐵𝐵𝑜𝑜 (𝑃𝑃�𝑈𝑈′ ) 𝐵𝐵𝑜𝑜 (𝑃𝑃�𝑈𝑈 )
Production Prediction.pdf (1.37 MiB) view on ChemRxiv download file