Walkthrough: Calculate Vp Using Fluid Substitution in Python
Fluid substitution is an important tool of the seismic rock physics analysis (e.g., AVO, 4D analysis), this method gives the geoscientist some insight on identification of fluid type that contained in reservoir. Fluid substitution is usually calculated using Gassmann’s equation which it has been discussed by numerous famous researchers on geoscience such as Mike Batzle, Zhijin Wang, James G. Berryman, and even the co-founder of Hampson-Russel (Brian Russel) has analyzed the formula of the Gassmann fluid Substitution.
Basically, the main target of fluid substitution is to calculate the seismic velocities (Vp and Vs) and density at particular condition of the reservoir, for example, pressure, temperature, porosity, mineral type and water salinity. The formula is also considering the content of fluid of the reservoir such as 80% of oil and 20% of water, or 90% water and 10% oil (the scenario after Enhanced Oil Recovery is being done). To estimate the value of seismic velocities, it could be calculated by using rock moduli (K and µ) and density with the assumption the rock is isotropic. These can be formulated as:
Whereas:
K is a moduli of incompressibility (Pascal)
µ is a shear moduli (Pascal)
𝜌 is a density of the rock (g/cm3)
Gassmann’s Fluid Substitution
Gassmann’s equation describes the relation between the bulk moduli of a rock with its pore, frame, and fluid properties that are contained with the assumption the rock’s material property is isotropy and homogeneous. The equation works best at low frequency (i.e. seismic data).
How to Do It
To perform Gassmann’s Fluid Substitution, we need to estimate bulk moduli of the frame, matrix and pore fluid.
To calculate the moduli of the frame, usually Vp, Vs and 𝜌 are given at initial state. These can be used to calculate Ksat and µsat. The formula comes from the modification of Vp and Vs formula that have been stated above. And we convert it into python function
def msat(vs, rho):
return rho*(vs**2)def ksat1(vp, vs, rho):
Ksat1 = ((vp**2)*rho) - 4/3*(msat(vs, rho))
return np.round(Ksat1, 3)
One of the methods to calculate bulk modulus of the matrix can be calculated using Voigt-Reuss-Hill average which averaging the Voigt upper bound and the Reuss lower bound. This can be expressed as:
Then this formula we transform into python function
def km(por, fq, fc, kq, kc):
kv = ((fq)*kq) + (((1-por)*fc)*kc)
kr = ((fq)/kq) + (((1-por)*fc)/kc))-1
km = (kv + kr)/2
return np.round(km, 3)Reuss lower bound can be also used to calculated the moduli of fluid mixtures such as:
And then we back to transform this formula into Python function
def kfl(sw, khc, kw):
kfl = ((sw/kw) + ((1-sw)/khc))**-1
return np.round(kfl, 3)where Kfl is the effective bulk modulus of the fluid mixture, K denotes the bulk moduli of the individual gas and fluid phases, and S represents their saturations.
When all is set and done, we could perform the Gassmann’s Fluid substitution. It can be explained into the following equation:
def ksat2(ksat1, k0, kfl1, kfl2, por):
a = (ksat1/(k0 - ksat1)) - (kfl1/(por*(k0 - kfl1))) +(kfl2/(por*(k0 - kfl2)))
ksat2 = (a*(k0))/(1+a)
return np.round(ksat2, 3)where Ksat is the effective bulk modulus of the rock with pore fluid, Kmineral is the bulk modulus of mineral material making up rock, Kfl is the effective bulk modulus of pore fluid, Ø is the porosity, and µsat is the effective shear modulus of rock with pore fluid. The numbering ((1) and (2)) Describing the condition of the rocks with different fluid properties.
The next step is to calculated the density at fluid 2:
def rhofluid(sw, so, rhoo, rhow):
rho1 = ((1-so) * rhow) + (so *rhoo)
rho2 = ((sw) * rhow) + ((1 - sw) * rhoo)
rho1 = np.round(rho1, 3)
rho2 = np.round(rho2, 3)
return rho1, rho2def rhosat(rho, por, rhow, sw, so, rhoo):
rho1, rho2 = rhofluid(sw, so, rhoo, rhow)
rhosat = rho+(por*(rho2 - rho1))
return np.round(rhosat, 3)
And finally, we could compute the velocity of the substituted fluid with the first formula.
def vp(ksat2, vs, rhosat):
vp = ((ksat2 + (4/3*(msat(vs, rhosat)))) / (rhosat))**0.5
return np.round(vp, 3)Why don’t we give it a try in a simple case study?
Imagine a sandy reservoir with the following elastic properties :
𝜌b= 2.13 g/cm3, Vp = 3.349 km/s & Vs = 1.836 km/s
with mineral composition consisted of 70% quartz and 30% calcite.The reservoir was found oil bearing, with
Ø= 30%, and So = 70%.
Elastic properties of the fluid were also derived whereas :
𝜌water = 1 g/cm3 and 𝜌oil= 0.8 g/cm3.
Some additional information:
Kquartz = 36 GPa | Kcalcite = 75 GPa | Koil = 1.6 GPa | Kwater = 2.38 GPa
Calculate a new Vp with a new condition where the reservoir are completely filled by water (Sw = 100%)
Vp = 3.349
Vs = 1.836Sw = 1
So = 0.7Por = 0.3Fq = 0.7
Fc = 0.3
Kq = 36
Kc = 75
Ko = 1.6
Kw = 2.83Rho = 2.13
Rhoo = 0.8
Rhow = 1
So, first of all we calculate The Ksat and µsat of the reservoir,
Ksat1 = ksat1(Vp, Vs, Rho)
Next, we estimate the bulk moduli of the fluid before and after fluid substitution,
Kfl1 = kfl((1-So), Ko, Kw)
Kfl2 = kfl(Sw, Ko, Kw)
After that, we evaluate the Bulk modulus of the Mineral that contained in the rocks. This can be estimate from Hill formula, this formula include Voight-Reuss Bounds calculation too,
K0 = km(Por, Fq, Fc, Kq, Kc)
After all of the required parameters have been calculated, then we can proceed to Gassmann’s Fluid Substitution,
Ksat2 = ksat2(Ksat1, K0, Kfl1, Kfl2, 0.3)
To calculate the new Vp, we need to calculate the new saturated density first,
Rhosat = rhosat(Rho, Por, Rhow, Sw, So, Rhoo)
Finally, we could recompute the fluid substituted velocities,
Vp = vp(Ksat2, Vs, Rhosat)
And the new Vp after fluid substitution is 3.44 km/s. If you interest with this subject, you may can play with the algorithm through our github page.
To make it more interesting we turn the calculation above into a software. We introduce you,Reflow!. It’s built using PyQt5 and CSS. The video below shown how Reflow works to calculate Vp using Gassmann’s Fluid Substitution
