<jupyter_start><jupyter_text>Volumetrics: HCIP calculationWe'll implement the volumetric equation:$$ V = A \times T \times G \times \phi \times N\!\!:\!\!G \times S_\mathrm{O} \times \frac{1}{B_\mathrm{O}} $$ Gross rock volume $$ \mathrm{GRV} = A \times T $$<jupyter_code>thick = 80            # metres
area = 20000 * 30000  # metres

grv = thick * area
grv<jupyter_output><empty_output><jupyter_text>Wouldn't it be cool if we could carry units around with our calculations? With [`pint`](https://pint.readthedocs.io/en/latest/index.html), we can!<jupyter_code>import pint

unit = pint.UnitRegistry()

thick = 80 * unit.m
area = 20000 * unit.m * 30000 * unit.m
grv = thick * area
grv<jupyter_output><empty_output><jupyter_text>Make a function that computes the GRV.<jupyter_code>def calculate_grv(thickness, area):
    return thickness * area<jupyter_output><empty_output><jupyter_text>Now we can just call this function, instead of remembering the equation. (Admittedly, the equation is rather easy to remember in this case!)<jupyter_code>calculate_grv(thick, area)<jupyter_output><empty_output><jupyter_text>It works!Now we need to compensate for the prospect not being a flat slab of rock &mdash; using the geometric factor. We will implement the equations implied by this diagram:<jupyter_code>top = input("What shape is the prospect? ")
height = 100 * unit.m
ratio = thick / height

if top == 'round':
    g = -0.6 * ratio + 1

if top == 'flat':
    g = -0.3 * ratio + 1

if top == 'slab':
    g = 1

g<jupyter_output><empty_output><jupyter_text>ExerciseTurn this one into a function.<jupyter_code>def geometric_factor(thick, height, top='slab'):

    # Your code here
    
    return g
def geometric_factor(thick, height, top='slab'):

    ratio = thick / height

    if top == 'round':
        g = -0.6 * ratio + 1
    elif top == 'flat':
        g = -0.3 * ratio + 1
    else:
        g = 1

    return g
geometric_factor(thick, height=100*unit.m, top='round')
grv *= g<jupyter_output><empty_output><jupyter_text>What if we have multiple prospects?<jupyter_code>import numpy as np

thicknesses = np.array([10, 25, 15, 5, 100]) * unit.m
heights = np.array([75, 100, 20, 100, 200]) * unit.m
geometric_factor(thicknesses, heights, top='round')<jupyter_output><empty_output><jupyter_text>It works!  HC pore volumeWe need:- net:gross &mdash; the ratio of reservoir-quality rock thickness to the total thickness of the interval.- porosity- $S_\mathrm{O}$ &mdash; the oil saturation, or proportion of oil to total pore fluid.<jupyter_code>netg = 0.5   # fraction
por = 0.24   # fraction
s_o = 0.8    # fraction

hcpv = netg * por * s_o
hcpv<jupyter_output><empty_output><jupyter_text>We'll leave that as a fraction for now.  EXERCISE- Turn this into a function by rearranging the following lines of code:<jupyter_code>"""A function to compute the hydrocarbon pore volume."""
return hcpv
hcpv = netg * por * s_o
def calculate_hcpv(netg, por, s_o):
# Put your code here:


def calculate_hcpv(netg, por, s_o):
    """A function to compute the hydrocarbon pore volume."""
    hcpv = netg * por * s_o
    return hcpv<jupyter_output><empty_output><jupyter_text>Formation volume factorOil shrinks when we produce it, especially if it has high GOR. The FVF, or $B_O$, is the ratio of a reservoir barrel to a stock-tank barrel (25 deg C and 1 atm). Typically the FVF is between 1 (heavy oil) and 1.7 (high GOR).<jupyter_code>fvf = 1.1<jupyter_output><empty_output><jupyter_text>EXERCISEFor gas, $B_\mathrm{G}$ is $0.35 Z T / P$, where $Z$ is the correction factor, or gas compressibility factor. $T$ should be in kelvin and $P$ in kPa. $Z$ is usually between 0.8 and 1.2, but it can be as low as 0.3 and as high as 2.0.Can you write a function to calculate $B_\mathrm{G}$?<jupyter_code>def calculate_Bg(     ):  # Add the arguments.
    """Write a docstring."""
    
    
    return         # Don't forget to return something!
def calculate_Bg(Z=1, T=273.15, P=101.325):
    """
    Compute B_G from correction factor Z, temperature T (K),
    and pressure P (kPa).
    """
    Bg = 0.35 * Z * T / P
    return Bg
calculate_Bg(T=293*unit.K, P=1000*unit.kPa)<jupyter_output><empty_output><jupyter_text>Put it all togetherNow we have the components of the volumetric equation:<jupyter_code>hcip = grv * hcpv / fvf
hcip<jupyter_output><empty_output><jupyter_text>Pint can convert to other units, e.g. Imperial barrels, for us.<jupyter_code>hcip.to('imperial_barrel')<jupyter_output><empty_output><jupyter_text>An Imperial barrel is 43 gallons ([Wikipedia](https://en.wikipedia.org/wiki/Barrel_(unit)), whereas an oil barrel is only 42 gallons. [For more on conversion to bbl, BOE, etc.](https://en.wikipedia.org/wiki/Barrel_of_oil_equivalent).So let's define a custom unit:<jupyter_code>unit.define('oil_barrel = 42 gallon = bbl')
hcip.to('bbl')<jupyter_output><empty_output><jupyter_text>EXERCISECan you write a function to compute the volume (i.e. the HCIP), given all the inputs?Try to use the functions you have already written.<jupyter_code># Put your code here.


    
# Solution.
def calculate_hcip(thickness, area, height, top, netg, por, s_o, fvf):
    grv = calculate_grv(thickness, area)
    g = geometric_factor(thickness, height, top)
    grv *= g
    hcpv = calculate_hcpv(netg, por, s_o)
    return grv * hcpv / fvf
calculate_hcip(thick, area, height, top, netg, por, s_o, fvf)<jupyter_output><empty_output><jupyter_text>Monte Carlo simulationWe can easily draw randomly from distributions of properties:- Normal: https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.normal.html- Uniform: https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.uniform.html- Lognormal: https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.lognormal.html<jupyter_code>import numpy as np

netg = np.random.normal(loc=0.5, scale=0.1, size=100)
por = np.random.normal(loc=0.15, scale=0.025, size=100)
import matplotlib.pyplot as plt

_ = plt.hist(por)
import seaborn as sns

sns.distplot(por, rug=True)
hcpv = calculate_hcpv(netg, por, s_o)
hcpv
_ = plt.hist(hcpv, bins=20)
hcpv.mean()<jupyter_output><empty_output><jupyter_text>The histogram looks a bit ragged, but this is probably because of the relatively small number of samples.  EXERCISE1. Compute HCIP with these distributions. Make a histogram of the result in millions of barrels.1. How does the histogram look if you take 1000 samples instead of 100?1. Make distributions for some of the other properties, like thickness and FVF.1. Maybe our functions should check that we don't get unreasonable values, like negative numbers, or decimal fractions over 1.0 Try to implement this if you have time.<jupyter_code>hcip = calculate_hcip(thick, area, height, top, netg, por, s_o, fvf)
_ = plt.hist(hcip.to('Mbbl'), bins=20)<jupyter_output>/Users/matt/opt/miniconda3/envs/geocomp/lib/python3.8/site-packages/numpy/core/_asarray.py:136: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
  return array(a, dtype, copy=False, order=order, subok=True)<jupyter_text>Using scipy's distributions`scipy` has some distribution functions too. While they're a bit more complicated to use than NumPy's, they are more rigorous. For example, they allow us to fit to data.Following [this post](https://mikulskibartosz.name/monte-carlo-simulation-in-python-d63f0cfcdf6f)<jupyter_code>from scipy.stats import norm

mean = 3
standard_deviation = 2
normal_distribution = norm(loc=mean, scale=standard_deviation)

x = np.linspace(-6, 12, 200)

_, ax = plt.subplots()
ax.plot(x, normal_distribution.pdf(x), '-', lw=2)
plt.title(f'Normal distribution with mean {mean} and standard deviation {standard_deviation}')
plt.show()
number_of_simulations = 20

normal_distribution.rvs(number_of_simulations)<jupyter_output><empty_output><jupyter_text>Let's fit a log-normal distribution to some porosity data.<jupyter_code>import pandas as pd

uid = "1QcSw_xRAYgJzD9HsIXNjmS7o4Zb6qkRBgIWhmp4f2mI"
url = f"https://docs.google.com/spreadsheets/d/{uid}/export?format=csv"
df = pd.read_csv(url)
df.head()
for name, group in df.groupby('Gross environment'):
    plt.scatter(group.Porosity, np.log10(group.Permeability), label=name)
plt.legend()
from scipy.stats import lognorm

deltaic = df.loc[df['Gross environment']=='Deltaic']
s, loc, scale = lognorm.fit(deltaic.Porosity)
s, loc, scale
x = np.linspace(0, 40, 200)

_, ax = plt.subplots()
plt.hist(df.Porosity, bins=80, range=(0, 80))
ax2 = ax.twinx()
ax2.plot(x, lognorm.pdf(x, loc=loc, scale=scale, s=s), 'r')
plt.show()<jupyter_output><empty_output><jupyter_text>Now we can use this in our calculation:<jupyter_code>hcip = calculate_hcip(thick, area, height, top, netg, por, s_o, fvf).magnitude

sns.distplot(hcip)
plt.title('100 realizations')
plt.xlabel('$\mathrm{HCIP\ in\ billion\ m^3}$')
plt.show()
print(f"P10  is {np.percentile(hcip, 10):.2e} Sm³")
print(f"P50  is {np.median(hcip):.2e} Sm³")
print(f"mean is {np.mean(hcip):.2e} Sm³")
print(f"P90  is {np.percentile(hcip, 90):.2e} Sm³")<jupyter_output>P10  is 1.73e+09 Sm³
P50  is 2.45e+09 Sm³
mean is 2.49e+09 Sm³
P90  is 3.30e+09 Sm³<jupyter_text>Compute on a DataFrameSuppose we have a spreadsheet of prospect data:<jupyter_code>import pandas as pd

uid = "1P2JxXG_jLZ0vx8BlFvm0hD6sBBZH2zU8tk9T-SI27mE"
url = f"https://docs.google.com/spreadsheets/d/{uid}/export?format=csv"

df = pd.read_csv(url)
df.head()<jupyter_output><empty_output><jupyter_text>We have to map the row contents to our HCIP function:<jupyter_code>names = {
    'thickness': 'Thick [m]',
    'area': 'Area [km2]',
    'netg': 'N:G',
    'por': 'phi',
    's_o': 'So',
    'fvf': 'Bo'
}

def hcip_row(row):
    params = {k: row[v] for k, v in names.items()}
    params['height'] = 1e6
    params['top'] = 'slab'
    hcip = calculate_hcip(**params) * row['GeomFactor']
    return hcip
df['HCIP'] = df.apply(hcip_row, axis=1)

df.head()<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup