Consider a two-aquifer system that contains one inhomogeneity. Inside the inhomogeneity the transmissivity of the top aquifer is much lower and the transmissivity of the bottom aquifer is much higher than outside the inhomogeneity. Aquifer properties are given in Table 1 and the inhomogeneity data is given in table 2. There is a uniform gradient of 0.002 in Southeastern direction.
| Layer | $k$ (m/d) | $z_b$ (m) | $z_t$ | $c$ (days) |
|---|---|---|---|---|
| Aquifer 0 | 10 | 0 | 20 | |
| Leaky Layer 1 | -10 | 0 | 4000 | |
| Aquifer 1 | 20 | -30 | 10 |
| Layer | $k$ (m/d) | $z_b$ (m) | $z_t$ | $c$ (days) |
|---|---|---|---|---|
| Aquifer 0 | 2 | 0 | 20 | - |
| Leaky Layer 1 | - | -10 | 0 | 500 |
| Aquifer 1 | 80 | -30 | -10 | - |
A layout of the nodes of the inhomogeneity are shown in Fig. 1 (inhomogeneity 2 will be added later on). A well is located in the top aquifer inside inhomogeneity 1 (the black dot).
In [1]:
%matplotlib inline
from timml import *
from pylab import *
figsize = (8, 8)
In [2]:
ml = ModelMaq(kaq=[10, 20], z=[20, 0, -10, -30], c=[4000])
xy1 = [(0, 600), (-100, 400), (-100, 200), (100, 100), (300, 100), (500, 100),
(700, 300), (700, 500), (600, 700), (400, 700), (200, 600)]
p1 = PolygonInhomMaq(ml, xy=xy1,
kaq=[2, 80], z=[20, 0, -10, -30], c=[500],
topboundary='conf', order=3, ndeg=2)
rf = Constant(ml, xr=1000, yr=0, hr=40)
uf = Uflow(ml, slope=0.002, angle=-45)
w = Well(ml, xw=400, yw=400, Qw=500, rw=0.2, layers=0)
ml.solve()
In [3]:
print('Leakage factor of the background aquifer is:', ml.aq.lab)
print('Leakage factor of the inhomogeneity is:', p1.lab)
In [4]:
ml.contour(win=[-200, 800, 0, 800], ngr=50, layers=[0, 1], levels=50, labels=1, decimals=2, figsize=figsize)
In [5]:
ml.contour(win=[-1200, 1800, -1000, 1800], ngr=50, layers=[0, 1], levels=50, labels=1, decimals=2, figsize=figsize)
In [6]:
htop = ml.headalongline(linspace(101, 499, 100), 100 + 0.001 * ones(100))
hbot = ml.headalongline(linspace(101, 499, 100), 100 - 0.001 * ones(100))
figure()
plot(linspace(101, 499, 100), htop[0])
plot(linspace(101, 499, 100), hbot[0])
Out[6]:
In [7]:
qtop = zeros(100)
qbot = zeros(100)
layer = 1
x = linspace(101, 499, 100)
for i in range(100):
qx, qy = ml.disvec(x[i], 100 + 0.001)
qtop[i] = qy[layer]
qx, qy = ml.disvec(x[i], 100 - 0.001)
qbot[i] = qy[layer]
figure()
plot(x, qtop)
plot(x, qbot)
Out[7]:
In [8]:
ml.plot(win=[-200, 800, 0, 800], orientation='both', figsize=figsize)
w.plotcapzone(hstepmax=50, nt=20, zstart=10, tmax=20 * 365.25, orientation='both')
Change the elevation of the bottom of aquifer 1 from -30 to -20 inside inhomogeneity 1 (you need to recreate the model). Make a contour plot with a cross-section below it and start some pathlines from $x=-200$, $y=700$. (note that the cross-section shows the elevation layers in the background aquifer, not the inhomogeneity). Note that the pathlines jump when they enter and exit inhomogeneity 1. This is caused by the jump in the base. It meets all continuity conditions and is an approximation of the smooth change in elevation that occurs over a distance of approximately one aquifer thickness from the boundary.
In [9]:
ml = ModelMaq(kaq=[10, 20], z=[20, 0, -10, -30], c=[4000])
xy1 = [(0, 600), (-100, 400), (-100, 200), (100, 100), (300, 100), (500, 100),
(700, 300), (700, 500), (600, 700), (400, 700), (200, 600)]
p1 = PolygonInhomMaq(ml, xy=xy1,
kaq=[2, 80], z=[20, 0, -10, -40], c=[500],
topboundary='conf', order=5, ndeg=3)
rf = Constant(ml, xr=1000, yr=0, hr=40)
uf = Uflow(ml, slope=0.002, angle=-45)
w = Well(ml, xw=400, yw=400, Qw=500, rw=0.2, layers=0)
ml.solve()
In [10]:
ml.plot(win=[-200, 800, 0, 800], orientation='both', figsize=figsize)
ml.tracelines(-200 * np.ones(2), 700 * np.ones(2), [-25, -15], hstepmax=25, orientation='both')
A second inhomogeneity is added, which shares part of its boundary with the first inhomogeneity. The aquifer properties for the inhomogeneity are provided in table 3. Inside this second inhomogeneity, the transmissivity of both the bottom aquifer and the resistance of the leaky layer are reduced. The input is now somewhat complicated. First the data of the two inhomgeneities is entered. Second, analytic elements are placed along the boundary of the inhomogeneity with MakeInhomPolySide.
This routine places line elements along a string of points, but it requires that the
aquifer data is the same on the left and right sides of the line. Hence, for this case we need to break the boundary up
in three sections: One section with the background aquifer on one side and inhom1 on the other, one section
with the background aquifer on one side and inhom2 on the other, and one section with inhom1 on one side
and inhom2 on the other. The input file is a bit longer
| Layer | $k$ (m/d) | $z_b$ (m) | $z_t$ | $c$ (days) |
|---|---|---|---|---|
| Aquifer 0 | 2 | -20 | 0 | - |
| Leaky Layer 1 | - | -40 | -20 | 50 |
| Aquifer 1 | 8 | -80 | -40 | - |
In [11]:
ml = ModelMaq(kaq=[10, 20], z=[20, 0, -10, -30], c=[4000])
xy1 = [(0, 600), (-100, 400), (-100, 200), (100, 100), (300, 100), (500, 100),
(700, 300), (700, 500), (600, 700), (400, 700), (200, 600)]
p1 = PolygonInhomMaq(ml, xy=xy1,
kaq=[2, 80], z=[20, 0, -10, -30], c=[500],
topboundary='conf', order=4, ndeg=2)
xy2 = [(0, 600), (200, 600), (400, 700), (400, 900), (200, 1100), (0, 1000), (-100, 800)]
p2 = PolygonInhomMaq(ml, xy=xy2,
kaq=[2, 8], z=[20, 0, -10, -30], c=[50],
topboundary='conf', order=4, ndeg=2)
rf = Constant(ml, xr=1000, yr=0, hr=40)
uf = Uflow(ml, slope=0.002, angle=-45)
w = Well(ml, xw=400, yw=400, Qw=500, rw=0.2, layers=0)
ml.solve()
In [12]:
ml.contour(win=[-200, 1000, 0, 1200], ngr=50, layers=[0, 1],
levels=20, figsize=figsize)
In [13]:
xy1[::-1]
Out[13]:
In [14]:
ml = ModelMaq(kaq=[10, 20], z=[20, 0, -10, -30], c=[4000])
p1 = PolygonInhomMaq(ml, xy=xy1[::-1],
kaq=[2, 80], z=[20, 0, -10, -30], c=[500],
topboundary='conf', order=3, ndeg=2)
rf = Constant(ml, xr=1000, yr=0, hr=40)
uf = Uflow(ml, slope=0.002, angle=-45)
w = Well(ml, xw=400, yw=400, Qw=500, rw=0.2, layers=0)
ml.solve()
In [15]:
ml.contour(win=[-200, 800, 0, 800], ngr=50, layers=[0, 1], levels=50, labels=1, decimals=2, figsize=figsize)
In [ ]: