Beggs and Brill Pressure Drop Calculation: Step by Step

By Saad Iqbal

Picture a gas-lifted well teetering on the edge of liquid loading. The vendor nodal-analysis software is down, and the production engineer needs an answer this afternoon: add gas, swap the choke, or pull the string. The 1973 Beggs and Brill correlation is the tool that gets that answer by hand, in any pipe, at any angle, whether the flow is segregated, slugging, or fully mixed. It has survived fifty years of flashier multiphase models because it is simple enough to code in an afternoon and robust enough to trust in a horizontal gathering line or a near-vertical tubing string alike.

This tutorial walks the full calculation end to end — the flow-pattern map, the liquid holdup, the inclination correction, the friction term — with one fully worked numeric example you can check your own spreadsheet against, then turns it into a Python function you never have to run by hand again.

What You Need Before You Start

The correlation works on one pipe segment at a time — a joint of tubing, a flowline run, one increment of a wellbore traverse. You need the segment geometry and the in-situ (flowing) fluid properties at that segment’s average pressure and temperature. Here is the worked example used throughout this tutorial:

Input	Symbol	Value
Pipe inside diameter	d	3.0 in (0.250 ft)
Inclination from horizontal	θ	90° (vertical, uphill)
Segment length	ΔZ	1,000 ft
Superficial liquid velocity	Vsl	3.0 ft/s
Superficial gas velocity	Vsg	6.0 ft/s
Liquid density	ρL	50.0 lb/ft³
Gas density (at 1,500 psia, 140°F, z=0.85, γg=0.75)	ρg	5.96 lb/ft³
Gas–liquid surface tension	σ	20 dyn/cm
Liquid viscosity	μL	5.0 cp
Gas viscosity	μg	0.02 cp

All of this needs to be in field units — ft, ft/s, lb/ft³, cp. The correlation’s coefficients were regressed in these units, and mixing in SI numbers partway through is the single most common way to get a nonsense answer.

Step 1: Convert Rates to Velocities and Find the No-Slip Holdup

Add the superficial velocities to get the mixture velocity, then take the liquid’s share of it — this is the no-slip liquid holdup, λ_L, the input liquid fraction before the correlation accounts for gas slipping past the liquid:

V_m = V_sl + V_sg 
 λ_L = V_sl / V_m 
 N_Fr = V_m² / (g × d), with g = 32.2 ft/s²

Plugging in the example: V_m = 3.0 + 6.0 = 9.0 ft/s; λ_L = 3.0 / 9.0 = 0.3333; N_Fr = 9.0² / (32.2 × 0.250) = 10.062. The Froude number N_Fr is the backbone of the next step — it measures how gravity-dominated versus inertia-dominated the flow is, which is exactly what separates a smoothly distributed flow from a lurching slug.

Step 2: Determine the Flow Pattern

Beggs and Brill map the flow into four regimes — segregated, transition, intermittent, and distributed — using four boundary curves that depend only on λ_L:

• L1 = 316 × λ_L^0.302
• L2 = 0.0009252 × λ_L^−2.4684
• L3 = 0.10 × λ_L^−1.4516
• L4 = 0.5 × λ_L^−6.738

With λ_L = 0.3333: L1 = 226.78, L2 = 0.0139, L3 = 0.4927, L4 ≈ 820.0. The regime follows from where N_Fr falls relative to these lines: segregated flow sits below L1/L2, distributed flow sits above L1 (for λ_L < 0.4) or L4, intermittent flow sits in between, and a narrow transition band sits between L2 and L3. Here, N_Fr = 10.062 falls between L3 (0.4927) and L1 (226.78) with λ_L < 0.4, so this segment is intermittent (slug) flow — exactly the regime you’d expect from a moderate-rate gas-lifted well.

Step 3: Calculate Liquid Holdup and Apply the Inclination Correction

First, the horizontal holdup: E_L(0) = a × λ_L^b / N_Fr^c, never allowed to fall below λ_L itself.

Regime	a	b	c
Segregated	0.980	0.4846	0.0868
Intermittent	0.845	0.5351	0.0173
Distributed	1.065	0.5824	0.0609

For intermittent flow: E_L(0) = 0.845 × 0.3333^0.5351 / 10.062^0.0173 = 0.4510. Next, the liquid velocity number N_Lv = 1.938 × V_sl × (ρ_L/σ)^0.25 = 7.311, which feeds the inclination factor C for uphill intermittent flow (coefficients d=2.96, e=0.305, f=−0.4473, g=0.0978):

C = (1−λ_L) × ln(2.96 × λ_L^0.305 × N_Lv^−0.4473 × N_Fr^0.0978) = 0.0574 (floored at 0 if negative — it isn’t here)

The angle correction is ψ = 1 + C × [sin(1.8θ) − ⅓sin³(1.8θ)]. At θ=90°, sin(162°)=0.309, so ψ = 1 + 0.0574 × (0.309−0.0098) = 1.0172. The true, slip-corrected holdup is then E_L(θ) = ψ × E_L(0) = 0.4588 — noticeably higher than the no-slip value of 0.3333, because gas slips upward past the slower liquid and leaves more liquid behind per foot of pipe. That gap between λ_L and E_L(θ) is the entire reason this correlation exists.

Step 4: Compute the Hydrostatic (Elevation) Pressure Component

The slip-corrected mixture density carries the weight of the column: ρ_m = ρ_L × E_L(θ) + ρ_g × (1−E_L(θ)) = 50.0 × 0.4588 + 5.96 × 0.5412 = 26.16 lb/ft³. The hydrostatic pressure drop over the segment is ΔP_el = ρ_m × ΔZ / 144 (for a deviated well, replace ΔZ with the true vertical rise, L × sinθ):

ΔP_el = 26.16 × 1,000 / 144 = 181.7 psi over the 1,000 ft segment.

Step 5: Compute the Frictional Pressure Component

Friction needs a two-phase friction factor, built by scaling a no-slip friction factor by a ratio that depends on y = λ_L / E_L(θ)² = 0.3333/0.4588² = 1.584. Since y > 1.2:

S = ln(y) / [−0.0523 + 3.182·ln(y) − 0.8725·ln(y)² + 0.01853·ln(y)⁴] = 0.4598 / 1.2273 = 0.3747, so f_tp/f_ns = e^S = 1.4545.

The no-slip density and viscosity use λ_L, not E_L(θ): ρ_ns = 50.0×0.3333 + 5.96×0.6667 = 20.64 lb/ft³; μ_ns = 5.0×0.3333 + 0.02×0.6667 = 1.68 cp. That gives a no-slip Reynolds number N_Re,ns = 1,488 × ρ_ns × V_m × d / μ_ns = 41,133 — comfortably turbulent. Using the Drew–Koo–McAdams smooth-pipe correlation that the original Beggs and Brill paper used, f_ns = 0.0056 + 0.5/N_Re,ns^0.32 = 0.02229 (Fanning basis), so f_tp = 0.02229 × 1.4545 = 0.03242.

The frictional pressure drop (Fanning form) is ΔP_f = 2×f_tp×ρ_ns×V_m²×ΔZ / (g_c×d×144) = 93.5 psi over the segment.

Step 6: Sum the Total and Automate It in Python

Total pressure drop is simply the two components added together: ΔP = 181.7 + 93.5 = 275.2 psi over 1,000 ft, a gradient of 0.275 psi/ft. Chain enough of these segments up the tubing — recomputing PVT properties at each new average pressure — and you have a full pressure traverse, which is exactly the VLP curve half of the nodal analysis plot that tells you where a well will actually operate.

Here is the whole calculation as a single Python function, using only the standard library for one segment (swap in Python with pandas and NumPy to vectorize it across hundreds of depth steps and dozens of wells at once):

import math

def beggs_brill_dp(d_ft, theta_deg, dZ_ft, Vsl, Vsg, rhoL, rhog, sigma, muL, mug):
    g, gc = 32.2, 32.2
    Vm = Vsl + Vsg
    lamL = Vsl / Vm
    NFr = Vm**2 / (g * d_ft)

    L1 = 316 * lamL**0.302
    L2 = 0.0009252 * lamL**-2.4684
    L3 = 0.10 * lamL**-1.4516
    L4 = 0.5 * lamL**-6.738

    if (lamL < 0.01 and NFr < L1) or (lamL >= 0.01 and NFr < L2):
        pattern = "segregated"
    elif lamL >= 0.01 and L2 <= NFr <= L3:
        pattern = "transition"
    elif (lamL < 0.4 and L3 < NFr <= L1) or (lamL >= 0.4 and L3 < NFr <= L4):
        pattern = "intermittent"
    else:
        pattern = "distributed"

    coeff = {"segregated": (0.98, 0.4846, 0.0868),
             "intermittent": (0.845, 0.5351, 0.0173),
             "distributed": (1.065, 0.5824, 0.0609)}
    incl = {"segregated": (0.011, -3.768, 3.539, -1.614),
            "intermittent": (2.96, 0.305, -0.4473, 0.0978),
            "distributed": None}

    def EL0_of(p):
        a, b, c = coeff[p]
        return max(a * lamL**b / NFr**c, lamL)

    def C_of(p):
        if incl[p] is None:
            return 0.0
        dc, e, f, gg = incl[p]
        return max((1 - lamL) * math.log(dc * lamL**e * NFr**gg * NLv**f), 0.0)

    NLv = 1.938 * Vsl * (rhoL / sigma)**0.25

    if pattern == "transition":
        A = (L3 - NFr) / (L3 - L2)
        EL0 = A * EL0_of("segregated") + (1 - A) * EL0_of("intermittent")
        C = A * C_of("segregated") + (1 - A) * C_of("intermittent")
    else:
        EL0 = EL0_of(pattern)
        C = C_of(pattern)

    theta = math.radians(theta_deg)
    psi = 1 + C * (math.sin(1.8 * theta) - (1/3) * math.sin(1.8 * theta)**3)
    ELtheta = psi * EL0

    rho_s = rhoL * ELtheta + rhog * (1 - ELtheta)
    dP_elev = rho_s * dZ_ft / 144.0

    y = lamL / ELtheta**2
    if 1 < y < 1.2:
        S = math.log(2.2 * y - 1.2)
    else:
        lny = math.log(y)
        S = lny / (-0.0523 + 3.182*lny - 0.8725*lny**2 + 0.01853*lny**4)

    rho_ns = rhoL * lamL + rhog * (1 - lamL)
    mu_ns = muL * lamL + mug * (1 - lamL)
    NRe_ns = 1488 * rho_ns * Vm * d_ft / mu_ns
    f_ns = 0.0056 + 0.5 / NRe_ns**0.32
    f_tp = f_ns * math.exp(S)

    dP_fric = 2 * f_tp * rho_ns * Vm**2 * dZ_ft / (gc * d_ft * 144.0)

    return {"pattern": pattern, "EL_theta": ELtheta,
            "dP_elevation_psi": dP_elev, "dP_friction_psi": dP_fric,
            "dP_total_psi": dP_elev + dP_fric}

result = beggs_brill_dp(d_ft=0.250, theta_deg=90, dZ_ft=1000,
                         Vsl=3.0, Vsg=6.0, rhoL=50.0, rhog=5.96,
                         sigma=20.0, muL=5.0, mug=0.02)
print(result)
# {'pattern': 'intermittent', 'EL_theta': 0.4588, 'dP_elevation_psi': 181.7,
#  'dP_friction_psi': 93.5, 'dP_total_psi': 275.2}

Running it reproduces the hand calculation exactly. For a daily job across your whole well list, drop the per-well inputs into a pandas DataFrame, apply the function row by row (or vectorize the math with NumPy), and push the output pressure traverses straight into Power BI or Looker Studio. If you'd rather not run Python at all, n8n or Power Automate can call the same logic as a scripted step inside a scheduled workflow that pulls rates from your production database and writes the gradient back.

How to Verify the Result and Common Pitfalls

Expected result: 275.2 psi over the 1,000 ft segment, split 181.7 psi hydrostatic and 93.5 psi frictional. Sanity check it two ways: first, the no-slip-only hydrostatic gradient (using ρ_ns = 20.64 lb/ft³ instead of the slip-corrected 26.16 lb/ft³) would give only 143.3 psi — since gas always slips past liquid in upward flow, the true holdup E_L(θ) must come out above the no-slip λ_L, and it does (0.4588 > 0.3333). If your code ever returns E_L(θ) below λ_L for uphill flow, you have a sign or floor-check bug. Second, where you have vendor software or a second correlation like Hagedorn–Brown available, cross-check: the two should typically agree within 10–15% in bubble or slug flow.

• Unit mixing — the L1–L4 and E_L(0) coefficients are fit to field units (ft, ft/s, lb/ft³); feeding in meters or kg/m³ partway through silently produces a wrong flow pattern.
• Skipping the E_L(0) ≥ λ_L floor — easy to miss right at the segregated/intermittent boundary, where the raw power-law formula can dip below the physically required minimum.
• Using uphill coefficients on a downhill run — a return bend or undulating flowline needs the separate downhill coefficient set (d=4.70, e=−0.3692, f=0.1244, g=−0.5056), not just a negated angle.
• Confusing λ_L and E_L(θ) — the no-slip input fraction and the true, slip-corrected holdup are never interchangeable; the gap between them is the entire physical point of the correlation.

Once the wellbore traverse is solid, the natural next stop is checking what happens at the surface: see our walkthrough of choke size calculation with the Gilbert equation for how the flowing wellhead pressure computed here sets the choke's performance, or the Turner equation tutorial for checking whether this same well is at risk of liquid loading.