<jupyter_start><jupyter_text>**Unit 9 Dry-Gas Reservoirs**<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
!git clone https://www.github.com/yohanesnuwara/reservoir-engineering<jupyter_output>Cloning into 'reservoir-engineering'...
warning: redirecting to https://github.com/yohanesnuwara/reservoir-engineering.git/
remote: Enumerating objects: 18, done.[K
remote: Counting objects: 100% (18/18), done.[K
remote: Compressing objects: 100% (18/18), done.[K
remote: Total 927 (delta 9), reused 0 (delta 0), pack-reused 909
Receiving objects: 100% (927/927), 12.07 MiB | 39.76 MiB/s, done.
Resolving deltas: 100% (404/404), done.<jupyter_text>**Functions (Chapter 9.2 Volumetrics and Recovery Factors)**<jupyter_code>def recovery_factor_drygas_waterdrive(swi, Bgi, sgr, Bga):
  "recovery factor of dry gas reservoir with waterdrive"
  # gas FVF at abandonment pressure (RB/scf)
  # sgr: residual gas saturation
  rf = (((1 - swi) / Bgi) - (sgr / Bga)) / ((1 - swi) / Bgi)
  return(rf)

def recovery_factor_drygas_volumetric(Bgi, Bga):
  "recovery factor of dry gas reservoir with NO waterdrive"  
  # gas FVF at abandonment pressure (RB/scf)
  # sgr: residual gas saturation
  rf = 1 - (Bgi / Bga)
  return(rf)<jupyter_output><empty_output><jupyter_text>**Chapter 9.3 Gas Equivalence to Condensates and Water**It's already mentioned in book that **1 STB of condensate is equivalent with 700 to 750 scf of gas** measured at standard conditions (`T = 519.67 R; 15.6 C; 60 F` and `P = 14.7 psi; 101.325 kPa`). I'm interested to find out the **API gravity** of condensates, by re-arranging Equation 9.5, 9.6, and 9.7)<jupyter_code>"Gas equivalence to condensate"

Tsc = 519.57 # standard temp in Rankine
psc = 14.7 # standard pressure in psi
R = 10.732
rho_w = 350.16 # water density in lbm/STB
V1 = 700 # the smallest value of gas equivalence (scf)
V2 = 750 # the largest value of gas equivalence (scf)

# solve for api from Equation 9.5, 9.6, 9.7
from scipy.optimize import fsolve

def f(api):
  so = 141.5 / (api + 131.5) # so: specific gravity of oil (dimensionless)
  Mo = 5854 / (api - 8.811) # molecular weight of oil
  n = (rho_w * so) / Mo
  # real gas equation: V = nRT/p
  f1 = ((n * R * Tsc) / psc) - V1
  return(f1)

api1 = fsolve(f, 40)
api1 = np.float64(api1)

def f(api):
  so = 141.5 / (api + 131.5) # so: specific gravity of oil (dimensionless)
  Mo = 5854 / (api - 8.811) # molecular weight of oil
  n = (rho_w * so) / Mo
  # real gas equation: V = nRT/p
  f2 = ((n * R * Tsc) / psc) - V2
  return(f2)

api2 = fsolve(f, 40)
api2 = np.float64(api2)

print("API gravity of condensate is between:", api1, "API and:", api2, "API")<jupyter_output>API gravity of condensate is between: 47.93319100231802 API and: 51.57941144067795 API<jupyter_text>In most references, the value of API for condesates is greater than this value range (50 to 80 API). For example: **100 STB of condensate with API gravity of 65**, is therefore equivalent to:<jupyter_code>def condensate_to_gas_equivalence(api, stb):
  "Derivation from real gas equation"
  Tsc = 519.57 # standard temp in Rankine
  psc = 14.7 # standard pressure in psi
  R = 10.732
  rho_w = 350.16 # water density in lbm/STB

  so = 141.5 / (api + 131.5) # so: specific gravity of oil (dimensionless)
  Mo = 5854 / (api - 8.811) # molecular weight of oil
  n = (rho_w * so) / Mo  

  V1stb = ((n * R * Tsc) / psc)
  V = V1stb * stb
  return(V1stb, V)

stb = 100; api = 65

V1stb, V = condensate_to_gas_equivalence(api, stb)
print("1 STB of condensate is equivalent to:", V1stb, "scf of gas.", stb, "STB of condensate is equivalent to:", V, "scf of gas")<jupyter_output>1 STB of condensate is equivalent to: 918.0501371833875 scf of gas. 100 STB of condensate is equivalent to: 91805.01371833876 scf of gas<jupyter_text>The gas equivalence is also simplified into a more general form (Eq 9.5)$$V = 132849\frac{\gamma}{M}$$$M$ is molecular weight * for condensate/oil, $M_o = \frac{42.43\gamma_o}{1.008-\gamma_o}$ with $\gamma_o$ is oil specific gravity; $\gamma_o = \frac{141.5}{api + 131.5}$ and $api$ is the API gravity* for water, $M_w = 18$ and $\gamma_w = 1$ **100 STB of water**, is therefore equivalent to:<jupyter_code>stb = 100
V1stb = 132849 * (1 / 18)
V = stb * V1stb
print(stb, "STB of water is equivalent to:", V, "scf of gas")<jupyter_output>100 STB of water is equivalent to: 738050.0 scf of gas<jupyter_text>**Chapter 9.4 Material-Balance Analysis of Dry-Gas Reservoirs**General equation of material balance in dry gas reservoirs:$$F = Gfgi \cdot Eg + \Delta W + Gfgi \cdot Bgi \cdot Efw$$$Gfgi \cdot Bgi \cdot Efw = 0$ because gas normally has compressibility higher than water and rock. In **no waterdrive**, aquifer influx $\Delta W = 0$. The material balance becomes:<jupyter_code>def matbalance_drygas_volumetric(Gp, Bg, Bgi, Gfgi):
  # no waterdrive
  right_term = Gp * Bg
  left_term = Gfgi * (Bg - Bgi)
  return(right_term, left_term) # right term = left term (balance)<jupyter_output><empty_output><jupyter_text>$B_g = \frac{zTp_s}{pT_s}$ from Chapter 2 is subtituted into the material balance, it becomes Eq 9.15: $$\frac{p}{z} = \frac{pi}{zi} - \frac{pi \cdot Gp}{zi \cdot Gfgi}$$ Re-arranging the equation gives **estimation of OGIP**:<jupyter_code>def ogip_matbalance_drygas_volumetric(Gp, p, z, pi, zi):
  # no waterdrive
  Gfgi = Gp / (1 - ((p * zi) / (pi * z)))
  return(Gfgi)<jupyter_output><empty_output><jupyter_text>Example 9.1 Determination of OGIP by Material-Balance<jupyter_code>"Task. determine OGIP"

# known: volumetric (no waterdrive) dry gas reservoir
pi = 3000 # initial pressure, psi
temp = 190 # formation temperature, F
p_current = 2876 # current pressure after production, psia
Gp = 384 # gas production, MMscf
pc = 671 # critical pressure, psia
Tc = 358 # critical temperature, Rankine

# conversion
Gp = Gp * 1E+06 # from MMscf to scf

# determine gas compressibility factor z using Dranchuk-Abou Kassem EOS rather than using chart
import os, sys
sys.path.append('/content/reservoir-engineering/Unit 2 Review of Rock and Fluid Properties/functions')
from dranchuk_aboukassem import dranchuk

# convert temperature from F to Rankine
temp_conv = temp + 459.67

P_pr_initial = pi / pc
P_pr_current = p_current / pc
T_pr = temp_conv / Tc

rho_pr_initial, z_initial = dranchuk(T_pr, P_pr_initial)
rho_pr_current, z_current = dranchuk(T_pr, P_pr_current)
print("Initial gas compressibility factor:", z_initial) 
print("Current gas compressibility factor:", z_current) 

# calculate OGIP
Gfgi = ogip_matbalance_drygas_volumetric(Gp, p_current, z_current, pi, z_initial)
print("Original Gas in Place:", Gfgi / 1E+09, "Bscf")<jupyter_output>Initial gas compressibility factor: 0.9069024150713222
Current gas compressibility factor: 0.9039322509456235
Original Gas in Place: 10.056746912653262 Bscf<jupyter_text>**Chapter 9.5 Material-Balance Plots to Determine OGIP** Plot 9.1: $G_p B_g$ vs $B_g - B_gi$* **Straight line** indicates gas is volumetric (waterdrive is negligible). Intercept $b = 0$ and slope $m = Gfgi$* **Curving upward** indicates waterdrive is present ($\Delta W$ is not negligible, See Chapter 9.6)* **Curving downward** indicates that $Efg$ is not negligible and formation compressibility $ct$ is comparable with gas compressibility (See Chapter 9.7). Plot 9.2: $\frac{p}{z}$ vs $G_p$* **Straight line** indicates gas is volumetric (waterdrive is negligible). Intercept $b = Gfgi$ and slope $m = \frac{p_i}{z_i \cdot Gfgi}$* **Curving upward** indicates waterdrive is present ($\Delta W$ is not negligible, See Chapter 9.6)* **Curving downward** indicates that $Efg$ is not negligible and formation compressibility $ct$ is comparable with gas compressibility (See Chapter 9.7).  Example 9.2 Determination of OGIP by Material-Balance Plot<jupyter_code>"Task. determine OGIP"

# known Table 9.1 
pressure = np.array([3000, 2876, 2824, 2755, 2688, 2570, 2435, 2226, 2122, 1866])
Gp = np.array([0, 384, 550, 788, 1002, 1445, 1899, 2670, 3113, 3982])

# properties similar to Example 9.1
temp = 190 + 459.67 # F to Rankine
Tc = 358
pc = 671

# calculate gas z factor and FVF
z_arr = []
Bg_arr = []
for i in range(len(pressure)):
  P_pr = pressure[i] / pc
  T_pr = temp / Tc
  rho_pr, z = dranchuk(T_pr, P_pr)
  Bg = 0.0282793 * z * temp / pressure[i] # Eq 2.2, temp in Rankine, p in psia
  z_arr.append(float(z))
  Bg_arr.append(float(Bg))

"Using Plot 9.1"
Eg = np.array(Bg_arr) - Bg_arr[0]
Gp_Bg = Gp * np.array(Bg_arr)

# regression
sys.path.append('/content/reservoir-engineering/Unit 7 Introduction to Well-Test Analysis/functions')
from regression import regression
c1, m1 = regression(Eg, Gp_Bg)

# straight line has slope m
Gfgi_plot1 = np.abs(m1)
print("OGIP calculated using Plot 9.1:", Gfgi_plot1 / 1E+03, "Bscf")

"Using Plot 9.2"
p_z = pressure / z_arr

# regression
sys.path.append('/content/reservoir-engineering/Unit 7 Introduction to Well-Test Analysis/functions')
from regression import regression
c2, m2 = regression(Gp, p_z)

# straight line has equation of: y=mx+c. Gfgi is the line intercept when y=0
Gfgi_plot2 = c2 / np.abs(m2)
print("OGIP calculated using Plot 9.2:", Gfgi_plot2 / 1E+03, "Bscf")

# Material-Balance Plots
plt.figure(figsize=(15, 5))

plt.subplot(1, 2, 1)
plt.plot(Eg, Gp_Bg, '.')
# plot regression
F_reg = m1 * Eg + c1
plt.plot(Eg, F_reg, '--')
plt.title('$F$ vs $Eg$ Plot', size=20)
plt.ylabel('$G_p B_g$ (million res cubic ft)'); plt.xlabel('$B_g-B_{gi}$ (res cubic ft/scf)')
plt.xlim(xmin=0); plt.ylim(ymin=0)

plt.subplot(1, 2, 2)
plt.plot(Gp, p_z, '.')
# plot regression
p_z_reg = m2 * Gp + c2
plt.plot(Gp, p_z_reg, '--')
plt.title('$p/z$ Plot', size=20)
plt.xlabel('$G_p$ (MMscf)'); plt.ylabel('$p/z$ (psia)')
plt.xlim(xmin=0)<jupyter_output>OGIP calculated using Plot 9.1: 10.748033224262944 Bscf
OGIP calculated using Plot 9.2: 10.758397571436433 Bscf<jupyter_text>**Chapter 9.6 Waterdrive Gas Reservoirs** **Plot 9.3**: $\frac{F}{E_g + Bgi \cdot Efw}$ vs. $G_p$ or $t$ (Waterdrive Diagnostic Plot)* If the line is **horizontal**, indicates **volumetric gas reservoir** (no waterdrive), equal to $Gfgi$.* If the line is **sloping**, the gas reservoir has **waterdrive** **Plot 9.4**: $\frac{F}{E_g}$ vs. $\sum (\Delta p \cdot We_d)$ (Havlena and Odeh)$We_d$ is calculated using **aquifer infux** in Unit 8. The slope is constant, equals to **aquifer constant, $k'$**. Intercept equals to $Gfgi$ **Plot 9.5**: $F$ vs. $E_c$ (Sills)Where $E_c = E_g + Bgi \cdot Efw + \frac{2(c_f + c_w)Bgi \sum (\Delta p \cdot We_d)}{1-Sw_i} \cdot (\frac{aquiferthickness}{reservoirthickness})$$m = Gfgi$ and $intercept = 0$  Example 9.3 Verification of Volumetric Nature of a Gas Reservoir<jupyter_code>"Task. use Plot 9.3 to verify that reservoir in Example 9.2 is volumetric (no waterdrive)"

axisx = (np.array(Gp))[1:]
F = Gp_Bg
axisy = (np.array(F / Eg))[1:]

# regression
sys.path.append('/content/reservoir-engineering/Unit 7 Introduction to Well-Test Analysis/functions')
from regression import regression
c3, m3 = regression(axisx, axisy)
y_reg = m3 * axisx + c3

# waterdrive diagnostic plot
plt.plot(axisx, axisy, '.')
plt.plot(axisx, y_reg, '--')
plt.title('$F/E_g$ vs $G_p$ Plot', size=20)
plt.xlabel('$G_p$ (MMscf)'); plt.ylabel('$F/E_g$ (MMscf)')
plt.xlim(xmin=0); plt.ylim(8000, 11000)

print("Slope:", m3, "is relatively horizontal")
print("Intercept is OGIP:", c3 / 1E+03, "Bscf")<jupyter_output>/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in true_divide
  """<jupyter_text>Therefore, the reservoir is most likely to be **volumetric (NO waterdrive)**  Example 9.4 Analysis of a Waterdrive Gas Reservoir<jupyter_code># load data
dataprod = pd.read_csv('/content/reservoir-engineering/Unit 9 Dry-Gas Reservoirs/data/Table 9.4-Reservoir Data and Documented Pressures and Production Example 9.4.csv')

# visualize data Gp, Np, Wp in one plot
time = np.array(dataprod.t); Gp = np.array(dataprod.Gp); Np = np.array(dataprod.Np); Wp = np.array(dataprod.Wp)

fig = plt.figure()
host = fig.add_subplot(111)

par1 = host.twinx()
par2 = host.twinx()

host.set_xlabel("Time (years)")
host.set_ylabel("Water Production (STB)")
par1.set_ylabel("Gas Production (STB)")
par2.set_ylabel("Oil Production (STB)")

color1 = plt.cm.viridis(0)
color2 = plt.cm.viridis(0.5)
color3 = plt.cm.viridis(.9)

p1, = host.plot(time, Wp, color=color1,label="Water production")
p2, = par1.plot(time, Gp, color=color2, label="Gas production")
p3, = par2.plot(time, Np, color=color3, label="Oil production")

lns = [p1, p2, p3]
host.legend(handles=lns, loc='best')

# right, left, top, bottom
par2.spines['right'].set_position(('outward', 60))      
pressure = np.array(dataprod.p); z = np.array(dataprod.z)
p_z = pressure / z

# conversion
Gp = np.array(dataprod.Gp) * 1E+06 # mmscf to scf
Bg = np.array(dataprod.Bg) * (1 / (1E+03)) # res ft3/mscf to resft3/scf
F = Gp * Bg # result in res ft3
Eg = Bg - Bg[0] # resft3/scf
F_Eg = F / Eg # res ft3 * (scf/res ft3), result in scf
F_Eg = F_Eg / (1E+09) # convert scf to Bscf

plt.figure(figsize=(13,5))

plt.subplot(1,2,1)
plt.plot((Gp / 1E+06), p_z, '.') # plot x in mmscf
plt.title('p/z vs $G_p$ plot', size=15, pad=10)
plt.xlabel('$G_p$ (MMscf)'); plt.ylabel('$p/z$ (psia)')

plt.subplot(1,2,2)
plt.plot((Gp / 1E+06), F_Eg, '.')
plt.title('$F/E_g$ vs $G_p$ plot', size=15, pad=10)
plt.xlabel('$G_p$ (MMscf)'); plt.ylabel('$F/E_g$ (Bscf)')

plt.tight_layout(pad=3.0) # set distance between two plots
# regression to estimate OGIP
c4, m4 = regression(Gp, p_z)
Gfgi_pz = c4 / np.abs(m4)
print("Original Gas in Place from p/z plot:", Gfgi_pz / 1E+12, "Tscf") # convert scf to Tscf

# the first row contains NaN, so filter it out
Gp_new = Gp[1:]; F_Eg_new = F_Eg[1:]
c5, m5 = regression(Gp_new, F_Eg_new)

y_reg4 = c4 + m4 * Gp
y_reg5 = c5 + m5 * Gp

plt.figure(figsize=(13,5))

plt.subplot(1,2,1)
plt.plot((Gp / 1E+06), p_z, '.')
plt.plot((Gp / 1E+06), y_reg4)
plt.title('p/z vs $G_p$ plot', size=15, pad=10)
plt.xlabel('$G_p$ (MMscf)'); plt.ylabel('$p/z$ (psia)')

plt.subplot(1,2,2)
plt.plot((Gp / 1E+06), F_Eg, '.')
plt.plot((Gp / 1E+06), y_reg5)
plt.title('$F/E_g$ vs $G_p$ plot', size=15, pad=10)
plt.xlabel('$G_p$ (MMscf)'); plt.ylabel('$F/E_g$ (Bscf)')

plt.tight_layout(pad=3.0) # set distance between two plots<jupyter_output>Original Gas in Place from p/z plot: 1.3384478020004833 Tscf<jupyter_text>From p/z plot, the OGIP estimated is 1.33 Tscf. However, the **waterdrive diagnostic plot** of $F/E_g$ vs $G_p$ indicates a sloping upward, meaning there's waterdrive in the reservoir, therefore the OGIP estimated from p/z plot is **misleading**. OGIP calculated by Sills is 231 Bscf.After knowing that the reservoir has waterdrive, calculate the aquifer influx ($W_e$) (Eq 9.18)Check units first:$$W_e=B_wW_p+F-G_{fgi}E_g-G_{fgi}B_{gi}E_{fw}$$$$W_e=(\frac{RB}{STB} STB) + (res ft^3) - (Bscf \cdot \frac{res ft^3}{scf}) - (Bscf \cdot \frac{res ft^3}{scf} \cdot)$$$G_{fgi}$ should be converted to $scf$<jupyter_code>Gfgi = 231 * 1E+09 # 231 from Sills, bscf to scf<jupyter_output><empty_output><jupyter_text>In Table 9.5, $F$ is in unit million res bbl. It's wrong, **should be million res $ft^3$** instead of million  res bbl.<jupyter_code>F_conv = F / 1E+06 # from ft3 to million ft3<jupyter_output><empty_output><jupyter_text>If we stick to the book, we will use this number instead.<jupyter_code># known
cf = 3.5E-6 # psi^-1
cw = 2.9E-6 # psi^-1
swi = 0.21

Bw = np.array(dataprod.Bw)
pi = pressure[0]
Efw = ((cf + cw * swi) / (1 - swi)) * (pi - pressure)

Bgi = Bg[0]
Bw_Wp = Bw * Wp
Gfgi_Eg = Gfgi * Eg
Gfgi_Bgi_Efw = Gfgi * Bgi * Efw
We = (Bw * Wp) + F - (Gfgi * Eg) - (Gfgi * Bgi * Efw)
We
result = pd.DataFrame({"Years": time, "Pressure": pressure, "F(million bbl)": F_conv, "Eg(RB/Mscf)": Eg * 1E+03, "F/Eg(Bscf)": F_Eg, "z": z, "We(res bbl)": We, "p/z(psia)": p_z})
print("Table 9.5-Calculated Parameters for Example 9.4")
result
# visualize We
plt.plot(time, We, '.-')
plt.ylim(ymin=0)
plt.xlabel("Time (years)"); plt.ylabel("Aquifer influx $W_e$ (res bbl)")
plt.title("Aquifer Influx over Production Years")<jupyter_output><empty_output><jupyter_text>***<jupyter_code>"Plot more than one graphs in one plot matplotlib"

import matplotlib.pyplot as plt 

x = [1, 2, 3, 4, 5, 6, 7]
y1 = [20, 34, 56, 10, 70, 80, 90]
y2 = [2, 3, 6, 10, 1, 3, 6]
y3 = [100, 200, 300, 200, 100, 90, 80]

fig = plt.figure()
host = fig.add_subplot(111)

par1 = host.twinx()
par2 = host.twinx()

# host.set_xlim(0, 2)
# host.set_ylim(0, 2)
# par1.set_ylim(0, 4)
# par2.set_ylim(1, 65)

host.set_xlabel("Distance")
host.set_ylabel("Density")
par1.set_ylabel("Temperature")
par2.set_ylabel("Velocity")

color1 = plt.cm.viridis(0)
color2 = plt.cm.viridis(0.5)
color3 = plt.cm.viridis(.9)

# p1, = host.plot([0, 1, 2], [0, 1, 2], color=color1,label="Density")
# p2, = par1.plot([0, 1, 2], [0, 3, 2], color=color2, label="Temperature")
# p3, = par2.plot([0, 1, 2], [50, 30, 15], color=color3, label="Velocity")

p1, = host.plot(x, y1, color=color1,label="Density")
p2, = par1.plot(x, y2, color=color2, label="Temperature")
p3, = par2.plot(x, y3, color=color3, label="Velocity")

lns = [p1, p2, p3]
host.legend(handles=lns, loc='best')

# right, left, top, bottom
par2.spines['right'].set_position(('outward', 60))      
# no x-ticks                 
par2.xaxis.set_ticks([])
# Sometimes handy, same for xaxis
par2.yaxis.set_ticks_position('right')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15077681591469358, 'fraction_numerical': 0.0583041907928682, 'mean_word_length': 4.594611930724824, 'pattern_counts': {'":': 8, '<': 53, '<?xml version=': 0, '>': 53, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '354987', 'license': 'MIT', 'n_tokens_mistral': 7763, 'n_tokens_neox': 6952, 'n_words': 2007, 'path': 'Unit 9 Dry-Gas Reservoirs/notebook/9_examples.ipynb', 'repo_name': 'vijaykumar1588/reservoir-engineering', 'size': '17443'}