tNavigator Reservoir Simulation Emulator

Build a .DATA file for a 3D black oil model:
- Grid: 20x20x5, dx=dy=200ft, dz varies 10-30ft
- Porosity: 0.15 avg, Perm: 50 md
- Oil API: 35, GOR: 600 scf/STB
- Initial pressure: 4000 psia at 8000ft datum
- One producer, BHP limit 1500 psia

Example .DATA output (partial):

RUNSPEC
TITLE
  Black Oil Model — Educational Example

DIMENS
  20 20 5 /

OIL
WATER
GAS
DISGAS

FIELD

START
  1 JAN 2025 /

TABDIMS
  1 1 40 20 /

WELLDIMS
  2 5 1 2 /

-- ============================================
GRID

DX
  2000*200.0 /
DY
  2000*200.0 /
DZ
  400*10.0 400*15.0 400*20.0 400*25.0 400*30.0 /

TOPS
  400*8000.0 /

PORO
  2000*0.15 /

PERMX
  2000*50.0 /
PERMY
  2000*50.0 /
PERMZ
  2000*5.0 /

Material Balance (Havlena-Odeh)

Generalized form for oil reservoirs with gas cap and water influx:

F = N * [Eo + m*Eg + Efw] + We

Where:

• F = total fluid withdrawal (underground voidage)
• N = OOIP (STB)
• Eo = oil expansion term
• Eg = gas cap expansion term
• Efw = connate water and pore compressibility term
• We = cumulative water influx
• m = ratio of initial gas cap to oil zone volume

See references/equations.md for full definitions.

Workflow:

1. User provides: Np, Gp, Wp (cumulative production), PVT data at multiple pressures
2. Calculate F and Et at each pressure step
3. Plot F vs Et — slope = N (OOIP)
4. Assess drive mechanism from plot shape

Arps Decline Curve Analysis

Three decline models:

Workflow:

1. User provides production history (time, rate pairs)
2. Fit exponential first (log q vs t — look for straight line)
3. If not straight, fit hyperbolic (optimize b and Di)
4. Report: qi, Di, b, EUR, remaining reserves at economic limit

Sensitivity: Pwf vs Recovery Factor

Vary bottom-hole flowing pressure and compute recovery factor using material balance or decline analysis:

1. Define Pwf range (e.g., 500 to 3000 psia in 250 psi steps)
2. At each Pwf, compute IPR rate, then cumulative production to abandonment
3. Report RF = Np/N at each Pwf
4. Output table and recommended operating Pwf

Module 2 — PVT Designer

Capabilities

• Z-factor: Standing-Katz (chart), Hall-Yarborough, Dranchuk-Abou-Kassem (DAK)
• Bubble point pressure: Standing, Vasquez-Beggs
• Oil FVF (Bo): Standing correlation
• Gas FVF (Bg): from Z-factor and P-T
• Water FVF (Bw): McCain correlation
• Solution GOR (Rs): Standing correlation
• Oil viscosity: Beggs-Robinson (dead oil), Vasquez-Beggs (live oil)
• Gas viscosity: Lee-Gonzalez-Eakin
• Equations of state: Peng-Robinson, SRK (for compositional work)

Z-Factor Calculation (DAK Correlation)

The Dranchuk-Abou-Kassem correlation solves iteratively:

Input: Tpr (pseudo-reduced temperature), Ppr (pseudo-reduced pressure)
Output: Z-factor

Pseudo-critical properties from gas gravity (Standing correlations):

Tpc = 168 + 325*gamma_g - 12.5*gamma_g^2  (R)
Ppc = 677 + 15*gamma_g - 37.5*gamma_g^2   (psia)
Tpr = T / Tpc
Ppr = P / Ppc

DAK uses 11 coefficients (A1-A11) — see references/equations.md.

Workflow:

1. User provides: T (F or R), P (psia), gamma_g (gas gravity, air=1)
2. Calculate Tpc, Ppc, Tpr, Ppr
3. Iterate DAK equation for Z
4. Report Z, Bg, gas density

Bubble Point Pressure (Standing)

Pb = 18.2 * ((Rs/gamma_g)^0.83 * 10^(0.00091*T - 0.0125*API) - 1.4)

Workflow:

1. User provides: Rs (scf/STB), gamma_g, T (F), API gravity
2. Calculate Pb
3. Also calculate Bo at Pb using Standing FVF correlation

Oil Formation Volume Factor (Standing)

Below bubble point:

Bo = 0.9759 + 0.00012 * (Rs*(gamma_g/gamma_o)^0.5 + 1.25*T)^1.2

Solution GOR (Standing)

Rs = gamma_g * (P/18.2 * 10^(0.0125*API - 0.00091*T))^1.2048

Oil Viscosity

Dead oil (Beggs-Robinson):

mu_od = 10^(10^(3.0324 - 0.02023*API) * T^(-1.163)) - 1

Live oil (Beggs-Robinson):

a = 10.715 * (Rs + 100)^(-0.515)
b = 5.44 * (Rs + 150)^(-0.338)
mu_o = a * mu_od^b

Gas Viscosity (Lee-Gonzalez-Eakin)

mu_g = K * exp(X * rho_g^Y) * 1e-4   (cp)
K = (9.4 + 0.02*M)*T^1.5 / (209 + 19*M + T)
X = 3.5 + 986/T + 0.01*M
Y = 2.4 - 0.2*X

Where M = gas molecular weight, T in Rankine, rho_g in g/cm3.

Example PVT Output

## PVT Properties at P = 3000 psia, T = 200F

| Property | Value | Units |
|----------|-------|-------|
| Z-factor | 0.842 | — |
| Bubble point (Pb) | 2850 | psia |
| Bo | 1.354 | RB/STB |
| Bg | 0.00523 | RB/scf |
| Rs | 580 | scf/STB |
| Oil viscosity (mu_o) | 0.65 | cp |
| Gas viscosity (mu_g) | 0.018 | cp |
| Oil density | 46.2 | lb/ft3 |

**Note:** Above bubble point — oil is undersaturated. Bo calculated
using above-Pb compressibility correction. Values are from Standing
and Beggs-Robinson correlations (field units).

Module 3 — Automated Assisted History Matching (AHM)

Capabilities

• Latin Hypercube Sampling (LHS) for parameter space exploration
• Particle Swarm Optimization (PSO) for history matching
• Mismatch objective functions (RMS, normalized RMS)
• Simplified proxy model for quick screening

Parameter Space

Typical matching parameters with default ranges:

LHS Workflow

1. User defines: parameters, ranges, number of samples (N)
2. Generate N samples using Latin Hypercube (stratified random)
3. For each sample: run proxy model or simplified material balance
4. Compute mismatch (objective function) against observed data
5. Rank samples by mismatch, report top 10

Objective Function

RMS = sqrt( (1/n) * SUM[(q_sim_i - q_obs_i)^2] )
NRMS = RMS / (q_obs_max - q_obs_min)

Good match: NRMS < 0.10 (10% normalized error).

PSO Workflow

1. Initialize swarm (20-50 particles) from LHS best results
2. Each particle = set of matching parameters
3. Evaluate objective function for each particle
4. Update velocity and position using personal best and global best
5. Iterate 50-200 generations
6. Report: best-fit parameters, convergence history, final mismatch

PSO update equations:

v_i(t+1) = w*v_i(t) + c1*r1*(pbest_i - x_i) + c2*r2*(gbest - x_i)
x_i(t+1) = x_i(t) + v_i(t+1)

Where w=0.7 (inertia), c1=c2=1.5 (cognitive/social), r1,r2 ~ U(0,1).

Example AHM Output

## History Match Results — Well PRD-01

### Best-Fit Parameters
| Parameter | Value | Unit | Range |
|-----------|-------|------|-------|
| kh | 125 | md | 1-500 |
| phi | 0.18 | frac | 0.05-0.30 |
| Skin | 2.5 | — | -3 to 20 |
| Aquifer Jaq | 35 | STB/D/psi | 0-100 |

NRMS = 0.062 (good match)
Iterations: 150, Swarm size: 30

**Note:** History matching is non-unique. Multiple parameter sets may
yield similar mismatch. Always validate against additional data (BHP,
water cut, GOR) before using for forecasting.

Module 4 — Well Designer and VFP Tables

Capabilities

• VFP table generation (tNav VFPPROD keyword format)
• Beggs and Brill multiphase flow correlation
• Tubing pressure drop calculations
• IPR curves (Vogel for two-phase, Fetkovich for gas wells)
• Nodal analysis (IPR-VFP intersection)
• Sensitivity: tubing ID, GLR, water cut

Vogel IPR (Two-Phase Oil Wells)

Below bubble point:

q/qmax = 1 - 0.2*(Pwf/Pr) - 0.8*(Pwf/Pr)^2

Where:

• q = oil rate at Pwf (STB/D)
• qmax = AOF (absolute open flow) at Pwf=0
• Pr = average reservoir pressure (psia)
• Pwf = flowing bottomhole pressure (psia)

Determine qmax from one test point (q_test, Pwf_test):

qmax = q_test / (1 - 0.2*(Pwf_test/Pr) - 0.8*(Pwf_test/Pr)^2)

Fetkovich IPR (Gas Wells)

q = C * (Pr^2 - Pwf^2)^n

Where C = deliverability coefficient, n = exponent (0.5-1.0).

Beggs and Brill Pressure Drop

Simplified approach for tubing pressure gradient:

(dP/dL) = (rho_m * g * sin(theta) + f * rho_m * v_m^2 / (2*d)) / (1 - rho_m*v_m*v_sg/P)

Flow pattern map (Beggs-Brill):

1. Segregated (stratified, annular)
2. Intermittent (plug, slug)
3. Distributed (bubble, mist)

Determine flow pattern from Froude number and liquid content (lambda_L), then compute liquid holdup correction.

Workflow:

1. User provides: tubing ID, depth, surface/BH temperature, fluid properties
2. Divide tubing into segments (50-100 ft each)
3. March from bottom to top (or top to bottom), computing dP per segment
4. Sum pressure drops for total tubing loss
5. Build VFP table: Pwh vs q at various GLR and water cut

VFP Table Format (tNav VFPPROD)

VFPPROD
-- Table  Datum   Rate    WFR     GFR     THP     ALQ
--  No.   Depth   Type    Type    Type    Type    Type
     1    8000.0  OIL     WCT     GOR     THP     ' '  /
-- Rate values (STB/D)
  100 500 1000 2000 3000 5000 /
-- THP values (psia)
  100 200 300 400 500 /
-- WCT values (fraction)
  0.0 0.2 0.5 0.8 /
-- GOR values (scf/STB)
  200 400 600 800 1000 /
-- ALQ values
  0 /
-- BHP data (psia) — one value per (Rate, THP, WCT, GOR, ALQ) combination
  1200 1350 1500 ...
/

Nodal Analysis

Find the operating point where IPR and VFP (tubing performance) intersect:

1. Generate IPR curve: q vs Pwf (Vogel or Fetkovich)
2. Generate VFP curve: q vs Pwf (Beggs-Brill from wellhead to bottomhole)
3. Plot both on same axes (Pwf vs q)
4. Intersection = operating point (q_op, Pwf_op)
5. Sensitivity: overlay curves at different tubing sizes, GLR, Pwh

Example Nodal Analysis Output

## Nodal Analysis — Well PRD-01

### Operating Conditions
| Parameter | Value | Unit |
|-----------|-------|------|
| Reservoir pressure (Pr) | 3500 | psia |
| Wellhead pressure (Pwh) | 200 | psia |
| Tubing ID | 2.992 | inches |
| Depth (TVD) | 8000 | ft |
| GOR | 600 | scf/STB |
| Water cut | 0.20 | fraction |

### Operating Point
| Parameter | Value | Unit |
|-----------|-------|------|
| Oil rate (q_op) | 1850 | STB/D |
| BHP (Pwf_op) | 2200 | psia |
| Drawdown | 1300 | psi |

### Tubing Size Sensitivity
| Tubing ID (in) | q_op (STB/D) | Pwf_op (psia) |
|-----------------|--------------|----------------|
| 2.441 | 1520 | 2450 |
| 2.992 | 1850 | 2200 |
| 3.500 | 2050 | 2050 |

**Note:** Increasing tubing from 2.441" to 2.992" gains 330 STB/D.
Further increase to 3.5" gains only 200 STB/D — diminishing returns.
Validate against actual well test data before sizing.

Module 5 — Geology and Property Modeling (Petrophysics)

Capabilities

• Kriging (simple, ordinary) for spatial interpolation
• Sequential Gaussian Simulation (SGS) concepts
• Synthetic 2D grid from porosity/permeability distributions
• Archie's equation for water saturation
• Permeability-porosity transforms (Kozeny-Carman, empirical)

Archie's Equation

Sw = (a / (phi^m * Rt) * Rw)^(1/n)

Where:

• Sw = water saturation (fraction)
• phi = porosity (fraction)
• Rt = true formation resistivity (ohm-m)
• Rw = formation water resistivity (ohm-m)
• a = tortuosity factor (typ. 0.62-1.0)
• m = cementation exponent (typ. 1.8-2.5)
• n = saturation exponent (typ. 1.8-2.2)

Typical Archie parameters:

Permeability-Porosity Transforms

Kozeny-Carman:

k = (phi^3 * d^2) / (72 * tau^2 * (1-phi)^2)

Where d = grain diameter (cm), tau = tortuosity.

Empirical (semi-log):

log(k) = a + b * phi

Fit a and b from core data (log-linear regression).

Kriging (Simple and Ordinary)

Kriging provides the Best Linear Unbiased Estimate (BLUE) of a spatial variable.

Variogram model (spherical):

gamma(h) = C0 + C * [1.5*(h/a) - 0.5*(h/a)^3]   for h <= a
gamma(h) = C0 + C                                  for h > a

Where C0 = nugget, C = sill - nugget, a = range, h = lag distance.

Simple kriging:

Z*(x0) = mu + SUM[lambda_i * (Z(xi) - mu)]

Weights lambda_i from solving the kriging system using the variogram.

Workflow:

1. User provides: known data points (x, y, value), variogram parameters
2. Build variogram matrix for known points
3. For each unknown grid node, solve kriging system for weights
4. Estimate value and kriging variance at each node
5. Output: interpolated grid + variance grid

SGS (Sequential Gaussian Simulation)

Produces multiple equiprobable realizations:

1. Transform data to normal distribution
2. Define a random path through all grid nodes
3. At each node: krige from nearby data + previously simulated nodes
4. Draw from conditional distribution (kriged mean + kriged variance)
5. Back-transform to original units
6. Repeat for multiple realizations

Note: Full SGS requires significant computation. This skill provides the methodology and can generate simplified 2D examples.

Example Petrophysics Output

## Water Saturation — Archie's Equation

### Input
| Parameter | Value |
|-----------|-------|
| Porosity (phi) | 0.18 |
| Rt (formation resistivity) | 25 ohm-m |
| Rw (water resistivity) | 0.04 ohm-m |
| a (tortuosity) | 0.62 |
| m (cementation) | 2.15 |
| n (saturation) | 2.0 |

### Result
| Property | Value |
|----------|-------|
| Sw | 0.23 (23%) |
| So (1-Sw) | 0.77 (77%) |

Hydrocarbon saturation of 77% indicates a productive zone.
Validate with capillary pressure data and log analysis.

Workflow Summary

Step 1 — Identify Module

Map user request to one of the 5 modules:

• .DATA file, material balance, decline curve -> Module 1
• PVT, Z-factor, Bo, Rs, viscosity, EOS -> Module 2
• History matching, LHS, PSO, parameter fitting -> Module 3
• VFP, IPR, nodal analysis, tubing, wellhead -> Module 4
• Kriging, Archie, petrophysics, grid, porosity-perm -> Module 5

Step 2 — Gather Inputs

Collect required parameters from the user. Provide sensible defaults (with field units) if user omits values. Always state assumptions.

Step 3 — Calculate

Use Python (stdlib only) or analytical formulas. Reference references/python_example.py for implementations. Reference references/equations.md for equation details.

Step 4 — Output

Format as:

1. Input summary table — echo back all parameters and assumptions
2. Results table — computed values with units
3. Narrative — interpretation, trends, recommendations
4. Caveats — correlation validity ranges, assumptions, limitations

Step 5 — .DATA Snippet (if applicable)

If the calculation feeds into a simulation model, provide the relevant Eclipse/tNav keyword block ready for copy-paste.

Output Format

All outputs follow this pattern:

## [Calculation Title]

### Input Parameters
| Parameter | Value | Unit | Source |
|-----------|-------|------|--------|
| ... | ... | ... | User / Default |

### Results
| Property | Value | Unit |
|----------|-------|------|
| ... | ... | ... |

**Summary:** [1-3 sentence interpretation of results]

**Caveats:**
- [Correlation validity range]
- [Assumptions made]
- [When to use full simulation instead]

Error Handling

Units Reference

All correlations use field units by default:

If user provides SI inputs, convert to field units before calculation, then present results in both unit systems.

Implementation Notes

• Use Python (stdlib only) for calculations — math, statistics, random, json, csv
• Reference implementations in references/python_example.py
• Equation derivations and variable definitions in references/equations.md
• For .DATA file generation, follow Eclipse 100/tNavigator keyword syntax
• All iterative solvers (Z-factor, PSO) should have max iteration limits
• Correlation validity ranges should be checked and warned about
• This is an educational tool — always recommend validation against lab data and full simulation for production decisions