Within Danomics Flows you can add a Python block, such as shown below:
To utilize the Danomics Petrophysics library you will need to import it. This can be done as shown above with the statement on line 2:
from src.main.python import petrophysicsAlternatively, the import can be performed for only specific functions. For example, to import the function for Vclay from Gamma Ray and Vclay from Neutron Density you’d use a statement like:
from src.main.python.petrophysics import calculate_nd_clay, calculate_gr_clayRegardless of the import style it could then be used in the execute block. For the former, calling the function would look like:
vclay_gr = petrophysics.calculate_gr_clay(gr, gr_clean, gr_clay)
vclay_nd = petrophysics.calculate_nd_clay(rhob, nphi, rhob_matrix, rhob_clean, rhob_clay, nphi_matrix, nphi_clean, nphi_clay)Similarly, in the latter the same commands would be:
vclay_gr = calculate_gr_clay(gr, gr_clean, gr_clay)
vclay_nd = calculate_nd_clay(rhob, nphi, rhob_matrix, rhob_clean, rhob_clay, nphi_matrix, nphi_clean, nphi_clay)Notes on usage
Remember that Danomics treats well log curves as numpy arrays, and these functions take advantage of that. This significantly increases the speed and for many operations removes the need to use for loops to iterate through arrays. However, it also requires that you ensure all your arrays are the same length – and although this is the case by default in Danomics if you work with array slices care should be taken. If you use arrays of different lengths as arguments then you’ll receive assertion errors.
Many of these functions take both curves and parameters. The curves, as mentioned above, are arrays. The params you will likely use (at least in this version) will likely be scalar values. This is fine and the library accommodates this. If you do have an array of parameter values, this will work as well, as long as that array has the same length as the input curves.
The code is given below.
Code
"""
Date: May 17th, 2021
Copyright Danomics LLC, 2021. All rights reserved.
Provided for use with Danomics software only.
No authorized for distribution or publication.
@author: CameronSnow
"""
import numpy as np
#################################
### Helper Conversions ###
#################################
def degf2degc(value):
"""
Convert value from degrees Farenheit to Celsius.
"""
return (value-32)*5/9
def mpsi2gpa(value):
"""
Convert value from degrees milli-pound force per square inch to Gigapascale.
"""
return value * 6.89475729
def psi2ppg(value, depth_curve):
"""
Convert value from psi to pounds per gallon.
"""
return (value / depth_curve)/0.052
def ppg2psi(value, depth_curve):
"""
Convert value from pounds per gallon to psi.
"""
return (value*depth_curve)*0.052
def fs2kms(value):
"""
Convert value from feet/second to km/second.
"""
return value * 0.000305
def kms2fs(value):
"""
Convert value from km/second to feet/second.
"""
return value * 3280.84
def api2gcc(value):
"""
Convert value from API gravity to g/cm3.
"""
return 141.5/(value+131.5)
def psi2mpa(value):
"""
Convert value from psi to megapascale.
"""
return value*6.894754*0.001
def mpa2psi(value):
"""
Convert value from megapascale to psi.
"""
return value/6.894754/0.001
def scfbbl2ll(value):
"""
Convert value from scf/barrel to liters/liter.
"""
return value/5.6145
def usf2fs(value):
"""
Convert value from microseconds/foot to feet/second.
"""
return 1/value*1000000
#################################
### Clay Volume Functions ###
#################################
def calculate_single_clay(curve, clean, clay):
"""
Calculate clay volume by a single clay method
Parameters
----------
curve : ARRAY
Curve to use.
clean : FLOAT
Clean point for curve.
clay : FLOAT
Clay point for curve.
Returns
-------
vclay : ARRAY
Clay Volume.
"""
if np.isscalar(clay):
temp = np.zeros(len(curve))
temp[:] = clay
clay = temp
if np.isscalar(clean):
temp = np.zeros(len(curve))
temp[:] = clean
clean = temp
assert len(clean) == len(curve), "Input clean array must be of same length as curve."
assert len(clay) == len(curve), "Input clay array must be of same length as curve."
vclay = (curve - clean) / (clay-clean)
vclay = np.minimum(np.maximum(vclay, 0), 1)
return vclay
def calculate_gr_clay(gr, gr_clean, gr_clay):
"""
Calculate clay volume by gamma ray method
Parameters
----------
gr : ARRAY
Gamma Ray curve.
gr_clean : FLOAT
Clean value for gamma ray curve.
gr_clay : FLOAT
Clay value for gamma ray curve.
Returns
-------
vcl_gr : ARRAY
Clay Volume.
"""
vcl_gr = calculate_single_clay(gr, gr_clean, gr_clay)
return vcl_gr
def calculate_gr_clay_by_method(gr, gr_clean, gr_clay, method):
"""
Calculate clay volume from gamma ray by a specific method.
Parameters
----------
gr : ARRAY
Gamma Ray curve.
gr_clean : TYPE
Clean value for gamma ray curve.
gr_clay : FLOAT
Clay value for gamma ray curve.
method : STRING
Options are 'larionov_young', 'larionov_old', 'clavier', 'steiber'.
Returns
-------
vcl : ARRAY
Clay Volume.
"""
if method == 'larionov_young':
vcl = 0.083*(2**(3.7*((gr - gr_clean) / (gr_clay - gr_clean)))-1)
elif method == 'larionov_old':
vcl =0.33*(2**(2*((gr - gr_clean) / (gr_clay - gr_clean)))-1)
elif method == 'clavier':
vcl = 1.7-(3.38-(((gr - gr_clean) / (gr_clay - gr_clean))+0.7)**2)**0.5
elif method == 'steiber':
vcl = 0.5*((gr - gr_clean) / (gr_clay - gr_clean))/(1.5-((gr - gr_clean) / (gr_clay - gr_clean)))
vcl = np.minimum(np.maximum(vcl, 0), 1)
return vcl
def calculate_sp_clay(sp, sp_clean, sp_clay):
"""
Calculate clay volume by gamma ray method
Parameters
----------
sp : ARRAY
SP curve.
sp_clean : FLOAT
Clean value for SP curve.
sp_clay : FLOAT
Clay value for SP curve.
Returns
-------
vcl_sp : ARRAY
Clay Volume.
"""
vcl_sp = calculate_single_clay(sp, sp_clean, sp_clay)
return vcl_sp
def calculate_nphi_clay(nphi, nphi_clean, nphi_clay):
"""
Calculate clay volume by neutron porosity method
Parameters
----------
nphi : ARRAY
Neutron curve.
nphi_clean : FLOAT
Neutron curve clean value.
nphi_clay : FLOAT
Neutron curve clay value.
Returns
-------
vcl_nphi : ARRAY
Clay Volume.
"""
if np.isscalar(nphi_clay):
temp = np.zeros(len(nphi))
temp[:] = nphi_clay
nphi_clay = temp
if np.isscalar(nphi_clean):
temp = np.zeros(len(nphi))
temp[:] = nphi_clean
nphi_clean = temp
assert len(nphi_clean) == len(nphi), "Input clean array must be of same length as curve."
assert len(nphi_clay) == len(nphi), "Input clay array must be of same length as curve."
vcl_nphi = np.sqrt((nphi/nphi_clay)*((nphi-nphi_clean)/(nphi_clay-nphi_clean)))
vcl_nphi = np.minimum(np.maximum(vcl_nphi, 0), 1)
return vcl_nphi
def calculate_dual_clay(c1, c2, c1cl1, c1cl2, c1clay, c2cl1, c2cl2, c2clay):
"""
Calculate clay volume by a dual clay method.
Parameters
----------
c1 : ARRAY
Curve one.
c2 : ARRAY
Curve two.
c1cl1 : FLOAT
Matrix point of curve one.
c1cl2 : FLOAT
Clean point of curve one.
c1clay : FLOAT
Clay point of curve one.
c2cl1 : FLOAT
Matrix point of curve two.
c2cl2 : FLOAT
Clean point of curve two.
c2clay : FLOAT
CLay point of curve two.
Returns
-------
vclay : ARRAY
Clay volume.
"""
assert len(c1)== len(c2), "Input curves must be of same length."
if np.isscalar(c1cl1):
temp = np.zeros(len(c1))
temp[:] = c1cl1
c1cl1 = temp
if np.isscalar(c1cl2):
temp = np.zeros(len(c1))
temp[:] = c1cl2
c1cl2 = temp
if np.isscalar(c1clay):
temp = np.zeros(len(c1))
temp[:] = c1clay
c1clay = temp
if np.isscalar(c2cl1):
temp = np.zeros(len(c1))
temp[:] = c2cl1
c2cl1 = temp
if np.isscalar(c2cl2):
temp = np.zeros(len(c1))
temp[:] = c2cl2
c2cl2 = temp
if np.isscalar(c2clay):
temp = np.zeros(len(c1))
temp[:] = c2clay
c2clay = temp
assert len(c1cl1) == len(c1), "Input curve one,clean one array must be of same length as curve."
assert len(c1cl2) == len(c1), "Input curve one,clean two array must be of same length as curve."
assert len(c1clay) == len(c1), "Input curve one,clay array must be of same length as curve."
assert len(c2cl1) == len(c1), "Input curve two,clean one array must be of same length as curve."
assert len(c2cl2) == len(c1), "Input curve two,clean two array must be of same length as curve."
assert len(c2clay) == len(c1), "Input curve two,clay array must be of same length as curve."
vclay = (((c1cl2 - c1cl1)*(c2 - c2cl1) - (c1 - c1cl1)*(c2cl2 - c2cl1)) /
((c1cl2 - c1cl1)*(c2clay - c2cl1)-(c1clay - c1cl1)*(c2cl2 - c2cl1)))
vclay = np.minimum(np.maximum(vclay, 0),1)
return vclay
def calculate_nd_clay(rhob, nphi, rhob_matrix, rhob_clean, rhob_clay, nphi_matrix, nphi_clean, nphi_clay):
"""
Calculate vclay by neutron density method
Parameters
----------
rhob : ARRAY
Bulk density curve.
nphi : ARRAY
Neutron curve.
rhob_matrix : FLOAT
Matrix density value.
rhob_clean : FLOAT
RhoB clean one point.
rhob_clay : FLOAT
Clay density value.
nphi_matrix : FLOAT
Neutron matrix value.
nphi_clean : FLOAT
Neutron clean one point.
nphi_clay : FLOAT
Clay neutron value.
Returns
-------
vcl_nd : ARRAY
Clay volume.
"""
vcl_nd = calculate_dual_clay(rhob, nphi, rhob_matrix, rhob_clean, rhob_clay, nphi_matrix, nphi_clean, nphi_clay)
return vcl_nd
def calculate_sd_clay(rhob, dt, rhob_matrix, rhob_clean, rhob_clay, dt_matrix, dt_clean, dt_clay):
"""
Calculate vclay by sonic density method
Parameters
----------
rhob : ARRAY
Bulk density curve.
dt : ARRAY
Sonic curve.
rhob_matrix : FLOAT
Matrix density value.
rhob_clean : FLOAT
RhoB clean one point.
rhob_clay : FLOAT
Clay density value.
dt_matrix : FLOAT
Sonic matrix value.
dt_clean : FLOAT
Sonic clean one point.
dt_clay : FLOAT
Clay sonic value.
Returns
-------
vcl_sd : ARRAY
Clay volume.
"""
vcl_sd = calculate_dual_clay(rhob, dt, rhob_matrix, rhob_clean, rhob_clay, dt_matrix, dt_clean, dt_clay)
return vcl_sd
def calculate_ns_clay(nphi, dt, nphi_matrix, nphi_clean, nphi_clay, dt_matrix, dt_clean, dt_clay):
"""
Calculate vclay by neutron sonic method
Parameters
----------
nphi : ARRAY
Neutron curve.
dt : ARRAY
Sonic curve.
nphi_matrix : FLOAT
Matrix neutron value.
nphi_clean : FLOAT
Neutron clean one point.
nphi_clay : FLOAT
Clay neutron value.
dt_matrix : FLOAT
Sonic matrix value.
dt_clean : FLOAT
Sonic clean one point.
dt_clay : FLOAT
Clay sonic value.
Returns
-------
vcl_ns : ARRAY
Clay volume.
"""
vcl_ns = calculate_dual_clay(nphi, dt, nphi_matrix, nphi_clean, nphi_clay, dt_matrix, dt_clean, dt_clay)
return vcl_ns
def calculate_vwcl(*args):
"""
Calculate the minimum from multiple input clay volume curves
Parameters
----------
*args : ARRAY
Two or more clay volume curves.
Returns
-------
vwcl : ARRAY
Volume wet clay.
"""
if len(args) == 1:
vwcl = args
return vwcl
else:
for i in range(len(args)-1):
vwcl = np.minimum(args[i], args[i+1])
return vwcl
#################################
### TOC Analysis Functions ###
#################################
def calculate_schmoker_toc(rhob, schmoker_a=154.497, schmoker_b=57.261):
"""
Calculate TOC wt. % via Schmoker method
Parameters
----------
rhob : ARRAY
Bulk density curve.
schmoker_a : FLOAT, optional
Schmoker a parameter. The default is 154.497.
schmoker_b : FLOAT, optional
Schmoker b parameter. The default is 57.261.
Returns
-------
toc : Array
Total Organic Carbon (wt %).
"""
if np.isscalar(schmoker_a):
temp = np.zeros(len(rhob))
temp[:] = schmoker_a
schmoker_a = temp
if np.isscalar(schmoker_b):
temp = np.zeros(len(rhob))
temp[:] = schmoker_b
schmoker_b = temp
assert len(schmoker_a) == len(rhob), "Input schmoker_a array must be of same length as rhob curve."
assert len(schmoker_b) == len(rhob), "Input schmoker_b array must be of same length as rhob curve."
toc = ((schmoker_a/rhob-schmoker_b)*.01)
toc = np.minimum(np.maximum(toc, 0),1)
return toc
def calculate_lom(vre):
"""
Calculate level of Maturity from vitrinite reflectance
Parameters
----------
vre : FLOAT
Vitrinite Reflectance.
Returns
-------
FLOAT
Level of Maturity.
"""
return (-8.3180734+(53.367852*vre)-(45.319273*vre**2)+(13.422786*vre**3)) / (1+(1.0895416*vre)-(1.4275785*vre**2)+(0.54231042*vre**3))
def calculate_mod_passey_toc(resistivity, porosity, vre, rwno_clay):
"""
Calculate TOC wt% via modified Passey method
Parameters
----------
resistivity : ARRAY
Deep resistivity curve.
porosity : ARRAY
Total porosity curve.
vre : FLOAT
Vitrinite reflectance.
rwno_clay : FLOAT
Resistivity water in non-organic clay interval.
Returns
-------
toc : ARRAY
Total Organic Carbon (wt %).
"""
assert len(resistivity) == len(porosity), "Resistivity and porosity arrays must be of same length"
if np.isscalar(vre):
temp = np.zeros(len(porosity))
temp[:] = vre
vre = temp
if np.isscalar(rwno_clay):
temp = np.zeros(len(porosity))
temp[:] = rwno_clay
rwno_clay = temp
assert len(vre) == len(porosity), "Input vre array must be of same length as curves."
assert len(rwno_clay) == len(porosity), "Input rwno_clay array must be of same length as curves."
toc_lom = calculate_lom(vre)
toc_ro = rwno_clay / porosity
toc = ((np.log10(resistivity) - np.log10(toc_ro)) * 10 ** (2.297-0.1688*toc_lom))/100
toc = np.minimum(np.maximum(toc, 0),1)
return toc
def calculate_faust_toc(dt, resistivity, depth_curve, fm=0.21):
"""
Calculate TOC wt % by Faust method
Parameters
----------
dt : ARRAY
Sonic curve.
resistivity : ARRAY
Deep Resistivity curve.
depth_curve : ARRAY
Depth curve.
fm : FLOAT, optional
Faust Modifier. The default is 0.21.
Returns
-------
toc : ARRAY
Total Organic Carbon (wt %).
"""
assert len(depth_curve) == len(resistivity) == len(dt), "All curves must be of same length."
if np.isscalar(fm):
temp = np.zeros(len(depth_curve))
temp[:] = fm
fm = temp
assert len(depth_curve) == len(fm), "Input fm array must be same length as input curves."
dt_faust = 513.347/(depth_curve*resistivity)**(1/6)
first = fm * ((dt - 56) / (189 - 56))
last = (dt_faust - 56) / (189 - 56)
toc = first - last
toc = np.minimum(np.maximum(toc, 0),1)
return toc
def calculate_passey_toc(resistivity, curve, resistivity_baseline, curve_baseline, vre, method):
"""
Calculate TOC wt% by designated Passey overlay method
Parameters
----------
resistivity : ARRAY
Deep Resistivity curve.
curve : ARRAY
Either a sonic, density, or neutron curve.
resistivity_baseline : FLOAT
Baseline value for resistivity in an inorganic shale interval.
curve_baseline : FLOAT
Baseline value for selected curve in an inorganic shale interval.
vre : FLOAT
Vitrinite reflectance.
method : STRING
Options are 'dt', 'rhob', or 'nphi'. Should correspond to curve entered.
Returns
-------
toc : ARRAY
Total Organic Carbon (wt %).
"""
assert len(resistivity) == len(curve), "All curves must be of same length."
assert method in ['dt', 'rhob', 'nphi'], "Method not recognized."
if np.isscalar(resistivity_baseline):
temp = np.zeros(len(curve))
temp[:] = resistivity_baseline
resistivity_baseline = temp
if np.isscalar(curve_baseline):
temp = np.zeros(len(curve))
temp[:] = curve_baseline
curve_baseline = temp
if np.isscalar(vre):
temp = np.zeros(len(curve))
temp[:] = vre
vre = temp
assert len(curve) == len(resistivity_baseline), "Input resistivity_baseline array must be same length as input curves."
assert len(curve) == len(curve_baseline), "Input curve_baseline array must be same length as input curves."
assert len(curve) == len(vre), "Input vre array must be same length as input curves."
toc_lom = calculate_lom(vre)
if method == 'dt':
xlogr = np.log10(resistivity/resistivity_baseline)+0.02 *(curve - curve_baseline)
elif method == 'rhob':
xlogr = np.log10(resistivity/resistivity_baseline)+2.5 *(curve - curve_baseline)
elif method == 'nphi':
xlogr = np.log10(resistivity/resistivity_baseline)+4.0 *(curve - curve_baseline)
toc = (xlogr*10**(2.297-0.1688*toc_lom))/100
toc = np.minimum(np.maximum(toc, 0),1)
return toc
def calculate_volume_kerogen(toc, ktoc = 0.80, method='slb', matrix_density = 2.65, dens_kerogen = 1.26):
"""
Convert TOC in wt. % to kerogen vol. %.
Parameters
----------
toc : ARRAY
Total organic carbon curve.
ktoc : FLOAT, optional
Kerogen to TOC ratio. The default is 0.80.
method : STRING, optional
Conversion method. Options are 'slb', 'danomics', 'standard'. The default is 'slb'.
matrix_density : FLOAT, optional
Matrix density. The default is 2.65.
dens_kerogen : FLOAT, optional
Kergogen density. The default is 1.26.
Returns
-------
v_ker : TYPE
DESCRIPTION.
"""
assert method in ['slb', 'danomics', 'standard'], "Method not recognized."
if np.isscalar(ktoc):
temp = np.zeros(len(toc))
temp[:] = ktoc
ktoc = temp
if np.isscalar(matrix_density):
temp = np.zeros(len(toc))
temp[:] = matrix_density
matrix_density = temp
if np.isscalar(dens_kerogen):
temp = np.zeros(len(toc))
temp[:] = dens_kerogen
dens_kerogen = temp
assert len(toc) == len(ktoc), "Input ktoc array must be of same length as toc curve."
assert len(toc) == len(matrix_density), "Input matrix_density array must be of same length as toc curve."
assert len(toc) == len(dens_kerogen), "Input dens_kerogen array must be of same length as toc curve."
if method == 'slb':
v_ker = (toc * ktoc * matrix_density / dens_kerogen)
elif method == 'danomics':
v_ker = ((toc/ktoc)/dens_kerogen)/(((toc/ktoc)/dens_kerogen) + ((1-(toc/ktoc))/matrix_density))
elif method == "standard":
v_ker = (toc * matrix_density / dens_kerogen)
return v_ker
#################################
### Porosity Functions ###
#################################
def calculate_density_porosity(rhob, rho_matrix=2.65, rho_fluid=1, phi_max=0.30):
"""
Calculate density porosity
Parameters
----------
rhob : ARRAY
Bulk Density Curve.
rho_matrix : FLOAT, optional
Matrix density. The default is 2.65.
rho_fluid : FLOAT, optional
Fluid density. The default is 1.
phi_max : Float, optional
Maximum allowed porosity. The default is 0.30.
Returns
-------
phit : ARRAY
Total Porosity.
"""
if np.isscalar(rho_matrix):
temp = np.zeros(len(rhob))
temp[:] = rho_matrix
rho_matrix = temp
if np.isscalar(rho_fluid):
temp = np.zeros(len(rhob))
temp[:] = rho_fluid
rho_fluid = temp
if np.isscalar(phi_max):
temp = np.zeros(len(rhob))
temp[:] = phi_max
phi_max = temp
assert len(rhob) == len(rho_matrix), "Input rho_matrix array must be of same length as rhob curve."
assert len(rhob) == len(rho_fluid), "Input rho_fluid array must be of same length as rhob curve."
assert len(rhob) == len(phi_max), "Input phi_max array must be of same length as rhob curve."
phit = (rho_matrix - rhob) / (rho_matrix - rho_fluid)
phit = np.maximum(np.minimum(phit, phi_max), 0.0001)
return phit
def calculate_nd_porosity(rhob, nphi, rho_matrix=2.65, rho_fluid=1, phi_max=0.30):
"""
Calculate neutron-density crossplot porosity
Parameters
----------
rhob : ARRAY
Bulk density curve.
nphi : ARRAY
Neutron porosity curve.
rho_matrix : FLOAT, optional
Matrix density. The default is 2.65.
rho_fluid : FLOAT, optional
Fluid density. The default is 1.
phi_max : TYPE, optional
Maximum allowed porosity. The default is 0.30.
Returns
-------
phit : ARRAY
Total porosity.
"""
assert len(rhob) == len(nphi), "Rhob and Nphi curves must be of same length."
if np.isscalar(rho_matrix):
temp = np.zeros(len(rhob))
temp[:] = rho_matrix
rho_matrix = temp
if np.isscalar(rho_fluid):
temp = np.zeros(len(rhob))
temp[:] = rho_fluid
rho_fluid = temp
if np.isscalar(phi_max):
temp = np.zeros(len(rhob))
temp[:] = phi_max
phi_max = temp
assert len(rhob) == len(rho_matrix), "Input rho_matrix array must be of same length as rhob curve."
assert len(rhob) == len(rho_fluid), "Input rho_fluid array must be of same length as rhob curve."
assert len(rhob) == len(phi_max), "Input phi_max array must be of same length as rhob curve."
dphi = (rho_matrix - rhob) / (rho_matrix - rho_fluid)
dphi = np.maximum(np.minimum(dphi, phi_max), 0.0001)
phit = (nphi + dphi)/2
phit = np.maximum(np.minimum(phit, phi_max), 0.0001)
return phit
def calculate_sonic_porosity(dt, dt_matrix=56, dt_fluid=189, dt_shale=120, phi_max = 0.30):
"""
Calculate sonic porosity.
Parameters
----------
dt : ARRAY
Compressional sonic curve.
dt_matrix : ARRAY, optional
Sonic matrix value. The default is 56.
dt_fluid : FLOAT, optional
Sonic fluid value. The default is 189.
dt_shale : FLOAT, optional
Sonic value in shale. The default is 120.
phi_max : FLOAT, optional
Maximum allowed porosity. The default is 0.30.
Returns
-------
phit : ARRAY
Total porosity.
"""
if np.isscalar(dt_matrix):
temp = np.zeros(len(dt))
temp[:] = dt_matrix
dt_matrix = temp
if np.isscalar(dt_fluid):
temp = np.zeros(len(dt))
temp[:] = dt_fluid
dt_fluid = temp
if np.isscalar(dt_shale):
temp = np.zeros(len(dt))
temp[:] = dt_shale
dt_shale = temp
if np.isscalar(phi_max):
temp = np.zeros(len(dt))
temp[:] = phi_max
phi_max = temp
assert len(dt) == len(dt_matrix), "Input dt_matrix array must be of same length as dt curve."
assert len(dt) == len(dt_fluid), "Input dt_fluid array must be of same length as dt curve."
assert len(dt) == len(dt_shale), "Input dt_shale array must be of same length as dt curve."
assert len(dt) == len(phi_max), "Input phi_max array must be of same length as dt curve."
phit = ((dt - dt_matrix) / (dt_fluid - dt_matrix) / np.maximum(1, dt_shale/100))
phit = np.maximum(np.minimum(phit, phi_max), 0.0001)
return phit
def calculate_deterministic_grain_density(rhob, nphi, vclay, rho_matrix=2.65, rho_fluid=1, phi_shale = 0.10):
"""
Calculate a deterministic grain density.
Parameters
----------
rhob : ARRAY
DESCRIPTION.
nphi : ARRAY
DESCRIPTION.
vclay : ARRAY
DESCRIPTION.
rho_matrix : FLOAT, optional
DESCRIPTION. The default is 2.65.
rho_fluid : FLOAT, optional
DESCRIPTION. The default is 1.
phi_shale : FLOAT, optional
DESCRIPTION. The default is 0.10.
Returns
-------
grain_density : ARRAY
Deterministic grain_density.
"""
assert len(rhob) == len(nphi) == len(vclay), "All input curves must be of same length."
if np.isscalar(rho_matrix):
temp = np.zeros(len(rhob))
temp[:] = rho_matrix
rho_matrix = temp
if np.isscalar(rho_fluid):
temp = np.zeros(len(rhob))
temp[:] = rho_fluid
rho_fluid = temp
if np.isscalar(phi_shale):
temp = np.zeros(len(rhob))
temp[:] = phi_shale
phi_shale = temp
assert len(rhob) == len(rho_matrix), "Input rho_matrix array must be of same length as rhob curve."
assert len(rhob) == len(rho_fluid), "Input rho_fluid array must be of same length as rhob curve."
assert len(rhob) == len(phi_shale), "Input phi_max array must be of same length as rhob curve."
dphi = (rho_matrix - rhob) / (rho_matrix - rho_fluid)
phix = np.sqrt((nphi**2 + dphi**2) / 2) - (vclay * phi_shale)
grain_density = (rhob-phix)/(1-phix)
return grain_density
def calculate_grain_density_porosity(rhob, grain_density, rho_fluid=1, phi_max=0.30):
"""
Calculate density porosity using grain density for rho_matrix.
Parameters
----------
rhob : ARRAY
Bulk Density Curve.
grain_density : ARRAY
Grain Density Curve.
rho_fluid : FLOAT, optional
Fluid density. The default is 1.
phi_max : Float, optional
Maximum allowed porosity. The default is 0.30.
Returns
-------
phit : ARRAY
Total Porosity.
"""
assert len(rhob) == len(grain_density), "Rhob and Grain Density curves must be of same length."
if np.isscalar(rho_fluid):
temp = np.zeros(len(rhob))
temp[:] = rho_fluid
rho_fluid = temp
if np.isscalar(phi_max):
temp = np.zeros(len(rhob))
temp[:] = phi_max
phi_max = temp
assert len(rhob) == len(rho_fluid), "Input rho_fluid array must be of same length as rhob curve."
assert len(rhob) == len(phi_max), "Input phi_max array must be of same length as rhob curve."
phit = (grain_density - rhob) / (grain_density - rho_fluid)
phit = np.maximum(np.minimum(phit, phi_max), 0.0001)
return phit
def calculate_phie_porosity(phit, vclay, phi_shale=0.10, phi_max = 0.30):
"""
Clay correct total porosity to effective porosity.
Parameters
----------
phit : ARRAY
Total porosity.
vclay : ARRAY
Clay volume.
phi_shale : FLOAT, optional
Porosity fraction of shale. The default is 0.10.
phi_max : FLOAT, optional
Maximum allowed porosity. The default is 0.30.
Returns
-------
phie : ARRAY
Effective Porosity.
"""
assert len(phit) == len(vclay), "Input curves must be of same length."
if np.isscalar(phi_shale):
temp = np.zeros(len(phit))
temp[:] = phi_shale
phi_shale = temp
if np.isscalar(phi_max):
temp = np.zeros(len(phit))
temp[:] = phi_max
phi_max = temp
assert len(phit) == len(phi_max), "Input phi_max array must be of same length as phit curve."
assert len(phit) == len(phi_shale), "Input phi_shale array must be of same length as phit curve."
phie = phit - (vclay * phi_shale)
phie = np.minimum(np.maximum(phie, 0.0001), phi_max)
return phie
#################################
### Water Saturation Functions###
#################################
def calculate_temperature_by_gradient(depth_curve, surf_temp = 55, gradient = 1.1):
"""
Calculate Temperature Curve from geothermal gradient
Parameters
----------
depth_curve : ARRAY
Depth curve in feet.
surf_temp : FLOAT, optional
Surface temperature at well location. The default is 55.
gradient : FLOAT, optional
Geothermal gradient in degrees Farenheit per 100 ft. The default is 1.1.
Returns
-------
temperature : ARRAY
TEMPERATURE curve in degrees Farenheit.
"""
if np.isscalar(surf_temp):
temp = np.zeros(len(depth_curve))
temp[:] = surf_temp
surf_temp = temp
if np.isscalar(gradient):
temp = np.zeros(len(depth_curve))
temp[:] = gradient
gradient = temp
assert len(depth_curve) == len(surf_temp), "Input surf_temp array must be of same length as depth curve."
assert len(depth_curve) == len(gradient), "Input gradient array must be of same length as depth curve."
temperature = (surf_temp + depth_curve * gradient/100)
return temperature
def calculate_rw_ft(temperature, rw_surf, surf_temp):
"""
Calculate the formation water resistivity at formation temperature.
Parameters
----------
rw_surf : FLOAT
Formation water resistivity at reference surface temperature.
surf_temp : FLOAT
Reference surface temperature for Rw.
temperature : ARRAY
Temperature curve.
Returns
-------
rw_ft : ARRAY
Formation water resistivity curve.
"""
if np.isscalar(surf_temp):
temp = np.zeros(len(temperature))
temp[:] = surf_temp
surf_temp = temp
if np.isscalar(rw_surf):
temp = np.zeros(len(temperature))
temp[:] = rw_surf
rw_surf = temp
assert len(temperature) == len(surf_temp), "Input surf_temp array must be of same length as temperature curve."
assert len(temperature) == len(rw_surf), "Input rw_surf array must be of same length as temperature curve."
rw_ft = (400000 / temperature / ((400000 / surf_temp) / ((rw_surf)**1.14)))**0.88
return rw_ft
def calculate_archie_sw(resistivity, porosity, rw_ft, a=1, m=2, n=2, use_variable_m = False):
"""
Calculate water saturation by Archie's equation.
Parameters
----------
resistivity : ARRAY
Deep resistivity curve.
porosity : ARRAY
Total or effective porosity curve.
rw_ft : ARRAY
Formation water resitivity at formation temperature.
a : FLOAT, optional
Tortuosity factor. The default is 1.
m : FLOAT, optional
Cementation exponent. The default is 2.
n : FLOAT, optional
Saturation exponent. The default is 2.
use_variable_m : BOOL, optional
Use a variable 'm' based on Shell method for cementation exponent. The default is False.
Returns
-------
ARRAY
Water Saturation Curve.
"""
assert len(resistivity) == len(porosity) == len(rw_ft), "All input curves must be of same length."
if np.isscalar(a):
temp = np.zeros(len(porosity))
temp[:] = a
a = temp
if np.isscalar(m):
temp = np.zeros(len(porosity))
temp[:] = m
m = temp
if np.isscalar(n):
temp = np.zeros(len(porosity))
temp[:] = n
n = temp
assert len(porosity) == len(a), "Input a array must have same length as input curves."
assert len(porosity) == len(m), "Input m array must have same length as input curves."
assert len(porosity) == len(n), "Input n array must have same length as input curves."
if use_variable_m is True:
m_curve = 1.87 + (1.9 * porosity)
else:
m_curve = m
sw = ((a * rw_ft/(resistivity * porosity**m_curve))**(1/n))
sw = np.minimum(np.maximum(sw, 0),1)
return sw
def calculate_simandoux_sw(resistivity, porosity, vclay, rw_ft, a=1, m=2, n=2, use_variable_m = False, Rt_clay=4):
"""
Calculate water saturation via Simandoux method.
Parameters
----------
resistivity : ARRAY
Deep resistivity curve.
porosity : ARRAY
Total or effective porosity curve.
vclay : ARRAY
Clay volume curve.
rw_ft: Array
Formation water resitivity at formation temperature.
a : FLOAT, optional
Tortuosity factor. The default is 1.
m : FLOAT, optional
Cementation exponent. The default is 2.
n : FLOAT, optional
Saturation exponent. The default is 2.
use_variable_m : BOOL, optional
Use a variable 'm' based on Shell method for cementation exponent. The default is False.
Rt_clay : FLOAT, optional
Clay resistivity. The default is 4.
Returns
-------
ARRAY
Water Saturation Curve.
"""
assert len(resistivity) == len(porosity) == len(rw_ft) == len(vclay), "All input curves must be of same length."
if np.isscalar(a):
temp = np.zeros(len(porosity))
temp[:] = a
a = temp
if np.isscalar(m):
temp = np.zeros(len(porosity))
temp[:] = m
m = temp
if np.isscalar(n):
temp = np.zeros(len(porosity))
temp[:] = n
n = temp
if np.isscalar(Rt_clay):
temp = np.zeros(len(porosity))
temp[:] = Rt_clay
Rt_clay = temp
assert len(porosity) == len(a), "Input a array must have same length as input curves."
assert len(porosity) == len(m), "Input m array must have same length as input curves."
assert len(porosity) == len(n), "Input n array must have same length as input curves."
assert len(porosity) == len(Rt_clay), "Input Rt_clay array must have same length as input curves."
if use_variable_m is True:
m_curve = 1.87 + (1.9 * porosity)
else:
m_curve = m
sw = (a * rw_ft/(2*porosity**m_curve)) *(np.sqrt(((vclay/Rt_clay)**2) + ((4*porosity**m_curve)/(a * rw_ft* resistivity))) - vclay/Rt_clay)
sw = np.minimum(np.maximum(sw, 0),1)
return sw
def calculate_mod_simandoux_sw(resistivity, porosity, vclay, rw_ft, a=1, m=2, n=2, use_variable_m = False, Rt_clay=4):
"""
Calculate water saturation via Modified Simandoux method.
Parameters
----------
resistivity : ARRAY
Deep resistivity curve.
porosity : ARRAY
Total or effective porosity curve.
vclay : ARRAY
Clay volume curve.
rw_ft: Array
Formation water resitivity at formation temperature.
a : FLOAT, optional
Tortuosity factor. The default is 1.
m : FLOAT, optional
Cementation exponent. The default is 2.
n : FLOAT, optional
Saturation exponent. The default is 2.
use_variable_m : BOOL, optional
Use a variable 'm' based on Shell method for cementation exponent. The default is False.
Rt_clay : FLOAT, optional
Clay resistivity. The default is 4.
Returns
-------
ARRAY
Water Saturation Curve.
"""
assert len(resistivity) == len(porosity) == len(rw_ft) == len(vclay), "All input curves must be of same length."
if np.isscalar(a):
temp = np.zeros(len(porosity))
temp[:] = a
a = temp
if np.isscalar(m):
temp = np.zeros(len(porosity))
temp[:] = m
m = temp
if np.isscalar(n):
temp = np.zeros(len(porosity))
temp[:] = n
n = temp
if np.isscalar(Rt_clay):
temp = np.zeros(len(porosity))
temp[:] = Rt_clay
Rt_clay = temp
assert len(porosity) == len(a), "Input a array must have same length as input curves."
assert len(porosity) == len(m), "Input m array must have same length as input curves."
assert len(porosity) == len(n), "Input n array must have same length as input curves."
assert len(porosity) == len(Rt_clay), "Input Rt_clay array must have same length as input curves."
if use_variable_m is True:
m_curve = 1.87 + (1.9 * porosity)
else:
m_curve = np.zeros(len(porosity))
m_curve[:] = m ###Need to pay attention to this when working with arrays for params
sw = np.empty(len(porosity))
sw[:] = np.nan
for i in range(len(porosity)):
if vclay[i] < 1:
sw[i] = ((((vclay[i]/Rt_clay[i])**2 + 4*porosity[i]**m_curve[i]/(a[i]* rw_ft[i]*(1-vclay[i])*resistivity[i]))**(1/n[i])) - vclay[i]/Rt_clay[i])/(2*porosity[i]**m_curve[i]/(a[i]*rw_ft[i]*(1-vclay[i])))
else:
sw[i] =1
sw = np.minimum(np.maximum(sw, 0),1)
return sw
def calculate_indonesian_sw(resistivity, porosity, vclay, rw_ft, a=1, m=2, n=2, use_variable_m = False, Rt_clay=4):
"""
Calculate water saturation via Indonesian method.
Parameters
----------
resistivity : ARRAY
Deep resistivity curve.
porosity : ARRAY
Total or effective porosity curve.
vclay : ARRAY
Clay volume curve.
rw_ft: Array
Formation water resitivity at formation temperature.
a : FLOAT, optional
Tortuosity factor. The default is 1.
m : FLOAT, optional
Cementation exponent. The default is 2.
n : FLOAT, optional
Saturation exponent. The default is 2.
use_variable_m : BOOL, optional
Use a variable 'm' based on Shell method for cementation exponent. The default is False.
Rt_clay : FLOAT, optional
Clay resistivity. The default is 4.
Returns
-------
ARRAY
Water Saturation Curve.
"""
assert len(resistivity) == len(porosity) == len(rw_ft) == len(vclay), "All input curves must be of same length."
if np.isscalar(a):
temp = np.zeros(len(porosity))
temp[:] = a
a = temp
if np.isscalar(m):
temp = np.zeros(len(porosity))
temp[:] = m
m = temp
if np.isscalar(n):
temp = np.zeros(len(porosity))
temp[:] = n
n = temp
if np.isscalar(Rt_clay):
temp = np.zeros(len(porosity))
temp[:] = Rt_clay
Rt_clay = temp
assert len(porosity) == len(a), "Input a array must have same length as input curves."
assert len(porosity) == len(m), "Input m array must have same length as input curves."
assert len(porosity) == len(n), "Input n array must have same length as input curves."
assert len(porosity) == len(Rt_clay), "Input Rt_clay array must have same length as input curves."
if use_variable_m is True:
m_curve = 1.87 + (1.9 * porosity)
else:
m_curve = m
sw = (((np.sqrt( (vclay**(2-vclay))/Rt_clay ) + np.sqrt( ((porosity**m_curve)/rw_ft) ))**2)*resistivity)**(-1/n)
sw = np.minimum(np.maximum(sw, 0),1)
return sw
def calculate_dual_water_sw(resistivity, porosity, vclay, rw_ft, rwb_ft, a=1, m=2, n=2, use_variable_m = False, Rt_clay=4, phi_shale=0.10):
"""
Calculate water saturation via Dual Water method.
Parameters
----------
resistivity : ARRAY
Deep resistivity curve.
porosity : ARRAY
Total or effective porosity curve.
vclay : ARRAY, optional
Clay volume curve
rw_ft: Array
Formation water resitivity at formation temperature.
rwb_ft: Array
Formation water resitivity at formation temperature.
a : FLOAT, optional
Tortuosity factor. The default is 1.
m : FLOAT, optional
Cementation exponent. The default is 2.
n : FLOAT, optional
Saturation exponent. The default is 2.
use_variable_m : BOOL, optional
Use a variable 'm' based on Shell method for cementation exponent. The default is False.
Rt_clay : FLOAT, optional
Clay resistivity. The default is 4.
phi_shale : FLOAT, optional
Shale porosity fraction. The default is 0.10.
Returns
-------
ARRAY
Water Saturation Curve.
"""
assert len(resistivity) == len(porosity) == len(rw_ft) == len(rwb_ft) == len(vclay), "All input curves must be of same length."
if np.isscalar(a):
temp = np.zeros(len(porosity))
temp[:] = a
a = temp
if np.isscalar(m):
temp = np.zeros(len(porosity))
temp[:] = m
m = temp
if np.isscalar(n):
temp = np.zeros(len(porosity))
temp[:] = n
n = temp
if np.isscalar(Rt_clay):
temp = np.zeros(len(porosity))
temp[:] = Rt_clay
Rt_clay = temp
if np.isscalar(phi_shale):
temp = np.zeros(len(porosity))
temp[:] = phi_shale
phi_shale = temp
assert len(porosity) == len(a), "Input a array must have same length as input curves."
assert len(porosity) == len(m), "Input m array must have same length as input curves."
assert len(porosity) == len(n), "Input n array must have same length as input curves."
assert len(porosity) == len(Rt_clay), "Input Rt_clay array must have same length as input curves."
assert len(porosity) == len(phi_shale), "Input phi_shale array must have same length as input curves."
if use_variable_m is True:
m_curve = 1.87 + (1.9 * porosity)
else:
m_curve = m
sw = (((a**2*rw_ft**2) / ((porosity**m_curve) * ((vclay * phi_shale/porosity) * a*rw_ft + (1-(vclay * phi_shale/porosity)) * a*rwb_ft)))/resistivity)**(1/n)
sw = np.minimum(np.maximum(sw, 0),1)
return sw
def calculate_waxman_smits_sw(resistivity, porosity, rw_ft, temperature, a=1, m=2, n=2, use_variable_m = False, Qv=0.4):
"""
Parameters
----------
resistivity : ARRAY
Deep resistivity curve.
porosity : ARRAY
Total or effective porosity curve.
rw_ft: ARRAY
Formation water resitivity at formation temperature.
temperature: ARRAY
Temperature curve in degrees Farenheit.
a : FLOAT, optional
Tortuosity factor. The default is 1.
m : FLOAT, optional
Cementation exponent. The default is 2.
n : FLOAT, optional
Saturation exponent. The default is 2.
use_variable_m : BOOL, optional
Use a variable 'm' based on Shell method for cementation exponent. The default is False.
Qv : FLOAT, optional
Qv parameter for Waxman Smits. The default is 0.4.
Returns
-------
ARRAY
Water Saturation Curve.
"""
assert len(resistivity) == len(porosity) == len(rw_ft), "All input curves must be of same length."
if np.isscalar(a):
temp = np.zeros(len(porosity))
temp[:] = a
a = temp
if np.isscalar(m):
temp = np.zeros(len(porosity))
temp[:] = m
m = temp
if np.isscalar(n):
temp = np.zeros(len(porosity))
temp[:] = n
n = temp
if np.isscalar(Qv):
temp = np.zeros(len(porosity))
temp[:] = Qv
Qv = temp
assert len(porosity) == len(a), "Input a array must have same length as input curves."
assert len(porosity) == len(m), "Input m array must have same length as input curves."
assert len(porosity) == len(n), "Input n array must have same length as input curves."
assert len(porosity) == len(Qv), "Input Qv array must have same length as input curves."
if use_variable_m is True:
m_curve = 1.87 + (1.9 * porosity)
else:
m_curve = m
sw_start = 0.50
m_star = (m_curve)+ 0.5*(0.258*(Qv*porosity/(1-porosity))+0.2*(1-np.exp(-16.4*(Qv*porosity/(1-porosity)))))
b_curve = (-1.28 + 0.225 * degf2degc(temperature) - 0.0004059 * degf2degc(temperature)**2)/(1 + (rw_ft**1.23)*(0.045*degf2degc(temperature)-0.027))
for i in range(7):
sw_start = ((a*rw_ft/resistivity) /( 1+b_curve*Qv*rw_ft/sw_start) / (porosity**m_star))**(1/n)
sw = sw_start
sw = np.minimum(np.maximum(sw, 0),1)
return sw
def calculate_BVW(sw, phit, phie):
"""
Parameters
----------
sw : ARRAY
Water saturation curve.
phit : ARRAY
Total porosity curve.
phie : ARRAY
Effective porosity curve.
Returns
-------
ARRAY
Total Bulk Volume Water.
ARRAY
Effective Bulk Volume Water.
"""
assert len(sw) == len(phit) == len(phie), "All input curves must be of same length."
bvwt = sw*phit
bvwe = sw*phie
return bvwt, bvwe
def calculate_SWT_SWE(bvwt, bvwe, phit, phie, vclay, phi_shale = 0.10, Sw_Type="total"):
"""
Parameters
----------
bvwt : ARRAY
Total Bulk Volume Water.
bvwe : ARRAY
Effective Bulk Volume Water.
phit : ARRAY
Total porosity curve.
phie : ARRAY
Effective porosity curve.
vclay : ARRAY
Clay volume curve.
phi_shale : FLOAT, optional
Shale porosity fraction. The default is 0.10.
Sw_Type : STRING, optional
Type of saturation calculated Options are total or effective. The default is "total".
Returns
-------
SwT : TYPE
DESCRIPTION.
SwE : TYPE
DESCRIPTION.
"""
assert len(bvwe) == len(bvwt) == len(phit) == len(phie) == len(vclay), "All input curves must be of same length."
if np.isscalar(phi_shale):
temp = np.zeros(len(phit))
temp[:] = phi_shale
phi_shale = temp
assert len(phit) == len(phi_shale), "Input phi_shale array must have same length as input curves."
Sw = bvwt/phit
if Sw_Type == "total":
SwT = Sw
SwE = np.minimum(np.maximum((bvwt-(vclay * phi_shale))/phie,0),1)
else:
SwE = Sw
SwT = np.minimum(np.maximum( (bvwe+(vclay * phi_shale))/phit ,0),1)
return SwT, SwE
######################################
### Pore Pressure Functions ###
######################################
def calculate_hydrostatic_pressure(depth_curve, hpg=0.460):
"""
Calculate the hydrostatic pressure gradient
Parameters
----------
depth_curve : ARRAY
Depth curve in feet.
hpg : FLOAT, optional
Hydrostatic pressure gradient in psi/ft. The default is 0.460.
Returns
-------
hp : ARRAY
Hydrostatic pressure in psi.
"""
if np.isscalar(hpg):
temp = np.zeros(len(depth_curve))
temp[:] = hpg
hpg = temp
assert len(depth_curve) == len(hpg), "Input hpg array must have same length as depth curve."
hp = depth_curve * hpg
return hp
def calculate_overburden_pressure(depth_curve, phit, rho_matrix = 2.65, rho_fluid = 1.0, lpg=1.1):
"""
Calculate the overburden pressure gradient using the porosity.
Rho_matrix can be a float or grain density array.
LPG is used for intervals without no porosity measurement available.
Parameters
----------
depth_curve : ARRAY
Depth curve.
phit : ARRAY
Total porosity curve.
rho_matrix : FLOAT, optional
Matrix density. The default is 2.65.
rho_fluid : FLOAT, optional
Fluid density. The default is 1.0.
lpg : FLOAT, optional
Lithostatic pressure gradient. The default is 1.1.
Returns
-------
obp : ARRAY
Overburden pressure in psi.
"""
assert len(phit) == len(depth_curve), "All input curves must be of same length."
if np.isscalar(rho_matrix):
temp = np.zeros(len(depth_curve))
temp[:] = rho_matrix
rho_matrix = temp
if np.isscalar(rho_fluid):
temp = np.zeros(len(depth_curve))
temp[:] = rho_fluid
rho_fluid = temp
if np.isscalar(lpg):
temp = np.zeros(len(depth_curve))
temp[:] = lpg
lpg = temp
step = np.zeros(len(depth_curve))
for i in range(1,len(depth_curve)):
step[i] = (depth_curve[i] - depth_curve[i-1])
step[0] = step[1]
obp = 0.433 *( (1-phit)*rho_matrix*step + phit*rho_fluid*step)
obp = np.nan_to_num(obp, nan=lpg)
for i in range(1,len(obp)):
obp[i] = np.nansum([obp[i],obp[i-1]])
return obp
def calculate_eaton_sonic_porepressure(depth_curve, dt, vclay, obp, hp,
vclay_cutoff = 0.35,
method="eaton",
dtn=90,
dt_exp = 3.0,
dto = 50.0,
nctl= -6.0 ):
"""
Calculate pore pressure via eaton or modified eaton methodology from sonic log.
Parameters
----------
depth_curve : ARRAY
Depth curve.
dt : ARRAY
Compressional sonic curve in usec/foot.
vclay : ARRAY
Clay volume curve.
obp : ARRAY
Overburden curve.
hp : ARRAY
Hydrostatic pressure curve.
vclay_cutoff : FLOAT, optional
Cutoff value to determine if an interval is clay for analysis. The default is 0.35.
method : STRING, optional
Options are 'eaton' or 'mod_eaton'. The default is "eaton".
dtn : FLOAT, optional
Normal value of sonic transit time. The default is 90.
dt_exp : FLOAT, optional
Eaton's exponent. The default is 3.0.
dto : FLOAT, optional
DTo value in modified eaton NCTL equation. The default is 50.0.
nctl: FLOAT, optional
'b' exponent in modified eaton NCTL equation. The default is -6.0.
Returns
-------
pp : ARRAY
Pore pressure in psi/ft.
"""
assert len(depth_curve) == len(dt) == len(vclay) == len(obp) == len(hp), "All input curves must be of same length."
assert method in ['eaton', 'mod_eaton'], "Input method not recognized."
if np.isscalar(vclay_cutoff):
temp = np.zeros(len(depth_curve))
temp[:] = vclay_cutoff
vclay_cutoff = temp
if np.isscalar(dtn):
temp = np.zeros(len(depth_curve))
temp[:] = dtn
dtn = temp
if np.isscalar(dt_exp):
temp = np.zeros(len(depth_curve))
temp[:] = dt_exp
dt_exp = temp
if np.isscalar(dto):
temp = np.zeros(len(depth_curve))
temp[:] = dto
dto = temp
if np.isscalar(nctl):
temp = np.zeros(len(depth_curve))
temp[:] = nctl
nctl= temp
assert len(depth_curve) == len(vclay_cutoff), "Vclay cutoff array must be same length as curves."
assert len(depth_curve) == len(dtn), "dtnarray must be same length as curves."
assert len(depth_curve) == len(dt_exp), "dt_exp array must be same length as curves."
assert len(depth_curve) == len(dto), "dto array must be same length as curves."
assert len(depth_curve) == len(nctl), "nctlarray must be same length as curves."
if method == "eaton":
dt_nctl = dtn
pp = np.zeros(len(depth_curve))
for i in range(len(vclay)):
if vclay[i] > vclay_cutoff:
pp[i] = (obp[i] - ((obp[i] - hp[i]) * (dt_nctl / dt[i])**dt_exp))
else:
pp[i] = hp[i]
elif method == "mod_eaton":
dt_nctl = (dto + (185-dto) * np.exp((nctl/10000*depth_curve/3.28)))
pp = np.zeros(len(depth_curve))
for i in range(len(vclay)):
if vclay[i] > vclay_cutoff:
pp[i] = (obp[i] - ((obp[i] - hp[i]) * (dt_nctl[i] / dt[i])**dt_exp))
else:
pp[i] = hp[i]
pp = np.maximum(hp, pp)
return pp
def calculate_eaton_resistivity_porepressure(depth_curve, resistivity, vclay, obp, hp,
vclay_cutoff = 0.35,
method = 'eaton',
res_n = 3.0,
res_exp = 1.2,
res_a = 1,
nctl = 0.34):
"""
Calculate pore pressure via eaton or modified eaton methodology from resistivity log.
Parameters
----------
depth_curve : ARRAY
Depth curve.
resistivity : ARRAY
Deep resistivity curve.
vclay : ARRAY
Clay volume curve.
obp : ARRAY
Overburden curve.
hp : ARRAY
Hydrostatic pressure curve.
vclay_cutoff : FLOAT, optional
Cutoff value to determine if an interval is clay for analysis. The default is 0.35.
method : STRING, optional
Options are 'eaton' or 'mod_eaton'. The default is "eaton".
res_n : FLOAT, optional
Normal value of resistivity. The default is 3.0.
res_exp : FLOAT, optional
Eaton exponent. The default is 1.2.
res_a : FLOAT, optional
Initial resistivity value in NCTL for modified eaton. The default is 1.
nctl : FLOAT, optional
C exponent in modified eaton NCTL.. The default is 0.34.
Returns
-------
pp : TYPE
DESCRIPTION.
"""
assert len(depth_curve) == len(resistivity) == len(vclay) == len(obp) == len(hp), "All input curves must be of same length."
assert method in ['eaton', 'mod_eaton'], "Input method not recognized."
if np.isscalar(vclay_cutoff):
temp = np.zeros(len(depth_curve))
temp[:] = vclay_cutoff
vclay_cutoff = temp
if np.isscalar(res_n):
temp = np.zeros(len(depth_curve))
temp[:] = res_n
res_n = temp
if np.isscalar(res_exp):
temp = np.zeros(len(depth_curve))
temp[:] = res_exp
res_exp = temp
if np.isscalar(res_a):
temp = np.zeros(len(depth_curve))
temp[:] = res_a
res_a = temp
if np.isscalar(nctl):
temp = np.zeros(len(depth_curve))
temp[:] = nctl
nctl = temp
assert len(depth_curve) == len(vclay_cutoff), "Vclay cutoff array must be same length as curves."
assert len(depth_curve) == len(res_n), "res_n array must be same length as curves."
assert len(depth_curve) == len(res_exp), "res_exp array must be same length as curves."
assert len(depth_curve) == len(res_a), "res_a array must be same length as curves."
assert len(depth_curve) == len(nctl), "nctl array must be same length as curves."
if method == "eaton":
res_nctl = res_n
pp = np.zeros(len(depth_curve))
for i in range(len(vclay)):
if vclay[i] > vclay_cutoff:
pp[i] = (obp[i] - ((obp[i] - hp[i]) * (resistivity[i] / res_nctl)**res_exp))
else:
pp[i] = hp[i]
elif method == "mod_eaton":
res_nctl = (res_a*np.exp((nctl/10000*depth_curve/3.28)))
pp = np.zeros(len(depth_curve))
for i in range(len(vclay)):
if vclay[i] > vclay_cutoff:
pp[i] = (obp[i] - ((obp[i] - hp[i]) * (resistivity[i] / res_nctl[i])**res_exp))
else:
pp[i] = hp[i]
pp = np.maximum(hp, pp)
return pp
def calculate_overpressure(pp, hp):
"""
Calculate overpressure as difference between pore pressure and hydrostatic pressure.
Parameters
----------
pp : ARRAY
Pore pressure curve.
hp : ARRAY
Hydrostatic pressure curve.
Returns
-------
overpressure : ARRAY
Overpressure curve.
"""
assert len(pp) == len(hp), "Input curves must be of same length."
overpressure = np.subtract(pp,hp)
return overpressure
######################################
### Shear Log Modelling Functions ###
######################################
def calculate_vp(dt):
"""
Calculate P-wave velocity.
Parameters
----------
dt : ARRAY
Compressional sonic log.
Returns
-------
ARRAY
P-wave velocity.
"""
return fs2kms(usf2fs(dt))
def calculate_vs(dts):
"""
Calculate S-wave velocity
Parameters
----------
dts : ARRAY
Shear sonic log.
Returns
-------
ARRAY
S-wave velocity.
"""
return fs2kms(usf2fs(dts))
def calculate_gc_ls_shear(vp):
"""
Calculate Shear log via Greenberg-Castagna limestone end-member model.
Parameters
----------
vp : ARRAY
P-wave velocity.
Returns
-------
ARRAY
Synthetic shear log.
"""
return usf2fs(kms2fs(-0.05508 * vp**2 + 1.01677 * vp - 1.03049))
def calculate_gc_dol_shear(vp):
"""
Calculate Shear log via Greenberg-Castagna dolomite end-member model.
Parameters
----------
vp : ARRAY
P-wave velocity.
Returns
-------
ARRAY
Synthetic shear log.
"""
return usf2fs(kms2fs(0.58321 * vp - 0.07775))
def calculate_gc_shale_shear(vp):
"""
Calculate Shear log via Greenberg-Castagna shale end-member model.
Parameters
----------
vp : ARRAY
P-wave velocity.
Returns
-------
ARRAY
Synthetic shear log.
"""
return usf2fs(kms2fs(0.7969 * vp - 0.86735))
def calculate_gc_mudrock_shear(vp):
"""
Calculate Shear log via Greenberg-Castagna mudrock line end-member model.
Parameters
----------
vp : ARRAY
P-wave velocity.
Returns
-------
ARRAY
Synthetic shear log.
"""
return usf2fs(kms2fs(0.8621 * vp - 1.1724))
def calculate_gc_ss_shear(vp):
"""
Calculate Shear log via Greenberg-Castagna sandstone end-member model.
Parameters
----------
vp : ARRAY
P-wave velocity.
Returns
-------
ARRAY
Synthetic shear log.
"""
return usf2fs(kms2fs(0.80416 * vp - 0.85588))
######################################
### Geomechanics Functions ###
######################################
def calculate_youngs_modulus(rhob, dt, dts):
"""
Calculate the Young's Modulus.
Parameters
----------
rhob : ARRAY
Bulk density curve.
dt : ARRAY
Compressional sonic curve.
dts : ARRAY
Shear sonic curve.
Returns
-------
youngs_modulus : ARRAY
Young's modulus.
"""
assert len(rhob) == len(dt) == len(dts), "All input curves must be of same length"
shear_modulus = 13474.45 * (rhob / dts**2)
bulk_modulus = 13474.45 * (rhob / (dt**2)) - (4/3) * shear_modulus
youngs_modulus = (9 * shear_modulus * bulk_modulus) / (shear_modulus + 3 * bulk_modulus)
return youngs_modulus
def calculate_bulk_modulus(rhob, dt, dts):
"""
Calculate the bulk modulus.
Parameters
----------
rhob : ARRAY
Bulk density curve.
dt : ARRAY
Compressional sonic curve.
dts : ARRAY
Shear sonic curve.
Returns
-------
bulk_modulus : ARRAY
Bulk modulus.
"""
assert len(rhob) == len(dt) == len(dts), "All input curves must be of same length"
shear_modulus = 13474.45 * (rhob / dts**2)
bulk_modulus = 13474.45 * (rhob / (dt**2)) - (4/3) * shear_modulus
return bulk_modulus
def calculate_shear_modulus(rhob, dts):
"""
Calculate the shear modulus
Parameters
----------
rhob : ARRAY
Bulk density curve.
dts : ARRAY
Shear sonic curve.
Returns
-------
shear_modulus : ARRAY
Shear Modulus.
"""
assert len(rhob) == len(dts), "All input curves must be of same length"
shear_modulus = 13474.45 * (rhob / dts**2)
return shear_modulus
def calculate_poissons_ratio(rhob, dt, dts):
"""
Calculate poissons ratio.
Parameters
----------
rhob : ARRAY
Bulk density curve.
dt : ARRAY
Compressional sonic curve.
dts : ARRAY
Shear sonic curve.
Returns
-------
poissons_ratio : ARRAY
Poisson's ratio.
"""
assert len(rhob) == len(dt) == len(dts), "All input curves must be of same length"
shear_modulus = 13474.45 * (rhob / dts**2)
bulk_modulus = 13474.45 * (rhob / (dt**2)) - (4/3) * shear_modulus
poissons_ratio = (3 * bulk_modulus - 2 * shear_modulus) / (6 * bulk_modulus + 2 * shear_modulus)
return poissons_ratio
