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
 37from mpl_tools.place import freshfig
 38import numpy as np
 39from scipy.ndimage import uniform_filter as smooth
 40
 41from minires import ResSim
 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)
 50logK = 3 * smooth(smooth(rng.standard_normal((40, 40))))
 51
 52def make(cls=ResSim):
 53    model = cls(Lx=1, Ly=1, Nx=40, Ny=40, K=np.exp(logK))
 54    X, Y = model.mesh
 55    model.active = ((X - .5) / .45)**2 + ((Y - .5) / .3)**2 <= 1
 56    # The boundary cells: active, with a neighbour that is not (or is off-grid)
 57    act = np.pad(model.active, 1, constant_values=False)
 58    boundary = model.active & ~(act[:-2, 1:-1] & act[2:, 1:-1] & act[1:-1, :-2] & act[1:-1, 2:])
 59    contact = boundary & (X < .25)  # the western end of the ellipse is the aquifer's
 60    model.wells = [dict(name="Aq", xy=np.column_stack([X[contact], Y[contact]]),
 61                        aquifer=True, bhp=p_i),
 62                   dict(name="Prd", xy=[.85, .5], bhp=p_bh, rw=1e-3)]
 63    return model
 64
 65## The Fetkovich aquifer: its pressure falls with the cumulative influx
 66W_ei = .25  # the capacity (about half the pore volume)
 67
 68class Fetkovich(ResSim):
 69    def well_controls(self, S, P, k):
 70        ctrl = super().well_controls(S, P, k)
 71        aq = self.wells.group == 0                          # the aquifer's completions
 72        W_e = self.wells.actual_rates[aq, :k].sum() * dt    # cumulative influx so far
 73        ctrl["bhp"][aq] = p_i * (1 - W_e / W_ei)
 74        return ctrl
 75
 76## Simulate
 77infinite, finite = make(), make(Fetkovich)
 78S0 = np.zeros(infinite.Nxy)
 79SS_inf, PP_inf = infinite.sim(dt, nSteps, S0, pbar=False)
 80SS_fin, PP_fin = finite.sim(dt, nSteps, S0, pbar=False)
 81q_inf, q_fin = infinite.wells.rates_by_well[0], finite.wells.rates_by_well[0]
 82p_aq = finite.wells.actual_bhp[0]  # the aquifer's, as controlled
 83
 84# The aquifer supplies the producer exactly (incompressible), and its water shows up
 85assert np.allclose(infinite.wells.rates_by_well.sum(0), 0)
 86assert np.allclose(finite.wells.rates_by_well.sum(0), 0)
 87assert SS_inf[-1].max() > .5 and (q_inf > 0).all()
 88# The finite aquifer depletes: exponentially, on the time scale W_ei / (J p_i),
 89# where J is the aquifer's initial productivity, q_fin[0] / (p_i - p_bh)
 90assert q_fin[-1] < .2 * q_fin[0] and (np.diff(p_aq) < 0).all()
 91tau = W_ei / (q_fin[0] / (p_i - p_bh))
 92assert np.isclose(np.log(q_fin[0] / q_fin[-1]) / (tt[-1] - tt[0]), 1 / tau, rtol=.25)
 93
 94## Plot
 95fig, axs = freshfig("Aquifer", ncols=3, figsize=(12, 3.8))
 96kws: dict = dict(finalize=False, labels=False, colorbar=False)
 97k = nSteps // 3
 98kws["wells"] = dict(exclude=["Aq"])  # the contact is drawn as a stroke instead
 99xy_aq = infinite.wells.xy[infinite.wells.group == 0]
100for ax, model, S, title in [(axs[0], infinite, SS_inf[k], f"const. pressure), t = {k*dt:.2f}"),
101                            (axs[1], finite, SS_fin[-1], f"Fetkovich), t = {nSteps*dt:.2f}")]:
102    model.plt_field(ax, S, "oil", title="Oil sat. (" + title, **kws)
103    model.plt_faces(ax, xy_aq, label="Aquifer")
104axs[0].legend(loc="lower right", fontsize="small")
105ax = axs[2]
106ax.plot(tt, q_inf, label="Influx, constant pressure")
107ax.plot(tt, q_fin, label="Influx, Fetkovich")
108ax.plot(tt, q_fin[0] * np.exp(-(tt - tt[0]) / tau), "k--", lw=1,
109        label=r"$q_0 \, e^{-t / \tau}$, $\tau = W_{ei} / (J p_i)$")
110ax.set(xlabel="Time", ylabel="Aquifer influx rate", ylim=(0, None))
111ax.legend(loc="upper right", fontsize="small")
112ax2 = ax.twinx()
113ax2.plot(tt, p_aq, "C1:", label="Aquifer pressure (Fetkovich)")
114ax2.set(ylabel="Aquifer pressure", ylim=(0, None))
115ax2.legend(loc="center right", fontsize="small")
116fig.tight_layout()
117
118# Regression values, checked by `tests/test_examples.py`.
119__digest__ = dict(sat_inf   = SS_inf[-1][infinite.active.ravel()],
120                  sat_fin   = SS_fin[-1][finite.active.ravel()],
121                  q_inf     = q_inf,
122                  q_fin     = q_fin,
123                  p_aq      = p_aq)
124
125if __name__ == "__main__":
126    show()