Fractional Flow vs Mobility Ratio

Level 2

Chapter 11: Reservoir Simulation

descriptionProblem

Plot the fractional flow curve for three different oil viscosities: 1 cp (light oil), 5 cp (medium), and 50 cp (heavy oil). How does the shock front saturation change? What does this mean for waterflood recovery efficiency in heavy oil reservoirs?

---

A waterflood on OML 58 sweeps oil toward the producer only as fast as the fractional-flow curve $f_w(S_w)$ allows. That curve (and therefore how much oil you recover before water breaks through) is controlled by the mobility ratio, here parameterised as the oil/water viscosity ratio $M = \mu_o / \mu_w$ . A favourable (low- $M$ ) flood pushes a piston-like front with a high front saturation and clean sweep; an unfavourable (high- $M$ ) flood fingers water through early, leaving oil behind.

You'll quantify this with the embedded corey_relperm and fractional_flow helpers (do not re-derive them). Saturations live on a fixed grid SW = np.linspace(SWI, 1 - SOR, NPTS) with SWI = 0.2, SOR = 0.2. Water viscosity is held at MU_W_CP = 0.5; you vary the oil viscosity so that mu_o_cp = M * MU_W_CP.

Your tasks:

Write fw_curves(M_list) returning a dict mapping each M in M_list

to its fractional-flow array fractional_flow(SW, SWI, SOR, M*MU_W_CP, MU_W_CP) evaluated on the shared SW grid. Each value is a length-NPTS array.

Write welge_front(Sw, fw, Swi) returning the Welge front saturation

Sw_f: the saturation on the grid that maximises the secant slope fw / (Sw - Swi) from the connate-water point. (The first grid point has Sw == Swi, so guard the division and ignore that point; return the Sw value at the arg-max as a float.)

Build front_by_M, a dict mapping each M in M_LIST = [2, 10, 50] to its

front saturation welge_front(SW, fw_curves(...)[M], SWI).

Plot all three f_w(S_w) curves on one axes (one ax.plot line per M),

with an x-label, a y-label, a title, and a legend. Return the Figure as fig. The grader inspects the axes only, never pixels.

> Think about it: as $M$ climbs from 2 to 50 the front saturation $S_{w,f}$ > drops and the fractional flow at the front falls; water reaches the > producer at a lower saturation, so less oil is displaced ahead of it. That is > exactly why an unfavourable mobility ratio (heavy oil, hot water cut) wrecks > waterflood sweep efficiency.

lightbulbHints (0/4)

Stuck? Reveal hints one at a time — they progress from nudge to near-solution.

codeYour solution

main.py

import numpy as np
import matplotlib.pyplot as plt

# ── Embedded relperm + fractional-flow helpers (do not edit) ─────────────
def corey_relperm(Sw, Swi, Sor, kro_max=0.8, krw_max=0.3, no=2.0, nw=3.0):
    Sn = np.clip((Sw - Swi) / (1 - Swi - Sor), 0, 1)
    return kro_max * (1 - Sn) ** no, krw_max * Sn ** nw

def fractional_flow(Sw, Swi, Sor, mu_o, mu_w, **kw):
    kro, krw = corey_relperm(Sw, Swi, Sor, **kw)
    return (krw / mu_w) / (krw / mu_w + kro / mu_o + 1e-30)

# ── OML 58 waterflood: fixed saturation grid + viscosity sweep ───────────
SWI = 0.2                 # connate water saturation (v/v)
SOR = 0.2                 # residual oil saturation (v/v)
NPTS = 501
SW = np.linspace(SWI, 1 - SOR, NPTS)   # shared Sw grid (sw, v/v)
MU_W_CP = 0.5             # water viscosity held fixed (cp)
M_LIST = [2, 10, 50]      # oil/water viscosity ratios M = mu_o / mu_w

def fw_curves(M_list):
    """Dict mapping each M -> fractional_flow on the shared SW grid.

For each M the oil viscosity is mu_o_cp = M * MU_W_CP; water stays MU_W_CP.
    """
    # TODO: return {M: fractional_flow(SW, SWI, SOR, M * MU_W_CP, MU_W_CP) for M in M_list}
    pass

def welge_front(Sw, fw, Swi):
    """Welge front saturation: the Sw that maximises fw / (Sw - Swi).

The first grid point has Sw == Swi (zero denominator) -- guard it and
    ignore that point. Return the Sw value at the arg-max as a float.
    """
    # TODO: slope = fw / (Sw - Swi + 1e-12)
    # TODO: valid = Sw > Swi + 1e-9        # drop the connate-water point
    # TODO: idx = int(np.argmax(np.where(valid, slope, -np.inf)))
    # TODO: return float(Sw[idx])
    pass

# TODO: front_by_M = {M: welge_front(SW, fw_curves([M])[M], SWI) for M in M_LIST}
front_by_M = None

def plot_fw_curves(M_list):
    """One f_w(Sw) line per M, with labels, a title, and a legend; return fig."""
    # TODO: fig, ax = plt.subplots(figsize=(7, 5))
    # TODO: curves = fw_curves(M_list)
    # TODO: for M in M_list: ax.plot(SW, curves[M], label=f"M = {M}")
    # TODO: ax.set_xlabel(...); ax.set_ylabel(...); ax.set_title(...); ax.legend()
    # TODO: return fig
    pass

fig = plot_fw_curves(M_LIST)
print("front saturations Sw_f by M:", front_by_M)
plt.show()

visibilityReveal reference solutionexpand_more

Try solving it yourself first — the hints walk you through it. The solution below is one valid approach; yours may differ and still be correct.

import numpy as np
import matplotlib.pyplot as plt


# ── Embedded relperm + fractional-flow helpers (do not edit) ─────────────
def corey_relperm(Sw, Swi, Sor, kro_max=0.8, krw_max=0.3, no=2.0, nw=3.0):
    Sn = np.clip((Sw - Swi) / (1 - Swi - Sor), 0, 1)
    return kro_max * (1 - Sn) ** no, krw_max * Sn ** nw


def fractional_flow(Sw, Swi, Sor, mu_o, mu_w, **kw):
    kro, krw = corey_relperm(Sw, Swi, Sor, **kw)
    return (krw / mu_w) / (krw / mu_w + kro / mu_o + 1e-30)


# ── OML 58 waterflood: fixed saturation grid + viscosity sweep ───────────
SWI = 0.2                 # connate water saturation (v/v)
SOR = 0.2                 # residual oil saturation (v/v)
NPTS = 501
SW = np.linspace(SWI, 1 - SOR, NPTS)   # shared Sw grid (sw, v/v)
MU_W_CP = 0.5             # water viscosity held fixed (cp)
M_LIST = [2, 10, 50]      # oil/water viscosity ratios M = mu_o / mu_w


def fw_curves(M_list):
    """Dict mapping each M -> fractional_flow on the shared SW grid.

    For each M the oil viscosity is mu_o_cp = M * MU_W_CP; water stays MU_W_CP.
    """
    return {M: fractional_flow(SW, SWI, SOR, M * MU_W_CP, MU_W_CP) for M in M_list}


def welge_front(Sw, fw, Swi):
    """Welge front saturation: the Sw that maximises fw / (Sw - Swi).

    The first grid point has Sw == Swi (zero denominator) -- guard it and
    ignore that point. Return the Sw value at the arg-max as a float.
    """
    Sw = np.asarray(Sw, dtype=float)
    fw = np.asarray(fw, dtype=float)
    slope = fw / (Sw - Swi + 1e-12)
    valid = Sw > Swi + 1e-9                     # drop the connate-water point
    idx = int(np.argmax(np.where(valid, slope, -np.inf)))
    return float(Sw[idx])


front_by_M = {M: welge_front(SW, fw_curves([M])[M], SWI) for M in M_LIST}


def plot_fw_curves(M_list):
    """One f_w(Sw) line per M, with labels, a title, and a legend; return fig."""
    fig, ax = plt.subplots(figsize=(7, 5))
    curves = fw_curves(M_list)
    for M in M_list:
        ax.plot(SW, curves[M], lw=1.8, label=f"M = {M}")
    ax.set_xlabel("Water saturation Sw (v/v)")
    ax.set_ylabel("Fractional flow of water fw")
    ax.set_title("OML 58 fractional flow vs mobility ratio M")
    ax.legend(loc="upper left")
    ax.grid(True, alpha=0.25)
    return fig


fig = plot_fw_curves(M_LIST)
print("front saturations Sw_f by M:", front_by_M)
plt.show()

lockCopying code is a Full Access feature.

arrow_back

11.5 Time-Step Sensitivity of the Implicit Simulator

11.7 Analytical Validation - Simulator vs the Ei Line-Source

arrow_forward