SPE-188418-MS

A Recovery Curve Method Based Workflow for Reserves Estimation of

Naturally Fractured Reservoirs – A Case Study

Justin Brand F., Technical University of Denmark; Georg M. Mittermeir, Heinemann Consulting GmbH; Gabor F.
Heinemann, Heinemann Oil GmbH

Copyright 2017, Society of Petroleum Engineers

This paper was prepared for presentation at the Abu Dhabi International Petroleum Exhibition & Conference held in Abu Dhabi, UAE, 13-16 November 2017.

This paper was selected for presentation by an SPE program committee following review of information contained in an abstract submitted by the author(s). Contents
of the paper have not been reviewed by the Society of Petroleum Engineers and are subject to correction by the author(s). The material does not necessarily reflect
any position of the Society of Petroleum Engineers, its officers, or members. Electronic reproduction, distribution, or storage of any part of this paper without the written
consent of the Society of Petroleum Engineers is prohibited. Permission to reproduce in print is restricted to an abstract of not more than 300 words; illustrations may
not be copied. The abstract must contain conspicuous acknowledgment of SPE copyright.

Abstract
Naturally fractured carbonate reservoirs (NFRs) account for a majority of the world’s proven oil (~60%)
and gas reserves (~40%). Many operators either significantly over-predict production and reserves, or non-
optimize NFRs lacking suitable methods for modelling the fracture-matrix fluid interactions. Even complex
functional shape factor based formulations fail. Only numerical fine-scale single-matrix-block (SMB)
calculations capture the matrix-fracture fluid exchange in its complexity. A workflow and a successful case
study (exceeding 80 years history) are presented.
Heinemann (2004) proposed a concept (Recovery Curve Method, RCM) utilizing predefined recovery
curves instead of functional descriptions to model the fracture-matrix interaction. Meanwhile, the necessary
methods and tools have been researched, developed and proven by field applications.
Our iterative multi-step workflow reduces uncertainty gradually. At first, matrix oil recovery is
investigated on SMB models (Steiner and Mittermeir (2017)). The resulting recovery curves are verified
and tuned with Mittermeir’s (2015) Dual Porosity Material Balance method. The link to full field simulation
is established by a numerical dual porosity column model. Finally, the full field is calculated and matched
by the RCM.
Based on a first conceptual reservoir model, burdened with uncertainty, the matrix and fracture rock
system is examined gradually. The reservoir volume, time frame and complexity covered is stepwise
extended. Using the finely gridded SMB, the physics of matrix-fracture interaction are investigated under
different drive mechanisms including single phase expansion, solution gas drive, capillary imbibition,
gravity drainage, viscous and diffusive forces. Effects of these mechanisms can be studied in an isolated
or combined way, including its time dependency. The resulting recovery curves give a good estimate of
reservoir behavior under different operating conditions. Dynamic data, such as pressure, production and
phase contact history are matched subsequently by dual porosity material balance calculations. Results are
water influx requirements (analytical aquifer models) and verified recovery curves. The latter are used to
estimate the missing or highly uncertain data used to build the conceptual model by reversing the SMB
calculation procedure. Consistency of the obtained improved data set is verified by a numerical dual porosity
column model. It resembles the twin-barrels of the material balance model, however accounts for vertical
2 SPE-188418-MS

transmissibility distribution. Individual well performance is matched in the full field model, which in
addition to the material balance is used for forecasting.
The recovery curve based workflow provides a comprehensive concept and a consistent closed loop
modelling of NFR exploitation. The shortcomings of the commonly applied shape factor concept (Kazemi
et al. (1976)), which are well documented in the literature, are eliminated by omitting an analytical transfer
function. Instead, diligently evaluated recovery curves are used directly to model the complex matrix-
fracture interaction under different recovery mechanisms at each stage of the field life.

Introduction
Objective, content and scope
The objective of this paper is to present a new workflow for the dynamic modeling of naturally fractured
reservoirs (NFRs). This workflow is based on the recovery-curve method introduced by Heinemann and
Mittermeir (2004, 2015). This paper discusses differences between a classical modeling approach and this
new variant, as well as demonstrating the superiority of the latter.
The field model utilized in this paper is derived from a full-field study. The parametrized grid model,
created with a geological modelling package, was history matched utilizing the classical Kazemi approach.
For want of another option, it is assumed that the matched model represents the true history, since it does
not represent a true handicap. Thus, a partial lack of this data set can be successfully circumvented in order
to apply the new methodology.
Finally, the dynamic-modeling workflow, using recovery curves, is applied. This consists of the following
steps, which are also depicted schematically in Figure 1:

• The recovery curves are calculated utilizing single matrix block (SMB) calculations;

• Dual-porosity material balance (MB) is calculated utilizing the previously calculated recovery
curves, before the phase contacts are matched by tuning the curves;
• SMB simulation is utilized to match the final recovery curves from the MB calculation, as well as
to acquire information about the true values of the associated matrix block height and the relative
permeability curves;
• A dual-porosity column model is simulated based on the twin-barrel model of the MB, creating
the link between the MB and full-field simulation;
• The tuned recovery curves are then utilized for a full field simulation in which the recovery-curve
method is applied.

Figure 1—Schematic representation of recovery curve based workflow for dynamic modelling of NFRs
SPE-188418-MS 3

Modelling Fracture-Matrix Fluid Transfer

All industrially -utilized methods for dynamic modeling of NFRs are based on the dual continuum approach
(DCM). At every point in space within a DCM, two values exist for each property: one for the fracture
and another for the matrix continuum. Initially, the fissures are filled with oil, and the oil from the matrix
is produced in the fissures according to the fluid/rock expansion and solution-gas drive. This is a purely
pressure-dependent process. When the water-oil or gas-oil interface reaches the matrix block, production
is mainly due to capillary and gravity effects. The difference between the DCM implementations is in their
calculation of the fracture-matrix fluid transfer. In current industry practice, the most commonly utilized
methods for calculating the matrix-fracture mass transfer are the Kazemi et al. (1969) formulation with
either the Gilman and Kazemi (1988) or the Quandalle and Sabathier (1989) gravity-drainage model.

The Reiss et al. approach

In this approach, Reiss et al. (1973) assumed that during NFR production, the phases are completely
segregated in the fissures, and the matrix block is entirely surrounded by water or gas. The matrix-fracture
exchange was represented by time-dependent source or sink functions, also known as recovery functions
or recovery curves. These functions were derived from laboratory experiments and/or from numerical
simulations of the oil-recovery mechanism in a finely gridded single-matrix-block model.
Mattax and Kyte (1962) performed water imbibition experiments. They demonstrated that in situations
in which water imbibition is the dominant force in oil displacement from the matrix, the cumulative oil
recovery by imbibition from a piece of rock surrounded by water can generally be approximated via an
exponential function, such as the Aronofsky et al. (1958) equation:
(1)
Rossen (1977) used the same approach during the semi-implicit handling of the source terms. Both Reiss
et al. and Rossen assumed that the oil recovery from a matrix block is mainly due to capillary and gravity
effects, and that it is only a function of the elapsed time after water (or gas) comes into contact with the
matrix block via the surrounding fractures. The effects of compressibility and viscous flow in the fracture
on matrix-fracture transfer were neglected by these authors, since they did not manage to derive transfer
functions that are only time-dependent and valid at each stage of the field exploitation. The professional
community did not acknowledge Reiss et al.’s assertion that this simplification is not overly restrictive to
usual field conditions. The Mattax and Kyte (1962) experiments-based, exponential recovery curves were
subsequently applied primarily to capillary imbibition. Examples of additional authors who are adopting
this approach are de Swaan (1978), Kazemi et al. (1992), Di Donato et al. (2007) and Lu et al. (2008).

The Kazemi et al. approach

Kazemi et al. (1976) extended the Warren-Root (1963) shape-factor approach, which is valid for single-
phase flows, to multiphase cases. The resultant transfer equation, in its general form, is valid for oil, water
and free gas:

(2)

where Vcell is the bulk volume of a single simulation cell; p stands for oil, water and free gas; ka is the
apparent matrix permeability; μ the viscosity; B the formation-volume factor; kr the relative permeability
and Φ the phase potential in the fracture (f) and in the matrix (m). The shape factor, σ is a characteristic value
of the matrix block, and it will be calculated based on the size and form of the individual matrix blocks.
4 SPE-188418-MS

The most complete version is the generalized Kazemi-Gilman-ElSharkawy (1992) shape factor as derived
by Heinemann and Mittermeir (2012).
Kazemi presumes that the driving force can be expressed as the difference between two discrete potential
values: one is valid at the center of the matrix block and the second is valid for every point within the
surrounding fractures. In this case, Eq. (2) is applicable. However, this is not the case for the gravitational
and viscous driving mechanisms. The difference between the average potential values is 0, and Eq. (2) is
consequently not directly applicable. To overcome this difficulty, the difference in the potential gradients
in the matrix and the fracture must be considered. With this correction in the water-oil two phase flow, the
potential differences for the water (index p = w) is expressed as:
(3)
where κ= 1/4. The term PLwo represents the potential differences originating from differences between
the water-pressure gradients in the fracture and those in the oil pressure in the matrix. Multiple authors
have made numerous supplementary attempts to make the transfer function more appropriate, without
mentionable success.

The Heinemann and Mittermeir approach

Heinemann and Mittermeir (2004, 2015, 2016) combined the two aforementioned concepts with the
underlying aim of separating the matrix-to-fracture and the fracture-to-matrix flow calculations. The Kazemi
et al. (1976) approach and one similar to that of Reiss et al. (1973) are used for the first and second of these
calculations, respectively. According to Reiss, the rate at which the fracture injects the displacing fluids
into the matrix is given by a time-dependent function. However, according to Heinemann and Mittermeir, it
should be a function of the matrix-recovery factor that is corrected to the actual pressure. The rate at which
the displaced fluid is expelled from the matrix must be calculated by the Kazemi (i.e. Warren-Root) type
transfer equation. Therefore, the basic assumption is that the forces of compressions and the capillary and
gravitational forces act independently, and their results will be superposed.
The recovery from capillary imbibition and gravitational displacement at constant pressure can be
calculated on SMB models (and also measured in laboratories), resulting in recovery curves. One example
is shown in Figure 7.15. The additional oil production during time-step Δt is estimated using the increment
of the recovery factor during the actual time-step. The following equation is the difference between the
recovery factor at time tνj + βΔj +1t and the beginning of the time-step tνj (Figure 7.15.), where β is the time-
scaling factor:
(4)
If the pressure has already dropped but the validity of Darcy’s equation is still assumed, then ΔE can be
scaled to the actual pressure and saturations:

(5)

where the index b relates to the reference pressure of the recovery curve. This type of correction was
already used by Gharsalla (2015) and Heinemann and Mittermeir (2016). The same equation is valid for
gas displacement by changing the superscript w to g. The volume of the displacing agent (water or gas) that
the fracture injects into the matrix for one standard initial oil volume is then as follows:
(6)
SPE-188418-MS 5

Figure 2—Estimation of the time and recovery factor on a gas recovery curve at the beginning of the time-step

Recovery-Curve Workflow for Numerical Simulation

When the Kazemi approach is adopted, the workflow is simply the traditional history-matching workflow,
and when the Recovery Curve Method is applied, a more complex workflow can be implemented. An
efficient workflow is proposed that consists of SMB simulation, dual-porosity MB and full-field simulation
utilizing the recovery-curve method. By using all these methods and tools, it is possible to match an observed
phase contact history. At the end of the workflow, there is a history-matched full field and MB model and
an improved understanding of the matrix-block properties.
The dynamic modeling workflow, utilizing recovery curves, consists of the following tasks:

• The first step involves the calculation of the recovery curves using SMB simulation. In this step, it is
necessary to acquire an understanding of the possible matrix block properties, especially the matrix
block height. After the estimation of a geologically reasonable matrix block height, distribution-
lumped gas and water recovery curves can be created.
• The next step is a dual-porosity MB calculation using the previously outlined recovery curves. With
this MB formulation, it is possible to analyze the matrix-fracture mass transfer. Tuned recovery
curves are obtained by matching the observed phase contacts.
• In a next step, SMB simulation is used again to match the scaled recovery curves from the MB
calculation. This is done to gain information about the true values of the matrix block height and
relative permeability curves, amongst other values.
• To ensure the existence of a link between the MB and the full-field model, a dual-porosity column
model is created based on the twin-barrel concept of the MB.
• Subsequently, those tuned recovery curves can be utilized for a full-field simulation using the
recovery curve method. In this process, a further tuning of the recovery curves is performed to
match the phase contact, pressure and production history. When acceptable matches are achieved
for those quantities, water and gas recovery curves are then derived that correctly describe the
matrix-fracture mass transfer.

The Field Case

The aforementioned workflow will now be applied to a real field case.
6 SPE-188418-MS

Background
It is difficult to identify a real field case on which the modeling workflow, with all its related concepts,
methodologies and tools, could be demonstrated in a generally valid extension. Fractured reservoirs differ
markedly in various matrix and associated fracture properties (e.g. Aguilera 1980 and Nelson 1992), whether
this be in internal connectivity (dual porosity or dual permeability), in fluid and rock properties (black
oil or compositional description, and wettability) or in depletion mechanisms (expansion, water and gas
drives, modelled either with or without oil resaturation). In most cases, the observational data related to such
reservoirs is incomplete and heavily burdened by uncertainty, with the available history covering only a
limited period of the reservoir’s overall field life. Furthermore, field data are generally treated confidentially,
which often leads to insufficient access to data to support real-world applications of new approaches. For
this reason, we will call our reservoir ABC-Field.

Field Description
The ABC-Field is a Type II NFR, and it is assumed to be of dual porosity and dual permeability. Type II
NFRs have a low matrix porosity and permeability (Allan and Sun, 2003). The matrix provides some storage
capacity, and the fractures provide the fluid pathways. The ABC-Field-model version is derived from a full-
field study. The original model was modified to protect the proprietary rights of the field owner. However,
the synthetic model still exhibits all of the original’s essential characteristics.
The ABC-Field’s original-oil-in-place (OOIP) is estimated to be 220 million sm3 with 95% stored within
the matrix. The current field-recovery factor is approximately 30%. The average reservoir permeability,
porosity and net-to-gross ratio are listed in Table 1 (matrix-block heights and shape-factor distributions will
be discussed later on).

Table 1—Summary of rock properties

property mean min. max. unit

matrix 2.23 1.00 13.94 millidarcy

X fracture 750.00 750.00 750.00 millidarcy

matrix 2.23 1.00 13.94 millidarcy

Y fracture 750.00 750.00 750.00 millidarcy

matrix 1.57 0.01 13.70 millidarcy

permeability Z fracture 100.00 100.00 100.00 millidarcy

matrix 17.79 1.00 34.63 %

porosity fracture 0.45 0.45 0.45 %

matrix 0.96 0.05 1.00

net/gross fracture 1.00 1.00 1.00

The case study’s production history recounts nearly all properties and difficulties that could arise when
dealing with a typical NFR in the Middle East. The field has been in production operation for more than
80 years (Figure 3), and the reservoir is modeled with a black-oil fluid description. With a standard density
of 739 kg/m³ or 46.1 °API, the oil can be classified as light crude oil, and is slightly under-saturated. The
standard densities for water and gas are 1094 kg/m³ and 0.788 kg/m³, respectively. The bubble-point pressure
is 155 bara, while the initial HC-weighted reservoir pressure is 168.9 bara. The reservoir is a two-phase oil-
water system (with dissolved gas) at initial conditions. During the depletion, a secondary gas cap develops,
prompting its reclassification as a three-phase gas-oil-water system. The anticlinal structure (Figure 4) is
bounded at all sides by an edge aquifer of considerable strengths. Thus, all possible drive mechanisms, such
SPE-188418-MS 7

as expansion, water and gas displacement as well as oil resaturation of the matrix, take place in certain parts
of the reservoir. Therefore, capturing the matrix-fracture mass transfer in its complexity is crucial.

Figure 3—ABC-Field production and average pressure history

Figure 4—Top depth of the ABC-Field with well locations.

The input data for the capillary pressures and the relative permeabilities can be seen in Figure 5 and Figure
6. In initializing the model, which involved calculating the initial pressure and the saturation distributions,
the drainage oil-water capillary curve is applied. The imbibition curve is utilized while simulating the oil
recovery, and the capillary pressure is set to 0 for the fracture analyses.

Figure 5—Water-oil capillary pressures and relative permeabilities

8 SPE-188418-MS

Figure 6—Gas-oil capillary pressures and relative permeabilities

Application of the Recovery-Curve Workflow for Numerical Simulation

Single Matrix Block Simulation elaboration on initial recovery and resaturation curves
Here, an SMB analysis tool, which was initially developed by Pirker et al. (2007) and integrated into
Heinemann’s (2017) research simulator H5, is utilized. The matrix-block model is discretized into grid cells.
The matrix block is homogenous; however, it can be anisotropic, having different permeability values in
x-, yand/or z-directions. Figure 7 depicts an SMB model used in this work. It allows for the calculation of
the recovery under different boundary conditions, considering capillary, gravitational and viscous forces as
well as the oil resaturation after its displacement by water or gas. An in-depth discussion of this model is
given by Steiner and Mittermeir (2017).

Figure 7—The quarter of an SMB model (a) fully surrounded by ractures and (b) bounded by vertical fractures only.

Figure 8—Shape-factor distribution for lumped single-matrix-block simulation

SPE-188418-MS 9

During the long production history of the ABC-Field, periods with an increasing oil-rim thickness have
also been observed. This means that, for some matrix blocks where the oil was already displaced, the
saturation boundary condition changes from either water or gas to the oil phase, and a resaturation of the
matrix block with oil occurs. This mass transfer is modeled with an oil resaturation curve.
Lumped recovery curves can be calculated to consider different matrix blocks inside a simulation cell.
A distribution of shape factors and matrix-block heights must be estimated. The applied shape-factor
distributions can be seen in Figure 8. Associated matrix-block heights are in the range of 2.5 to 16 m.
It should be noted that applying average values and conducting a single SMB simulation does not lead
to the same results as lumping calculated recovery curves. The resulting initially calculated recovery and
resaturation curves are presented in Figure 9.

Figure 9—Recovery curves from single-matrix-block calculations

for gas and water drive and associated oil-resaturation curves

Dual Porosity Material Balance – Tuning of recovery and resaturation curves by phase contact
matching
The application of the dual-porosity MB method, as described by Mittermeir (2015), is essential for
understanding the matrix oil-recovery mechanisms. The method is based on the recognition that the
performance of water and gas displacement from matrix blocks can be depicted by plotting recovery factors
against time. These recovery curves determine the matrix-fracture mass transfer. The reservoir’s pressure
change depends on the original fluids in place and the strength of the aquifer. Therefore, a close relationship
exists between the recovery curves and the aquifer parameters, the matrix-fracture oil transfer and the
observed reservoir state (in terms of pressure, position of phase contacts, water cut and GOR).
Within this same publication the approach is further explained as follows. The reservoir model for the
dual-porosity MB calculation is built by two columns: one representing the fracture and another, the matrix
continuum. The columns are divided into tranches (slices of unit thickness), starting from the initial OWC
and progressing upwards. The vertical distributions of the extensive values (for example, OOIP and pore
volume) are the sum of the horizontal tranches throughout the reservoir. Summing up the fluid contents
of the tranches results in two columns of barrels. One column represents the matrix and the second, the
fracture continuum. As a consequence of the vertical resolution, the distribution of the fluids and the position
of both the initial OWC and the oil-gas contact (OGC) are known. Intensive values, such as Swi (initial
water saturation), are the averaged sum of them. Figure 10 schematically illustrates the transformation
process from a 3D-populated grid model to the twin barrels of the dual-porosity MB. Furthermore, the most
10 SPE-188418-MS

important aspects of the latter are also illustrated including behaviors, such as: fluid production and water
encroachment through the fracture network only, and, as well as, the different recovery mechanisms.

Figure 10—Seven- and eight-blade 16-mm PDC bit designs considered in the analysis.

Dual-porosity MB can be run either on its own or embedded in the workflow promoted herein. Within
the workflow, the objectives are to verify the previously elaborated SMB-derived recovery and resaturation
curves and to subsequently tune them by matching the observed phase-contact movements.
As a starting point for the demonstrated phase-contact matching process, two recovery curves (water and
gas drive) and two oil-resaturation curves (resaturation after water flooding and gas flooding) are required.
For each recovery and associated resaturation curve, two factors can be used to influence the phase-contact
movements:

• Recovery-scaling factors for adjusting the ultimate recovery, and

• Timescale factors for influencing the pace of oil recovery.

The appropriate factors are determined in an iterative process until the observed and calculated phase-
contact movements match each other to a reasonable extent. The water influx from the outside aquifer and
thus the pressure support were simultaneously regulated using the target pressure method (Mittermeir et
al., 2004). This method ensures that the calculated pressure always matches the historical pressure. At the
end of such a run, an analytical aquifer model that reproduces this pressure history within close limits will
be automatically calculated.
For this case study, an acceptable history match for phase contact movements was achieved after eight
calculation runs. At this point, the aquifer model was changed from the target pressure method to the
Fetkovich analytical model. Seven additional iteration cycles were performed to fine-tune the phase contact
match.
Two sets of scaling factors for the case setups eight and 15 can be seen in Table 2. Recovery scaling
factors for water drive and gas drive are abbreviated to SRECw and SRECg, respectively, while TREC
and TRES represent the timescale factor for the recovery curve and the resaturation curves, respectively.
A visual comparison of the scaled (setup 15) and the original input recovery curve can be seen in Figure
11 and Figure 12.
SPE-188418-MS 11

Table 2—Scaling factors and applied aquifer model for the iterative phase-contact match.

# SRECw SRECg SRESw SRESg TRECw TRECg TRESw TRESg AquiModel

8 1.75 1.00 1.75 1.00 1.00 0.10 1.00 0.10 PTARGET

15 1.65 1.15 1.65 1.15 0.75 0.10 0.75 0.10 FETKOV

Figure 11—Comparison of initial and scaled (#15) water drive RC

Figure 12—Comparison of initial and scaled (#15) gas drive RC

Figure 13 depicts the calculated phase contacts as solid lines and the phase-contact history as circles. The
displayed heights are those above the initial oil-water contact. The initial oil-water contact depth is at 945
m subsea. Figure 14 illustrates the belonging pressure, oil rate, water cut and gas-oil ratio. Each quantity is
valid at the field level. The scattered data is the historical data, and the solid lines are the calculated results.
The deviations in GOR are a direct consequence of the PVT data. As long as no secondary gas cap is formed,
the production GOR must be equal to the solution GOR. To paraphrase more generally, when considering
the production GOR, two cases must be examined: first, the production GOR cannot be less than the solution
GOR; second, it cannot be greater than the solution GOR in the absence of free gas in the fracture.
12 SPE-188418-MS

Figure 13—Best fit of the phase-contact history with Setup 15

Figure 14—Production and pressure as calculated with Setup 15

The objective of these iterative steps was to successfully match the oil-rim thickness at the end of
the production history. Naturally, deviations from the measured values are possible. The phase contacts,
especially the oil-water contact, do not form a perfect horizontal plane in the reservoir.

Single Matrix Block Simulation – Determination of apparent matrix properties

The single matrix block analysis described above were performed to obtain a first educated guess of the
matrix oil recovery based primarily on data distributions derived from the geological model and PVT and
SCAL data. No dynamic, measured data, such as observed rates, pressures and phase contacts, have been
considered. In the subsequent dual porosity MB, those data were accounted for and thus the previous "static"
recovery curves became "dynamic" ones. The next step in the workflow will now be to improve the quality
of the geological model.
Based on the results of the previous phase contact matching procedure, recovery curves can be calculated,
using H5’s SMB tool, that match the scaled recovery curves of the dual porosity MB. Multiple combinations
of parameters can be manipulated to achieve a desired match of the recovery curves. Therefore, the
determination of the matrix properties requires an understanding of the reservoir and knowledge of which
parameters are relatively certain or uncertain. This section on the determination of SMB parameters is not
intended to provide a single correct solution. Instead, its purpose is to demonstrate that it is possible to
derive multiple sets of parameters by matching the recovery curves. Therefore, this part of the workflow is
suitable for conducting a sensitivity analysis and an uncertainty estimation.
SPE-188418-MS 13

Possible paths for matching the water drive recovery curve include adjusting the matrix permeability,
the matrix-block height and the oil-water capillary pressure curve. Manipulating the matrix permeability
by multipliers produces the expected linear relationship between permeability value and time to achieve a
certain recovery value. However, this does not influence the potential difference between the matrix and
fractures in the individual phases. Therefore, as long as permeability is greater than 0, this factor does not
influence the ultimate theoretical recovery. It must be noted that very low permeability values limit the
practical application of this methodology, and do not provide reasonable ultimate recovery estimates.
Both the matrix-block height and the oil-water capillary pressure curve have a more complex impact on
the recovery curve by influencing both the nature and the speed of the ultimate recovery. In this case, the
phase-potential difference is not only scaled by a constant factor, but the fundamental driving forces are
also manipulated (in this context, gravitational and capillary forces).
For water drive, the scaled recovery curves can be reproduced twofold with SMB calculations—either
by setting the permeability to 4 md and the matrix-block height to 14 m or by only adjusting the oil-water
imbibition capillary pressure. A graphical comparison of both these two possibilities and the un-scaled and
scaled water recovery curves is presented in Figure 15.
While matching the gas-drive recovery curve, only the period that was also used in the MB calculation
was considered (approximately 30,000 days). The results of the matching process can be seen in Figure 16.
The match has been achieved by adjusting the permeability to 0.35 md and the matrix-block height to 18 m.

Figure 15—Determination of matrix-block parameters (water drive)

14 SPE-188418-MS

Figure 16—Determination of matrix-block parameters (gas drive)

Investigation of numerical column model

The twin-barrel model, which is illustrated in Figure 10 and which was created for the dual-porosity MB
calculation, can also be utilized as a dual-continuum grid model (1D column). In this case, the block
transmissibilities must also be determined based on the parametrized (geo-cellular) model. The sum of all
vertical transmissibilities—now being connected in parallel in the column model—must be summed up. The
fracture-matrix transmissibilities must be built as an oil-in-place weighted average, and the column model
can be initialized in either equilibrium or in non-equilibrium states. In the second case, the OOIP will be
equal to that of the geocellular and MB models, although some initial instabilities will have to be endured.
In the actual field project, the column is built from 2x454 cells of a uniform thickness of 1 m. The
equilibrium initialization was used, and a single well was introduced, representing the entire field’s
production. This single well is placed within the column model’s middle 20 cells, and then opened.
Furthermore, the recovery curve scaling factors evaluated at the MB calculation were used. The target
pressure and phase method (TPPM), which was introduced by Abrahem et al. (2010), was used with the
possibility of opening pseudo-perforations at the top and the bottom of the fracture column. This ensures
an exact pressure and that the WC and GOR match.
As illustrated in Figure 17, the OWC movement is exactly the same as in the MB calculations; however,
the OGC calculated with the column model is generally in a deeper position. The main difference between
the MB and column models is the pressure utilized in each approach. The MB operates with the average
pressure, which is equal for the fracture and all matrix tranches. However, in the column model, there is
a pressure gradient in both the matrix and the fracture columns. The pressure difference between the top
and the bottom is 38 bara in the column. This results from the associated hydrostatic pressure gradient on
a large reservoir thickness.
SPE-188418-MS 15

Figure 17—Comparison between phase-contact movements as calculated by dual-porosity MB and numerical 1D DCM

Full-field simulation using recovery curves

The full-field example was also calculated using the recovery-curve method. Calculating an entire field
history in a full-field simulation is much more time-consuming than the previous steps of the workflow.
Therefore, the generation of recovery curves using SMB and the scaling of these curves using MB and
column models is a necessary prerequisite to starting this simulation.
In this setup, the wells are only controlled by their net oil-production rates. In contrast to the MB or the
column model, the full-field model is not only a dual-porosity model but also one of dual permeability. This
means that matrix cells are communicating, and a possible water or gas displacement of the oil can occur.
The same recovery curves as those for the dual-porosity MB calculation are used. The RC and timescale
factors have also been taken from the MB analysis. If necessary, the position of the phase contacts of the
full-field model could be influenced by adjusting the values of the recovery and timescale factors. When
the same factors as those in the MB calculation are used, the resulting calculated phase contacts can be
seen in Figure 18.
Again, as for the column model, TPPM was applied to achieve a satisfactory history match for pressure
and produced phases. Water cut and gas/oil ratio matches are shown in Figure 19.

Figure 18—Oil-gas-contact and oil-water-contact history matches with full-field simulation using RC
16 SPE-188418-MS

Figure 19—WC and GOR match with full-field simulation using RC

Summary and Conclusions

For the first time, the recovery-curve-based workflow provides a closed-loop numerical modelling approach
for the development planning and associated recovery factor estimation of NFRs. This numerical modelling
approach can be carried out at any time of the field’s life to perform the associated investigations of the
reservoir’s dual porosity behavior.
At any point in the field’s life, the analyses will always begin with the first workflow step of analyzing
the SMBs’ behaviors under possible recovery mechanisms, which can be performed in any number of
combinations. On the basis of a conceptual reservoir configuration, the concurrence of the fracture and
matrix block system in space and time can be examined at any level of complexity. Material balance,
column or cross-section model calculations and single well modeling make it possible to adjust the initial
starting point of a field’s analyses. Moreover, together with practical observations, they make it possible to
improve the reservoir’s associated development concept. At this point, it is possible to return to the SMB
investigation, matching fundamental properties, such as matrix block height, rock wettability and capillary
functions, before starting the loop again. All investigations can now be undertaken utilizing the same tool
and based on a common data repository. None of these features had been previously possible utilizing any
classical modeling approach.
Until now, this workflow was only exercised on a few field cases (three real-world applications and
associated data sets). For future research initiatives, supplementary efforts will be necessary to test and
generalize the recovery curve method.

Acknowledgements
The authors thank Zoltán E. Heinemann, emeritus professor at the Mining University of Leoben, for his
encouragement and suggestions, and the PHDG association for supporting this work.

Nomenclature
B formation volume factor, L3/L3, STB/res bbl
E efficiency/recovery factor
g specific gas amount, scaled to unit oil in place
ka apparent matrix permeability, L², md
kr relative permeability
p pressure, m/L³, psia
Pc capillary pressure, m/L³, psia
PL potential difference due to different pressure gradients in matrix and fracture, m/L³, psia
SPE-188418-MS 17

q production rate or flow rate, L³/t, STB/day

t time, t, day
S saturation
V volume, L³, res bbl
w specific water amount, scaled to unit oil in place

Greek symbols:
β time scaling factor
Δ difference operator
κ multiplier
λ exponential recovery constant, 1/t, 1/day
μ viscosity, m/Lt, cp
Φp phase potential, m/L³, psia
σ shape factor, 1/L2, 1/ft2

Subscripts:
a apparent
b bubblepoint
f fracture
g gas phase
j time point index
m matrix
o oil phase
p phase
R recovery
w water phase

Greek Subscripts:
V apparent

Superscripts:
g gas
r reference
w water

Conversion Factors

141.5/(131.5+°API) kg/m3
bar x 1.0* E+05 = Pa
bbl × 1.589873 E-01 =m3
cp × 1.0* E-03 = Pa.s
day × 8.64* E+04 =s
ft × 3.048* E-01 =m
18 SPE-188418-MS

ft³ × 2.831685 E-02 = m3

lbm × 4.535924 E-01 = kg
lbm/ft3 × 1.601846 E+01 = kg/m³
md × 9.869233 E-04 = m²
psi × 6.894757 E+00 = kPa
* Conversion factor is exact

References
Abrahem F. A., Heinemann Z. E. and Mittermeir G. M. 2010. A New Computer Assisted History Matching Method. Paper
SPE 130426-MS presented at the SPE EUROPEC/EAGE Annual Conference and Exhibition held in Barcelona, Spain,
14-17 June. http://dx.doi.org/10.2118/130426-MS
Aguilera, R. 1980. Naturally Fractured Reservoirs, 1st edition. p.8. Tulsa, Oklahoma: Penn Well Books
Allan, J. and Sun, S. Q. 2003. Controls on Recovery Factor in Fractured Reservoirs: Lessons Learned from 100 Fractured
Fields. Paper SPE 84590-MS presented at the SPE ATCE held in Denver, Colorado, U.S.A, 5-8 October. http://
dx.doi.org/10.2118/84590-MS.
Aronofsky, J. S., Masse, L., Natanson, S. G. 1958. A model for the mechanism of oil recovery from the porous matrix
due to water invasion in fractured reservoirs. Trans. AIME 213, 17–19. (SPE-932-G).
De Swaan, A. 1978. Theory of Waterflooding in Fractured Reservoirs. SPE J. 18 (2): 117–122. SPE-5892-PA. http://
dx.doi.org/10.2118/5892-PA
Di Donato, G., Lu, H., Tavassoli, Z., Blunt, M. J. 2007. Multirate-transfer dual-porosity modeling of gravity drainage and
imbibition. SPE J. 12 (1), 77–88. http://dx.doi.org/10.2118/93144-PA.
Gharsalla M. M. 2015. Application of the Recovery Curve Method to Petroleum Reservoir Material Balance. PhD Thesis,
Montanuniversität Leoben, Leoben, Austria (April 2015).
Gilman, J. R., Kazemi, H.. 1988. Improved calculations for viscous and gravity displacement in matrix blocks in dual-
porosity simulators. SPE J. Pet. Technol. 40 (1) http://dx.doi.org/10.2118/16010-PA.
Heinemann, Z. E. 2004. Using Recovery Curves in Modeling Natural Fractured Hydrocarbon Reservoirs. Proposal for a
PhD Research Project at the Montanuniversität Leoben, Austria (March 2004)
Heinemann, Z. E. and Mittermeir, G. M. 2012. Derivation of the Kazemi-Gilman-Elsharkawy Generalized Dual Porosity
Shape Factor. Transport in Porous Media, 91 (1) 123–132. http://dx.doi.org/10.1007/s11242-011-9836-4.
Heinemann, Z. E. and Mittermeir, G. M. 2016. Generally Applicable Method For Calculation Of The Matrix-fracture
Fluid Transfer Rates. Presented at the SPE Europec, Vienna, Austria, 30 May - 2 June. SPE-180121-MS. http://
dx.doi.org/10.2118/180121-MS
Heinemann, Z. E. 2017. H5 Reservoir Simulator Manual, Version 2017.2, Heinemann Consulting GmbH, Leoben, Austria.
Kazemi H., Merrill L. S., Porterfield K. L., Zeman P. R. 1976. Numerical Simulation of Water-Oil Flow in Naturally
Fractured Reservoirs, SPE. J. 16 (6) 317–326. http://dx.doi.org/10.2118/5719-PA
Kazemi, H., Gilman, J. R., Elsharkawy, A. M. 1992. Analytical and numerical solution of oil recovery from fractured
reservoirs using empirical transfer functions. SPE Res. Eng. 7 (2), 219–227. http://dx.doi.org/10.2118/19849-PA.
Lu, H., Donato, G., Blunt, M. J. 2008. General transfer functions for multiphase flow in fractured reservoirs. SPE J. 13
(3), 289–297. http://dx.doi.org/10.2118/102542-PA.
Mattax, C. C., Kyte, J. R. 1962. Imbibition oil recovery from fractured, water-drive reservoirs. SPE J. 2 (2), 177–184.
http://dx.doi.org/10.2118/187-PA.
Mittermeir, G. M. 2015. Material-Balance Method for Dual-Porosity reservoirs With Recovery Curved to
Model the Matrix/Fracture Transfer, SPE Reservoir Evaluation & Engineering, 18 (2), 171–186. http://
dx.doi.org/10.2118/174082-PA.
Mittermeir, G. M., Pichelbauer, J. and Heinemann, Z. E. 2004. Automated Determination of Aquifer Properties from Field
Production Data, Presented at the 9th European Conference on Mathematics of Oil Recovery (ECMOR IX), Cannes,
France, 30 August-2 September
Nelson, R. A. 1992. An Approach to Evaluating Fractured Reservoirs. J Pet Technol 34 (9): 2167–2170. SPE-10331-PA.
http://dx.doi.org/10.2118/10331-PA
Quandalle, P., Sabathier, J. C., 1989. Typical features of a multipurpose reservoir simulator. SPE Res. Eng. 4 (4), 475–480.
http://dx.doi.org/10.2118/16007-PA.
SPE-188418-MS 19

Pirker, B. Mittermeir, G. M., Heinemann, Z. E., 2007. Numerically Derived Type Curves for Assessing Matrix Recovery
Factors. Presented at In: Proceedings of the EUROPEC/EAGE Conference and Exhibition, London, U.K., 11-14 June.
SPE-107074-MS. http://dx.doi.org/10.2118/107074-MS.
Reiss L. H., Bossie Codreanu D., Lefebvre du Prey E. J. 1973. Flow in Fissured Reservoirs, Paper SPE 4343-MS presented
at the Second Annual European Meeting of AIME, London, England, April 2-3. https://dx.doi.org/10.2118/4343-MS
Rossen, R. H., Shen, E. L. 1987. Simulation of gas/oil drainage and Water/oil imbibition in naturally-fractured reservoirs.
SPE Reserv. Eng. 4 (4), 464–470. http://dx.doi.org/10.2118/16982-PA.
Steiner C. F. and Mittermeir G. M. 2017. Applicability of the shape-factor concept for naturally fractured reservoirs
and an alternative approach, Journal of Petroleum Science and Engineering, 154, 60–75, http://dx.doi.org/10.1016/
j.petrol.2017.04.009.
Warren, J. E., Root, P. J. 1963. The behavior of naturally fractured reservoirs. SPE J. 3 (3), 245–255. http://
dx.doi.org/10.2118/426-PA.