examples.aquifer

  1"""An aquifer: water beyond the boundary, feeding the reservoir at its own pressure.
  2
  3Beyond some stretch of the reservoir's boundary there may lie water-bearing
  4rock -- an *aquifer* -- at a pressure of its own, which then feeds the
  5reservoir cells touching it at a rate proportional to the pressure difference
  6across the boundary face. That is the law of a BHP-controlled well, so an
  7aquifer *is* one: a "well" completed in every contact cell, with the aquifer
  8pressure for `bhp` and, for `WI`, the transmissibility of the boundary face(s)
  9of the cell -- `minires.wells.aquifer_WI`, which the record key
 10`aquifer=True` applies for you. Nothing else of the model is involved: the
 11influx enters the equations like a well's, anchors the pressure (so that even
 12an incompressible reservoir needs no injector), and is reported like a well's
 13(`minires.wells.Wells.actual_rates`).
 14
 15Here an elliptic reservoir (ref `examples.inactive_cells`) has an aquifer
 16along its western end, and a single producer to the east, on BHP control,
 17so that the rate is left to follow the pressures. Two aquifer models:
 18
 19- **Constant pressure**, i.e. an infinite aquifer: the influx settles to the
 20  rate the mobilities admit, and the front sweeps the oil towards the producer.
 21- **Fetkovich**, i.e. a finite one: its pressure falls in proportion to the
 22  cumulative influx, $ p_\\mathrm{aq} = p_i (1 - W_e / W_{ei}) $, where the
 23  capacity $ W_{ei} $ is the aquifer's compressibility times its volume times
 24  $ p_i $. The influx (and hence the production) then declines exponentially,
 25  with time constant $ W_{ei} / (J p_i) $, $ J $ being the aquifer's total
 26  productivity. This is closed-loop control, implemented as a
 27  `minires.ResSim.well_controls` override, as in `examples.well_control`.
 28
 29In the figure: the oil saturation under the constant-pressure aquifer as the
 30front advances (left), and under the Fetkovich aquifer at the end, the sweep
 31having stalled as it depleted (middle) -- the aquifer's contact stroked along
 32the boundary faces (`minires.plotting.Plot2D.plt_faces`) in place of its
 33ring of well markers (hidden by `wells=dict(exclude=...)`); and the influx
 34rates and aquifer pressure of the two models (right).
 35"""
 36
 37import matplotlib.pyplot as plt
 38import numpy as np
 39
 40from minires import Grid2D, ResSim
 41from minires.geostat import gaussian_fields
 42from minires.plotting import show
 43
 44rng = np.random.default_rng(3)  # Reproducibility (the values are regression tested)
 45
 46## An elliptic reservoir, incompressible, heterogeneous
 47p_i, p_bh = 1., 0.  # initial aquifer pressure, and the producer's
 48dt, nSteps = .02, 120
 49tt = dt * np.arange(1, nSteps + 1)
 50mesh = Grid2D(Lx=1, Ly=1, Nx=40, Ny=40).mesh
 51logK = gaussian_fields(mesh, r=.15, rng=rng)[0].reshape(mesh[0].shape)
 52# logK = 3 * smooth(smooth(rng.standard_normal((40, 40))))  # the old way
 53
 54def make(cls=ResSim):
 55    model = cls(Lx=1, Ly=1, Nx=40, Ny=40, K=np.exp(logK))
 56    X, Y = model.mesh
 57    model.active = ((X - .5) / .45)**2 + ((Y - .5) / .3)**2 <= 1
 58    # The boundary cells: active, with a neighbour that is not (or is off-grid)
 59    act = np.pad(model.active, 1, constant_values=False)
 60    boundary = model.active & ~(act[:-2, 1:-1] & act[2:, 1:-1] & act[1:-1, :-2] & act[1:-1, 2:])
 61    contact = boundary & (X < .25)  # the western end of the ellipse is the aquifer's
 62    model.wells = [dict(name="Aq", xy=np.column_stack([X[contact], Y[contact]]),
 63                        aquifer=True, bhp=p_i),
 64                   dict(name="Prd", xy=[.85, .5], bhp=p_bh, rw=1e-3)]
 65    return model
 66
 67## The Fetkovich aquifer: its pressure falls with the cumulative influx
 68W_ei = .25  # the capacity (about half the pore volume)
 69
 70class Fetkovich(ResSim):
 71    def well_controls(self, S, P, k):
 72        ctrl = super().well_controls(S, P, k)
 73        aq = self.wells.group == 0                          # the aquifer's completions
 74        W_e = self.wells.actual_rates[aq, :k].sum() * dt    # cumulative influx so far
 75        ctrl["bhp"][aq] = p_i * (1 - W_e / W_ei)
 76        return ctrl
 77
 78## Simulate
 79infinite, finite = make(), make(Fetkovich)
 80S0 = np.zeros(infinite.Nxy)
 81SS_inf, PP_inf = infinite.sim(dt, nSteps, S0, pbar=False)
 82SS_fin, PP_fin = finite.sim(dt, nSteps, S0, pbar=False)
 83q_inf, q_fin = infinite.wells.rates_by_well[0], finite.wells.rates_by_well[0]
 84p_aq = finite.wells.actual_bhp[0]  # the aquifer's, as controlled
 85
 86# The aquifer supplies the producer exactly (incompressible), and its water shows up
 87assert np.allclose(infinite.wells.rates_by_well.sum(0), 0)
 88assert np.allclose(finite.wells.rates_by_well.sum(0), 0)
 89assert SS_inf[-1].max() > .5 and (q_inf > 0).all()
 90# The finite aquifer depletes: exponentially, on the time scale W_ei / (J p_i),
 91# where J is the aquifer's initial productivity, q_fin[0] / (p_i - p_bh)
 92# (approximately: J follows the mobilities, which change as the water advances)
 93assert q_fin[-1] < .2 * q_fin[0] and (np.diff(p_aq) < 0).all()
 94tau = W_ei / (q_fin[0] / (p_i - p_bh))
 95assert np.isclose(np.log(q_fin[0] / q_fin[-1]) / (tt[-1] - tt[0]), 1 / tau, rtol=.3)
 96
 97## Plot
 98fig, axs = plt.subplots(num="Aquifer", clear=True, ncols=3, figsize=(12, 3.8))
 99kws: dict = dict(finalize=False, labels=False, colorbar=False)
100k = nSteps // 3
101kws["wells"] = dict(exclude=["Aq"])  # the contact is drawn as a stroke instead
102xy_aq = infinite.wells.xy[infinite.wells.group == 0]
103for ax, model, S, title in [(axs[0], infinite, SS_inf[k], f"const. pressure), t = {k*dt:.2f}"),
104                            (axs[1], finite, SS_fin[-1], f"Fetkovich), t = {nSteps*dt:.2f}")]:
105    model.plt_field(ax, S, "oil", title="Oil sat. (" + title, **kws)
106    model.plt_faces(ax, xy_aq, label="Aquifer")
107axs[0].legend(loc="lower right", fontsize="small")
108ax = axs[2]
109ax.plot(tt, q_inf, label="Influx, constant pressure")
110ax.plot(tt, q_fin, label="Influx, Fetkovich")
111ax.plot(tt, q_fin[0] * np.exp(-(tt - tt[0]) / tau), "k--", lw=1,
112        label=r"$q_0 \, e^{-t / \tau}$, $\tau = W_{ei} / (J p_i)$")
113ax.set(xlabel="Time", ylabel="Aquifer influx rate", ylim=(0, None))
114ax.legend(loc="upper right", fontsize="small")
115ax2 = ax.twinx()
116ax2.plot(tt, p_aq, "C1:", label="Aquifer pressure (Fetkovich)")
117ax2.set(ylabel="Aquifer pressure", ylim=(0, None))
118ax2.legend(loc="center right", fontsize="small")
119fig.tight_layout()
120
121# Regression values, checked by `tests/test_examples.py`.
122__digest__ = dict(sat_inf   = SS_inf[-1][infinite.active.ravel()],
123                  sat_fin   = SS_fin[-1][finite.active.ravel()],
124                  q_inf     = q_inf,
125                  q_fin     = q_fin,
126                  p_aq      = p_aq)
127
128if __name__ == "__main__":
129    show()