Nodal Analysis: Find Your Well’s Operating Point

By Saad Iqbal

A well on a 2-7/8″ string is choked back to 150 psi at the wellhead, and the production supervisor wants to know one thing: will bumping the choke actually move more oil, or is the tubing already the bottleneck? Pull up the decline curve and it won’t tell you. Pull up the reservoir model and it won’t tell you either. The only way to answer it honestly is to put the inflow and the outflow on the same chart and see where they cross. That’s nodal analysis, and it’s the single most useful fifteen minutes of engineering you can spend on a well that’s “producing fine, I guess.”

Nodal analysis sounds like software-vendor jargon, but the idea underneath it is simple: pick a point in the system (usually bottomhole), write an equation for how much the reservoir can deliver to that point, write a second equation for how much pressure the wellbore needs at that point to lift fluid to surface, and find where the two agree. That intersection — not the reservoir’s absolute open flow, not the tubing’s theoretical capacity — is what the well will actually do. This tutorial builds both curves for an undersaturated oil well from first principles, finds the operating point with a few lines of Python, and shows how tubing size shifts it.

Prerequisites: what you need before you start

• Python 3 with NumPy installed (a plain calculator and graph paper work too, just slower)
• Static reservoir pressure Pr, productivity index J, and bubble point pressure Pb from a well test or PVT report
• Tubing ID, perforation (or tubing shoe) depth, and fluid properties (API gravity, viscosity)
• A fixed flowing wellhead pressure (THP) — set by your separator or choke
• The worked example below uses: Pr = 3,500 psia, J = 1.2 STB/d/psi, Pb = 1,800 psia, depth = 8,000 ft, 2-7/8″ tubing (2.441″ ID), 35° API oil, 2 cp viscosity, THP = 150 psia

One scoping decision up front: this example assumes the well stays above its bubble point throughout the tubing string — single-phase liquid flow, no free gas. That’s common for a newer undersaturated reservoir and it lets us use exact single-phase hydraulics instead of an empirical multiphase correlation. If your well is already producing below Pb, the IPR becomes Vogel’s equation instead of a straight line, and the VLP needs a multiphase correlation like Hagedorn-Brown or Beggs-Brill — see our Vogel IPR walkthrough for that inflow case.

Step 1: Build the IPR curve

Above the bubble point, inflow is governed by Darcy’s law and the relationship between rate and flowing bottomhole pressure (Pwf) is a straight line, defined entirely by the productivity index J:

q = J × (Pr − Pwf), or equivalently Pwf = Pr − q/J

With Pr = 3,500 psia and J = 1.2 STB/d/psi, a Pwf of 3,000 psia implies q = 1.2 × (3,500 − 3,000) = 600 STB/d. The straight line is only valid while Pwf stays above the bubble point — once Pwf drops below Pb, gas comes out of solution, relative permeability to oil drops, and the curve bends over (that’s exactly what Vogel’s equation corrects for). Here that sets a hard ceiling: q_max = J × (Pr − Pb) = 1.2 × (3,500 − 1,800) = 2,040 STB/d. Anything the VLP calculation suggests above that rate isn’t physically valid for this IPR.

Step 2: Build the VLP curve

The vertical lift performance (VLP, or outflow) curve answers a different question: for a given surface rate, how much pressure does the tubing need at the bottom to deliver that rate to a fixed wellhead pressure? For single-phase liquid flow, that pressure has exactly two components — the hydrostatic column and friction — and both have exact, non-empirical formulas.

Hydrostatic (elevation) pressure: ΔP_elev = ρ × h / 144, where ρ is fluid density in lb/ft³ and h is vertical depth in feet (144 converts lb/ft² to psi). For 35° API oil, specific gravity γo = 141.5 / (131.5 + API) = 141.5/166.5 = 0.850, so ρ = 0.850 × 62.4 = 53.0 lb/ft³. At 8,000 ft: ΔP_elev = 53.0 × 8,000 / 144 = 2,946 psi — and notice this term doesn’t depend on rate at all.

Friction pressure (Darcy-Weisbach): ΔP_fric = f × (L/D) × (ρ × v²) / (2 × gc) / 144, where f is the Darcy friction factor, L/D is length over diameter, v is fluid velocity, and gc = 32.174 lbm·ft/(lbf·s²). Rather than iterate the Colebrook-White equation for f, the Swamee-Jain approximation gets within about 1% without iteration:

f = 0.25 / (log10(eps/(3.7*D) + 5.74/Re**0.9))**2   # turbulent, Re > 2300
f = 64 / Re                                          # laminar, Re <= 2300

where eps is absolute pipe roughness (0.0006 in for fairly new tubing) and Re = ρvD/μ. Unlike the coiled tubing friction job we walked through for a CT string, there's no need to iterate here — velocity comes straight from the trial rate, so f can be computed directly for each point on the curve.

Running this across a range of rates (full code in Step 3) shows friction is nearly irrelevant at low rates and grows fast at high ones — exactly what the physics predicts, since friction pressure scales with velocity squared:

Rate (STB/d)	Elevation ΔP (psi)	Friction ΔP (psi)	VLP Pwf (psia)
200	2,946	1.6	3,098
600	2,946	10.3	3,106
1,000	2,946	25.1	3,121
1,400	2,946	45.4	3,142
1,800	2,946	70.9	3,167
2,000	2,946	85.6	3,182

VLP Pwf = THP + ΔP_elev + ΔP_fric. Because the hydrostatic term dominates, the VLP curve is almost flat — it rises only about 85 psi from 200 to 2,000 STB/d. That flatness matters: it means this well's performance is controlled almost entirely by the reservoir (the IPR), not the tubing. That's a useful diagnostic on its own, before you even find the intersection.

Step 3: Solve for the operating point

The operating point is the rate where IPR_Pwf(q) = VLP_Pwf(q). Since IPR is a straight line falling with rate and VLP is a curve rising with rate, there's exactly one crossing in the valid range — a textbook case for bisection or scipy.optimize.brentq. Here's the complete, runnable version, including the friction-factor logic from Step 2:

import math

# --- Inputs ---
Pr, Pb, J = 3500.0, 1800.0, 1.20          # psia, psia, STB/d/psi
depth, d_in = 8000.0, 2.441               # ft, in
API, mu_cp, THP = 35.0, 2.0, 150.0        # deg API, cp, psia
eps_in = 0.0006                           # in, pipe roughness

gamma_o = 141.5/(131.5+API)
rho = gamma_o*62.4                        # lb/ft3
mu = mu_cp*6.72e-4                        # lb/(ft.s)
d_ft = d_in/12.0
A = math.pi/4*d_ft**2
eps_ft = eps_in/12.0
gc = 32.174

def pwf_ipr(q):
    return Pr - q/J

def pwf_vlp(q):
    v = (q*5.615/86400.0)/A               # bbl/d -> ft/s
    Re = rho*v*d_ft/mu if v > 0 else 0
    if Re <= 2300:
        f = 64.0/max(Re, 1e-6)
    else:
        rel_rough = eps_ft/d_ft
        f = 0.25/(math.log10(rel_rough/3.7 + 5.74/Re**0.9))**2
    dp_elev = rho*depth/144.0
    dp_fric = f*(depth/d_ft)*(rho*v**2)/(2*gc)/144.0
    return THP + dp_elev + dp_fric

def mismatch(q):
    return pwf_ipr(q) - pwf_vlp(q)

# Bisection between a tiny rate and the bubble-point-limited max rate
q_max = J*(Pr-Pb)
lo, hi = 1.0, q_max - 1
for _ in range(100):
    mid = (lo+hi)/2
    if mismatch(lo)*mismatch(mid) <= 0:
        hi = mid
    else:
        lo = mid

q_op = (lo+hi)/2
print(f"Operating point: q = {q_op:.0f} STB/d, Pwf = {pwf_ipr(q_op):.0f} psia")
# Operating point: q = 476 STB/d, Pwf = 3103 psia

For this well, the answer is 476 STB/d at a flowing bottomhole pressure of about 3,103 psia — well inside the 2,040 STB/d bubble-point ceiling, so the single-phase assumption holds. That's the number the choke-bumping question actually depends on, and it's not something you'd have gotten from the IPR or the VLP in isolation.

Step 4: Run the sensitivity that actually answers the field question

Now go back to the original question: does bumping the choke move more oil? Lowering THP shifts the whole VLP curve down, which — since VLP is so flat here — barely moves the intersection, because the IPR is nearly vertical relative to it near the crossing point. The bigger lever is tubing size. Re-running pwf_vlp with a smaller 1.995" ID (2-3/8" tubing) versus a larger 3.958" ID (3-1/2" tubing), holding everything else fixed, shows how much friction — and therefore the operating point — actually moves:

• 2-3/8" tubing (1.995" ID): friction roughly triples at the same rate versus the base case, nudging the operating point down
• 2-7/8" tubing (2.441" ID, base case): operating point ≈ 476 STB/d
• 3-1/2" tubing (3.958" ID): friction becomes nearly negligible, operating point edges up slightly — the well is reservoir-limited either way

The honest answer for this well: it's reservoir-limited, not tubing-limited, so a bigger choke or bigger tubing buys almost nothing. That's a conclusion worth having before anyone spends money on a workover. For wells with a higher GOR, a longer string, or a smaller completion, the balance tips the other way — which is exactly why you build both curves instead of assuming.

To automate this for a whole portfolio, wrap the solver above in a loop over well files (Python plus pandas handles that cleanly — read Pr, J, Pb, depth and tubing size from a spreadsheet, solve for q_op per well, flag anything within 10% of its bubble-point ceiling), and feed the results into a dashboard with Power BI or Looker Studio so engineers can see which wells are reservoir-limited versus tubing-limited at a glance. If your reservoir has a free gas phase throughout, our integrated asset modelling roundup covers where full multiphase nodal tools like PROSPER and PIPESIM fit in versus a hand-rolled script like this one.

How to verify it worked

Three sanity checks before you trust the number:

1. Zero-rate check: at q = 0, VLP_Pwf should equal THP + the full hydrostatic column (no friction). Here that's 150 + 2,946 = 3,096 psia — confirm your code returns that.
2. Sign check: IPR must be decreasing in q and VLP increasing (or flat); if VLP is decreasing somewhere, a unit got mixed up (bbl/d vs ft³/s is the usual culprit).
3. Validity check: the solved Pwf must stay above Pb (3,103 > 1,800 here ✓). If it doesn't, the single-phase model no longer applies and you need Vogel plus a multiphase VLP correlation instead.

Common pitfalls

• Mixing rate units mid-calculation. bbl/d, ft³/s, and gal/min all show up in petroleum hydraulics — convert once, at the top of the function, and never again mid-formula.
• Forgetting the bubble-point ceiling. A linear IPR extrapolated past Pb will hand you a rate that looks plausible and is physically wrong.
• Using a single average friction factor across the whole rate range. f changes with Reynolds number, which changes with rate — compute it fresh at every trial point, as the code above does.

Nodal analysis is one of those techniques that feels like overkill for a single well and becomes essential the moment you're triaging twenty of them. Once the solver works for one well, it works for the field — and that's where the real time savings show up.