From this pull request: https://github.com/agile-geoscience/bruges/pull/79
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import lasio
%cd ../bruges
import bruges.rockphysics.rockphysicsmodels as brp
import bruges.rockphysics.fluidsub as bfl
from bruges.rockphysics import vrh
%cd ../geophysical_notes
C:\Users\ag19324\GoogleDrive\PYTHON\_GITHUB\bruges C:\Users\ag19324\GoogleDrive\PYTHON\_GITHUB\geophysical_notes
Reproducing the results from the example at page 260 of Mavko et al., The Rock Physics Handbook (2009):
from IPython.display import Image
Image(filename='Mavko_2009_pag260.png', width=600)
Note that the bulk and shear elastic moduli ($K$ and $\mu$ from the book's example above, K
and G
in the code block below) are in GPa, density is in g/cc, pressure in MPa.
When calculating the final velocities with:
$$ V_p = \sqrt{ \frac{K + 4/3 \mu}{rho}} ; V_s = \sqrt{ \frac{\mu}{rho}} $$remember that with these input units the output will be in km/s. To get velocities in m/s just use Pa for the moduli and kg/m3 for density; i.e. K = K * 1e9
and rho = rho * 1e3
.
K0,G0,D0 = 36.6,45,2.65
sigma =20
phi = 0.3
Khm,Ghm = brp.hertz_mindlin(K0, G0, sigma, phi_c=0.36, Cn=9)
Kd,Gd = brp.soft_sand(K0, G0, phi, sigma, phi_c=0.36, Cn=9)
rho = D0*(1-phi)
vp_soft = np.sqrt((Kd+4/3*Gd)/rho)
vs_soft = np.sqrt(Gd/rho)
print('Hertz-Mindlin pure quartz at effective pressure={} MPa, phi_c={:.2f}: K={:.2f} GPa, G={:.2f} GPa'.format(sigma,.36,Khm,Ghm))
print('Soft sand model at porosity={:.2f}: K={:.2f} GPa, G={:.2f} GPa'.format(phi,Kd,Gd))
print('Soft sand model velocities at porosity={:.2f}: Vp={:.2f} km/s, Vs={:.2f} km/s, rho={:.3f} g/cc'.format(phi,vp_soft,vs_soft,rho))
Hertz-Mindlin pure quartz at effective pressure=20 MPa, phi_c=0.36: K=2.05 GPa, G=3.02 GPa Soft sand model at porosity=0.30: K=3.05 GPa, G=3.99 GPa Soft sand model velocities at porosity=0.30: Vp=2.12 km/s, Vs=1.47 km/s, rho=1.855 g/cc
The differences from the book's example (3.05 GPa instead of 3 for the dry bulk modulus, 2.12 km/s instead of 2.05 for Vp, etc.) are probably just a consequence of rounding off floats.
Reproducing the results from figure 5.4.7 at page 261 of The Rock Physics Handbook. Will not bother with Raymer-Hunt-Gardner and intermediate stiff sand models. Use different linewidths to better differentiate between Stiff and Soft Sand models. Remember that the P-wave modulus is $K + 4/3 \mu = \rho {V_p}^2$.
Image(filename='Mavko_2009_fig547.png', width=800)
phi = np.linspace(0,0.36,200)
Kd_soft,Gd_soft = brp.soft_sand(K0, G0, phi, sigma, phi_c=0.36, Cn=9)
Kd_stiff,Gd_stiff = brp.stiff_sand(K0, G0, phi, sigma, phi_c=0.36, Cn=9)
rho = D0*(1-phi)
vp_soft = np.sqrt((Kd_soft+4/3*Gd_soft)/rho)
vs_soft = np.sqrt(Gd_soft/rho)
vp_stiff = np.sqrt((Kd_stiff+4/3*Gd_stiff)/rho)
vs_stiff = np.sqrt(Gd_stiff/rho)
f,ax=plt.subplots(1,2,figsize=(10,8))
ax[0].plot(phi,Kd_soft+4/3*Gd_soft, '-k', lw=1, label='P-wave modulus Soft Sand')
ax[0].plot(phi,Kd_stiff+4/3*Gd_soft, '-k', lw=3, label='P-wave modulus Stiff Sand')
ax[0].plot(phi,Gd_soft, '--k', lw=1, label='Shear Modulus Soft Sand')
ax[0].plot(phi,Gd_stiff, '--k', lw=3, label='Shear Modulus Stiff Sand')
ax[1].plot(phi,vp_soft, '-k', lw=1, label='Vp Soft Sand')
ax[1].plot(phi,vp_stiff, '-k', lw=3, label='Vp Stiff Sand')
ax[1].plot(phi,vs_soft, '--k', lw=1, label='Vs Soft Sand')
ax[1].plot(phi,vs_stiff, '--k', lw=3, label='Vs Stiff Sand')
ax[0].set_ylim(0,100)
ax[1].set_ylim(0,6.5)
ax[0].set_ylabel('Modulus (GPa)')
ax[1].set_ylabel('Velocity (km/s)')
ax[1].yaxis.tick_right()
ax[1].yaxis.set_label_position('right')
for aa in ax:
aa.set_xlim(0,0.4)
aa.legend(fontsize='medium')
aa.grid()
aa.set_xlabel('Porosity')
f.suptitle('Mavko et al. (2009), Fig. 5.4.7, p.261');
Reproducing the results from Dvorkin & Nur, Elasticity of High-Porosity Sandstones: Theory for Two North Sea Data Sets, Geophysics (1996), figure 5. The three subplots only differ for the plot of the Oseberg samples at different confining pressures, but it is only relevent for this exercise comparing the plots of the theoretical models which are the same for the three plots.
Image(filename='Dvorkin-Nur_1996_fig5.png', width=800)
K0,G0,D0 = 38,44,2.65
phi = np.linspace(0.15,0.4,200)
rho = D0*(1-phi)
Kc, Gc = 36.6, 45 # quartz cement
Kd,Gd = brp.contact_cement(K0, G0, phi, Kc=Kc, Gc=Gc)
vp_cc_qtz = np.sqrt((Kd+4/3*Gd)/rho)
vs_cc_qtz = np.sqrt(Gd/rho)
Kc, Gc = 21, 7 # clay cement
Kd,Gd = brp.contact_cement(K0, G0, phi, Kc=Kc, Gc=Gc)
vp_cc_cl = np.sqrt((Kd+4/3*Gd)/rho)
vs_cc_cl = np.sqrt(Gd/rho)
f,ax=plt.subplots(figsize=(5,6))
ax.plot(phi,vp_cc_qtz, 'k-', label='contact cement model, cement=qtz')
ax.plot(phi,vs_cc_qtz, 'k--', label='contact cement model, cement=qtz')
ax.plot(phi,vp_cc_cl, 'r-', label='contact cement model, cement=clay')
ax.plot(phi,vs_cc_cl, 'r--', label='contact cement model, cement=clay')
ax.set_xlim(0,0.4)
ax.set_ylim(0,5)
ax.set_xlabel('Porosity')
ax.set_ylabel('Velocities (km/s)')
ax.legend(fontsize='medium')
ax.grid()
ax.set_title('Dvorkin & Nur (1996), fig.5');
Reproducing the results from Vernik & Kachanov, Modeling Elastic Properties of Siliciclastic Rocks, Geophysics (2010), figure 4.
Image(filename='Vernik_2010_fig4.png', width=400)
K0,G0 = 35.7,33
phi = np.linspace(0,.5,200)
Kd_10, Gd_10 = brp.vernik_consol_sand(K0,G0,phi,sigma=10)
Kd_20, Gd_20 = brp.vernik_consol_sand(K0,G0,phi,sigma=20)
Kd_50, Gd_50 = brp.vernik_consol_sand(K0,G0,phi,sigma=50)
f,ax=plt.subplots(figsize=(6,6))
ax.plot(phi,Kd_10+4/3*Gd_10,'-g', label='M-modulus 10Mpa')
ax.plot(phi,Kd_20+4/3*Gd_20,'-b', label='M-modulus 20Mpa')
ax.plot(phi,Kd_50+4/3*Gd_50,'-r', label='M-modulus 50Mpa')
ax.plot(phi,Gd_10,'--g', label='G-modulus 10Mpa')
ax.plot(phi,Gd_20,'--b', label='G-modulus 20Mpa')
ax.plot(phi,Gd_50,'--r', label='G-modulus 50Mpa')
ax.set_xlim(0,0.4)
ax.set_ylim(0,80)
ax.set_xlabel('Porosity')
ax.set_ylabel('Dry moduli (GPa)')
ax.legend(fontsize='medium')
ax.grid()
ax.set_title('Vernik & Kachanov (2010), fig.4\nConsolidated Sandstone Model (CSM)');
Reproducing the results from Vernik & Kachanov, Modeling Elastic Properties of Siliciclastic Rocks, Geophysics (2010), figure 5.
Image(filename='Vernik_2010_fig5.png', width=400)
K0,G0,D0 = 35.7,33,2.65
phi = np.linspace(0,.4,200)
rho = D0*(1-phi)
Kd, Gd = brp.vernik_sand_diagenesis(K0,G0,phi,sigma=50,phi_c=0.4,phi_con=0.3)
vp = np.sqrt((Kd + 4/3*Gd)/rho)
vs = np.sqrt(Gd/rho)
f,ax=plt.subplots(figsize=(6,6))
ax.plot(phi,vp,'-r', label='Vp')
ax.plot(phi,vs,'-g', label='Vs')
ax.set_xlim(0,0.4)
ax.set_ylim(0,6)
ax.set_xlabel('Porosity')
ax.set_ylabel('Dry velocities (km/s)')
ax.legend(fontsize='medium')
ax.grid()
ax.set_title('Vernik & Kachanov (2010), fig.5\nSand Diagenesis Model (SDM)');
Reproducing the results from Vernik & Kachanov, Modeling Elastic Properties of Siliciclastic Rocks, Geophysics (2010), figure 11.
Image(filename='Vernik_2010_fig11.png', width=400)
vv = np.linspace(0.2,0.8,6)
phi = np.linspace(0.01,0.4,200)
f,ax=plt.subplots(figsize=(6,6))
for i in vv:
vp, vs, rho = brp.vernik_shale(i, phi)
ax.plot(phi,vp,'-k', lw=3)
ax.plot(phi,vs,'-k', lw=1)
ax.set_xlim(0,0.4)
ax.set_ylim(0,5e3)
ax.set_xlabel('Porosity')
ax.set_ylabel('velocities (m/s)')
ax.grid()
ax.set_title('Vernik & Kachanov (2010), fig.11\nShale Model (thick=Vp, thin=Vs)');
Reproducing the results from Vernik & Kachanov, Modeling Elastic Properties of Siliciclastic Rocks, Geophysics (2010), figure 9. Only the "consolidated sand model" lines will be plotted. The Soft Sand Model 2 as it is called in the library, can be found at eq.14 in the paper:
$$ M_{dry} = M_0 \Big( 1 + p \frac{\phi}{1-phi} \Big)^{-1} $$$$ G_{dry} = G_0 \Big( 1 + q \frac{\phi}{1-phi} \Big)^{-1} $$Image(filename='Vernik_2010_fig9.png', width=400)
K0,G0,D0 = 35.6,33,2.65
phi= np.linspace(0,0.4,200)
poreshape=[4,6,10,18,45]
Kf, Df = 3.20, 1.05
f,ax=plt.subplots(figsize=(6,6))
for mm in poreshape:
Kd, Gd = brp.vernik_soft_sand_2(K0,G0,phi,p=mm,q=mm)
vp,*_ = bfl.vels(Kd,Gd,K0,D0,Kf,Df,phi)
ax.plot(phi,vp/1e3, label='SSM2, pore shape={}'.format(mm))
ax.set_xlim(0,0.4)
ax.set_ylim(1.5,5.5)
ax.set_xlabel('Porosity')
ax.set_ylabel('Velocities (km/s)')
ax.legend(fontsize='medium')
ax.grid()
ax.set_title('Vernik & Kachanov (2010), fig.9\nSoft Sand Model 2');
A small example on how to create Vp, Vs and density logs from scratch, using a rock physics model:
phi
), shale and sand volumes (vsh
, vqz
), water saturation (sw
).res
, identified by the interval top
-base
).phi + vsh + vqz = 1
....continues below...
#--- elastic moduli for minerals and fluids
#--- bulk and shear moduli in GPa, density in g/cc
Dsh, Ksh, Gsh = 2.4, 15, 5
Dqz, Kqz, Gqz = 2.6, 37, 44
Db, Kb = 1.1, 2.2
Do, Ko = 0.8, 1.5
Dg, Kg = 0.2, 0.06
hc='gas'
if hc=='oil':
Dhc, Khc = Do, Ko
else:
Dhc, Khc = Dg, Kg
#--- creates dummy depth scale and empty petrophysical logs
ns=100
top,base=30,50
z=np.arange(ns)
phi,vsh,vqz = (np.zeros(ns) for _ in range(3))
sw = np.ones(ns)
#--- cuts out a reservoir section between top and base
#--- and assigns higher porosity and lower Sw
res = (z>=top) & (z<base)
phi[res] = 0.3
sw[res] = 0.1
ns_res = np.count_nonzero(res)
#--- adds random variations to Phi and Sw
phi[~res] += np.random.rand(ns-ns_res)*.05
phi[res] += np.random.rand(ns_res)*.2
sw[res] -= np.random.rand(ns_res)*.1
#--- fills in mineralogical fraction volumes
vqz = np.random.rand(ns)*.15*(1-phi)
vqz[res] += .4
vsh = 1-phi-vqz
...
bruges.rockphysics.vrh
) for bulk (K0
) and shear (G0
) elastic moduli...continues below...
#--- normalize mineral volumes
vshN = vsh/(1-phi)
vqzN = vqz/(1-phi)
#--- computes rock matrix elastic moduli via Voigt-Reuss-Hill average
K0 = vrh(Ksh,Kqz,vshN)
G0 = vrh(Gsh,Gqz,vshN)
D0 = vshN*Dsh+vqzN*Dqz
#--- computes total fluid elastic moduli via Reuss average
f = np.array([sw,1-sw]).T
k = np.resize(np.array([Kb,Khc]),np.shape(f))
Kf = 1./np.sum(f/k,axis=1)
Df = sw*Db+(1-sw)*Dhc
...
bruges.rockphysics.rpmodels.softsand
) to calculate dry rock properties then the final saturated velocities using Gassmann's equation (bruges.rockphysics.rpmodels.vels
)#--- use soft sand model to create velocities and density logs
Kdry, Gdry = brp.soft_sand(K0, G0, phi, sigma)
vp,vs,rho,_= bfl.vels(Kdry,Gdry,K0,D0,Kf,Df,phi)
#--- assemble everything into a Pandas DataFrame
temp={'VQZ':vqz,'VSH':vsh,'PHIE':phi,'SWE':sw,'VP':vp,'VS':vs,'RHO':rho}
ww0 = pd.DataFrame(data=temp, index=z)
ww0.index.rename('DEPTH', inplace=True)
Using the describe()
method of Pandas, the final velocities and density averaged over the entire sequence can be inspected:
ww0[['VP','VS','RHO']].describe()
VP | VS | RHO | |
---|---|---|---|
count | 100.000000 | 100.000000 | 100.000000 |
mean | 2733.720103 | 1419.932775 | 2.236246 |
std | 585.542515 | 201.555743 | 0.302957 |
min | 1401.999246 | 1008.755166 | 1.432287 |
25% | 2751.875933 | 1285.647087 | 2.348763 |
50% | 2961.610931 | 1450.999296 | 2.381535 |
75% | 3094.514107 | 1576.694784 | 2.395229 |
max | 3408.842794 | 1836.666584 | 2.428962 |
Pandas also allows to do a rough summary plot of this synthetic dataset with just one line of code:
ww0.plot(subplots=True,figsize=(10,8), grid=True);
Create a second dataset using the Stiff Sand rock physics model and inspect the differences with describe()
only in the reservoir section:
#--- use stiff sand model to create velocities and density logs
Kdry, Gdry = brp.stiff_sand(K0, G0, phi, sigma=20)
vp,vs,rho,_= bfl.vels(Kdry,Gdry,K0,D0,Kf,Df,phi)
#--- assemble everything into a Pandas DataFrame
temp={'VQZ':vqz,'VSH':vsh,'PHIE':phi,'SWE':sw,'VP':vp,'VS':vs,'RHO':rho}
ww1 = pd.DataFrame(data=temp, index=z)
ww1.index.rename('DEPTH', inplace=True)
print('--- Soft Sand model:')
print(ww0[['VP','VS','RHO']][res].describe())
print('--- Stiff Sand model:')
print(ww1[['VP','VS','RHO']][res].describe())
--- Soft Sand model: VP VS RHO count 20.000000 20.000000 20.000000 mean 1606.703959 1102.885585 1.641811 std 96.090800 46.658509 0.106798 min 1401.999246 1008.755166 1.432287 25% 1550.694607 1074.909727 1.572600 50% 1614.409488 1106.890248 1.653390 75% 1668.569112 1128.120028 1.708809 max 1773.161433 1191.871493 1.818819 --- Stiff Sand model: VP VS RHO count 17.000000 18.000000 20.000000 mean 1805.115830 1148.714934 1.641811 std 400.956907 344.454117 0.106798 min 902.517688 47.020888 1.432287 25% 1671.135283 1118.325304 1.572600 50% 1815.639874 1208.411846 1.653390 75% 2022.282639 1337.622796 1.708809 max 2431.872158 1551.560785 1.818819