PEMODELAN AVO DAN SUBSTITUSI FLUIDA: ALGORITMA DAN

PEMROGRAMANNYA DALAM MATLAB

Hilfan Khairy, Maman Hermana
Geosciences and Petroleum Engineering Department
Block 16, Level 3, Universiti Teknologi PETRONAS
Tronoh, 31750, Malaysia
Corresponding author:
hilfan.khairy@petronas.com.my

Abstrak
Respon seismik dari suatu medium ditentukan oleh banyak faktor salah satunya adalah
kandungan fluida di dalam medium tersebut. Dalam lithologi yang sama, ketika kandungan
fluida berubah maka respon seismik ikut berubah terhadap jarak offset. Fenomena ini menjadi
dasar bagi karakterisasi fluida reservoar berdasarkan analisis amplitudo terhadap offset (AVO).
Dalam analisis AVO, perubahan properti medium dapat dikaitkan dengan properti seismik
berdasarkan persamaan Gassmann. Dalam tulisan ini akan dikupas secara detil bagaimana
penggunaan persamaan Gassman dalam pemodelan AVO beserta penulisan algoritmanya yang
ditulis dalam program MATLAB. Dua aproksimasi terhadap persamaan Zoeppritz yaitu Aki
Richard dan Shuey akan dipakai sebagai contoh penggambaran pemodelan AVO. Diharapkan
melalui tulisan ini, pembaca bisa mendapat gambaran bagaimana melakukan pemodelan
substitusi fluida dan AVO. Tentusaja pengembangan algoritma perlu terus dilakukan agar dapat
digunakan langsung untuk melakukan analisa AVO pada data seismik sesunguhnya.

Pendahuluan
Tujuan utama investigasi seismik adalah untuk mendeteksi jenis fluida, mengarakterisasi
reservoar dan memonitor proses recovery. Pengetahuan mengenai hal ini sangat penting terutama
dalam tahapan produksi hidrokarbon. Tentu saja, akan terjadi perubahan parameter seismik pada
saat pergantian jenis fluida atau saat terjadi perubahan saturasinya seperti pada proses water
flooding. Disamping itu, tekanan yang keluar selama proses produksi juga akan menginduksi
terjadinya perubahan parameter seismik. Oleh karenanya, sejauh mana perubahan parameter
tersebut bisa dideteksi adalah kunci untuk memonitor kondisi reservoar.Manfaat lainnya, analisa
perubahan tersebut bisa dimanfaatkan dalam evaluasi fisibilitas survei seismik 4 dimensi.
Salah satu atribut seismik yang bisa digunakan untuk memonitor perubahan jenis fluida adalah
perubahan amplitudo terhadap offset, biasa dikenali sebagai AVO (amplitude versus offset).
Kontras densitas antara minyak, gas dan air akan menghasilkan perbedaan respon impedansi
akustik, bergantung kepada cepat-rambat gelombang seismik. Sejumlah penelitian menunjukkan
terjadi pengurangan kecepatan gelombang pada heavy oil atau tar-sand dengan bertambahnya
temperature (Tosaya et.al., 1987; Wang, 1988) danflooding CO
2
pada batuan tersaturasi minyak
(Wang and Nur, 1989, Hirsche etl.al, 1990). Demikian pula, perbedaan kecepatan gelombang
kompresi didapati pada batuan yang tersaturasi oleh jenis fluida yang berbeda seperti gas, air
ataupun minyak(Nur and Simon., 1969; Gregory., 1976, Khairy et.al., 2002).
Pemodelan amplitudo seismik terhadap offset telah dipopulerkan oleh Ostrander (1984) dan
menjadi semakin dikenal seiringkemajuan teknologi komputer. Dengan menggunakan persamaan
Zoepprits (1919) yang kemudian disederhanakan oleh beberapa peneliti seperti Aki-Richards
(1980) dan Shuey (1985), dapatlah dikalkulasi perubahan amplitudo gelombang seismik sejalan
dengan bertambahnya jarak offset.
Demikian pula halnya dengan pemodelan substitusi fluida (Fluid Replacement
Modeling/FRM).Menggunakan persamaan Gassmanns (1951), dapat dihitung perubahan
kecepatan gelombang seismik terhadap saturasi atau jenis fluida yang berbeda.Integrasi antara
pemodelan AVO dan substitusi fluida bisa menjadi landasan evaluasi seberapa besar atribut
seismik bisa digunakan.
Banyak peneliti yang telah mendiskusikan formulasi substitusi fluida menggunakan persamaan
Gassman (1951), berikut dengan limitasi dan kelebihannya (Smith et.al., 2003) ataupun
kaitannya dengan respon AVO (Li, et.al., 2007). Akan tetapi, jarang sekali dari para peneliti
mengelaborasinya dalam bentuk pemrograman numerik.Tentu saja ini sangat bermanfaat,
dimana dengan mudah dapat dievaluasi model AVO saat saturasi atau jenis fluida dimodelkan.
Dalam makalah ini, akan dielaborasi formulasi dan algoritma pemrograman menggunakan
MATLAB. Diharapkan, mahasiswa atau praktisi bisa memanfaatkannya sebagai sarana
pembelajaran.

Pemodelan AVO

Pemodelan AVO dilakukan berdasarkan formulasi analitik pentulan gelombang pada bidang
batas yang bergantung pada sudut datang (Zoeppritz., 1919). Akan tetapi, formulasi ini kurang
bisa memberikan gambaran yang jelas bagaimana kaitan antara perubahan amplitudo terhadap
parameter fisis seperti kecepatan dan densitas.Kemudian, beberapa aproksimasi dilakukan oleh
sejumlah peneliti untuk menyederhanakannya, sehingga keterkaitan antara amplitudo and
parameter fisis seperti kecepatan gelombang-P, gelombang-S, densitas dan poisson ratio. Dua
daripadanya yang akan dielaborasi pada makalah ini adalah Aki-Richards (1980) dan Shuey
(1985).

Aproksimasi Aki-Richards (1980) bisa dituliskan sebagai berikut:
( ) ( )
2 2 2
sin tan sin R A B C u u u u = + + (1)
dengan
1
2
p
p
V
A
V

| | A
A
= + |
|
\ .
(2)
2 2
2 2
1
4 2
2
p
s s s
p p s p
V
V V V
B
V V V V

| | A
A A
=
|
|
\ .

2
p
p
V
C
V
A
=

( ) / 2
i t
u u u = +

Bila aproksimasi Aki-Richards terkait dengan parameter Vp (kecepatan gelombang-P), Vs
(kecepatan gelombang-S) and (densitas), lain halnya dengan Shuey (1985) dimana ia
mengaitkannya dengan Vp, and o (poisson ratio). Persamaan ini dituliskan sebagai berikut
( ) ( )
2 2 2
sin tan sin R A B C u u u u = + +
dengan

1
2
p
p
V
A
V

| | A
A
= + |
|
\ .


( )
0 2
1
B AA
o
o
(
A
= + (
(



2
p
p
V
C
V
A
=


( )
1 2
2
o o
o
+
=



(4)
(5)
(3)
(6)
(7)
(8)
(9)
(10)
(11)




Sebagai catatatan, kedua pendekatan tersebut bisa dituliskan dalam bentuk yang serupa untuk
sudut datang yang kecil dan rasio Vp terhadap Vs adalah sebagai berikut:
( )
2
sin R A B u u = +

atau
( )
2
sin
p
R R G u u = + (14)
dimana
( )
2
p s
G R R = aproksimasiAki-Richards (15)
9
4
p
G R o
(
= A
(

aproksimasi Shuey (16)
1
2
p
p
p
V
R
V

| | A
A
= +
|
|
\ .

1
2
s
s
s
V
R
V

| | A A
= +
|
\ .

Terminologi R
p
and G biasa dikenal masing-masing sebagai reflektifitas gelombang-P dan
gradient.

Pemodelan Substitusi Fluida
Analisa mengenai bagaimana jenis fluida dan saturasinya dapat mengubah kecepatan gelombang
seismik telah diberikan oleh persamaan Gassman (1951).Formula ini menguantifikasi respon
kecepatan gelombang terhadap perubahan stiffness batuan sebagai akibat terjadinya perubahan
saturasi dan jenis fluida.
Dalam aplikasi praktis, meskipun telah banyak persamaan empiris yang mengaitkan kecepatan
gelombang terhadap saturasi, teori Gassman masih yang terbaik untuk diterapkan. Teori ini
menghubungkan saturasi fluida terhadap porositas, modulus fluida, modulus rangka (frame), dan
( )
2 1
0 0 0
0
1 2
2 1
1
/
/ /
p p
p p
A B B
V V
B
V V
o o o
o
o

A =

= +

A
=
A + A
(12)
(13)
(17)
(18)
modulus matrik penyusun batuan. Secara lengkap persamaan Gassmann dapat ditulis dalam
bentuk:
( )
2
*
*
*
2
1
1
m
sat
f m m
K
K
K K
K
K K K
|
|
| |

|
\ .
= +

+
(19)
denganK
sat
= modulus bulk batuan tersaturasi fluida, K
m
= modulus bulk mineral penyusunnya, K
f

= modulus bulk fluida tersaturasi, K
*
= modulus bulk rangka (frame)and | = porositas.
Bila telah didapatkan harga modulus bulk batuan tersaturasi, maka kecepatan gelombang dapat
dikalkulasi menggunakan persamaan:

4 / 3
sat sat
p
b
sat
s
b
K
V
V

+
=
=

dimana
sat
dan
b
masing-masing adalah modulus geser (shear) dan densitas bulk batuan.Disini,
diasumsikan nilai modulus geser tidak berubah selama pemodelan substitusi fluida.
Untuk melakukan pemodelan dengan baik, tentu saja sejumlah asumsi harus dipenuhi, yaitu:
a. Batuan bersifat homogen dan isotropis. Adanya variasi mineral dalam media harus mempunyai
perbedaan bulk modulus yang tidak besar.
b. Persamaan Gassmann sebaiknya diterapkan pada data seismik berfrekuensi rendah. Pada skala
wireline log, persamaan ini terkadang bisa atau tidak bisa diterapkan (Berrymann, 1999). Namun
untuk batu pasir unconsolidated, persamaan Gassmann masih berkelakuan baik karena frekuensi
sonic dan seismik ada dalam rentang yang sama (Wang, 2000).
c. Semua pori saling terkoneksi dan bisa mengalirkan fluida.
d. Pergerakan relatif antara fluida dan solid matrik diabaikan.

Formulasi
a. Pemodelan AVO 1D
Model refleksi gelombang pada dua bidang batas dapat digambarkan seperti pada Gambar
1.Lapisan paling atas adalah lapisan penutup seperti shale sedangkan dibawahnya adalah
(20)
(21)
reservoar batupasir. Batuan ini dapat berisi minyak, gas ataupun air dengan derajat saturasi
berbeda yang nanti akan dilakukan pemodelan menggunakan persamaan Gassmann seperti telah
diulas sebelumnya.Disebabkan efek fluida dapat memberikan respon cepat rambat gelombang
kompresi berbeda, maka diharapkan pemodelan AVO yang dilakukan bisa mengkarakterisasi
perubahan ini.






Gambar 1. Model Geologi dua lapis dalam pemodelan AVO.

Merujuk pada persamaan Aki-Richards dan Shuey, detail dari keduanya dapat diuraikan sebagai
berikut:

( )
2 1
2 1
/ 2
p p p
p p p
V V V
V V V
A
=
+


( )
2 1
2 1
/ 2
s s s
s s s
V V V
V V V
A
=
+


( )
2 1
2 1
/ 2


A
=
+


( ) ( )
2
1
/ 2 arcsin sin / 2
p
i t i i
p
V
V
u u u u u
| | ( | |
| = + = + ( |
|
|
(
\ . \ .

( )
( )
( )
( )
2 1
2 1
2 1 2 1
1
2 / 2 / 2
p p
p p
V V
A
V V


(

= + (
+ +
(


( ) ( )
2 1
2
2 1
1 / 2
o
B AA
o o
o o
| |

|
= +
|
+
\ .

(23)
(24)
(25)
Sand
V
p2
, V
s2
,
2


Shale
V
p1
, V
s1
,
1


u
i
u
t
(27)
(28)
(26)
( )
( )
2 1
2 1
1
2 / 2
p p
p p
V V
C
V V
(

= (
+
(



b. Properti Mineral Penyusun
Tanpa menghiraukan detil geometri mineral penyusun, modulus bulk batuan dapat diestimasi
melalui model batas (boundary model) seperti Voigt-Reuss-Hill (Hill, 1963). Jika pengamatan
sayatan tipis sampel batuan tersedia, maka ini bisa digunakan sebagai acuan untuk menentukan
komposisi volumetrik mineral penyusun.Dalam kondisi yang berbeda, komposisi mineral secara
sederhana bisa diestimasi melalaui pengukuran wireline log seperti log Gamma-ray. Melalui cara
ini dapat diperkirakan komposisi mineral quartz dan clay. Asumsi yang biasa diterapkan pada
perhitungan komposisi batuan penutup shale ialah tersusun dari 65% mineral clay dan sisanya
adalah jenis mineral lainnya (Pettijohn, 1975). Oleh karenanya, berikut adalah formulasi properti
matriks batuan:
65%
clay shale
V V = (V
sh
didapatkan dari pengukuran log Gamma-ray)
1
quartz clay
V V =

( ) ( )
1
1
2
clay quartz
m clay clay quartz quartz
clay quartz
V V
K V K V K
K K

(
| |
(
= + + +
|
|
(
\ .



m clay clay quartz quartz
V V = +

c. Properti Fluida
Estimasi properti fluida dalam makalah ini merujuk kepada hasil penelitian yang dilakukan oleh
Batzle-Wang (1992).Keterangan detil mengenai proses dan penurunannya bisa dilihat pada
makalah tersebut.
c.1Fluida Gas
Persamaan yang digunakan untuk mengkalkulasi densitas gas adalah:

( )
28.8
273.15
gas
o
GP
ZR T C
~
(
+


dimana
(29)
(31)
(32)
(33)
(34)

( ) ( )
3
4
0.03 0.00527 3.5 0.642 0.007 0.52
pr pr pr pr
Z T P T T E
(
= + + +
(



( ) ( )
1.2
2 2
0.109 3.85 exp 0.45 8 0.56 1/
pr
pr pr
pr
P
E T T
T


(
= +
`
(


)


( )
273.15
94.72 170.75
o
pr
T C
T
G
+
=
+


4.892 0.4048
pr
P
P
G
=

G is specific gravity gas (API), yaitu rasio densitas gas terhadap udara pada suhu 15.6
o
C dan
tekanan atmosfer. Biasanya nilai G ada pada rentang 0.56 sampai dengan 1.8 bergantung kepada
bilangan karbon.Z adalah factor kompressibilitas,R adalah konstanta gas bernilai 8.314, T adalah
temperature dan P adalah tekanan.
Modulus bulk gas dihitung melalui persamaan:

0
1
gas
pr
pr
T
P
K
P
Z
Z P
~
| |
c

|
|
c
\ .

dimana

( )
( )
( )
0 2
5.6 27.1
0.85 8.7exp 0.65 1
2
3.5
pr
pr
pr
P
P
P

(
= + + +

+
+


( ) ( )
3 2
2 2
0.2 1.2
0.03 0.00527 3.5 0.109 3.85
1 1
1.2 0.45 8 0.56 exp 0.45 8 0.56
pr pr
pr
pr pr
pr pr pr pr
Z
T T
P
P P
T T T T
| |
c
= + +
|
|
c
\ .
( | |
| | | |

(
|
+ +
| | ` `
| |
(
|
\ . \ .
) \ . )

c.2Fluida Air (Brine water)
Air adalah jenis fluida yang paling umum dan banyak dijumpai dalam batuan reservoar.Densitas
fluidanya bergantung kepada salinitas, tekanan dan temperature.Semakin tinggi salinitas,
semakin besar densitasnya. Batzle-Wang (1992) memformulasikannya sebagai berikut:
(35)
(36)
(37)
(38)
(39)
(40)
(41)
6
80 3 3300
0.668 0.44 10 300 2400
13 47
br w
T S
S S P PS T
P PS


+ ( | |
= + + + +
` ( |
+
\ . )


2 3 2
6
3 5 2 2
80 3.3 0.00175 489 2 0.016
1 10
1.3 10 0.333 0.002
w
T T T P TP T P
T P P TP

| | + + +
= +
|
|

\ .

Dengan
w
,
br
, S, P, T adalah densitas air murni, densitas brine, salinitas, tekanan dan
temperature.
Sedangkan kecepatan gelombang pada fluida brine berhubungan dengan temperatur, tekanan dan
salinitas sebagai berikut:
( )
( )
2 5 3 2
1.5 2 2
1170 9.6 0.055 8.5 10 2.6 0.0029 0.0476
780 10 0.16 820
br w
V V S T T x T P TP P
S P P S

= + + + +
+


4 3
0 0
i j
w ij
i j
V w T P
= =
=

dengan koefisien w
ij
diuraikan pada Tabel-1.
Tabel 1.Harga koefisien wuntuk menghitung kecepatan fluida air
w
00
=1402.85 w
02
= 3.437 x 10
-3
w
10
=4.871 w
12
= 1.739 x 10
-4
w
20
=-0.04783 w
22
= -2.135 x 10
-6
w
30
=1.487x10
-4
w
32
= -1.455 x 10
-8
w
40
=-2.197 x 10
-7
w
42
= 5.230 x 10
-11
w
01
=1.524 w
03
= -1.197 x 10
-5
w
11
=-0.0111 w
13
= -1.628 x 10
-6

w
21
=2.747 x 10
-4
w
23
= 1.237 x 10
-8
w
31
=-6.503 x 10
-7
w
33
= 1.327 x 10
-10
w
41
=7.987 x 10
-10
w
43
= -4.614 x 10
-13

Setelah diperoleh kecepatan gelombang melalui fluida brine dan densitasnya, maka modulus
Bulk terkait dihitung dengan persamaan:

2
br br br
K V =
c.3Fluida Minyak
Densitas minyak bervasiasi dari sangat ringan (light oil) sampai berat (heavy oil) seperti halnya
tar, bergantung kepada bilangan carbon.Minyak ringan yang terkompresi dapat melarutkan gas
(42)
(43)
(44)
(45)
(46)
didalamnya sehingga dapat mengubah densitas dan modulus kompresibilitas fluida.
Ketergantungan tekanan dan juga temperatur dijabarkan dalam betuk persamaan polinomial oleh
Batzle-Wang (1992) sebagai berikut:

( )( )
( )
( )
2
3 7 4
0 0
1.175
4
0.00277 1.71 10 1.15 3.49 10
0.972 3.81 17.78 10
oil
P P P
T

+ +
=
+ +


( )
1/ 2
1/ 2
1 0
0
0
2096 3.7 4.64 0.0115 4.12 1.08 1 1
2.6
oil
V T P TP

| |
(
= + +
|
(

\ .

atau dalam bentuk API

( ) ( )
1/2
1/2
15450 77.1 3.7 4.64 0.0115 0.36 1
oil
V API T P API TP

= + + +
Dimana T, P,
oil
, dan
0
berturut-turut adalah temperatur, tekanan, densitas minyak dan densitas
minyak referensi yang diukur pada tekanan atmosfir dan suhu 15.6
o
C.
Pada kasus dimana minyak tersaturasi oleh gas (live oil), densitas minyak referensi (
o
) pada
persamaan (47) dan (48) harus diubah menjadi:

( )
0
0
0.0012
G
ob
GR
B

+
=
dan

( )
0
0
1 0.001
ol
G
B R

=
+

dengan

1.175
1/ 2
0
0
0.972 0.00038 2.4 17.8
G
G
B R T

(
| |
( = + + +
|
(
\ .


dimana R
G
adalah gas to oil ratio (GOR).
Setelah didapatkan nilai densitas minyak dan kecepatan gelombangnya, parameter modulus bulk
dapat dihitung kembali menggunakan Persamaan (46).
d. Percampuran Fluida
Setelah diperoleh modulus bulk fluida (gas, minyak dan brine), percampuran antara ketiganya
dapat dikalkulasi menggunakan model Reuss (1929) sebagai berikut:
(47)
(48)
(49)
(50)
(51)
(52)

1
g
w o
f br oil gas
S
S S
K K K K
= + +
Sementara itu, densitas efektif dihitung dengan persamaan:

f w br g gas o oil
S S S = + +
Sebagai contoh, dalam kasus dua fasa antara minyak dan air, kedua persamaan di atas ditulis
ulang sebagai berikut:

( ) 1
1
w w
f br oil
S S
K K K

= +
( ) 1
f w br w oil
S S = +
e. Modulus Rangka
Terdapat sejumlah teknik untuk menghitung modulus ini yaitu dengan pengukuran langsung
pada sampel core,hubungan empiris antara modulus matriks batuan dan porositas (Geertsma dan
Smith., 1961) atau dengan menginversi persamaan Gassmann. Dalam aplikasi praktis, dimana
pemodelan dilakukan pada skala log, modulus ini diperoleh melalui manipulasi persamaan
Gassman sehingga diperoleh persamaan sebagai berikut:



dimana


( ) 1
b m f
| | = +
dengan
m
,
w
,
h
dan S
w
secara berurutan adalah densitas matriks, densitas air, densitas
hidrokarbon dan saturasi air. Kecepatan gelombang yang digunakan pada Persamaan (58)
diekstrak dari pengukuran log.
Terakhir, asumsi yang diterapkan pada pemodelan substitusi fluida, modulus rangka (frame)
dianggap konstan (tidak banyak berubah).

2 2
4
3
sat b p s
K V V
| |
=
|
\ .
(57)
(58)
(59)
(53)
(54)
(55)
(56)
*
1
1
m
sat m
fl
m sat
fl m
K
K K
K
K
K K
K K
|
|
|
|
| |
+
|
|
\ .
=
+
Algoritma
Algoritma integrasi pemodelan kedepan AVO dan substitusi fluida dapat digambarkan pada
bagan alir di bawah:























Input
b

(Pers. 59)
Properti
Fluida
|, Input S
w

Hitung
m
(Pers. 33)
Hitung K
sat
,
sat
(Pers. 58, 21 )
Hitung
f
(Pers. 54)
Input S
w

clay
,
quartz
Hitung
gas

(Pers. 34-38)
Gas
Hitung K
gas

(Pers. 39-41)
Brine
Hitung
br

(Pers. 42-43)
Hitung V
br

(Pers. 44-45)
Hitung K
br

(Pers. 46)
Minyak
P,R,T,G

Hitung
ol

(Pers. 50, 52)
Hitung V
ol

(Pers. 48)

ob

o

(Pers. 51, 52)

R
G
,G,P,
T,
0
Hitung K
oil

(Pers. 46)
Hitung K
f
(Pers. 53)
Input S
w

Input V
p
,
V
s

A
Hitung K
m
(Pers. 32)
K
quartz
, K
clay

P,S,T

Mulai
























Gambar 2. Alir algoritma pemodelan AVO


Konstan selama FRM
A
Hitung K
*

(Pers. 57 )
Hitung model
f
(Pers. 56 )
Model S
w

Hitung model K
f
(Pers. 55 )
Model S
w

Hitung model
b
(Pers. 59 )
Hitung model K
sat
(Pers. 19 )
Hitung model V
p
, V
s
(Pers. 20-21 )
Hitung AV
p
/ V
p
, AV
s
/ V
s
,
A / ,D,E,F (Pers. 23-29 )
FRM BLOCK
AVO BLOCK
Reservoir Rock
V
1p
, V
1s
,

1
, u
Cap Rock
Pemodelan
AVO
R (u) (Pers. 1 )
Aki-Richard
Shuey
R (u) (Pers. 6 )
R
pp
, Grad
(Pers. 14-18)
B

B
Krossplot
a. R(u) vs u
b. R
p
vs G

Selesai
Ubah jenis
fluida atau
saturasinya

Contoh
Dalam tulisan ini, kami mencoba menuliskan kode program dengan menggunakan MATLAB
yang dilengkapi dengan fasilitas GUI (graphical user interface) agar pengguna mudah untuk
melakukan pemodelan perubahan saturasi atau jenis fluida.Akan tetapi, dikarenakan enkripsi
binary low level file GUI, kode listing GUI tidak dapat disertakan. Pada lampiran, hanya
disertakan program inti perhitungan numerik substitusi fluida dan pemodelan kedepan AVO.
Terdapat dua buah input utama pada program yaitu, pemodelan substitusi fluida dan pemodelan
AVO. Hasil dari substitusi fluida akan digunakan sebagai masukan pada properti reservoar di
dalam pemodelan AVO. Karenanya, kita dapat melakukan sejumlah skenario perhitungan model
kecepatan sebagai respon dari perubahan fluida dan saturasinya (Gambar 1-3).Demikian pula
halnya pada pemodelan AVO (Gambar 4-6). Detil input untuk pemodelan substitusi fluida dan
AVO terdapat pada Tabel 2.

Tabel2.Uraian parameter yang digunakan pada program.
Nama Input Nilai Satuan Keterangan
Tekanan (P) 22 MPa FRM
Salinitas (S) 4500 ppm FRM
Temperatur (T) 160 Celcius FRM
Input awal sat. air (ISW) 30 % FRM
Volume shale (VSH) 35 % FRM
Kec. Gel-P reservoar (Vp) 3000 m/s FRM
Kec. Gel-S reservoar (Vs) 1800 m/s FRM
Input model sat. air (NSW) 70 % FRM
Porositas (Por) 20 % FRM
Oil gravity 55 Deg. API FRM
Gas to oil ratio (GOR) 150 L/L FRM
Gas gravity 1.2 API FRM
Kec. Gel-P shale (Vp) 2300 m/s AVO
Kec. Gel-S shale (Vs) 1400 m/s AVO
Densitas shale 1.9 g/cc AVO
Sudut datang 45 Derajat AVO







(a) (b)

(c)
Gambar 3. Contoh parameter input pada pemodelan substitusi fluida. Untuk memodelkan
saturasi penuh air (brine water), masukkan nilai NSW=100% pada untuk sembarang nilai pada
target fluid (a). Gambar (b) dan (c) berturut-turut adalah contoh kalkulasi substitusi fluida 30%
pada target fluida gas (NSM=70%) dan kalkulasi 20% pada fluida minyak (NSW=80%).



Gambar 4.Perhitungan pemodelan substitusi fluida dengan variasi saturasi air 0% sampai
dengan 100%. Pemodelan dilakukan untuk perpindahan fluida dari minyak menuju gas.Grafik
disebelah kanan adalah plot hasil kecepatan gelombang-P dan S terhadap saturasi air.


Gambar 5. Perhitungan pemodelan substitusi fluida dengan variasi saturasi air 0% sampai
dengan 100%. Pemodelan dilakukan untuk jenis fluida sama yaitu, minyak. Grafik disebelah
kanan adalah plot hasil kecepatan gelombang-P dan S terhadap saturasi air.

Gambar 6.Parameter pemodelan AVO. Parameter input bagi reservoar dilakukan pada blok FRM
yang memodelkan perubahan saturasi fluida minyak dengan saturasi air: 0%, 40%, 80% dan
100%. Grafik sebelah kanan adalah crossplot antara R
pp
dengan sudut datang menggunakan
persamaan Aki-Richard.

Gambar 7.Parameter pemodelan AVO. Parameter input bagi reservoar dilakukan pada blok FRM
yang memodelkan perubahan jenis fluida dari minyak menuju gas dengan saturasi air: 0%, 40%,
80% dan 100%. Grafik sebelah kanan adalah krosplot antara R
pp
dengan sudut dating
menggunakan persamaan Aki-Richard.
Sw=0%
Sw=40%
Sw=80%
Sw=100%
Sw=0%
Sw=80%
Sw=100%
Setelah memasukkan parameter input, user diminta untuk menentukan jenis fluida awal dan
akhir yang difasilitasi oleh menu pop-up dibagian atas. Dengan mengeksekusi perhitungan
substitusi fluida, maka proses pemodelan mulai dilakukan. Untuk perhitungan pemodelan AVO,
diperlukan data kecepatan dan densitas lapisan penutup (seal rock) dan reservoarnya.Properti
reservoar diperoleh dari perhitungan substitusi fluida sebelumnya. Sedangkan lapisan penutup
memiliki input yang terpisah. Pemodelan AVO dapat dipillih salah satu diantara aproksimasi
Shuey ataupun Aki-Richard (Gambar 8).

Gambar 8.Contoh pemodelan AVO menggunakan aproksimasi Shuey dan Aki-Richard.Kolom
bagian kanan bawah adalah koefisien perhitungan keduanya.
Atribut AVO lain yang dapat diplot adalah gradien terhadap offset. Plot kurva ini perlu dilakukan
untuk melihat jenis kelas pada AVO dan perubahannya bila terjadi variasi saturasi fluida selama
substitusi (Gambar 9).

Gambar 9. Plot gradient (G) terhadap intercept (R
pp
) pada variasi saturasi air berbeda: 0%, 50%
dan 100% (arah anak panah).
Shuey
Aki-Richard
Diskusi dan Kesimpulan
Formulasi, algoritma substitusi fluida dan pemodelan AVO telah diuraikan secara detil. Pada
prinsipnya pemodelan dan komputasi numerik bisa dilakukan dengan menggunakan bahasa
pemrograman apapun. Hanya saja, hal terpenting yang harus dilakukan adalah selalu memeriksa
konversi antar parameter dan konsistensi penggunaan satuan pada setiap tahap perhitungan.
Terdapat tiga bagian yang mesti digarisbawahi dalam perhitungan substitusi fluida yaitu properti
matriks, fluida dan rangka. Jika didapatkan analisa sayatan tipis, maka informasi detil volumetrik
dan jenis mineral harus diakomodir dalam modulus efektifnya. Sebagai tambahan, dalam
makalah ini densitas bulk batuan diturunkan melalui perhitungan porositas dan saturasi fluida
tersaturasi. Apabila digunakan data densitas dari log, maka proses perhitungan ini dapat
dihilangkan. Tentu saja, dalam menggunakan persaman Gassmann untuk kalkulasi substitusi
fluida harus terpenuhi dahulu asumsi-asumsi yang melekat padanya.

Pustaka
Aki, K and Richards, P. G., 1980., Quantitative seismology, Theory and methods: W.H. Freeman
& Co.
Berrymann, J.G., 1999.,Origin of Gassmanns equation., Geophysics, 64, 1627-1629.
Batzle, M and Wang, Z., 1992.,Seismic properties of pore fluids., Geophysics, 57, 1396-1404.
Gassmann, F., 1951., Uber die Elastizitat Poroser Medien: Vier. der Natur. Gesellschaft in
Zurich, 96, 1-23.
Geertsma, I., Smith, D.C., 1961.,Some aspect of elastic wave propagation in fluid saturated
porous rock, Geophysics, 26, 169-181.
Gregory, A.R, 1976.,Fluid saturation effects on dynamic elastic properties of sedimentary rocks.,
Geophysics, 41, 895-921.
Hill, R., 1963, Elastic properties of reinforced solids: Some theoretical principles., J. Mech.
Phys. Solids, 11, 357372
Hirsche, W.K., Wang, Z., Sedgwick, G., 1990.,Seismic monitoring of miscible/immiscible floods:
Part II,Model studies and field test., 60th Ann.Internat, Mtg. Soc. Expl. Geophysics.,
Expanded Abstract, 233-236.
Khairy, H., Nurhandoko, B.E.B andSantoso, D., 2002.,Karakterisasi gelombang seismik pada
reservoar batupasir dalam medium tersaturasi fluida., Proceeding of 27th HAGI annual
meeting, Indonesia.
Li, Y., Downton, J and Xu, Y., 2007., Practical aspect of AVO modeling., Leading Edge, March,
295-311.
Nur A and Simmons, ,D., 1969., The effect of saturation on velocity in low porosity rocks., Earth
Plan.Sci.Lett., 7., 183-193
Ostrader, W.J, 1984, plane wave reflection coefficients for gas sand at non normal angles of
incident., Geophysics, 49, 1637-1648.
Pettijohn, F. J., 1975.,Sedimentary Rocks., 3rd edition. Harper and Row Publishers, New York.
628.
Reuss, A., 1929.,Berechnung der Fliegrenze von Mischkristallen., Angew. Mathem. U Mech., 9,
49-58.
Shuey ., 1985., A simplification of the Zeoppritz equations., Geophysics 50, 609-614.
Smith, T.M., Sondelgeld, C.H and Rai, C.S., Gassmann fluid substitution: A tutorial.,
Geophysics, 68., 430-440
Tosaya, C., Nur, A., Vo-Thanh, D., 1987.,Laboratory seismic methods for remote monitoring of
thermal EOR., SPE Res.Eng., 2, 235-242
Wang Z., 1988.,Wave velocities in hydrocarbons and hydrocarbon saturated rocks, with
applications to EOR monitoring: Ph.D thesis., Stanford University
Wang Z., 2000.,The Gassmann equation revisited: Comparing laboratory data with Gassmanns
Predictions., Geophysics reprint series, V.3
Wang Z and Nur A., 1989.,Effect of CO2 flooding on the velocities in rocks with hydrocarbons.,
SPE Res.Eng., 4, 429-436
Zoeppritz, K., 1919, Erdbebenwellen VIIIB, On the reflection and propagation of seismic waves.,
Gottinger Nachrichten, I,66-84.
















Lampiran
Kode MATLAB untuk pemodelan AVO dan substitusi fluida
% Program code to calculate Fluid Replacement and AVO modeling
% Created by Hilfan Khairy
% Universiti Teknologi PETRONAS, April 2011
%
% %%%%%Input of the program are:%%%%%%
% ------FRM BLOCK--------
% Pressure (P) in MPa
% Salinity (S) in ppm
% Temperature (Ta) in deg C
% Initial water saturation (ISW) in fraction
% Volume of shale (VSH) in fraction
% Insitu Vp velocity (Vp) in m/s
% Insitu Vs velocity (Vs) in m/s
% Porosity (Por) in fraction
% Target or New water saturation (NSW) in fraction
% Gas gravity (G) in API
% Gas to oil ratio/GOR (Rg) in L/L
% Oil gravity (OD) in deg API
% -------AVO BLOCK-------
% Shale P-velocity (Vpseal) in m/s
% Shale S-velocity (Vsseal) in m/s
% Shale density (Rhoseal) in g/cc
% Incident angle (IncAngle)

clc;

%--------FRM------------
Kq=36.6*10^(9); %in Pa
Kc=20.9*10^(9); %in Pa
Rho_c=2.58; %in g/cc
Rho_q=2.65;
R=8.314;

VSH=0.35;
Por=0.2;
OD=55; %deg API
P=22;
Ta=160;
Rg=160; %Gas-Oil ratio
G=1.2;% API
S=4500;
ISW=0.3;
NSW=0.7; % to model full saturation, type NSW=100% for any initial fluid type
Vp=3000;
Vs=1800;
Influid=1; %initial fluid: (1) for oil and (2) for gas
Nefluid=1; %target fluid: (1) for oil and (2) for gas

T=Ta+273.15;
OD=141.5/(OD+131.5); % g/cc
Vc=0.65*VSH;
Vq=1-Vc;
Kvoigt=Vc*Kc+Vq*Kq;
Kreuss=1/(Vc/Kc+Vq/Kq);
Km=(Kvoigt+Kreuss)*0.5;
Rho_m=Vc*Rho_c+Vq*Rho_q;

%%% Calculating Frame Modulus

%water properties
S=S/1000000;
%----water coefficient------
w(1,1)=1402.85;
w(2,1)=4.871;
w(3,1)=-0.04783;
w(4,1)=1.487*10^(-4);
w(5,1)=-2.197*10^(-7);
w(1,2)=1.524;
w(2,2)=-0.0111;
w(3,2)=2.747*10^(-4);
w(4,2)=-6.503*10^(-7);
w(5,2)=7.987*10^(-10);
w(1,3)=3.437*10^(-3);
w(2,3)=1.739*10^(-4);
w(3,3)=-2.135*10^(-6);
w(4,3)=-1.455*10^(-8);
w(5,3)=5.23*10^(-11);
w(1,4)=-1.197*10^(-5);
w(2,4)=-1.628*10^(-6);
w(3,4)=1.237*10^(-8);
w(4,4)=1.327*10^(-10);
w(5,4)=-4.614*10^(-13);
%----------------------------

Rho_w=1+(-80*Ta-3.3*Ta*Ta+0.00175*Ta*Ta*Ta+489*P-2*Ta*P+0.016*Ta*Ta*P-
1.3*Ta*Ta*Ta*P*10^(-5)-0.333*P*P-0.002*Ta*P*P)*10^(-6);
Rho_br=Rho_w+S*(0.668+0.44*S+(300*P-2400*P*S+Ta*(80+3*Ta-3300*S-
13*P+47*P*S))*10^(-6)); %brine density g/cc
Vw=0;
for i=1:5
for j=1:4
Vw=Vw+w(i,j)*(Ta^(i-1))*P^(j-1);
end;
end;

V_br=Vw+S*(1170-9.6*Ta+0.055*Ta*Ta-8.5*Ta*Ta*Ta*10^(-5)+2.6*P-0.0029*Ta*P-
0.0476*P*P)+(780-10*P+0.16*P*P)*S^(1.5)-820*S*S; %brine velocity in m/s
K_br=Rho_br*(V_br^2)*10^3; %brine bulk modulus in Pa, density is changed into
kg/m3

if Influid==1
%oil properties
B=0.972+0.00038*(2.495*Rg*sqrt(G/OD)+Ta+17.8)^(1.175);
Rho_ol=OD/(B*(1+0.001*Rg));%in g/cc
Rho_ob=(OD+0.0012*G*Rg)/B;%in g/cc
Rho_oil=(Rho_ob+(0.00277*P-1.71*P*P*P*0.0000001)*(Rho_ob-1.15)*(Rho_ob-
1.15)+3.49*P*0.0001)/(0.972+3.81*0.0001*(Ta+17.78)^1.175); % oil density in
g/cc
V_oil=2096*sqrt(Rho_ol/(2.6-Rho_ol))-3.7*Ta+4.64*P+0.0115*(4.12*sqrt(-
1+1.08/Rho_ol)-1)*Ta*P; % oil velocity m/s

K_oil=Rho_oil*(V_oil*V_oil)*10^3; %oil bulk modulus in Pa

%mixing fluids
Rho_f=ISW*Rho_br+(1-ISW)*Rho_oil;

K_f=K_br*K_oil/(ISW*K_oil+(1-ISW)*K_br);
else
%gas properties
Tpr=T/(94.72+170.75*G);
Ppr=P/(4.892-0.4048*G);
E=0.109*(3.85-Tpr)*(3.85-Tpr)*exp(-(0.45+8*(0.56-
1/Tpr)^2)*(Ppr^(1.2))/Tpr);
Z=((0.03+0.00527*(3.5-Tpr)^3)*Ppr)+(0.642*Tpr-0.007*Tpr^(4)-0.52)+E;
Rho_gas=28.8*G*P/(Z*R*T); %gas density in g/cc

Gamma=0.85+5.6/(Ppr+2)+27.1/((Ppr+3.5)*(Ppr+3.5))-8.7*exp(-0.65*(Ppr+1));
last=-1.2*((Ppr^(0.2))/Tpr)*(0.45+8*(0.56-1/Tpr)^2)*exp(-(0.45+8*(0.56-
1/Tpr)^2)*(Ppr^(1.2))/Tpr);
DzDpr=0.03+0.00527*(3.5-Tpr)^3+0.109*(3.85-Tpr)^2*last;

K_gas=(P*10^(6))*Gamma/(1-(Ppr/Z)*DzDpr); %bulk modulus of gas in Pa

%mixing fluids
Rho_f=ISW*Rho_br+(1-ISW)*Rho_gas;

K_f=K_br*K_gas/(ISW*K_gas+(1-ISW)*K_br);
end;

Rho_b=Por*Rho_f+(1-Por)*Rho_m; % g/cc
IKS=Rho_b*(Vp*Vp-4*Vs*Vs/3)*10^(3); %initial Ksat in Pa will be changed into
GPa for display
IMS=Rho_b*Vs*Vs*10^(3); %initial Miu_sat in Pa
Ko=(IKS*(Por*Km/K_f+1-Por)-Km)/(Por*Km/K_f+IKS/Km-1-Por); %Frame bulk modulus
remain constant during modeling

%%Calculating New Fluid Modeling
if Nefluid==1
%oil properties
B=0.972+0.00038*(2.495*Rg*sqrt(G/OD)+Ta+17.8)^(1.175);
Rho_ol=OD/(B*(1+0.001*Rg));%in g/cc
Rho_ob=(OD+0.0012*G*Rg)/B;%in g/cc
Rho_oil=(Rho_ob+(0.00277*P-1.71*P*P*P*0.0000001)*(Rho_ob-1.15)*(Rho_ob-
1.15)+3.49*P*0.0001)/(0.972+3.81*0.0001*(Ta+17.78)^1.175); % oil density in
g/cc
V_oil=2096*sqrt(Rho_ol/(2.6-Rho_ol))-3.7*Ta+4.64*P+0.0115*(4.12*sqrt(-
1+1.08/Rho_ol)-1)*Ta*P; % oil velocity m/s

K_oil=Rho_oil*(V_oil*V_oil)*10^3; %oil bulk modulus in Pa

% New mixing fluids
NRho_f=NSW*Rho_br+(1-NSW)*Rho_oil;

NK_f=K_br*K_oil/(NSW*K_oil+(1-NSW)*K_br); %New fluid mixing modulus

else
%gas properties
Tpr=T/(94.72+170.75*G);
Ppr=P/(4.892-0.4048*G);
E=0.109*(3.85-Tpr)*(3.85-Tpr)*exp(-(0.45+8*(0.56-
1/Tpr)^2)*(Ppr^(1.2))/Tpr);
Z=((0.03+0.00527*(3.5-Tpr)^3)*Ppr)+(0.642*Tpr-0.007*Tpr^(4)-0.52)+E;
Rho_gas=28.8*G*P/(Z*R*T); %gas density in g/cc

Gamma=0.85+5.6/(Ppr+2)+27.1/((Ppr+3.5)*(Ppr+3.5))-8.7*exp(-0.65*(Ppr+1));
last=-1.2*((Ppr^(0.2))/Tpr)*(0.45+8*(0.56-1/Tpr)^2)*exp(-(0.45+8*(0.56-
1/Tpr)^2)*(Ppr^(1.2))/Tpr);
DzDpr=0.03+0.00527*(3.5-Tpr)^3+0.109*(3.85-Tpr)^2*last;

K_gas=(P*10^(6))*Gamma/(1-(Ppr/Z)*DzDpr); %bulk modulus of gas in Pa

% New mixing fluids
NRho_f=NSW*Rho_br+(1-NSW)*Rho_gas;

NK_f=K_br*K_gas/(NSW*K_gas+(1-NSW)*K_br);
end;

NRho_b=Por*NRho_f+(1-Por)*Rho_m; %g/cc

% Calculating new Ksat, new Vp and new Vs
NKS=Ko+(1-Ko/Km)*(1-Ko/Km)/((Por/NK_f)+(1-Por)/Km-Ko/(Km*Km));
NVp=sqrt((NKS+4*IMS/3)*0.001/NRho_b);
NVs=sqrt(IMS*0.001/NRho_b);

disp(['New Vp = ',num2str(NVp)]);
disp(['New Vs = ',num2str(NVs)]);


%-------AVO Modeling--------
Vpres=NVp;
Vsres=NVs;
Rhores=NRho_b*1000; %in kg/m3
IncAngle=45;
Choise=1 % (1) Shuey and (2) Aki-Richard (3) Plot Gradient-Intercept

Vpseal=2300;
Vsseal=1400;
Rhoseal=1.9*1000; %in kg/m3

tou1=((Vpseal/Vsseal)*(Vpseal/Vsseal)-
2)/(2*((Vpseal/Vsseal)*(Vpseal/Vsseal))-2);
tou2=((Vpres/Vsres)*(Vpres/Vsres)-2)/(2*((Vpres/Vsres)*(Vpres/Vsres))-2);
tou=(tou1+tou2)/2;
deltou=tou2-tou1;
delAlpha=(Vpres-Vpseal);
delBeta=(Vsres-Vsseal);
delRho=(Rhores-Rhoseal);
beta=(Vsres+Vsseal)/2;
alpha=(Vpres+Vpseal)/2;
rho=(Rhores+Rhoseal)/2;
TranAngle=asind((Vpres/Vpseal)*sind(IncAngle));
teta=(IncAngle+TranAngle)/2;
angle=[0:teta];

if Choise==1
%Shuey
Bo=(delAlpha/alpha)/(delAlpha/alpha+delRho/rho);
Ao=Bo-2*(1+Bo)*(1-2*tou)/(1-tou);
A=0.5*(delAlpha/alpha+delRho/rho);
B=A*Ao+deltou/((1-tou)*(1-tou));
C=0.5*delAlpha/alpha;
Rteta1=A+B*sind(angle).*sind(angle)+C*((tand(angle).*tand(angle))-
(sind(angle).*sind(angle)));

plot(angle,Rteta1,'*r');
xlabel('Angle(degree)');
ylabel('Rp');
hold on;
grid on;
ishold;
axis auto;

elseif Choise==2
%Aki-Richard
A=0.5*((delAlpha/alpha)+(delRho/rho));
B=(0.5*(delAlpha/alpha)-(4*beta*beta*delBeta/(alpha*alpha*beta)));
C=0.5*delAlpha/alpha;
Rteta2=A+B*sind(angle).*sind(angle)+C*((tand(angle).*tand(angle))-
(sind(angle).*sind(angle)));

plot(angle,Rteta2,'^b');
xlabel('Angle(degree)');
ylabel('Rp');
hold on;
grid on;
ishold;
axis auto;

elseif Choise==3
Rpp=0.5*(delAlpha/alpha+delRho/rho);
Rss=0.5*(delBeta/beta+delRho/rho);
Grad=Rpp-2*Rss;
plot(Rpp,Grad,'ok');
xlabel('Rpp');
ylabel('G');
axis([-1 1 -1 1]);
grid on;
end;

%% End of Program