3-D
streamline simulation is widely considered as one of the major
developments in
petroleum reservoir simulation and performance forecasting in the last
decade.
The technology has been rapidly assimilated by the industry for highly
detailed
flow simulation, reservoir management, history matching and uncertainty
assessment. Several commercial reservoir simulators have been developed
based
on the streamline ‘time of flight’ concept introduced by
Dr. Datta-Gupta and
King (1995). With the advancement in high resolution data acquisition
and
seismic technolgies, geologic models now routinely consist of
multimillion
cells. This resulted in a widening gap between geologic modeling, flow
simulation and uncertainty assessement. Streamline simulation has
effectively bridged
this gap. Our work developing streamline flow simulation and history
matching
of high resolution geologic models has been pivotal to these
developments.
·Fast Marching
Methods: Application to Unconventional Reservoirs
The
concept of depth of investigation is fundamental to well test analysis
and is
routinely used to design well tests and to understand the reservoir
volume
investigated. The depth of investigation can also be useful in
identifying new well
locations and planning, designing and optimizing hydraulic fractures in
unconventional reservoirs. One way to calculate the depth of
investigation is
through numerical simulation. The numerical approach is very general
and can be
applied to arbitrary reservoir and well conditions. The computation
time and
expenses, however, make the numerical simulation approach unfeasible
for
routine applications.
As
an alternative to numerical simulation, we can directly solve for the
propagation equation for the pressure ‘front’ defined as
the maximum of the
pressure response for an impulse source.Such
propagation equation can be derived using asymptotic
ray theory
which has been used extensively in electromagnetic and seismic wave
propagation
(Virieux et al., 1994). The asymptotic method draws upon an analogy
between the
propagating pressure ‘front’ and a propagating wave front,
and many of the
concepts such as rays and propagating fronts have their counterparts in
Petroleum Engineering (Datta-Gupta and King, 2007). Specifically,
a high frequency asymptotic solution of the diffusivity equation leads
to the
following equation for a propagating pressure ‘front’ for
an impulse source or
sink (Vasco et al., 2000; Datta-Gupta and King, 2007):
In
the above equation, is the
diffusivity (consisting of =
permeability, =
porosity, µ = fluid viscosity, ct
= total compressibility) which can be function of position andis
the propagation time of the pressure ‘front’ (also
called a ‘diffusive time of flight’). It simply tells us
that the pressure
‘front’ propagates in the reservoir with a velocity given
by the square root of
diffusivity. This propagation depends on reservoir and fluid properties
and is
independent of flow rate. Also, has
unit of square root of time which is consistent
with the scaling behavior of pressure diffusion. In fact, is
related to physical time through a simple
expression of the form where
the pre-factor c depends on the specific
flow geometry (Kim et al., 2009). For example, for radial flow, c
=
1/4, for linear flow, c = 1/2 and for spherical flow, c
=
1/6.
The
pressure ‘front’ equation is a form of the Eikonal equation
which appears in
many contexts such as elastic and electromagnetic wave propagation; its
properties are well developed in the literature (Sethian, 1999). Most
importantly, the Eikonal equation can be solved very efficiently (in
seconds
opposed to hours of simulation time for comparable problems) by a class
of
solution called Fast Marching Methods (FMM) (Sethian, 1999). Following
are
three examples of depth of investigation: 1) generalization to
heterogeneous
medium; 2) generalization to anisotropic permeability and unstructured
grid;
and 3) a horizontal well with multistage hydraulic fractures.
·Rapid Compositional Flow Simulation in Shale
Oil Reservoirs using Fast Marching
Method
Reliable performance assessment
of unconventional reservoirs
requires accurate modeling of inter-porosity flow
characteristics in hydraulic fractures, microfractures and reservoir
rock. In addition, phase behavior in nano-porous
rocks plays an important role in reservoir performance. High resolution
flow simulation incorporating the complex
underlying physics and detailed reservoir heterogeneity is
computationally expensive. Our proposed simulation approach is
based on a high frequency asymptotic solution of the diffusivity
equation in
heterogeneous and fractured reservoirs. The high frequency solution
leads to the Eikonal equation which is solved
for a ‘diffusive time of flight’ (DTOF) that governs the
propagation of the ‘pressure front’ in the reservoir. Our
approach consists of two decoupled steps: calculation of the DTOF using
the Fast Marching Method (FMM) and
fully-implicit compositional simulation using DTOF as a spatial
coordinate. The computational efficiency is
achieved by reducing the 3-D compositional flow equation into 1-D
equation using the DTOF as spatial coordinate,
leading to orders of magnitude faster computation over full 3-D
compositional simulation. The savings in
computation time increases significantly with grid refinement and for
high resolution models.
·Data-Driven Fracture Flow Diagnostics and
Characterization
Multistage
hydraulically fractured horizontal
wells provide an effective means to exploit unconventional reservoirs.
The
current industry practice in the interpretation of field response often
utilizes empirical decline curve analysis or pressure/rate transient
analysis
(PTA/RTA) for the characterization of these reservoirs and fractures.
These
analytical tools are based on many simplifying assumptions and cannot
provide a
detailed description of the evolving reservoir drainage volume from a
well, the
understanding of which is essential for unconventional reservoir and
fracture
assessment and optimization.We propose a novel
“data-driven” methodology for the production analysis
of shale gas and shale oil reservoirs. The approach is based on the
high
frequency asymptotic solution of the diffusivity equation in
heterogeneous
reservoirs. It allows us to determine the drainage volume from a well,
and the
instantaneous recovery ratio (IRR), which is defined as the ratio of
the
produced volume to the drainage volume, directly from the production
data. In
addition, in a manner analogous to the diagnostic plot in PTA, a new
w(τ) plot
has been proposed to provide better insight into the depletion
mechanisms and
the fracture geometry. The major advantages of this
current approach are the model free analysis
without presumptions of flow regimes, and a simple and intuitive
understanding
of the drainage volume and fracture conductivity. The results of the
analysis
are useful for well and hydraulic fracturing operation design
optimization and
matrix and fracture parameter estimation.
·Flow
Simulation of
Complex Fracture Systems with Unstructured Grids using Fast Marching
Method
With the current industry
practice of reduced cluster spacing
and increased fracturing proppant/fluid volume, the
hydraulic fracture treatment tends to generate more complex fracture
systems, where an unstructured computational
grid, instead of a Cartesian or corner point grid, is preferred to
accurately model the fracture geometry. With
unstructured grids, the reservoir performance is generally simulated
with finite volume simulation, for which one
major issue is the potentially heavy computational cost. A novel
approach has recently been introduced to provide a
rapid simulation of unconventional reservoirs, which first captures the
drainage volume during the transient pressure
propagation process using the Fast Marching Method (FMM) and then
rapidly solves fluid flow equation in an
equivalent 1D domain. However, this application is currently limited to
calculating the reservoir response with
Cartesian or corner-point grids.
·Multiscale Parameterization
and Model Calibration with Grid Connectivity Based Transforms (GCT)
The
motivation of this topic is to expand upon the current state of the art
in
parameterization by linear transformation for history matching. In
addition to
theoretical developments, an emphasis is placed on the development of
parameterization methods that are applicable to large reservoir models,
on par
with current industry standards and computational modeling
capabilities, and
also varying levels of prior geologic complexity.
Methods
of linear transform dependent on prior model assumptions have commonly
utilized
the parameter covariance matrix. The Karhunen-Loeve transform (KLT) or
principle component analysis (PCA) has been widely used [Gavalas et al.,
1976; Karhunen, 1947; Loéve, 1978; Reynolds
et al., 1996; Li
and Cirpka, 2006; Ma et al., 2008; Jafarpour and
McLaughlin,
2009]. In this approach the estimable property defined at any
grid cell
is represented as the linear expansion of the weighted eigenvectors of
the
property covariance matrix. The eigenvectors form the transform basis
and are
typically ranked by their corresponding eigenvalues, from largest to
smallest,
that are related to the variance contribution of each eigenvector to
the total
parameter variance. The expansion is optimal in the mean-squared-error
sense
among the linear class of transformations; therefore, the KLT
coefficients
present the fewest number of parameters that capture the maximum amount
of
variation for any low-rank approximation. Accordingly, the sorted
eigenvectors
convey larger to smaller scales of spatial parameter variation. A
limitation of
this approach is that in realistic problems with high-resolution models, the covariance can be unknown (or
uncertain),
resulting in misleading basis functions when prior model assumptions
are
incorrect [Jafarpour and
McLaughlin, 2009].
Eigendecomposition of the covariance
matrix, required for each parameter update, is also prohibitively
expensive for
high-resolution models relative to current computational capability.
Last, the
preservation of moments beyond second-order (i.e., covariance), as well
as
complex continuous geologic structures (e.g., channel systems), are not
guaranteed in the KL expansion during the parameter updates.
A
new model-independent basis constructed from grid-connectivity
information is
developed here. The parameterization is applicable to any grid
structure and
domain geometry. The development of the grid-connectivity-based
transform (GCT)
basis begins from first principles, merging discrete Fourier analysis
and
spectral graph theory. The basis functions represent the modal shapes
or
harmonics of the grid, are defined by a modal frequency, and converge
to the
discrete cosine transforms (DCT) for certain grid geometries and
boundary
assumptions; therefore, reservoir model calibration is performed in the
spectral domain. Using an adaptive multiscale workflow, the GCT
parameterization
is successfully applied for history matching of several synthetic
reservoir
models of varying geometry and geologic complexity and also of a field
case.
·A Hierarchical
Multi-scale Approach to History Matching
We
develop a hierarchical multiscale history matching
approach with two levels of hierarchy for inversion in parameter space
and
geometric space. The approach follows a sequence of calibrations from
global to
local parameters in coarsened and fine scales. At first, we identify
the ‘heavy
hitters’ in the large scale static and dynamic parameters and
then calibrate
them using an evolutionary algorithm. This global parameter
calibration,
matching and balancing field level energy, is followed by
streamline-assisted
multiscale inversion to match well by well production history with
local
parameter updates. Starting with a relatively coarse grid, we match the
production data at wells by dynamically refining the reservoir grid.
This
multiscale data integration offers not only benefits of enhancement in
computation efficiency and but also effective iterative minimization. The
approach has been successfully applied to a variety of field cases in
collaboration with operating companies.
·Robust
Streamline Tracing in Embedded Discrete Fracture Models
The streamline based technology
has proven to be effective
for various subsurface flow and transport modeling problems including
reservoir
simulation, model calibration and optimization. For naturally fractured
systems, current streamline models are well suited for dual porosity
single
permeability systems because streamlines exist need to be traced only
for the
fracture system. However, complications arise for dual porosity dual
permeability (DPDP) systems because streamlines need to be traced for
both
fracture and matrix systems. Also, the streamlines in the two systems
may
interact. We present a robust streamline tracing framework for use in
the DPDP
models via an embedded discrete fracture model (EDFM) framework.
·Flood Rate Optimization Using Streamlines
It
is well recognized that field-scale rate optimization problems often
involves
highly complex reservoir model, production and facility related
constraints and
geological uncertainty. Deployment of smart well completions with
inflow
control valves (ICV) to control production/injection rates for various
segments
along the wellbore further compounds the complexity of the problem. All
these
make optimal reservoir management via rate control difficult without
efficient
optimization algorithms. Alhuthali et
al.
(2007, 2008, 2010)
proposed a streamline-based approach to maximize waterflood sweep
efficiency.
The approach focuses on equalizing arrival times of the waterfront at
all
producers within selected sub-regions of a waterflood project. This
results in
delayed water breakthrough and reduces field water cut after water
breakthrough.
The optimization was performed under production and facility
constraints using
a sequential quadratic programming approach. A major advantage of the
streamline-based approach is the analytical computation of the
sensitivities of
the waterfront arrival times at the producers with respect to well
production/injection rates and also, analytic calculations of the
gradient and
Hessian of the objective function. This makes it computationally
efficient and
suitable for large field cases. It also accounts for geological
uncertainty via
a stochastic optimization framework based on a combination of the
expected
value and variance of a performance measure from multiple realizations.
We have
also studied the tradeoff between maximizing sweep efficiency and
accelerating
production (maximizing NPV) using an augmented objective function. In
addition
to waterflooding, we have applied this method to EOR processes like
polymer and
CO2 flooding.
·Diffuse Source Upscaling of Geologic Models
The
motivation for Diffuse Source (DS) upscaling is
to address the assumption of pressure equilibrium while coarsening. The
assumption that weakly connected and disconnected pay would have the
same
equilibrium pressure as the well-connected pay inside a coarse block is
error
prone. Existing flow based upscaling methods are based on Steady State
(SS)
concepts of flow and cannot distinguish between well-connected and
weakly
connected pay in the upscaling region. The current work looks at the
use of
diffusive time of flight and drainage volume concepts to build basis
functions
that effectively capture the subgrid heterogeneity that is
well-connected
especially in high contrast media. The application of so called Diffuse
Source
basis functions to set up local flow problems leads to non-uniform
pressure
gradients even for homogeneous systems. This forced us to revisit the
construction of transmissibility. A new construction of
transmissibility based
on Pseudo Steady State (PSS) flow concept is utilized which does not
require
the explicit specification of cell centers, instead the total PSS
pressure drop
across the volume is used as the reference pressure drop for the
transmissibility construction. The DS formulation is an extension of
PSS concepts
to high contrast media through the use of diffusive time of flight.
Figure
1
shows the pressure map comparison of the fine scale and a vertically
coarsened
case for the 7 million cell Amellago carbonate geologic model [Amellago
carbonate model provided by Dr. Sebastian Geiger and the
International
Centre for
Carbonate Reservoirs, Institute of Petroleum Engineering, Heriot-Watt
University]. The coarsening design is generated from the statistical
layering
algorithm in SWIFT giving a factor of 6 coarsening in this case and the
upscaling of properties is done by the Diffuse Source algorithm.
Figure 2
shows the water saturation maps for the same and in
Figure 3,
we have the quantitative comparison of cumulative oil volumes for
Steady State
and Diffuse Source upscaling algorithms. The field water cut comparison
is
given in
Figure 4.
Figure 1: Pressure comparison
between fine scale (left) and 1x1xN coarse models (right).
Figure 2: Water saturation
comparison between fine scale (left) and 1x1xN coarse models (right).
Figure 3: Field water cut vs. Pore Volume
Injected
Figure 4: Field oil preoduction vs. Pore Volume
Injected
