This is a little investigation on the reliability of the values given in two papers written by John Castagna et al. (1997, 1998 -- see references below).
In these papers I have found Vp, Vs and density for an AVO class 4 interface; I have always put them together with those found for the other 3 AVO classes in Hilterman, 2001 (see below for references) but then I realized there was no porosity assigned to the class 4 sand.
I started to look into it and look for the porosity that could be assigned to this sand. But then I realized that depending on the method there is no unique solutions, and that suggests that one of the inputs (either velocities or density) are wrong.
My suggestion it this: use the class 3 sand velocities and density in place of class 4 sand. With the same class 4 overburden shale we still get a class 4 response, and none of the absurdities that derive from trying to make sense of the parameters used by Castagna.
REFERENCES
THINGS THAT CAN BE DONE
The next cell block does the standard imports and defines a few other functions. This makes it a rather lenghty block which should be kept hidden.
The additional functions are:
gassmann
to implement Gassmann's equationgassmann_phi
to get porosity from elastic parameters using Gassmann's equationbulk
to get bulk modulus from Vp, Vs and densityshuey
to calculate the reflectivity variation with anglevels
to calculate velocities from density, bulk and shear moduliclassref
to plot AVO classeshertzmindlin
and softsand
to implement the Soft Sand rock physics model%matplotlib inline
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 100
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import ipywidgets as ipw
from IPython.display import display
def percdiff(start, end):
return (end-start)/start
def find_nearest(a, a0):
idx = np.abs(a - a0).argmin()
return idx, a[idx]
def bulk(vp, vs, rho):
# converts density to SI (kg/m3)
D = rho*1e3
K = D*vp**2 - 4/3*D*vs**2
return K/1e9
def gassmann(vp1, vs1, rho1, rhof1, kfluid1, rhof2, kfluid2, kmin, phi):
# convert density to kg/m3 and elastic moduli to Pa
d1 = rho1*1e3
df1 = rhof1*1e3
df2 = rhof2*1e3
k0 = kmin*1e9
kf1 = kfluid1*1e9
kf2 = kfluid2*1e9
d2 = d1 - phi * df1 + phi * df2
mu1 = d1 * vs1**2
k1 = d1 * vp1**2 - (4/3) * mu1
kd = (k1 * (phi*k0 / kf1 + 1-phi) - k0) / (phi*k0 / kf1 + k1 / k0 - 1-phi)
mu2 = mu1
with np.errstate(divide='ignore', invalid='ignore'):
k2 = kd + (1 - kd/k0 )**2 / (phi/kf2 + (1-phi)/k0 - kd/k0**2)
vp2 = np.sqrt((k2 + 4/3*mu2) / d2)
vs2 = np.sqrt(mu2 / d2)
return vp2, vs2, d2/1e3, k2/1e9, kd/1e9
def gassmann_phi(ksat1, ksat2, kf1, kf2, kmin):
a = (kmin-ksat1)*(kmin-ksat2)*(kf1-kf2)
b = (kmin-kf1)*(kmin-kf2)*(ksat1-ksat2)
return a / b
def hertzmindlin(K0, G0, sigma, phi_c=0.4, Cn=8.6, f=1):
sigma0 = sigma / 1e3 # converts pressure in same units as solid moduli (GPa)
pr0 = (3*K0-2*G0) / (6*K0+2*G0) # poisson's ratio of mineral mixture
Khm = (sigma0*(Cn**2*(1 - phi_c)**2*G0**2) / (
18*np.pi**2 * (1 - pr0)**2))**(1/3)
Ghm = ((2+3*f-pr0*(1+3*f)) / (5*(2-pr0))) * (
(sigma0 * (3 * Cn**2 * (1 - phi_c)**2 * G0**2) / (
2 * np.pi**2 * (1 - pr0)**2)))**(1/3)
return Khm, Ghm
def softsand(K0, G0, phi, sigma, phi_c=0.4, Cn=8.6, f=1):
Khm, Ghm = hertzmindlin(K0, G0, sigma, phi_c, Cn, f)
Kdry = -4/3 * Ghm + (((phi / phi_c) / (Khm + 4/3 * Ghm)) + (
(1 - phi / phi_c) / (K0 + 4/3 * Ghm)))**-1
gxx = Ghm / 6 * ((9 * Khm + 8 * Ghm) / (Khm + 2 * Ghm))
Gdry = -gxx + ((phi / phi_c) / (Ghm + gxx) + (
(1 - phi / phi_c) / (G0 + gxx)))**-1
return Kdry, Gdry
def vels(kdry, gdry, kmin, rho0, kfluid, rhof, phi):
# convert density to kg/m3 and elastic moduli to Pa
d0 = rho0*1e3
df = rhof*1e3
kd = kdry*1e9
gd = gdry*1e9
kf = kfluid*1e9
k0 = kmin*1e9
rho = d0 * (1 - phi) + df * phi
with np.errstate(divide='ignore', invalid='ignore'):
ksat = kd + (1 - kd/k0)**2 / (phi/kf + (1-phi)/k0 - kd/k0**2)
vp = np.sqrt((ksat+4/3*gd)/rho)
vs = np.sqrt(gd/rho)
return vp, vs, rho/1e3, ksat/1e9
def shuey(vp1, vs1, rho1, vp2, vs2, rho2, theta, approx=True, terms=False):
a = np.radians(theta)
dvp = vp2-vp1
dvs = vs2-vs1
drho = rho2-rho1
vp = np.mean([vp1, vp2], axis=0)
vs = np.mean([vs1, vs2], axis=0)
rho = np.mean([rho1, rho2], axis=0)
R0 = 0.5*(dvp/vp + drho/rho)
G = 0.5*(dvp/vp) - 2*(vs**2/vp**2)*(drho/rho+2*(dvs/vs))
F = 0.5*(dvp/vp)
# if angles is an array
if a.size>1:
R0 = R0.reshape(-1,1)
G = G.reshape(-1,1)
F = F.reshape(-1,1)
if approx:
R = R0 + G*np.sin(a)**2
else:
R = R0 + G*np.sin(a)**2 + F*(np.tan(a)**2-np.sin(a)**2)
if terms:
return R, R0, G
else:
return R
def classref(near=5, far=30, mx=.6, plot_brine=False, plot_colorzones=True):
tmp_shl = np.array([[3094, 1515, 2.40, 0],
[2643, 1167, 2.29, 0],
[2192, 818, 2.16, 0],
[3240, 1620, 2.34, 0]])
tmp_ssg = np.array([[4050, 2526, 2.21, .2],
[2781, 1665, 2.08, .25],
[1542, 901, 1.88, .33],
[1650, 1090, 2.07, .156]])
tmp_ssb = np.array([[4115, 2453, 2.32, .2],
[3048, 1595, 2.23, .25],
[2134, 860, 2.11, .33],
[2590, 1060, 2.21, .156]])
avocl = ['CLASS1', 'CLASS2', 'CLASS3', 'CLASS4']
logs = ['VP', 'VS', 'RHO', 'PHI']
shl = pd.DataFrame(tmp_shl, columns=logs, index=avocl)
ssg = pd.DataFrame(tmp_ssg, columns=logs, index=avocl)
ssb = pd.DataFrame(tmp_ssb, columns=logs, index=avocl)
opttxt = dict(weight='bold', ha='left', va='center')
mrkg = {'ms': 10, 'mew': 2, 'ls': 'none'}
mrkb = {'ms': 10, 'mew': 2, 'ls': 'none', 'mfc':'none'}
mrk_sel = {'marker': '*', 'mec': 'k', 'mfc': 'white', 'ms': 16, 'ls': 'none', 'mew': 2}
angs = np.array([near, far])
tmp = ['C0', 'C1', 'C2', 'C3']
cc = dict(zip(avocl, tmp))
tmp = ['s', 'P', 'v', '^']
mm = dict(zip(avocl, tmp))
f, ax = plt.subplots(constrained_layout=True)
ax.axhline(0, color='k', lw=3)
ax.axvline(0, color='k', lw=3)
for i, sh in shl.iterrows():
vpsh, vssh, dsh = sh['VP'], sh['VS'], sh['RHO']
vpb, vsb, db = ssb.loc[i, 'VP'], ssb.loc[i, 'VS'], ssb.loc[i, 'RHO']
vpg, vsg, dg = ssg.loc[i, 'VP'], ssg.loc[i, 'VS'], ssg.loc[i, 'RHO']
Ab, Ib, Gb = shuey(vpsh, vssh, dsh, vpb, vsb, db, angs, terms=True)
Ag, Ig, Gg = shuey(vpsh, vssh, dsh, vpg, vsg, dg, angs, terms=True)
ax.plot(Ig, Gg, fillstyle='full', label=sh.name, marker=mm[i], mfc=cc[i], mec=cc[i], **mrkg)
if plot_brine:
ax.plot(Ib, Gb, fillstyle='none', label=None, marker=mm[i], mec=cc[i], **mrkb)
ax.set_xlabel('Intercept')
ax.set_ylabel('Gradient')
ax.legend()
ax.set_xlim(-mx, mx)
ax.set_ylim(-mx, mx)
ax.set_aspect('equal', 'box')
ax.grid()
if plot_colorzones:
opt1 = dict(edgecolor='None', alpha=0.2)
cl1_area = patches.Rectangle((0.02, -1), .98, 1, facecolor=cc['CLASS1'], **opt1)
cl2_area = patches.Rectangle((-0.02, -1), .04, 2, facecolor=cc['CLASS2'], **opt1)
cl3_area = patches.Rectangle((-1, -1), .98, 1, facecolor=cc['CLASS3'], **opt1)
cl4_area = patches.Rectangle((-1, 0), .98, 1, facecolor=cc['CLASS4'], **opt1)
background = patches.Polygon([[-1, 1], [1, -1], [1, 1]], facecolor='w')
ax.add_patch(cl1_area)
ax.add_patch(cl2_area)
ax.add_patch(cl3_area)
ax.add_patch(cl4_area)
ax.add_patch(background)
opt2 = dict(ha='center', va='center', weight='bold', size='large')
ax.text(.15, -.3, 'Class 1', color=cc['CLASS1'], **opt2)
ax.text(0, -.25, 'Class 2/2p', color=cc['CLASS2'], **opt2)
ax.text(-.35, -.3, 'Class 3', color=cc['CLASS3'], **opt2)
ax.text(-.35, .15, 'Class 4', color=cc['CLASS4'], **opt2)
if plot_brine:
ax.set_title('Filled markers: gas, empty=brine')
return ax
The table below shows the elastic properties for the various AVO classes taken from the references above (Hilterman 2001, Castagna et al., 1997 and 1998):
facies | AVO | Vp | Vs | density | porosity |
---|---|---|---|---|---|
shale | 1 | 3094 | 1515 | 2.40 | |
gas sand | 1 | 4050 | 2526 | 2.21 | 0.20 |
brine sand | 1 | 4115 | 2543 | 2.32 | 0.20 |
shale | 2 | 2643 | 1167 | 2.29 | |
gas sand | 2 | 2781 | 1665 | 2.08 | 0.25 |
brine sand | 2 | 3048 | 1595 | 2.23 | 0.25 |
shale | 3 | 2192 | 818 | 2.16 | |
gas sand | 3 | 1543 | 901 | 1.88 | 0.33 |
brine sand | 3 | 2134 | 860 | 2.11 | 0.33 |
shale | 4 | 3240 | 1620 | 2.34 | |
gas sand | 4 | 1650 | 1090 | 2.07 | ? |
brine sand | 4 | 2590 | 1060 | 2.21 | ? |
Porosity is unknown for class 4.
First we build 3 tables detailing the elastic properties $V_p$, $V_s$ and $\rho$ plus porosity for shales, gas sands and brine sands for each AVO class:
tmp_shale = np.array([[3094,1515,2.40,0], [2643,1167,2.29,0], [2192,818,2.16,0], [3240,1620,2.34,0]])
tmp_sandg = np.array([[4050,2526,2.21,.2], [2781,1665,2.08,.25], [1542,901,1.88,.33], [1650,1090,2.07,0]])
tmp_sandb = np.array([[4115,2453,2.32,.2], [3048,1595,2.23,.25], [2134,860,2.11,.33], [2590,1060,2.21,0]])
avocl = ['CLASS1', 'CLASS2', 'CLASS3', 'CLASS4']
logs = ['VP', 'VS', 'RHO', 'PHI']
shale = pd.DataFrame(tmp_shale, columns=logs, index=avocl)
sandg = pd.DataFrame(tmp_sandg, columns=logs, index=avocl)
sandb = pd.DataFrame(tmp_sandb, columns=logs, index=avocl)
shale
VP | VS | RHO | PHI | |
---|---|---|---|---|
CLASS1 | 3094.0 | 1515.0 | 2.40 | 0.0 |
CLASS2 | 2643.0 | 1167.0 | 2.29 | 0.0 |
CLASS3 | 2192.0 | 818.0 | 2.16 | 0.0 |
CLASS4 | 3240.0 | 1620.0 | 2.34 | 0.0 |
sandg
VP | VS | RHO | PHI | |
---|---|---|---|---|
CLASS1 | 4050.0 | 2526.0 | 2.21 | 0.20 |
CLASS2 | 2781.0 | 1665.0 | 2.08 | 0.25 |
CLASS3 | 1542.0 | 901.0 | 1.88 | 0.33 |
CLASS4 | 1650.0 | 1090.0 | 2.07 | 0.00 |
sandb
VP | VS | RHO | PHI | |
---|---|---|---|---|
CLASS1 | 4115.0 | 2453.0 | 2.32 | 0.20 |
CLASS2 | 3048.0 | 1595.0 | 2.23 | 0.25 |
CLASS3 | 2134.0 | 860.0 | 2.11 | 0.33 |
CLASS4 | 2590.0 | 1060.0 | 2.21 | 0.00 |
We can derive porosity from given densities using the simple formula:
$$ \phi = \frac{\rho - \rho_m} {\rho_f- \rho_m} $$The returned values when compared to the given porosities are a bit off; I need to tweak brine and gas density to get porosities for classes 1-3 "almost" correct.
For example, starting from brine sands I initially set $\rho_{B} = 1$ for brine density but I ended up with $1.0135$ to match the given values.
Starting from gas sands instead I have to use an unusually high gas density ($\rho_{G} = 0.35$, usually this is around $0.1$) to get closer to the given values.
rho_m = 2.65 # mineral density
print('STARTING POINT: BRINE SAND DENSITIES')
for i, ss in sandb.iterrows():
rho_fluid = 1.0135
phi = (ss['RHO'] - rho_m)/(rho_fluid - rho_m)
itms = (ss.name, phi, ss.PHI)
print('{:s}: calculated={:.4f}, given={:.2f}'.format(*itms))
print('\nSTARTING POINT: GAS SAND DENSITIES')
for i, ss in sandg.iterrows():
rho_fluid = 0.35
phi = (ss['RHO'] - rho_m)/(rho_fluid - rho_m)
itms = (ss.name, phi, ss.PHI)
print('{:s}: calculated={:.4f}, given={:.2f}'.format(*itms))
STARTING POINT: BRINE SAND DENSITIES CLASS1: calculated=0.2016, given=0.20 CLASS2: calculated=0.2566, given=0.25 CLASS3: calculated=0.3300, given=0.33 CLASS4: calculated=0.2689, given=0.00 STARTING POINT: GAS SAND DENSITIES CLASS1: calculated=0.1913, given=0.20 CLASS2: calculated=0.2478, given=0.25 CLASS3: calculated=0.3348, given=0.33 CLASS4: calculated=0.2522, given=0.00
Instead of using fixed values for the fluid densities and assuming correct given porosities (at least for classes 1-3) and constant mineral density $\rho_{m} = 2.65 g/cm^2$, we can scan through realistic values for brine and gas density and select those values that return the given porosity. See next example for class 3:
i = 'CLASS3'
phi = sandb.loc[i, 'PHI']
rhoB = sandb.loc[i, 'RHO']
rhoG = sandg.loc[i, 'RHO']
rho_b = np.linspace(0.9, 1.2)
rho_g = np.linspace(0.05, 0.6)
phi_from_bri = (rhoB - rho_m)/(rho_b - rho_m)
phi_from_gas = (rhoG - rho_m)/(rho_g - rho_m)
rho_brine_exact = (rhoB-rho_m*(1-phi)) / phi
rho_gas_exact = (rhoG-rho_m*(1-phi)) / phi
mrkg = dict(marker='o', color='r', ms=12, mfc='none', ls='none')
mrkb = dict(marker='o', color='b', ms=12, mfc='none', ls='none')
plt.figure()
plt.plot(phi_from_gas, rho_g, '-r', label='from gas')
plt.plot(phi_from_bri, rho_b, '-b', label='from brine')
if i != 'CLASS4':
plt.axvline(sandb.loc[i, 'PHI'], ls='--', color='k')
plt.plot(phi, rho_brine_exact, label='{:.3f}'.format(rho_brine_exact), **mrkb)
plt.plot(phi, rho_gas_exact, label='{:.3f}'.format(rho_gas_exact), **mrkg)
plt.legend()
plt.xlabel('phi')
plt.ylabel('rho fluids')
plt.grid()
def phi_from_dens(avoclass, rho_m=2.65):
phi = sandb.loc[avoclass, 'PHI']
rhoB = sandb.loc[avoclass, 'RHO']
rhoG = sandg.loc[avoclass, 'RHO']
rho_b = np.linspace(0.9, 1.2)
rho_g = np.linspace(0.05, 0.6)
phi_from_bri = (rhoB - rho_m)/(rho_b - rho_m)
phi_from_gas = (rhoG - rho_m)/(rho_g - rho_m)
rho_brine_exact = (rhoB-rho_m*(1-phi)) / phi
rho_gas_exact = (rhoG-rho_m*(1-phi)) / phi
mrkg = dict(marker='o', color='r', ms=12, mfc='none', ls='none')
mrkb = dict(marker='o', color='b', ms=12, mfc='none', ls='none')
f, ax = plt.subplots(constrained_layout=True)
ax.plot(phi_from_gas, rho_g, '-r', label='from gas')
ax.plot(phi_from_bri, rho_b, '-b', label='from brine')
ax.axvline(sandb.loc[avoclass, 'PHI'], ls='--', color='k')
ax.plot(phi, rho_brine_exact, label='{:.3f} g/cc'.format(rho_brine_exact), **mrkb)
ax.plot(phi, rho_gas_exact, label='{:.3f} g/cc'.format(rho_gas_exact), **mrkg)
ax.legend(loc='lower right')
ax.set_xlabel('Porosity')
ax.set_ylabel('Fluid density (g/cc)')
ax.grid()
w_avoc = ipw.ToggleButtons(description='AVO Class', options=['CLASS1', 'CLASS2', 'CLASS3'])
w_rhom = ipw.FloatSlider(value=2.65, description='Matrix density', max=2.8, min=2.5, step=0.02)
ui0 = ipw.VBox([w_avoc, w_rhom])
tool0 = ipw.interactive_output(phi_from_dens, {'avoclass': w_avoc, 'rho_m': w_rhom})
display(ui0, tool0)
The fluid densities that I have found can be stored into dictionaries for later use:
rho_b = {}
rho_g = {}
for i, ss in sandb.iterrows():
if i != 'CLASS4':
rhoB = sandb.loc[i, 'RHO']
rhoG = sandg.loc[i, 'RHO']
rho_b[i] = (rhoB-rho_m*(1-sandb.loc[i, 'PHI']))/sandb.loc[i, 'PHI']
rho_g[i] = (rhoG-rho_m*(1-sandb.loc[i, 'PHI']))/sandb.loc[i, 'PHI']
print('{:s}: rho brine={:.3f}, rho gas={:.3f}'.format(i, rho_b[i], rho_g[i]))
And copy the values for class 3 to class 4 since they should occur at similar depths (but this is unclear from Castagna's papers, so we are hunting in the dark right now):
rho_b['CLASS4'] = rho_b['CLASS3']
rho_g['CLASS4'] = rho_g['CLASS3']
Obviously when using these class-specific fluid densities starting from either brine or gas values I will now get the same porosities. The following code is redundant but since this is not a paper I will put it here as well:
for i, ss in sandb.iterrows():
phi_from_bri = (ss['RHO'] - rho_m)/(rho_b[i] - rho_m)
phi_from_gas = (sandg.loc[i, 'RHO'] - rho_m)/(rho_g[i] - rho_m)
print('{:s}: given phi={:.2f}, phi from brine={:.2f}, phi from gas={:.2f}'.format(i, ss.PHI, phi_from_bri, phi_from_gas))
I now try a different approach using Gassmann's formula to derive porosity:
$$ \frac{k_{satB}} {k_m - k_{satG}} = \frac {k_d} {k_m - k_d} + \frac{k_B}{\phi \cdot (k_m - k_B)} $$and
$$ \frac{k_{satG}} {k_m - k_{satG}} = \frac {k_d} {k_m - k_d} + \frac{k_G}{\phi \cdot (k_m - k_G)} $$These two combined become:
$$ \phi = \frac{(k_m - k_{satB}) \cdot (k_m - k_{satG}) \cdot (k_B - k_G)}{(k_m - k_B) \cdot (k_m - k_G) \cdot (k_{satB} - k_{satG})}$$An example of how to code the above:
k_b = 2.2
k_g = 0.02
k_m = 37
i = 'CLASS3'
q, w, e, phi_given = sandg.loc[i].values
ksat_g = bulk(q, w, e)
q, w, e, _ = sandb.loc[i].values
ksat_b = bulk(q, w, e)
phi = gassmann_phi(ksat_b, ksat_g, 2.1, k_g, k_m)
print('Porosity: calculated={:.4f}, given={:.4f}'.format(phi, phi_given))
I also need to take into account that the fluids' bulk moduli are depth-dependant (and densities too): I should derive them systematically using Batzle-Wang equations but for now with simple trial and error I have derived the following brine bulk moduli:
Applying the same formula above to class 4:
i = 'CLASS4'
q, w, e, phi_given = sandg.loc[i].values
ksat_g = bulk(q, w, e)
q, w, e, _ = sandb.loc[i].values
ksat_b = bulk(q, w, e)
phi = gassmann_phi(ksat_b, ksat_g, k_b, k_g, k_m)
print('Porosity: calculated={:.4f}, given={:.4f}'.format(phi, phi_given))
The idea is to compute dry-rock bulk modulus for a range of porosities starting from brine and sand properties. Where they intersect, that's where we expect to find the 'correct' porosity. This works for classes 1-3 as expected, while for class 4 it obviously returns the same porosity I get by inverting Gassmann (0.163) which is quite clearly wrong.
Dry-rock bulk moduli are obviously much more prone to error if we use a denser fluid such as brine, that's why the blue lines below also reach negative (clearly incorrect) values.
def phi_finder(avoclass, rho_g, k_g, rho_w, k_w, k_m):
phi_array = np.linspace(.1, .4)
vpG, vsG, rhoG, phi_given = sandg.loc[avoclass]
vpB, vsB, rhoB, _ = sandb.loc[avoclass]
vp_0, vs_0, rho_0, ksat_0, kdry_0 = gassmann(vpG, vsG, rhoG, rho_g, k_g, rho_w, k_w, k_m, phi_array)
vp_1, vs_1, rho_1, ksat_1, kdry_1 = gassmann(vpB, vsB, rhoB, rho_w, k_w, rho_g, k_g, k_m, phi_array)
# find intersection and return porosity
idx, val = find_nearest(kdry_0,kdry_1)
phi_calc = phi_array[idx]
f, ax = plt.subplots(figsize=(7, 5), constrained_layout=True)
ax.plot(phi_array, kdry_0, '-r', lw=4, alpha=0.5, label='K_dry from gas')
ax.plot(phi_array, kdry_1, '-b', lw=4, alpha=0.5, label='K_dry from brine')
ax.axvline(phi_given, ls=':', color='k', label='Given phi={:.3f}'.format(phi_given))
ax.plot(phi_calc, val, marker='o', mec='k', mfc='none', ms=15, ls='none', label='Calculated phi={:.3f}'.format(phi_calc))
ax.set_xlabel('Porosity')
ax.set_ylabel(r'$K_{dry}$ (GPa)')
ax.legend(bbox_to_anchor=(1.01, 1))
ax.set_ylim(-5, 25)
ax.set_title(avoclass)
def reset_values(b):
w_avoc.value = 'CLASS3'
w_km.value = 37
w_rhog.value = 0.3
w_kg.value = 0.02
w_rhow.value = 1.0
w_kw.value = 2.2
sty = {'description_width': 'initial'}
w_avoc = ipw.ToggleButtons(value='CLASS3', description='AVO Class', options=['CLASS1', 'CLASS2', 'CLASS3', 'CLASS4'])
w_km = ipw.FloatSlider(value=37, description='K matrix (GPa)', max=40, min=30, step=1, style=sty)
w_rhog = ipw.FloatSlider(value=0.3, description='rho gas (g/cc)', max=0.4, min=0.1, step=0.05, style=sty)
w_kg = ipw.FloatSlider(value=0.02, description='K gas (GPa)', max=0.2, min=0.01, step=0.01, style=sty)
w_rhow = ipw.FloatSlider(value=1.0, description='rho water (g/cc)', max=1.1, min=0.9, step=0.05, style=sty)
w_kw = ipw.FloatSlider(value=2.2, description='K water (GPa)', max=3.5, min=2, step=0.05, style=sty)
reset_button = ipw.Button(description = "Reset")
reset_button.on_click(reset_values)
lbox = ipw.VBox([w_kg, w_rhog])
rbox = ipw.VBox([w_kw, w_rhow])
lrbox = ipw.HBox([lbox, rbox])
ui1 = ipw.VBox([w_avoc, w_km, lrbox, reset_button])
tool1 = ipw.interactive_output(phi_finder, {
'avoclass': w_avoc,
'rho_g': w_rhog,
'k_g': w_kg,
'rho_w': w_rhow,
'k_w': w_kw,
'k_m': w_km})
display(ui1, tool1)
The same concept above, shown in a single figure:
phi = np.linspace(.1, .4)
f, ax = plt.subplots(ncols=4, figsize=(8,3), sharey=True,constrained_layout=True)
for n in range(4):
i = sandb.index[n]
vpG, vsG, rhoG, _ = sandg.loc[i].values
vp_0, vs_0, rho_0, ksat_0, kdry_0 = gassmann(vpG, vsG, rhoG, rho_g[i], k_g, rho_b[i], k_b, k_m, phi)
vpB, vsB, rhoB, _ = sandb.loc[i].values
vp_1, vs_1, rho_1, ksat_1, kdry_1 = gassmann(vpB, vsB, rhoB, rho_b[i], k_b, rho_g[i], k_g, k_m, phi)
ax[n].plot(phi, kdry_0, '-r', label='from gas')
ax[n].plot(phi, kdry_1, '-b', label='from brine')
ax[n].axhline(0, ls='-', color='k', lw=2)
if i != 'CLASS4':
ax[n].axvline(sandb.loc[i, 'PHI'], ls='--', color='k')
ax[n].set_xlabel('phi')
ax[n].set_ylabel('kdry')
ax[n].grid()
ax[n].set_title(i)
So using Gassmann I get 16% porosity while using densities I get 26% porosity:
check_vpb = sandb.loc['CLASS4', 'VP']
check_vpg = sandg.loc['CLASS4', 'VP']
print(check_vpb-check_vpg)
print(percdiff(check_vpg, check_vpb))
check_vpb = sandb.loc['CLASS3', 'VP']
check_vpg = sandg.loc['CLASS3', 'VP']
print(check_vpb-check_vpg)
print(percdiff(check_vpg, check_vpb))
Below is a little snippet of code showing fluid replacement on class 4 starting from both brine and gas sand with the two different porosities; we observe greater fluid effect with 16% porosity which is a bit stange to say the least:
i = 'CLASS4'
vpB, vsB, rhoB, _ = sandb.loc[i].values
vpG, vsG, rhoG, _ = sandg.loc[i].values
print('\nBRINE --> GAS')
print('initial: vp={:5.0f}, vs={:5.0f}, rho={:5.2f}'.format(vpB, vsB, rhoB))
print('expected: vp={:5.0f}, vs={:5.0f}, rho={:5.2f}'.format(vpG, vsG, rhoG))
phi = 0.26
vp2, vs2, rho2, ksat2, kdry2 = gassmann(vpB, vsB, rhoB, rho_b[i], k_b, rho_g[i], k_g, k_m, phi)
print('gassmann (phi={:.2f}): vp={:5.0f}, vs={:5.0f}, rho={:5.2f}, ksat={:5.2f}, kdry={:4.2f}'.format(phi, vp2, vs2, rho2, ksat2, kdry2))
phi = 0.16
vp2, vs2, rho2, ksat2, kdry2 = gassmann(vpB, vsB, rhoB, rho_b[i], k_b, rho_g[i], k_g, k_m, phi)
print('gassmann (phi={:.2f}): vp={:5.0f}, vs={:5.0f}, rho={:5.2f}, ksat={:5.2f}, kdry={:4.2f}'.format(phi, vp2, vs2, rho2, ksat2, kdry2))
print('\nGAS --> BRINE')
print('initial: vp={:5.0f}, vs={:5.0f}, rho={:5.2f}'.format(vpG, vsG, rhoG))
print('expected: vp={:5.0f}, vs={:5.0f}, rho={:5.2f}'.format(vpB, vsB, rhoB))
phi = 0.26
vp2, vs2, rho2, ksat2, kdry = gassmann(vpG, vsG, rhoG, rho_g[i], k_g, rho_b[i], k_b, k_m, phi)
print('gassmann (phi={:.2f}): vp={:5.0f}, vs={:5.0f}, rho={:5.2f}, ksat={:5.2f}, kdry={:4.2f}'.format(phi, vp2, vs2, rho2, ksat2, kdry2))
phi = 0.16
vp2, vs2, rho2, ksat2, kdry = gassmann(vpG, vsG, rhoG, rho_g[i], k_g, rho_b[i], k_b, k_m, phi)
print('gassmann (phi={:.2f}): vp={:5.0f}, vs={:5.0f}, rho={:5.2f}, ksat={:5.2f}, kdry={:4.2f}'.format(phi, vp2, vs2, rho2, ksat2, kdry2))
This thingie shows the numerical results of Gassmann's equation on the various AVO class sands depending on porosity.
k_m = 37
k_w, rho_w = 2.2, 1.0
k_g, rho_g = 0.02, 0.32
def frm_check_pors(avoclass, phi):
vpG, vsG, rhoG, phi_given = sandg.loc[avoclass]
vpB, vsB, rhoB, _ = sandb.loc[avoclass]
vp_0, vs_0, rho_0, ksat_0, kdry_0 = gassmann(vpG, vsG, rhoG, rho_g, k_g, rho_w, k_w, k_m, phi)
vp_1, vs_1, rho_1, ksat_1, kdry_1 = gassmann(vpB, vsB, rhoB, rho_w, k_w, rho_g, k_g, k_m, phi)
print('Given porosity: {:.2f}'.format(phi_given))
print('Assigned porosity: {:.2f}'.format(phi))
print('\nBRINE --> GAS')
print('Input (brine): Vp={:.0f}, Vs={:.0f}, rho={:.2f}'.format(vpB, vsB, rhoB))
print('Reference (gas): Vp={:.0f}, Vs={:.0f}, rho={:.2f}'.format(vpG, vsG, rhoG))
print('Output FRM gas: Vp={:.0f}, Vs={:.0f}, rho={:.2f}'.format(vp_1, vs_1, rho_1))
print('gassmann qc: ksat={:5.2f}, kdry={:4.2f}'.format(ksat_1, kdry_1))
itms = (percdiff(vpB, vp_1), percdiff(vsB, vs_1), percdiff(rhoB, rho_1))
print('Difference: Vp={:+4.0%}, Vs={:+4.0%}, rho={:+4.0%}'.format(*itms))
print('\nGAS --> BRINE')
print('Input (gas): Vp={:.0f}, Vs={:.0f}, rho={:.2f}'.format(vpG, vsG, rhoG))
print('Reference (brine): Vp={:.0f}, Vs={:.0f}, rho={:.2f}'.format(vpB, vsB, rhoB))
print('Output FRM brine: Vp={:.0f}, Vs={:.0f}, rho={:.2f}'.format(vp_0, vs_0, rho_0))
print('gassmann qc: ksat={:5.2f}, kdry={:4.2f}'.format(ksat_0, kdry_0))
itms = (percdiff(vpG, vp_0), percdiff(vsG, vs_0), percdiff(rhoG, rho_0))
print('Difference: Vp={:+4.0%}, Vs={:+4.0%}, rho={:+4.0%}'.format(*itms))
ipw.interact(frm_check_pors,
avoclass=ipw.ToggleButtons(options=['CLASS1', 'CLASS2', 'CLASS3', 'CLASS4'], value='CLASS4'),
phi=ipw.FloatSlider(min=0.1, max=0.40, step=0.01, value=0.26))
What is the conclusion from all of this:
caprock_vp, caprock_vs, caprock_rho, _ = shale.loc['CLASS4']
sand_vp, sand_vs, sand_rho, _ = sandg.loc['CLASS3']
_, i_newcl4, g_newcl4 = shuey(caprock_vp, caprock_vs, caprock_rho, sand_vp, sand_vs, sand_rho, np.arange(30), terms=True)
ax = classref(plot_brine=False)
ax.plot(i_newcl4, g_newcl4, '^', mec='r', mfc='k', ms=14)
In the plot above, Castagna's class 4 is the red triangle, the "new" class 4 is the black triangle outlined in red.
Work in progress. The idea is to set up a Monte Carlo simulation to capture the possibile variations of a real sand using rock physics modeling but I only arrived to setting up an initial porosity distribution. See also the lightning talk notebook
Make 500 porosity samples following a normal distribution with:
$$ \mu = 0.2 $$$$ \sigma = 0.03 $$from scipy.stats import norm
ns = 500
phi_dist = norm.rvs(loc=0.2, scale=0.03, size=ns)
pd.Series(phi_dist).plot(kind='density', figsize=(5,3), xlabel='Porosity',
xlim=(phi_dist.min(), phi_dist.max()),
title='Porosities',
rot=45, color='k', lw=4);