U.S. patent application number 14/880404 was filed with the patent office on 2016-02-04 for estimating petrophysical parameters and invasion profile using joint induction and pressure data inversion approach.
The applicant listed for this patent is SCHLUMBERGER TECHNOLOGY CORPORATION. Invention is credited to Aria Abubakar, Tarek Habashy, Lin Liang, Michael Thambynayagam.
Application Number | 20160033673 14/880404 |
Document ID | / |
Family ID | 42337609 |
Filed Date | 2016-02-04 |
United States Patent
Application |
20160033673 |
Kind Code |
A1 |
Liang; Lin ; et al. |
February 4, 2016 |
ESTIMATING PETROPHYSICAL PARAMETERS AND INVASION PROFILE USING
JOINT INDUCTION AND PRESSURE DATA INVERSION APPROACH
Abstract
Methods and related systems are described relating to an
inversion approach for interpreting the geophysical electromagnetic
data. The inversion can be constrained by using a multiphase fluid
flow simulator (incorporating pressure data if available) which
simulates the fluid flow process and calculates the spatial
distribution of the water saturation and the salt concentration,
which are in turn transformed into the formation conductivity using
a resistivity-saturation formula. In this way, the inverted
invasion profile is consistent with the fluid flow physics and
moreover accounts for gravity segregation effects. Jointly with the
pressure data, the inversion estimates a parametric one-dimensional
distribution of permeability and porosity. The fluid flow volume is
directly inverted from the fluid-flow-constrained inversion of the
electromagnetic data. The approach is not limited by the
traditional interpretation of the formation test, which is based on
a single-phase model without taking into account invasion or
assuming that the fluid, for example mud-filtrate, has been cleaned
up from the formation testing zone. The joint inversion of the
electromagnetic and pressure data provides for a more reliable
interpretation of formation permeability. One advantage of the
approaches described herein, is its possible generalization to
three-dimensional geometries, for example dipping beds and highly
deviated wells.
Inventors: |
Liang; Lin; (Belmont,
MA) ; Abubakar; Aria; (Sugar Land, TX) ;
Habashy; Tarek; (Burlington, MA) ; Thambynayagam;
Michael; (Sugar Land, TX) |
|
Applicant: |
Name |
City |
State |
Country |
Type |
SCHLUMBERGER TECHNOLOGY CORPORATION |
Sugar Land |
TX |
US |
|
|
Family ID: |
42337609 |
Appl. No.: |
14/880404 |
Filed: |
October 12, 2015 |
Related U.S. Patent Documents
|
|
|
|
|
|
Application
Number |
Filing Date |
Patent Number |
|
|
12607708 |
Oct 28, 2009 |
9176252 |
|
|
14880404 |
|
|
|
|
61145658 |
Jan 19, 2009 |
|
|
|
Current U.S.
Class: |
702/7 |
Current CPC
Class: |
G01V 3/28 20130101; G01V
3/38 20130101; G01N 15/088 20130101; Y02A 90/30 20180101; Y02A
90/344 20180101 |
International
Class: |
G01V 3/38 20060101
G01V003/38; G01N 15/08 20060101 G01N015/08; G01V 3/28 20060101
G01V003/28 |
Claims
1. A method of estimating porosity for a subterranean formation of
interest including one or more geologic beds, the method
comprising: receiving electromagnetic survey data of the formation
of interest; receiving fluid flow characteristics within the
formation of interest; and estimating porosity within at least a
portion of the formation of interest based at least in part on the
electromagnetic survey data and the fluid flow characteristics,
wherein the estimation has a spatial resolution such that at least
two porosity values are estimated for each of the one or more
geologic beds.
2. A method according to claim 1 wherein the estimated porosity has
a spatial resolution such that at least two porosity values are
estimated within each predefined depth interval.
3. A method according to claim 1 further comprising directly
inverting the electromagnetic data and the fluid flow
characteristics, wherein the estimated porosity is based at least
in part on the direct inversion.
4. A method according to claim 3 further comprising estimating a
mud filtrate invasion rate based at least in part on the
inversion.
5. A method according to claim 3 further comprising estimating a
well completion fluid invasion rate based at least in part on the
inversion.
6. A method according to claim 3 further comprising estimating a
fluid injection rate based at least in part on the inversion.
7. A method according to claim 3 wherein the fluid flow
characteristics are generated from a fluid flow simulation.
8. A method according to claim 3 wherein the fluid flow
characteristics are generated from pressure transient measurements
made in the borehole.
9. A method according to claim 8 further comprising estimating
permeability within the formation based at least in part on the
inversion.
10. A method according to claim 3 further comprising generating a
quantitative predictive formation model for at least part of the
formation of interest based at least in part on the inversion.
11. A method according to claim 10 further comprising predicting
future fluid flow in at least part of the formation of interest,
the prediction being based at least in part on the generated
quantitative predictive formation model.
12. A method according to claim 1 wherein the electromagnetic
survey data is obtained using at least a first downhole
electromagnetic transceiver deployed in a borehole passing through
the formation of interest.
13. A system of estimating porosity for a subterranean formation of
interest including one or more geologic beds, the system comprising
a processing system adapted and configured to receive
electromagnetic survey data of the formation of interest and fluid
flow characteristics within the formation of interest and to
estimate porosity within at least a portion of the formation of
interest based at least in part on the electromagnetic survey data
and the fluid flow characteristics, wherein the estimation has a
spatial resolution such that at least two porosity values are
estimated for each of the one or more geologic beds.
14. A system according to claim 13 wherein the processing system is
further adapted and configured to directly invert the
electromagnetic data and the fluid flow characteristics, wherein
the estimated porosity is based at least in part on the direct
inversion.
15. A system according to claim 14 wherein the processing system is
further adapted and configured to generate a quantitative
predictive formation model for at least part of the formation of
interest based at least in part on the inversion.
16. A system according to claim 13 further comprising a downhole
electromagnetic transceiver deployable in a borehole passing
through the formation of interest.
Description
CROSS REFERENCE TO RELATED APPLICATION
[0001] This patent application is a divisional application of U.S.
patent application Ser. No. 12/607,708, filed Oct. 28, 2009, which
claims the benefit of U.S. Provisional Patent Application Ser. No.
61/145,658, filed Jan. 19, 2009, both of which are hereby
incorporated by reference in their entireties.
FIELD OF THE INVENTION
[0002] This invention relates broadly to the investigation of
geological formations traversed by a borehole. More particularly,
this invention relates to a method for integrating electromagnetics
and multi-phase fluid flow in order to provide a robust, physically
consistent interpretation for reservoir characterization.
BACKGROUND OF THE INVENTION
[0003] Induction logging measurements are sensitive to water
saturation and brine concentration in the rock pores. In an active
reservoir, the formation's petrophysical parameters can have a
strong imprint on the temporal and spatial distribution of water
saturation and salt concentration, which in turn can be transformed
into the distribution of formation conductivity using an
appropriate saturation-resistivity equation. This relationship
between induction measurements and the petrophysical parameters
offers an opportunity to integrate electromagnetics and multi-phase
fluid flow to provide a robust, physically consistent
interpretation for reservoir characterization.
[0004] Induction logging tools can be used to determine the
formation resistivity and invasion profile via a model-based
inversion approach. See, Habashy, T. M. and A. Abubakar, 2004, A
general framework for constraint minimization for the inversion of
electromagnetic measurements: Progress In Electromagnetics
Research, 46, 265-312 (hereinafter "Habashy and Abubakar"); Barber,
T., B. Anderson, A. Abubakar, T. Broussard, K. Chen, S.
Davydycheva, V. Druskin, T. Habashy, D. Homan, G. Minerbo, R.
Rosthal, R. Schlein, and H. Wang, 2004, Determining formation
resisitivty anisotropy in the presence of invasion, in SPE Annual
Technical Conference and Exhibition, paper 90526; Abubakar, A., T.
M. Habashy, V. Druskin, L. Knizhnerman, and S. Davydycheva, 2006, A
3d parametric inversion algorithm for triaxial induction data:
Geophysics, 71, G1-G9; and Abubakar, A., T. M. Habashy, V. Druskin,
and L. Knizhnerman, 2006, An enhanced gauss-newton inversion
algorithm using a dual-optimal grid approach: IEEE Transactions on
Geoscience and Remote Sensing, 44, 1419-1427. However, this
application is constrained by some simplifying assumptions, for
example, assuming a well penetrating a layer-cake formation with a
step-profile three-parameter (R.sub.xo, R.sub.t, D.sub.i) invasion
model. The quality and accuracy of the inversion results are
affected by the complexity of the reservoir configuration and the
actual invasion profile. In horizontal or highly deviated wells or
in anisotropic formations, invasion profiles become too complex to
be described by a simple invasion model. An alternative inversion
approach is to employ the so-called pixel-based inversion method.
See: Abubakar, A. and P. M. van den Berg, 2000, Nonlinear inversion
in electrode logging in a highly deviated formation with invasion
using an oblique coordinate system: IEEE Transactions on Geoscience
and Remote Sensing, 38, 25-38; Alumbaugh, D. and M. Wilt, 2001, A
numerical sensitivity study of three-dimensional imaging from a
single borehole: Petrophysics, 42, 19-31; Abubakar, A. and P. van
den Berg, 2002, Application of a non-orthogonal coordinate system
to inverse single-well electromagnetic logging problem:
Three-Dimensional Electromagnetics, eds. M. S. Zhdanov and P. E.
Wannamaker, 215-231; Abubakar, A. and T. Habashy, 2006,
Three-dimensional single-well imaging of the multi-array triaxial
induction logging data, in SEG Technical Program Expanded
Abstracts, volume 25, 411-415. A disadvantage of this pixel-based
inversion method is that the number of unknowns used to describe
the configuration is large and hence, a large number of unknown
model parameters often need to be inverted. This approach will only
provide a qualitative inverted resistivity map (an image) around
the well-bore.
[0005] On the other hand, the mud-filtrate invasion affects the
interpretation of well testing or wireline formation test. Although
some researchers studied and discussed the influence of invasion to
the formation test interpretation, practically, the interpretation
of formation test usually still employs a single-phase fluid flow
model regardless of the mud-filtrate invasion. The error associates
with this can be neglected only when the mobility of the
mud-filtrate is very close to that of the formation fluid and the
capillary pressure is negligible or the invaded mud-filtrate has
already been cleaned up from the test target zone by a drawdown,
which may be time-consuming. A more realistic invasion profile will
benefit both induction data interpretation and formation test
interpretation.
[0006] Some previous works have already exploited the benefit of
integrating the electromagnetic model with the multiphase fluid
flow model. In Semmelbeck, M. E. and S. A. Holditch, 1988, The
effects of mud-filtrate invasion on the interpretation of induction
logs: SPE (hereinafter "Semmelbeck and Holditch (1988)"), an
attempt is made to develop a model to simulate the mud-filtrate
invasion process, which deals with the mud-cake and the formation
simultaneously, and incorporated the capability to simulate the
salt transport to calculate the formation resistivity. This model
was then used to analyze the effect of the mud-cake permeability,
the formation permeability, the porosity and the overbalance
pressure on the induction logging data. In Tobola, D. P. and S. A.
Holditch, 1991, Determination of reservoir permeability from
repeated induction logging: SPE Formation Evaluation, 20-26
(hereinafter "Tobola and Holditch (1991)"), the approach described
in Semmelbeck and Holditch (1988) was applied to a low permeability
reservoir and determined the reservoir permeability by history
matching the change in the induction logging data over time. The
effect of the initial water saturation on the induction logging
data was also analyzed. In this approach, the mud-cake permeability
profile over time must be accurately simulated to properly
interpret the induction logging data. In Yao, C. Y. and S. A.
Holditch, 1996, Reservoir permeability estimation from time-lapse
log data: SPE Formation Evaluation, 69-74 (hereinafter "Yao and
Holditch (1996)"), the approach described in Semmelbeck and
Holditch (1988) was extended to history match both the time-lapse
resistivity and neutron logging data. This technique was used to
estimate the reservoir permeability, which was verified against
production data and core analysis.
[0007] The approaches described in Semmelbeck and Holditch (1988),
Tobola and Holditch (1991) and Yao and Holditch (1996) simulate the
mud-filtrate invasion assuming a mud-cake with a constant thickness
and a variable permeability, which need to be accurately determined
from other independent measurement data. See: Ferguson, C. K. and
J. A. Klotz, 1954, Filtration from mud during drilling: Trans.,
AIME 201, 29-42; and Williams, M. and G. E. Cannon, 1938,
Evaluation of filtration properties of drilling muds: The Oil
Weekly, 25-32. These approaches ignore gravity and diffusion
effects, and aim to estimate the permeability in tight formations.
However, this estimation heavily depends on the precision of
determining mud-cake properties, which is another difficult problem
to be resolved.
[0008] Semmelbeck, M. E., J. T. Dewan, and S. A. Holditch, 1995,
Invasion-based method for estimating permeability from logs: SPE
30581 presented at the 1995 SPE Annual Technical Conference and
Exhibition, Dallas, 22-25 describes an extension of their previous
work to integrate the fluid flow model with the Dewan's mudcake
growth model, which is an experimentally-verified model for
predicting the mud-cake thickness and the permeability during
static and dynamic filtration conditions. See: Dewan, J. T. and M.
E. Chenevert, 1993, Mudcake buildup and invasion in low
permeability formations: application to permeability determination
by measurement while drilling: SPWLA 34th Annual Logging Symposium,
13-16. This approach was applied to the new generation of array
induction logging data. From a linear covariance analysis, they
concluded that their method could provide reasonable estimates of
the permeability in low to moderate permeable formations, but the
method could not provide accurate estimates of the
saturation-dependent properties.
[0009] Ramakrishnan, T. S. and D. J. Wilkinson, 1997, Formation
producibility and fractional flow curves from radial resistivity
variation caused by drilling fluid invasion: Physics of Fluids, 9,
833-844, and Ramakrishnan, T. S. and D. J. Wilkinson, Water-cut and
fractional flow logs from array-induction measurements: SPE
Reservoir Evaluation & Engineering, 2, 85-94 discuss
establishing a mathematical model to estimate the fractional flow
and the relative permeability curves from the array induction
logging measurements. Based on the work described in the two
Ramakrishnan papers, Zeybek, M., T. Ramakrishnan, S. AI-Otaibi, S.
Salamy, and F. Kuchuk, 2001, Estimating multiphase flow properties
using pressure and flowline water-cut data from dual packer
formation tester interval tests and openhole array resistivity
measurements: SPE 71568 discusses a methodology for estimating the
horizontal and vertical layer permeabilities and the relative
permeabilities using array induction logging tool measurements,
pressure transient measurements and water-cut measurements from a
packer-probe wireline formation tester. This joint inversion was
done in a sequential mode. In this method, fractional flow
parameters (connate water saturation, residual water saturation,
maximum residual oil saturation, filtrate loss and pore size
distribution) are estimated by matching resistivity measurements.
These estimated fractional flow parameters are then input into a
numerical model for the sampling process to ultimately match the
observed and simulated water cut and pressure data. The final
outputs are the relative, the horizontal and vertical
permeabilities.
[0010] Alpak, F. O., T. M. Habashy, C. Torres-Verdin, and V.
Dussan, 2004, Joint inversion of pressure and time-lapse
electromagnetic logging measurements: Petrophysics, 45, 251-267
(hereinafter "Alpak, et al. (2004)"); Alpak, F. O., C.
Torres-Verdin, and T. M. Habashy, 2004, Joint inversion of pressure
and dc resistivity measurements acquired with in-situ permanent
sensors: a numerical study: Geophysics, 69, 1173-1191; Alpak, F.
O., C. Torres-Verdin, and T. M. Habashy 2006, Petrophysical
inversion of borehole array-induction logs: Part I--numerical
examples: Geophysics, 71, F101-F119; and Alpak, F. O., C.
Torres-Verdin, T. M. Habashy, and K. Sephernoori, 2004,
Simultaneous estimation of in-situ multiphase petrophysical
properties of rock formations from wireline formation tester and
induction logging measurements: SPE 90960 all discuss the
simultaneous estimation of in-situ multiphase petrophysical
properties of rock formations from wireline formation tester and
induction logging measurements. The papers also discuss extensive
work to assess the sensitivity of array induction measurements to
the presence of the water-based mud-filtrate invasion for various
combinations of petrophysical parameters and fluid properties.
Torres-Verdin, C., F. O. Alpak, and T. M. Habashy, 2006,
Petrophysical inversion of borehole array-induction log: Part
II--field data examples: Geophysics, 71, G261-G268; and Salazar, J.
M., C. Torres-Verdin, F. O. Alpak, T. M. Habashy, and J. D. Klein,
2006, Estimation of permeability from borehole array induction
measurements: Application to the petrophysical appraisal of tight
gas sands: Petrophysics, 47, 527-544 discuss applying method
described by Alpak to some low or moderate permeability gas
reservoirs for estimating the permeability and the porosity. This
method employs a modified UTCHEM code (INVADE) (see Wu, J., C.
Torres-Verdin, K. Sepehrnoori, and M. Proett, 2005, The influence
of water-base mud properties and petrophysical parameters on
mudcake growth, filtrate invasion, and formation pressure:
Petrophysics, 46, 14-32) to calculate the time-dependent
mud-filtrate loss rate, and takes the average of this rate as the
injection rate to run a two-phase three-component fluid flow model
for simulating the mud-filtrate invasion process, which
consequently affects the resistivity measurements. Therefore, the
accuracy of the inverted petrophysical parameters can depend on the
mud-filtrate loss calculation, which is still a difficult problem
at present. Angeles, R., J. Skolnakorn, F. Antonsen, A. Chandler,
and C. Torres-Verdin, 2008, Advantages in joint-inversions of
resistivity and formation-tester measurements: Examples from a
norwegian field: SPWLA 49th Annual Logging Symposium, Edinburgh,
Scotland discusses advancing previous work by using a reservoir
simulator dynamically-coupled to a mud-filtrate invasion model.
Hence they integrated the calculation of mud-filtrate invasion rate
into the inversion process.
[0011] Pereira, N., R. Altman, J. Rasmus, and J. Oliveira, 2008,
Estimation of permeability and permeability anisotropy in
horizontal wells through numerical simulation of mud filtrate
invasion: Rio Oil & Gas Expo and Conference, Rio de Janeiro.
discusses an approach to estimate formation permeability and
permeability anisotropy using dual inversion of time-lapse
azimuthal laterolog-while-drilling and near well-bore numerical
simulation of the water-based mud-filtrate invasion. This method
requires a priori bottom hole pressure while drilling and assumes a
known constant mud-cake permeability. Kuchuk, F., L. Zhan, S. M.
Ma, A. M. Al-Shahri, T. S. Ramakrishnan, B. Altundas, M. Zeybek, R.
Loubens, and N. Chugunov, 2008, Determination of in-situ two-phase
flow properties through downhole fluid movement monitoring: SPE
116068 presented at the 2008 SPE Annual Technical Conference and
Exhibition held in Denver, Colo., USA, 21-24 discussses method for
in-situ estimation of two-phase transport properties of porous
media using the time-lapse DC resistivity data, pressure and flow
rate data. The time-lapse DC resistivity are recorded using a
permanent downhole electrode resistivity array. The pressure and
flow rate data are obtained from the injection test regardless of
the mud-filtrate invasion. The integrated interpretation is based
on an inversion workflow, which employs a multi-layered integrated
flow/electrical numerical simulator to solve the fluid flow, the
salt transport and the DC electrical array response. A field
experiment conducted in a carbonate reservoir has been used to
verify the proposed method.
[0012] Most of the above works involve the simulation of the
mud-filtrate invasion process to obtain the invasion rates as
accurately as possible. However, the mud-filtrate invasion is a
dynamic and complex process. It starts with a high invasion rate,
which is mainly controlled by the overbalance pressure and the
formation permeability when a well is drilled. After the mud-cake
is formed, the invasion rate will tend to be controlled by the
mud-cake since the permeability of mud-cake is normally much lower
than that of the formation. There have been extensive works on the
simulation of the mud-filtrate invasion process, which although
possible, is still too complex for any practical application,
especially when considering the complicated drilling and reaming
conditions, such as the mud circulation rate, workover, mud-cake
scraped by the tools, and so on.
SUMMARY OF THE INVENTION
[0013] This method according to some embodiments of the invention
paper includes the general framework developed in Alpak et al.
(2004). However, the method of the invention propose to directly
invert the average mud-filtrate invasion rates for each flow unit
or each depth interval, in order to avoid using complicated,
time-consuming mud-cake growth and invasion models. Moreover, based
on the derived invasion profile, the joint inversion technique
according to some embodiments of the invention also provide a more
reliable interpretation of formation test compared to the
conventional methods. The method according to the invention can be
applied to support horizontal and deviated wells with dipping
formation layers.
[0014] According to some embodiments, a method of estimating fluid
invasion rate for portions of a subterranean formation of interest
surrounding a borehole is provided. The method includes receiving
electromagnetic survey data of the subterranean formation of
interest and receiving fluid flow characteristics relating to the
formation of interest. A direct inversion is performed of the
electromagnetic data and constrained by the fluid flow simulator
for a fluid invasion rate, and the fluid invasion rate is thereby
estimated based. The fluid can be borehole originating fluids such
as mud filtrate or completion fluids. The fluid flow simulator can
also be used to generate the pressure transient data made in the
borehole. If the pressure transient measurements are available to
be compared with the one generated from the fluid flow simulator
then permeability can be estimated. According to some embodiments
porosity is also estimated based on the inversion. A quantitative
predictive formation model can be generated based on the inversion.
The electromagnetic data can be obtained using various tools
operated in a single well. The same or similar workflow or process
can be applied to different configurations such as cross-well,
surface to borehole, borehole to surface and surface to
surface.
[0015] According to some embodiments, a method of estimating
porosity for a subterranean formation of interest including one or
more geologic beds is provided. The method includes receiving
electromagnetic survey data and fluid flow characteristics of the
formation of interest, and estimating porosity within at least a
portion of the formation of interest based at least in part on the
electromagnetic survey data and the fluid flow characteristics. The
estimation has a spatial resolution such that at least two porosity
values are estimated for each of the one or more geologic beds.
[0016] According to some embodiments, systems are provided for
estimating fluid invasion rate and for estimating porosity for
portions of a subterranean formation of interest surrounding a
borehole.
[0017] Additional objects and advantages of the invention will
become apparent to those skilled in the art upon reference to the
detailed description taken in conjunction with the provided
figures.
BRIEF DESCRIPTION OF THE FIGURES
[0018] FIG. 1 is a flow chart showing workflow steps of a joint
inversion, according to some embodiments;
[0019] FIGS. 2a and 2b are diagrams showing certain configurations,
according to some embodiments;
[0020] FIGS. 3a, 3b and 3c are plots illustrating the saturation
dependent functions used in certain described examples;
[0021] FIGS. 4a-4h, are vertical cross section diagrams show the
results for oil-water example shown in FIG. 2a;
[0022] FIGS. 5a and 5b are log-log plots of the pressure and the
pressure derivative for the build-up for the example shown in FIG.
2a;
[0023] FIG. 6 is a two-dimensional vertical cross-section of a
reservoir having an injector and producer well, according to some
embodiments;
[0024] FIG. 7 shows the source-receiver setup and associated
surface equipment according to the example shown in FIG. 6; and
[0025] FIGS. 8a-h are vertical cross-sections of spatial
distributions of the porosity, the conductivity, the water
saturation and the salt concentration, according to some
embodiments.
[0026] It will be recognized by the person of ordinary skill in the
art, given the benefit of this disclosure, that certain dimensions,
features, components, and the like in the figures may have been
enlarged, distorted or otherwise shown in a non-proportional or
non-conventional manner to facilitate a better understanding of the
technology disclosed herein.
DETAILED DESCRIPTION OF EXEMPLARY EMBODIMENTS
[0027] The following description provides exemplary embodiments
only, and is not intended to limit the scope, applicability, or
configuration of the disclosure. Rather, the following description
of the exemplary embodiments will provide those skilled in the art
with an enabling description for implementing one or more exemplary
embodiments. It being understood that various changes may be made
in the function and arrangement of elements without departing from
the spirit and scope of the invention as set forth in the appended
claims.
[0028] Specific details are given in the following description to
provide a thorough understanding of the embodiments. However, it
will be understood by one of ordinary skill in the art that the
embodiments may be practiced without these specific details. For
example, systems, processes, and other elements in the invention
may be shown as components in block diagram form in order not to
obscure the embodiments in unnecessary detail. In other instances,
well-known processes, structures, and techniques may be shown
without unnecessary detail in order to avoid obscuring the
embodiments. Further, like reference numbers and designations in
the various drawings indicated like elements.
[0029] Also, it is noted that individual embodiments may be
described as a process which is depicted as a flowchart, a flow
diagram, a data flow diagram, a structure diagram, or a block
diagram. Although a flowchart may describe the operations as a
sequential process, many of the operations can be performed in
parallel or concurrently. In addition, the order of the operations
may be re-arranged. A process may be terminated when its operations
are completed, but could have additional steps not discussed or
included in a figure. Furthermore, not all operations in any
particularly described process may occur in all embodiments. A
process may correspond to a method, a function, a procedure, a
subroutine, a subprogram, etc. When a process corresponds to a
function, its termination corresponds to a return of the function
to the calling function or the main function.
[0030] Furthermore, embodiments of the invention may be
implemented, at least in part, either manually or automatically.
Manual or automatic implementations may be executed, or at least
assisted, through the use of machines, hardware, software,
firmware, middleware, microcode, hardware description languages, or
any combination thereof. When implemented in software, firmware,
middleware or microcode, the program code or code segments to
perform the necessary tasks may be stored in a machine readable
medium. A processor(s) may perform the necessary tasks.
[0031] Fluid Flow Simulation.
[0032] Further detail relating to fluid flow simulation, according
to some embodiments, will now be provided. After a well is drilled
the mud-filtrate invades into the surrounding rock formation. A
three-phase, four-component (oil, water, gas, salt) model is
employed to describe the mud-filtrate invasion process. The
mud-filtrate invasion is simulated as an injection process by
solving the pressure diffusive equation derived from Darcy's Law
and the mass balance equation. Ignoring the chemical reactions and
diffusive/dispersive transport, the mass balance equation for the
fluid phase can be stated as follows (see Aziz and Settari
(1979)):
.differential. ( .rho. l .phi. S l ) .differential. t + .gradient.
( .rho. l v l ) = - q l , l = 1 , 2 , 3 , ( 1 ) ##EQU00001##
where the subscript l is the index of phase, which can be oil,
water or gas. The symbols .rho., v, .phi., q and S denote the fluid
density, fluid velocity, porosity, source term, and phase
saturation, respectively. The symbol .gradient. denotes the spatial
differentiation with respect to spatial position (x,y,z). Darcy's
law is used as the constitutive equation for the velocity vector,
i.e.,
v l = - k k r ; i .mu. l ( .gradient. .rho. l - .gamma. l
.gradient. D z ) , ( 2 ) ##EQU00002##
where k is the tensor absolute permeability, k.sub.r is the phase
relative permeability, .mu. is the phase viscosity, p is the phase
pressure, .gamma. is the phase specific gravity and D.sub.z is the
vertical depth below a reference level. We assume that the fluid
(oil and water) and rock compressibilities are defined by:
c fluid = - 1 V .differential. V .differential. p T = 1 .rho.
.differential. .rho. .rho. .differential. p T , ( 3 )
##EQU00003##
They are constant over the pressure range of interest in the
simulation. The gas compressibility is treated as a function of the
pressure-volume-temperature (PVT) properties of gas for each
time-step. Additional relations about capillary pressures and fluid
saturations are given by:
p.sub.c=p.sub.nw-p.sub.w, (5)
and
.SIGMA..sub.iS.sub.1=1 (6)
where the subscripts nw and w stand for the nonwetting and the
wetting phase.
[0033] Salt is the forth component which is simulated, according to
some embodiments. Assuming there is a high salt concentration
contrast between the injected fluid and the in-situ formation brine
and ignoring the diffusive/dispersive transport, the isothermal
salt mixing in the aqueous-phase is solved by using the convective
transport equation as follows:
.differential. ( .rho. w .phi. S w C w ) .differential. t +
.gradient. ( .rho. w v w C w ) = - C wi q i , ( 7 )
##EQU00004##
where C.sub.w, C.sub.wi and q.sub.i is the salt concentration of
formation brine, the salt concentration of invading fluid and the
fluid injection rate, respectively. In our study the above
equations are solved numerically using a finite-difference based
reservoir simulator in fully implicit, black-oil mode with brine
tracking option (e.g. ECLIPSE.TM.). The pressure transient in the
well/formation test can be simulated as well as the temporal and
spatial distribution of the phase saturation and the salt
concentration. In this description, focus is made on the scenario
where a single well penetrates layered hydrocarbon-bearing
formations. In the case of a vertical well, a cylindrical
axisymmetric model is used in the numerical simulation and a
logarithmic radial grid with very fine cells around the borehole is
used to honor the mud-filtrate invasion and accurately catch the
pressure transient phenomena in the well/formation test.
[0034] Electromagnetic Simulation.
[0035] Further detail relating to electromagnetic simulation,
according to some embodiments, will now be provided.
Electromagnetic responses measured by the array induction tools can
be numerically simulated by solving a frequency-domain
electromagnetic induction problem, which is described by the
following diffusive Maxwell's equations:
.gradient..times.E+i.omega..mu.H=0, (8)
.gradient..times.H-.sigma.E=-J, (9)
where E is the electric field vector, H is the magnetic field
vector and J is the external electric current source vector. The
symbols .sigma.=.sigma.(x,y,z) and .mu. denote the electrical
conductivity and the magnetic permeability (which is constant). The
symbol .omega. denotes the radial frequency and i= {square root
over (-1)}. By substituting equation (8) into equation (9) we
obtain:
.sigma..sup.-1.gradient..times..gradient..times.E+i.omega.E=-i.omega..si-
gma..sup.-1J (10)
[0036] A frequency-domain finite-difference algorithm such as
described in Druskin, V., L. Knizhnerman, and P. Lee, 1999, New
spectral lanczos decomposition method for induction modeling in
arbitrary 3d geometry: Geophysics, 64, 701-706 (herein after
"Druskin et al. (1999)"), incorporated herein by reference, is
employed to numerically solve equation (10) on a staggered Yee grid
(See, Yee, K. S., 1966, Numerical solution of initial boundary
value problems involving maxwell's equations in isotropic media:
IEEE Transactions on Antennas and Propagation, AP-14, 302-307).
This approach uses a spectral Lanczos decomposition method (SLDM)
with Krylov subspaces generated from the inverse powers of the
Maxwell operator. Compared to the standard SLDM method, this
so-called SLDMINV method is significantly more efficient. As shown
in Druskin et al. (1999), the convergence rate of the SLDMINV
method increases as the frequency decreases, this makes SLDMINV
particularly attractive for low-frequency electromagnetic
simulations.
[0037] Petrophysical Transform.
[0038] Further detail relating to petrophysical transforms,
according to some embodiments, will now be provided. According to
some embodimetns, the spatial distribution of water saturation and
salt concentration at logging times are computed by a multiphase
fluid flow simulator such as ECLIPSE.TM. from Schlumberger. These
are then transformed to a conductivity map using Archie's law (see,
Archie, G. E., 1942, The electrical resistivity log as an aid in
determining some reservoir characteristics: Petroleum Transactions
of the AIME, 146, 54-62):
.sigma. = .sigma. w .phi. m S w n a , ( 11 ) ##EQU00005##
where .sigma..sub.w is the brine conductivity and S.sub.w is the
water saturation. Archie's parameters m, n and a are obtained from
core measurements. The brine conductivity is a function of the salt
concentration and temperature and is given by the following
equation (see, Zhang, J.-H., Q. Hu, and Z.-H. Liu, 1999, Estimation
of true formation resistivity and water saturation with a
time-lapse induction logging method: The Log Analyst, 40,
138-148):
.sigma. w = [ ( 0.0123 + 3647.5 C w 0.0955 ) 82 1.8 T + 39 ] - 1 ,
( 12 ) ##EQU00006##
where C.sub.w is the salt concentration (ppm) and T is the
formation temperature (.degree. C.). In this example, The
multiphase fluid flow model is solved on a cylindrical grid while
solving the electromagnetic model on a Cartesian grid. A robust
homogenization approach described in Moskow, S., V. Druskin, T. M.
Habashy, P. Lee, and S. Davydycheva, 1999, A finite difference
scheme for elliptic equations with rough coefficients using a
cartesian grid nonconforming to interfaces: SIAM Journal on
Numerical Analysis, 36, 442-464, incorporated herein by reference,
is employed to map the conductivities from the three-dimensional
(3D) cylindrical grid to the 3D Cartesian grid. This feature gives
us more flexibility on constructing the model.
[0039] Inversion Algorithm.
[0040] Further detail relating to inversion algorithms, according
to some embodiments, will now be provided. In the inversion process
the electromagnetic data and pressure data (whenever available) are
simultaneously inverted. After spatial discretization, the
nonlinear inverse problem can be described by the following
operator notations:
M.sup.e= S.sup.e( m)+ .epsilon..sup.e, (13)
M.sup.p= S.sup.p( m)+ .epsilon..sup.p, (14)
where M.sup.e is the measured electromagnetic data while M.sup.p is
the measured pressure data. The components of these vectors are
given by
M.sup.e=[M( r.sub.i.sup.s, r.sub.j.sup.r,t.sub.l),i=1, . . . ,
N.sup.s;j=1, . . . , N.sup.r;l=1, . . . , N.sup.l]T, (15)
M.sub.p=[M( r.sub.i.sup.g,t.sub.j),i=1, . . . , N.sup.g;j=1, . . .
, N.sup.t].sup.T, (16)
where N.sup.s, N.sup.r and N.sup.l are the number of sources,
receivers and logging snapshots in the electromagnetic measurements
while N.sup.g and N.sup.t are the number of pressure gauges and
time steps in the pressure measurements. In equations (13) and
(14), S.sup.e and S.sup.p are the vector of simulated
electromagnetic and pressure data as functions of the vector model
parameter m. The parametric model according to this example is
layer by layer horizontal permeabilities, vertical permeabilities,
porosities and average mud-filtrate invasion rates denoted by
k.sub.h, k.sub.v, .phi. and q. The vector of model parameters is
defined as follows:
m=[k.sub.h;1,k.sub.h;2, . . . ,k.sub.h;L,k.sub.v;1,k.sub.v;2, . . .
,k.sub.v;L,.phi..sub.1,.phi..sub.2, . . .
,.phi..sub.L,q.sub.1,q.sub.2, . . . ,q.sub.L].sup.T (17)
where L is the number of layers in our one-dimensional (1D)
parametric model. In equations (13) and (14), .epsilon..sup.e and
.epsilon..sup.p denote the measurement errors in the
electromagnetic and pressure data.
[0041] The inverse problems are posed as a minimization problem of
the following cost function:
C.sub.n( m)=.parallel. W.sub.d.sup.e[ M.sup.e- S.sup.e(
m)].sup.2.parallel.+.parallel. W.sub.d.sup.p[ M.sup.p- S.sup.p(
m)].parallel..sup.2+.lamda..sub.n.parallel.ln( m)-ln(
m.sub.n).parallel..sup.2, (18)
where .lamda..sub.n is the regularization parameter at iteration n.
In the regularization term, the logarithmic functions are employed
since the unknown model parameters have different dimensions and
can have different orders of magnitude. The diagonal matrices
W.sub.d.sup.e and W.sub.d.sup.p are the data weighting matrices.
There are several options for choosing these data weighting
matrices. One particular choice is to account for all data points
on equal footing. Hence, the data weighting matrices can be defined
as follows:
W d ; i , i e = 1 N e M i e , i = 1 , , N e , ( 19 ) W d ; i , i p
= 1 N p M i p i = 1 , , N p ( 20 ) ##EQU00007##
where N.sup.e=N.sup.s.times.N.sup.r.times.N.sup.l and
N.sup.p=N.sup.g.times.N.sup.t. Substituting equations (19) and (20)
into equation (18) we arrive at
C n ( m _ ) = 1 N e i = 1 N e M i e - S i e ( m _ ) M i e 2 + 1 N p
i = 1 N p M i p - S i p ( m _ ) M i p 2 + .lamda. n l = 1 4 L ln (
m l ) - ln ( m n , l ) 2 , ( 21 ) ##EQU00008##
Another choice is to apply uniform weighting on all data points as
follows:
W d ; i , i e = 1 j = 1 N e M j e 2 , i = 1 , , N e , ( 22 ) W d ;
i , i p = 1 j = 1 N p M j p 2 , i = 1 , , N p , ( 23 )
##EQU00009##
Substituting equations (22) and (23) into equation (18) we arrive
at
C n ( m _ ) = i = 1 N e M i e - S i e ( m _ ) 2 i = 1 N e M i e 2 +
i = 1 N p M i p - S i p ( m _ ) 2 i = 1 N p M i p 2 + .lamda. n l =
1 4 L ln ( m l ) - ln ( m n , l ) 2 . ( 24 ) ##EQU00010##
[0042] The minimization process is done by using the Gauss-Newton
approach described in Habashy and Abubakar, incorporated herein by
reference. The gradient vector of the cost function in equation
(18) is given by
g i ( m _ ) = .gradient. m i C n ( m _ ) = j = 1 N e J j , i e ( m
_ ) W d ; j , j e W d ; j , j e [ M j e - S j e ( m _ ) ] + j = 1 N
p J j , i p ( m _ ) W d ; j , j p W d ; j , j p [ M j p - S j p ( m
_ ) ] + .lamda. n 1 m i [ ln ( m i ) - ln ( m n , i ) ] 2 , i = 1 ,
, 4 L ( 25 ) ##EQU00011##
Hence, the gradient vector at the n-th iteration is given by
g n , i = g i ( m _ n ) = j = 1 N e J j , i e ( m _ n ) W d ; j , j
e W d ; j , j e [ M j e - S j e ( m _ n ) ] + j = 1 N p J j , i P (
m _ n ) W d ; j , j p W d ; j , j p [ M j p - S j p ( m _ n ) ] . (
26 ) ##EQU00012##
The jacobian matrices are calculated using finite-difference
approximations as follows:
J j , i e ( m _ n ) = .differential. S j e ( m _ ) .differential. m
i m _ = m _ n .apprxeq. - S j e ( m _ n + .delta. i , k .DELTA. m _
n ) - S j e ( m _ n ) .DELTA. m n , i , ( 27 ) j = 1 , , N e , i ,
k = 1 , , 4 L and J j , i p ( m _ n ) = - .differential. S j p ( m
_ ) .differential. m i m _ = m _ n .apprxeq. - S j p ( m _ n +
.delta. i , k .DELTA. m _ n ) - S j p ( m _ n ) .DELTA. m n , i (
28 ) j = 1 , , N p , i , k = 1 , , 4 L ##EQU00013##
The Hessian matrix is given by
H i , k ( m _ ) = .gradient. m q .gradient. m i C n ( m _ )
.apprxeq. j = 1 N e J j , i e ( m _ ) W d ; j , j e W d ; j , j e J
j , k e ( m _ ) + j = 1 N p J j , i p ( m _ ) W d ; j , j p W d ; j
, j p J j , k p ( m _ ) + .lamda. n 1 m i 1 m k , i , k = 1 , , 4 L
( 29 ) ##EQU00014##
Hence, the Hessian matrix at the n-th iteration is given by
H i , k ( m _ n ) .apprxeq. j = 1 N e J j , i e ( m _ n ) W d ; j ,
j e W d ; j , j e J j , k e ( m _ n ) + j = 1 N p J j , i p ( m _ n
) W d ; j , j p W d ; j , j p J j , k p ( m _ n ) + .lamda. n 1 m n
, i 1 m n , k . ( 30 ) ##EQU00015##
In the Gauss-Newton minimization process, the search vector p.sub.n
can be constructed from the following linear equation:
H.sub.n p.sub.n=- g.sub.n. (31)
[0043] Equation (31) is solved via a Gauss elimination procedure.
The regularization parameter .lamda..sub.n is chosen as
follows:
.lamda. n = .beta. max { .tau. i } , if min { .tau. i } max { .tau.
i } < .beta. , i = 1 , , 4 L , ( 32 ) ##EQU00016##
where .beta. is a constant value determined from numerical
experiments and .tau. are the eigenvalues of the Hessian matrix
H.
[0044] We employ a line-search procedure as in Habashy and Abubakar
to enforce the reduction of cost function in each iteration as
follows:
m.sub.n+1= m.sub.n+v p.sub.n, (33)
where 0<v.sub.n.ltoreq.1 is the step-length. The step-length is
calculated by solving the minimization problem:
v n = arg min v { C ( m _ n + v p _ n ) } . ( 34 ) ##EQU00017##
[0045] In order to avoid expensive evaluation of the cost function
C( m), we adopt an algorithm (see Dennis Jr., J. E. and R. B.
Schnabel, 1983, Numerical methods for unconstrained optimization
and nonlinear equations: Prentice Hall, Inc., New Jersey) whereby
v.sub.n is selected such that the average rate of decrease from C(
m.sub.n) to C( m.sub.n+v.sub.n p.sub.n) is at least some prescribed
fraction, .alpha., of the initial rate of decrease at m.sub.n along
the direction p.sub.n, i.e.,
C( m.sub.n+v.sub.n p.sub.n).ltoreq.C(
m.sub.n)+.alpha.v.sub.n.delta.C.sub.n+1, (35)
where 0<.alpha.<1 is a fractional number which, in practice,
is set quite small (we use 10.sup.-4 herein) so that hardly more
than a decrease in function value is required. .delta.C.sub.n+1 is
the rate of decrease of C( m) at m.sub.n along the direction
p.sub.n and is given by:
.delta. C n + 1 = .differential. .differential. v C ( m _ n + v p _
n ) v = 0 = g _ T ( m _ n ) p _ n . ( 36 ) ##EQU00018##
We first employ the full Gauss-Newton search step in the
minimization process. If v.sub.n=1 fails to satisfy the criterion
given by equation (35), then we reduce v.sub.n to backtrack along
the direction p.sub.n until an acceptable next m.sub.n+1 is found.
Assume, at the (n+1)-th Gauss-Newton iteration, v.sub.n.sup.k is
the current step-length that does not satisfy equation (35), we
compute the next backtracking step-length, v.sub.n.sup.k+1, by
searching for the minimum of the cost function, f(v).ident.C(
m.sub.n+v p.sub.n), which can be approximated by a quadratic
expression. Thus, v.sub.n.sup.k+1, which is the minimum of f(v),
for k=0, 1, 2, . . . is given by:
v n k + 1 = - ( v n k ) 2 2 .delta. C n + 1 C ( m _ n + v n k p _ n
) - C ( m _ n ) - v n k .delta. C n + 1 . ( 37 ) ##EQU00019##
[0046] To impose a priori information of maximum and minimum bounds
on the unknown parameters, we constrained them using a non-linear
transformation of the form:
m = f ( c , m min , m max ) = m min + m max - m min c 2 + 1 c 2 , -
.infin. < c < + .infin. . ( 38 ) ##EQU00020##
[0047] The Gauss-Newton search step in m,
p.sub.n=m.sub.n+1-m.sub.n, is related to the step in c,
q.sub.n=c.sub.n+1-c.sub.n, by
p = m c q = 2 m max - m m max - m min ( m max - m ) ( m - m min ) q
. ( 39 ) ##EQU00021##
The inverted parameter m is then updated using the following
equation,
m n + 1 = m min + m max - m min .alpha. n 2 + ( m n - m min ) ( m
max - m n ) 3 .alpha. n 2 , where ( 40 ) .alpha. n = ( m n - m min
) ( m max - m n ) + 1 2 ( m max - m min ) v n p n . ( 41 )
##EQU00022##
[0048] FIG. 1 is a flow chart showing workflow steps of a joint
inversion, according to some embodiments. The initial guess of the
unknown parameters (permeability, porosity and mud-filtrate
invasion rate) 112 and additional information such as PVT, Kr, Pc,
formation geometry, mud-filtrate and formation brine concentration
et al. 110 are used for building a reservoir model. A multiphase
fluid flow simulator 114 is run on the reservoir model generated
from the data 110 and the initial parameters 112, thereby
generating spatial distributions 116 of water saturation and salt
concentration. A conductivity-saturation model 118 is used to
transform the spatial distributions 116 into a conductivity
distribution 120. Electromagnetic simulator 122 generates simulated
electromagnetic data 126 based on the input conductivity
distribution 120. Simulated pressure and flow rate data 124 is
generated by multiphase fluid flow simulator 114. The pressure and
flow rate data 124, electromagnetic data 126, measured pressure and
flow rate data 128, and measured electromagnetic data 130 are input
into cost function 132. Cost function 132 is combination of
differences of 124 to 128, and 126 to 130. In step 136, values are
updated for permeability, porosity and mud-filtrate invasion rate
112 until the cost function 132 becomes less than a predefined
tolerance, thereby yielding the inverted parameters and derived
distributions of water saturation and conductivity 134.
[0049] Test Cases.
[0050] Further detail relating to two test cases, according to some
embodiments, will now be provided. To demonstrate the workflow,
according to some embodiments, two synthetic examples were employed
consisting of a vertical well penetrating layered
hydrocarbon-bearing horizontal formations with sealing upper and
lower shoulder beds. For reasons of facilitating comparisons with
known results, the same examples are employed which can be referred
to in Alpak, et al. (2004).
[0051] In this example, it is assumed that there is no-flux
surrounding the boundaries. The water-based mud-filtration invasion
is taking place since the well is drilled. The induction logging
data are assumed to be recorded using an array induction tool,
which has one transmitter and an eight-receiver array, providing
fourteen pairs of apparent conductivity logs for different depths
of investigation. Immediately after the first induction logging, a
formation test is carried out using a packer-probe formation tester
such as with Schlumberger's Modular Formation Dynamics Tester (MDT)
tools. Three time-series of pressure transient measurements are
recorded at two probes and the packer interval. This is followed by
a second induction logging to sense the perturbed resistivity
around the packer interval caused by the formation test. This
process is referred to as a log-test-log mode. If the perturbation
caused by the formation test is not significant, the second
induction logging can be discarded, so the process is simplified to
the log-test mode. In these two modes there are two measurements,
namely the induction logging and pressure transient data, which are
then used for the joint inversion. According to some embodiments,
another mode is provided that only uses two induction logging data
collected at different times. This mode is referred to herein as a
time-lapse log mode. In a time-lapse log mode, the multiphase fluid
flow simulator is only used to constrain the inversion of the
induction logging data.
[0052] FIGS. 2a and 2b are diagrams showing certain configurations
(refer to Alpak, et al. (2004)), according to some embodiments. For
FIGS. 2a and 2b the drainage radii for wells 230 and 228
respectively, are both about 300 m. FIGS. 2a and 2b are
two-dimensional vertical cross-sections of the reservoir model
configuration. Impermeable shoulder beds 208 and 210 with known
homogeneous formation conductivities are superimposed on the top
and bottom of the reservoirs 212 and 214 respectively. The
permeable zones of interest are subject to the water-based
mud-filtrate invasion. In FIG. 2a, the reservoir 212 comprises a
single-layer TI-anisotropic formation 216 with oil-water system.
Mud filtrate invasion zones 232 for well 230, and 234 for well 228
are shown in FIGS. 2a and 2b respectively.
[0053] The case shown in FIG. 2a is a single-layer model with a
TI-anisotropic formation and oil-water system, while the case shown
in FIG. 2b is a three-layer model with isotropic formation and
gas-water system. In these examples, the fluid properties,
saturation dependent functions and the parameters in Archie's
equation are available from laboratory and core measurements. Both
the oleic and aqueous phases are slightly compressible while the
gaseous phase is compressible with its compressibility dependent on
pressure.
[0054] FIGS. 3a, 3b and 3c are plots illustrating the saturation
dependent functions used in certain described examples (refer to
Alpak, et al. (2004)). FIG. 3a shows the relative permeability for
the oil-water case shown in FIG. 2a, where the capillary pressure
is negligible. Curve 310 is the relative permeability of water and
curve 312 is the relative permeability of oil in the formation
layer 214. FIG. 3b shows the layer-by-layer relative permeability
for the gas-water case shown in FIG. 2b. Curve 314 shows the
relative permeability for gas in the upper and lower layers 218 and
222; curve 316 shows the relative permeability for water in the
upper and lower layers 218 and 222; curve 318 shows the relative
permeability of gas for the middle layer 220; and curve 320 shows
the relative permeability of water for the middle layer 220. FIG.
3c shows the layer-by-layer capillary pressure for gas-water case
shown in FIG. 2c. Curve 324 shows the capillary pressure for the
upper layer 118; curve 326 shows the capillary pressure for the
middle layer 220; and curve 328 shows the capillary pressure for
the lower layer 222.
[0055] For the oil-water case shown in FIG. 2a, the following
operation schedule is assumed: The well is subjected to the
mud-filtrate invasion after drilling. The first induction log is
recorded on the third day and then immediately followed by the
wireline formation test, which is composed of a draw-down at a
constant liquid rate of 4.76962 m.sup.3/d for 100 minutes followed
by a build-up for 200 minutes. In this example, the pressure
transient data is only used in the build-up period for the
inversion. Another induction log is then recorded to capture the
perturbation caused by the formation test. For the time-lapse log
mode, the second induction log is recorded at the tenth day.
[0056] In the fluid flow simulator we used a 31.times.6.times.30
grid with the finest radial cell size of about 0.0305 m, adjacent
to the borehole. The model is initialized with the irreducible
water saturation. Since the fluid and salt movement are convective
and the diffusion and salt dispersion can be neglected, the
mud-filtrate invasion process is calculated using an average
mud-filtrate invasion rate. It is assumed that the average invasion
rate is 0.0067 m.sup.3/d/m for the first three days. For the
time-lapse log mode, an invasion rate of 0.00577 m.sup.3/d/m
between the third and tenth day is also assumed. Other specific
model parameters are listed in Table 1.
TABLE-US-00001 TABLE 1 Petrophysical and fluid parameters and other
relevant information employed in oil-water example shows in FIG.
2a. (also, Refer to Alpak, et al. (2004)) Variable Value Formation
rock compressibility 7.25'10.sup.-8 [bar.sup.-1] Aqueous phase
viscosity 1.274 [cp] Aqueous phase density 1001.074 [kg/m.sup.3]
Aqueous phase formation volume factor 0.996 [rm.sup.3/sm.sup.3]
Aqueous phase compressibility 3.7'10.sup.-5 [bar.sup.-1] Oleic
phase viscosity 0.355 [cp] Oleic phase API density 42 [.degree.
API] Oleic phase density 815.564 [kg/m.sup.3] Oleic phase formation
volume factor 1.471 [rm.sup.3/sm.sup.3] Oleic phase compressibility
2.76'10.sup.-4 [bar.sup.-1] Formation pressure (refer to formation
top) 206.84 [bar] Wellbore radius 0.108 [m] Formation horizontal
permeability 100 [mD] Formation vertical permeability 20 [mD]
Formation porosity 0.25 [fraction] Formation temperature 220
[.degree. F.] Formation brine salinity 120000 [ppm] Mud-filtrate
salinity 5000 [ppm] a-constant in the Archie's equation 1.00
[dimensionless] m-cementation exponent in the Archie's 2.00
[dimensionless] equation n-water saturation exponent in the 2.00
[dimensionless] Archie's equation Mud conductivity 2631.58 [mS/m]
Upper and lower shoulder bed conductivities 1000.00 [mS/m] Logging
interval 0.61 [m]
[0057] In the gas-water three-layer example shown in FIG. 2b, a
similar operation schedule is used. However, the first induction
log is recorded at the first day, the draw-down for the formation
test lasts 100 minutes at a constant rate of 0.2385 m.sup.3/d, then
the tool is shut down for another 200 minutes for pressure
build-up. The second induction logging was not carried out since
the perturbation in formation resistivity caused by the formation
test is not significant. The fluid flow simulator uses a
91.times.6.times.30 grid with the finest radial cell size of about
0.0061 m. The model is initialized with the irreducible water
saturation. Different average invasion rates, permeabilities and
porosities are used for the three layers. Referring to FIG. 2b, for
the first-layer 218: k=6.83 mD, .phi.=0.139, q=0.42523 m.sup.3/d/m;
for the second-layer 220: k=40.98 mD, .phi.=0.137, q=0.42485
m.sup.3/d/m; and for the third-layer 222: k=14.33 mD, .phi.=0.125,
q=0.42548 m.sup.3/d/m. Other specific model parameters can be found
in Table 2. Some parameters are the same as in oil-water example of
FIG. 2a.
TABLE-US-00002 TABLE 2 Petrophysical and fluid parameters and other
relevant information employed in example shown in FIG. 2b. (Also
refer to Alpak, et al. (2004)) Variable Value Aqueous phase
viscosity 1.000 [cp] Aqueous phase density 1186.007 [kg/m.sup.3]
Aqueous phase formation volume factor 1.000 [rm.sup.3/sm.sup.3]
Aqueous phase compressibility 1.45'10.sup.-5 [bar.sup.-1] Gaseous
phase viscosity [f(p) for T = const.] 0.01087 [cp] @1 bar Gaseous
phase density [f(p) for T = const.] 0.977 [kg/m.sup.3] @1 bar
Gaseous phase formation volume factor [f(p) 0.782
[rm.sup.3/sm.sup.3] @1 bar for T = const.] Formation pressure
(refer to formation top) 6.62 [bar] Wellbore radius 0.107 [m]
Formation temperature 96 [.degree. F.] Formation brine salinity
200000 [ppm] Mud-filtrate salinity 2000 [ppm] Mud conductivity
387.60 [mS/m] Upper and lower shoulder bed conductivities 400.00
[mS/m]
[0058] Inversion Results.
[0059] In this section, the inversion results are described for the
synthetic induction logging data and pressure transient data, which
are generated by forward modelling based on pre-defined true
petrophysical parameters and invasion rates, with added zero-mean
Gauss random noise. FIGS. 4-6, show the results for oil-water
example shown in FIG. 2a. In FIG. 4a, diagram 410 shows the
vertical cross-section of spatial distributions of the
reconstructed water saturation for the first logging run for the
oil-water example shown in FIG. 2a. In FIG. 4b, diagram 412 shows
the vertical cross-section of spatial distributions of the
reconstructed water saturation for the second logging run for the
oil-water example shown in FIG. 2a. In FIG. 4c, diagram 414 shows
the vertical cross-section of spatial distributions of the
reconstructed salt concentration for the first logging run for the
oil-water example shown in FIG. 2a. In FIG. 4d, diagram 416 shows
the vertical cross-section of spatial distributions of the
reconstructed salt concentration for the second logging run for the
oil-water example shown in FIG. 2a. In FIG. 4e, diagram 418 shows
the vertical cross-section of spatial distributions of the
reconstructed formation conductivity for the first logging run for
the oil-water example shown in FIG. 2a. In FIG. 4f, diagram 420
shows the vertical cross-section of spatial distributions of the
reconstructed formation conductivity for the second logging run for
the oil-water example shown in FIG. 2a. In all cases there was
Gaussian random noise of up to 7%. A high conductivity annulus is
observed, which is caused by the invasion of the low salinity
water-based mud-filtrate into the hydrocarbon-bearing formation
with high salinity connate water. FIGS. 4b, 4d and 4f show oil
breaking through the high conductivity annulus caused by draw-down
in the formation test. In FIG. 4g, diagram 422 plots the difference
between the true conductivity and the reconstructed conductivity
distributions for the first logging run for the oil-water example
shown in FIG. 2a. In FIG. 4h, diagram 424 plots the difference
between the true conductivity and the reconstructed conductivity
distributions for the second logging run for the oil-water example
shown in FIG. 2a. It can be seen from FIGS. 4h and 4g that there is
good agreement between the true conductivity and reconstructed
conductivity distributions.
[0060] FIGS. 5a and 5b are log-log plots of the pressure and the
pressure derivative for the build-up for the example shown in FIG.
2a. In FIG. 5a plots results with 0% noise, and in FIG. 5b plots
results with 7% noise. In FIG. 5a, curve 510 plots the measured
pressure, curve 512 plots the measured derivative pressure, curve
514 plots the post-inversion pressure, and curve 516 plots the
post-inversion derivative pressure. In FIG. 5b, curve 520 plots the
measured pressure, curve 522 plots the measured derivative
pressure, curve 524 plots the post-inversion pressure, and curve
526 plots the post-inversion derivative pressure.
[0061] The inverted permeability, porosity and average mud-filtrate
invasion rate are listed together with the number of Gauss-Newton
iterations in Tables 3 and 4.
TABLE-US-00003 TABLE 3 Log-test-log inversion results for example
shown in FIG. 2a Log-test-log .phi. Num. of Add. Gauss k.sub.h
k.sub.v (frac- q itera- noise level (md) (md) tion) (m.sup.3/d/m)
tions Init. 40 40 0.12 0.00052 guess True 100 20 0.25 0.00670
Inverted 0% 100.00 20.00 0.250 0.00670 8 Inverted 1% 100.78 19.08
0.250 0.00670 6 Inverted 3% 102.43 17.02 0.250 0.00670 5 Inverted
5% 104.24 14.86 0.250 0.00672 7 Inverted 7% 105.99 13.54 0.250
0.00672 6 Analytic 0% 80.26 18.32
TABLE-US-00004 TABLE 4 Time-lapse log inversion results for example
shown in FIG. 2a Time-lapse log .phi. Num. of Add. Gauss k (frac-
q.sub.1 q.sub.2 itera- noise level (md) tion) (m.sup.3/d/m)
(m.sup.3/d/m) tions Init. N/A 0.12 0.001 0.00052 guess True N/A
0.25 0.00670 0.00577 Inverted 0% N/A 0.250 0.00670 0.00577 4
Inverted 1% N/A 0.250 0.00677 0.00562 3 Inverted 3% N/A 0.249
0.00721 0.00526 2 Inverted 5% N/A 0.249 0.00732 0.00495 2
[0062] The inversion results demonstrate that we can obtain
reliable horizontal permeability, porosity and average invasion
rate for up to 7% noise level. The inversion of the porosity is
very robust with little sensitivity to the noise level. This is
because the induction logging data is governed by the formation
conductivity, which is strongly affected by the porosity. With
increasing noise level, the quality of the inverted vertical
permeability quickly deteriorates. This is expected because the
sensitivity of pressure transient measurements to the vertical
permeability are relatively low and would be overwhelmed by noise.
The inversion of the horizontal permeability shows good accuracy.
In the time-lapse log mode, the formation permeability cannot be
inverted from induction measurements only and without the pressure
data since the invasion process is mainly governed by the mud-cake
and moreover we are using average invasion rates as the boundary
condition in the fluid flow simulation.
[0063] The horizontal and vertical permeabilities can also be
independently obtained from the formation test by using the
analytic straight-line computation, which is based on the
single-phase model without taking into account the invasion
profile. This approach can be applied on the pressure transient
data, which are generated using the fluid flow simulator for the
oil-water example shown in FIG. 2a. Only the build-up section of
the data are interpreted. But the build-up time is prolonged to
reach a stable radial flow status (FIGS. 5a and 5b show that the
radial flow has not formed yet). No noise was introduced in this
case to the synthetic pressure transient data. The analytical
interpretation results of the build-up section are also listed in
Table 3 for comparison. It can be observed that the joint inversion
produces more accurate estimation of formation permeabilities.
TABLE-US-00005 TABLE 5 Log-test inversion results for example shown
in FIG. 2b. .phi. Num. of Add. Gauss k (frac- q itera- noise level
(md) tion) (m.sup.3/d/m) tions Init. 12.39 0.18 0.26081 guess Layer
1 True 6.83 0.139 0.42523 Inverted 0% 6.83 0.139 0.42523 8 Inverted
5% 6.61 0.139 0.41693 10 Layer 2 True 40.98 0.137 0.42485 Inverted
0% 40.98 0.137 0.42486 8 Inverted 5% 47.15 0.139 0.47190 10 Layer 3
True 14.33 0.125 0.42548 Inverted 0% 14.33 0.125 0.42547 8 Inverted
5% 14.33 0.123 0.40225 10
[0064] Thus, according to some embodiments, a parametric inversion
algorithm is described to jointly invert the induction logging and
pressure transient measurements, which can be used for the
determination of the horizontal and vertical permeability, the
porosity and the average mud-filtrate invasion rate for each flow
unit or unit height. When using pressure measurements, this
algorithm can reliably invert the permeability without
incorporating a mud-cake growth and invasion model. Additionally,
the porosity can be inverted robustly either using the induction
log together with pressure measurement, or only by using induction
logs but constrained by the fluid flow simulator.
[0065] Traditionally, the mud-filtrate invasion profile can be
derived from the inversion of induction measurements, which are
based on a simplified step-profile invasion model. Inversion of the
induction logs constrained by the fluid flow simulator improves the
estimation of the mud-filtrate invasion profile, which becomes
consistent with the fluid flow physics. On the other hand, the
conventional interpretation of formation test is based on the
single-phase flow model without taking into account fluid invasion
or the need to pump a long time to clean up the contamination of
the mud-filtrate before the test procedure. The joint inversion of
induction logging and pressure measurements enable a more reliable
result. This technique can accommodate the existence of the
mud-filtrate in the region of formation test, and hence can save
costly rig time for clean up.
[0066] As a by-product, through the fluid-flow constrained
time-lapse induction log inversion, the temporal and spatial
distributions of the saturation and conductivity can be obtained,
which can be used to monitor the fluid front movement. According to
some embodiments, the inversion techniques described herein are
used to generate a quantitative predictive formation model that is
dynamic, rather than static, to perform fluid flow simulations.
[0067] According to some embodiments, the inversion procedures
described herein with respect to mud-filtrate invasion rates are
applied to estimate completion fluid invasion rates. For example
the invasion zones 232 and 234 in FIGS. 2a and 2b respectively,
according to these embodiments, represent completion fluid invasion
zones. According to some embodiments, a subsequent electromagnetic
survey can be conducted for the purpose of estimating the
completion fluid invasion rates. However, according to other
embodiments, a quantitative predictive formation model based on
earlier electromagnetic surveys for mud filtrate invasion can also
be used to directly predict the completion fluid invasion.
[0068] According to some embodiments, the result of the inversion
procedures is used to aid well clean up job design. A "clean up"
job can be performed following a drilling procedure and/or
completion procedure to reduce the effects that the drilling fluid
and/or completion fluid has on the producing formation. Based on
the invasion profile and the quantitative predictive formation
model derived from the inversion techniques described herein, a
fluid flow simulation is performed to find out how the clean up
procedure should be operated.
[0069] According to some embodiments, the described workflow can be
applied to deviated and horizontal wells with dipping layers and
cross-well applications. In these cases the cylindrical grid in the
fluid flow simulator should be replaced by the Cartesian grid.
[0070] In addition to the 1D parametric inversions of porosity for
each layer, according to some embodiments, the approach can also be
applied to invert a spatial distribution of porosity. FIGS. 6 and 7
illustrate an example used to demonstrate this capability. FIG. 6
is a two-dimensional vertical cross-section of a reservoir having
an injector and producer well, according to some embodiments. The
reservoir 614 is composed of three layers, 618, 620 and 622, with a
total thickness of 25 meters. The sealing upper and lower shoulder
beds 610 and 612 allow for no-flux surrounding the boundaries. The
permeability and porosity of each layer are different. The upper
layer 618 has k.sub.h=30 mD, the middle layer 620 has k.sub.h=100
mD, and the lower layer 622 has k.sub.h=70 mD. It is assumed that
the formations are TI-anisotropic with constant anisotropy ratio,
k.sub.v/k.sub.h=0.1. An injection well 630 and a production well
634 penetrate all the layers 610, 618, 620, 622 and 612. Both wells
630 and 634 are vertical and the separation is 80 meters. An
oil-water system is used and it is assumed that the fluid
properties, saturation dependent functions and coefficients used in
the Archie's equation are available from laboratory and core
measurements. For simplicity, it is assumed that the injection well
630 is operated under bottomhole pressure control, which is a
constant 241.32 bar. The production well 634 produces a constant
oil rate of 158.99 m.sup.3/d. After 5 days injection and
production, a crosswell electromagnetic logging is conducted. Water
flooding zone 632 is also shown.
[0071] FIG. 7 shows the source-receiver setup and associated
surface equipment according to the example shown in FIG. 6. In the
crosswell electromagnetic logging process, sources 702 are located
in the injection well 630, while receivers 708 are located in the
production well 634. Although only 14 source stations and 14
receiver stations are shown on the toolstrings in FIG. 7 for
simplicity, in this example there are 29 source stations 702 and 29
receiver stations 708. Thus there are 29.times.29 source-receiver
combinations. For example, each of the receivers 708 receiver
signals form source 704 and from source 706 as shown by the dashed
and solid lines in FIG. 7. The data are the vertical component of
the magnetic fields and the sources operate at 500 Hz. For the
fluid flow modelling, a 198.times.11.times.31 grid is used
according to some embodiments. The grid is uniform and fine at the
center region then logarithmically coarsened towards the boundaries
in the x and y directions, however the grid is completely uniform
in the z direction. The model is initialized with the irreducible
water saturation. For the electromagnetic modelling, an 84.times.84
grid is used where the inner inversion domain is uniform and the
outer domain is optimized automatically. It is assumed that the
conductivities of shoulder beds are 1.183 S/m. According to some
embodiments, the electromagnetic measurements are acquired using a
tool or tools such as Schlumberger's Array Induction Imager Tools
(AIT) in either a single borehole (such as FIG. 2a and/or 2b, or in
two boreholes as shown in FIGS. 6 and 7.
[0072] Also shown in FIG. 7 is related surface equipment, according
to some embodiments. Sources 702 are deployed in well 630 on a tool
string suspended from a wireline cable 764 from logging truck 762.
Logging truck 762 also is used to control the sources 702 and in
communication with a processing center 750. Similarly receivers 708
are deployed in well 634 on a tool string suspended from a wireline
cable 774 from logging truck 772. Measurements from receivers 708
are recorded and can be processed in one of the logging trucks 762
or 772 or can be transmitted directly to a processing center 750.
Processing center 750 includes one or more central processing units
744 for carrying out the inversion procedures as described herein,
as well as other processing. Processing center 750 also includes a
storage system 742, communications and input/output modules 740, a
user display 746 and a user input system 748. According to some
embodiments, processing center 750 can be included in one or both
of the logging trucks 762 and 772, or may be located in a location
remote from the wellsites. Also shown in FIG. 7 is a surface
electromagnetic array 780 for obtaining electromagnetic data. Thus
the electromagnetic data used in the inversion procedures described
herein can be obtained using various arrangements of equipment
including single well (as in FIGS. 2a and 2b), cross well, surface
to borehole, borehole to surface and surface to surface. Although
wireline logging trucks are shown deploying the sources and
receivers in FIG. 7, when deployed in a marine environment the
sources and receivers may be deployed from a platform or vessel as
well known in the art. Furthermore, according to other embodiments
the sources and/or receivers can be deployed on a drill collar
during part of a drilling operation. In the case of surface
deployed electromagnetic arrays in the marine environment,
according to some embodiments, surface array 780 can be implemented
using deep-towed electric dipole antenna such as used with
Controlled Source Electromagnetics (CSEM) surveys from WesterGeco,
a business unit of Schlumbeger.
[0073] Because the crosswell electromagnetic inversion is not very
sensitive to permeability, we assume the formation permeabilities
are acquired from other independent measurements. FIGS. 8a-h are
vertical cross-sections of spatial distributions of the porosity,
the conductivity, the water saturation and the salt concentration,
according to some embodiments. In FIG. 8a diagram 810 shows the
true distribution of porosity, while in FIG. 8b diagram 812 shows
the inverted distribution of porosity. In FIG. 8c, diagram 814
shows the true distribution for conductivity, while in FIG. 8d
diagram 816 shows the reconstructed distribution for conductivity.
In FIG. 8e, diagram 818 shows the true distribution for water
saturation while in FIG. 8f, diagram 820 shows the reconstructed
distribution for water saturation. In FIG. 8g, diagram 822 shows
the true distribution for salt concentration, while in FIG. 8h,
diagram 824 shows the reconstructed distribution for salt
concentration. In all cases the data is corrupted with 2% Gauss
random noise in the inversion. Compared to the true distribution,
also as shown in FIG. 8, the inverted porosity recovers the true
distribution to some extent. However, the reconstructed
conductivity and water saturation are very close to the true
distributions, and from which we can clearly see the fluid front,
which is the main interest for the fluid movement monitoring.
[0074] According to some embodiments, a second example was
conducted composed of 20 dipping layers, which were categorized
into 4 rock types. Each rock type had different properties such as
permeability, porosity and relative permeability. Similar to the
example shown in FIGS. 6 and 7, water was simultaneously injected
and oil produced in certain zones. It was assumed that the
injection well was operated under bottomhole pressure control and
the production well produces a constant oil rate. The same
inversion procedure was implemented with additional 2% Gauss random
noise added to the data. The spatial porosity distribution was
recovered to some extent, and the water front was clearly seen from
the reconstructed conductivity distribution.
[0075] Whereas many alterations and modifications of the present
invention will no doubt become apparent to a person of ordinary
skill in the art after having read the foregoing description, it is
to be understood that the particular embodiments shown and
described by way of illustration are in no way intended to be
considered limiting. Further, the invention has been described with
reference to particular preferred embodiments, but variations
within the spirit and scope of the invention will occur to those
skilled in the art. It is noted that the foregoing examples have
been provided merely for the purpose of explanation and are in no
way to be construed as limiting of the present invention. While the
present invention has been described with reference to exemplary
embodiments, it is understood that the words, which have been used
herein, are words of description and illustration, rather than
words of limitation. Changes may be made, within the purview of the
appended claims, as presently stated and as amended, without
departing from the scope and spirit of the present invention in its
aspects. Although the present invention has been described herein
with reference to particular means, materials and embodiments, the
present invention is not intended to be limited to the particulars
disclosed herein; rather, the present invention extends to all
functionally equivalent structures, methods and uses, such as are
within the scope of the appended claims.
* * * * *