Cross-Sectional Models#

This notebook demonstrates how to build and solve cross-sectional (1D) models in TimML, including examples with strip inhomogeneities and infinitely long line elements.

import matplotlib.pyplot as plt
import numpy as np

import timml as tml

plt.rcParams["figure.figsize"] = (4, 3)
plt.rcParams["figure.autolayout"] = True

Two-layer model with head-specified line-sink#

Two-layer aquifer bounded on top by a semi-confined layer. Head above the semi-confining layer is 5. Head line-sink located at \(x=0\) with head equal to 2, cutting through layer 0 only.

ml = tml.ModelMaq(
    kaq=[1, 2], z=[4, 3, 2, 1, 0], c=[1000, 1000], topboundary="semi", hstar=5
)
ls = tml.HeadLineSink1D(ml, xls=0, hls=2, layers=0)
ml.solve()

x = np.linspace(-200, 200, 101)
h = ml.headalongline(x, np.zeros_like(x))
plt.plot(x, h[0], label="layer 0")
plt.plot(x, h[1], label="layer 1")
plt.legend(loc="best")
plt.grid()

Number of elements, Number of equations: 2 , 1
.
.

solution complete
../_images/9e2fe725b996d9da64814c3c49c1c7584544c188977181fa031439d2edc7f600.png

1D inhomogeneity#

Three strips with semi-confined conditions on top of all three

ml = tml.ModelXsection(naq=2)
tml.XsectionMaq(
    ml,
    x1=-np.inf,
    x2=-50,
    kaq=[1, 2],
    z=[4, 3, 2, 1, 0],
    c=[1000, 1000],
    npor=0.3,
    topboundary="semi",
    hstar=5,
)
tml.XsectionMaq(
    ml,
    x1=-50,
    x2=50,
    kaq=[1, 2],
    z=[4, 3, 2, 1, 0],
    c=[1000, 1000],
    npor=0.3,
    topboundary="semi",
    hstar=4.5,
)
tml.XsectionMaq(
    ml,
    x1=50,
    x2=np.inf,
    kaq=[1, 2],
    z=[4, 3, 2, 1, 0],
    c=[1000, 1000],
    npor=0.3,
    topboundary="semi",
    hstar=4,
)
ml.solve()

Number of elements, Number of equations: 7 , 8
.
.
.
.
.
.
.

solution complete

fig, ax = plt.subplots(1, 1, figsize=(10, 3))
ml.plots.xsection(xy=[(-150, 0), (150, 0)], ax=ax, params=True);
../_images/96093d631c480579c6d0548a319d5ae0487b39eb6b19c1b6ce3cebebb04fac46.png 
x = np.linspace(-200, 200, 101)
h = ml.headalongline(x, np.zeros(101))
plt.plot(x, h[0], label="layer 0")
plt.plot(x, h[1], label="layer 1")
plt.xlabel("x (m)")
plt.ylabel("head (m)")
plt.legend(loc="best")
plt.grid()
../_images/7e10b75454017bd0f152887857150850b4cac52bc6c44914fa0b10b2131f2205.png 
ml.plots.vcontoursf1D(x1=-200, x2=200, nx=100, levels=20, color="C0", figsize=(10, 3));
../_images/c15c2c8692689d8835da20f2df71a10e10f9cd93be2c49267400a0f085450dfe.png

Three strips with semi-confined conditions at the top of the strip in the middle only. The head is specified in the strip on the left and in the strip on the right.

ml = tml.ModelXsection(naq=2)
tml.XsectionMaq(
    ml,
    x1=-np.inf,
    x2=-50,
    kaq=[1, 2],
    z=[3, 2, 1, 0],
    c=[1000],
    npor=0.3,
    topboundary="conf",
)
tml.XsectionMaq(
    ml,
    x1=-50,
    x2=50,
    kaq=[1, 2],
    z=[4, 3, 2, 1, 0],
    c=[1000, 1000],
    npor=0.3,
    topboundary="semi",
    hstar=4,
)
tml.XsectionMaq(
    ml,
    x1=50,
    x2=np.inf,
    kaq=[1, 2],
    z=[3, 2, 1, 0],
    c=[1000],
    npor=0.3,
    topboundary="conf",
)
rf1 = tml.Constant(ml, -100, 0, 5)
rf2 = tml.Constant(ml, 100, 0, 5)

ml.solve()

ml.plots.xsection(xy=[(-100, 0), (100, 0)]);

Number of elements, Number of equations: 7 , 10
.
.
.
.
.
.
.

solution complete
../_images/69ffd0e1bbf767e219d2f9d803413cfccab4d457202e36d8b68574a474b82f75.png 
x = np.linspace(-200, 200, 101)
h = ml.headalongline(x, np.zeros_like(x))
Qx, _ = ml.disvecalongline(x, np.zeros_like(x))

plt.figure(figsize=(10, 3))
plt.subplot(121)
plt.plot(x, h[0], label="layer 0")
plt.plot(x, h[1], label="layer 1")
plt.plot([-100, 100], [4, 4], "b.", label="fixed heads")
plt.xlabel("x (m)")
plt.ylabel("head (m)")
plt.legend(loc="best")
plt.grid()
plt.subplot(122)
plt.plot(x, Qx[0], label="layer 0")
plt.plot(x, Qx[1], label="layer 1")
plt.xlabel("x (m)")
plt.ylabel("$Q_x$ (m$^2$/d)")
plt.grid()
../_images/02b4f746024e9ebbff0cd7af7c0dd60f9e8c124b4ca2afc68b81409df34bdf53.png 
ml.plots.vcontoursf1D(x1=-200, x2=200, nx=100, levels=20, color="C0", figsize=(10, 3));
../_images/5824eb20cf093f5fcbf95fef97e144a2b9b9e590eb0b29af51379ddd7df7a7c5.png

Impermeable wall#

Flow from left to right in three-layer aquifer with impermeable wall in bottom 2 layers

# need ModelMaq here since Uflow requires a confined background aquifer
ml = tml.ModelMaq(kaq=[1, 2, 4], z=[5, 4, 3, 2, 1, 0], c=[5000, 1000])
uf = tml.Uflow(ml, 0.002, 0)
rf = tml.Constant(ml, 100, 0, 20)
ld1 = tml.ImpLineDoublet1D(ml, xld=0, layers=[0, 1])

ml.solve()

Number of elements, Number of equations: 3 , 3
.
.
.

solution complete

x = np.linspace(-100, 100, 101)
h = ml.headalongline(x, np.zeros_like(x))
Qx, _ = ml.disvecalongline(x, np.zeros_like(x))

plt.figure(figsize=(10, 3))
plt.subplot(121)
plt.title("head")
plt.plot(x, h[0], label="layer 0")
plt.plot(x, h[1], label="layer 1")
plt.plot(x, h[2], label="layer 2")
plt.xlabel("x (m)")
plt.ylabel("head (m)")
plt.legend(loc="best")
plt.grid()
plt.subplot(122)
plt.title("Qx")
plt.plot(x, Qx[0], label="layer 0")
plt.plot(x, Qx[1], label="layer 1")
plt.plot(x, Qx[2], label="layer 2")
plt.xlabel("x (m)")
plt.ylabel("$Q_x$ (m$^2$/d)")
plt.legend(loc="best")
plt.grid()
../_images/29d462f9e185770e2fb1ef9a08b119e36d21cd6dab0ea7985c5fbdf65a042848.png 
ax = ml.plots.vcontoursf1D(
    x1=-200, x2=200, nx=100, levels=20, color="C0", figsize=(10, 3)
)
ld1.plot(ax);  # plot wall
../_images/e1aeadcb2082523c76bd834b051fe14af75838e860a5bc719ed3c877cd159c2f.png

Infiltration#

Comparing solution with Xsection inhomogeneities to XsectionAreaSink solution.

ml = tml.ModelXsection(naq=2)
tml.XsectionMaq(
    ml,
    x1=-np.inf,
    x2=-50,
    kaq=[1, 2],
    z=[3, 2, 1, 0],
    c=[1000],
    npor=0.3,
    topboundary="conf",
)
tml.XsectionMaq(
    ml,
    x1=-50,
    x2=50,
    kaq=[1, 2],
    z=[3, 2, 1, 0],
    c=[1000],
    npor=0.3,
    topboundary="conf",
    N=0.001,
)
tml.XsectionMaq(
    ml,
    x1=50,
    x2=np.inf,
    kaq=[1, 2],
    z=[3, 2, 1, 0],
    c=[1000],
    npor=0.3,
    topboundary="conf",
)
tml.Constant(ml, -100, 0, 10)
tml.Constant(ml, 100, 0, 10)
ml.solve()

ml.plots.vcontoursf1D(x1=-100, x2=100, nx=100, levels=20, color="C0", figsize=(10, 3));

Number of elements, Number of equations: 7 , 10
.
.
.
.
.
.
.

solution complete
../_images/034393a090df3f1b5cdac766d6c6692e8e3e3bf0c4e21aed480c190b77bc9c06.png 
ml2 = tml.ModelMaq(kaq=[1, 2], z=[3, 2, 1, 0], c=[1000], topboundary="conf")
tml.XsectionAreaSink(ml2, -50, 50, 0.001)
tml.Constant(ml2, -100, 0, 10)
ml2.solve()
ml2.plots.vcontoursf1D(x1=-100, x2=100, nx=100, levels=20, color="C0", figsize=(10, 3));

Number of elements, Number of equations: 2 , 1
.
.

solution complete

/tmp/ipykernel_1004/3881863634.py:2: DeprecationWarning: XsectionAreaSink is only for testing purposes. It is recommended to add infiltration through XsectionMaq or Xsection3D and specifying 'N'.
  tml.XsectionAreaSink(ml2, -50, 50, 0.001)
../_images/ce6346bce090cd13607c21e39a649152e11815e7084d49ccf6630e8a604806b8.png 
x = np.linspace(-100, 100, 100)
plt.plot(x, ml.headalongline(x, 0)[0], "C0")
plt.plot(x, ml.headalongline(x, 0)[1], "C0")
plt.plot(x, ml2.headalongline(x, 0)[0], "--C1")
plt.plot(x, ml2.headalongline(x, 0)[1], "--C1")
plt.xlabel("x (m)")
plt.ylabel("head (m)")
plt.grid()
../_images/9403068b4c2bf7887c1293726bf3a613aae6ac65b6317f532cb674b4c5e574e9.png 
ml = tml.ModelXsection(naq=50)
tml.Xsection3D(ml, x1=-np.inf, x2=-5, kaq=1, z=np.arange(5, -0.1, -0.1), kzoverkh=0.1)
tml.Xsection3D(
    ml,
    x1=-5,
    x2=5,
    kaq=1,
    z=np.arange(5, -0.1, -0.1),
    kzoverkh=0.1,
    topboundary="semi",
    hstar=5.5,
    topres=3,
    topthick=0.3,
)
tml.Xsection3D(ml, x1=5, x2=np.inf, kaq=1, z=np.arange(5, -0.1, -0.1), kzoverkh=0.1)
rf1 = tml.Constant(ml, -100, 0, 5.7)
rf2 = tml.Constant(ml, 100, 0, 5.47)

ml.solve()

ml.plots.vcontoursf1D(x1=-20, x2=20, nx=100, levels=20, color="C0", figsize=(10, 3));

Number of elements, Number of equations: 7 , 202
.
.
.
.
.
.
.

solution complete
../_images/05ab30b14d9a52078b05c3df8b06fce9f62b31cbf996718ec309e80e5d12cc58.png 
ml = tml.ModelXsection(naq=5)
tml.Xsection3D(ml, x1=-np.inf, x2=-5, kaq=1, z=np.arange(5, -0.1, -1), kzoverkh=0.1)
tml.Xsection3D(
    ml,
    x1=-5,
    x2=5,
    kaq=1,
    z=np.arange(5, -0.1, -1),
    kzoverkh=0.1,
    topboundary="semi",
    hstar=5.5,
    topres=3,
    topthick=0.3,
)
tml.Xsection3D(ml, x1=5, x2=np.inf, kaq=1, z=np.arange(5, -0.1, -1), kzoverkh=0.1)
rf1 = tml.Constant(ml, -100, 0, 5.7)
rf2 = tml.Constant(ml, 100, 0, 5.47)

ml.solve()

ml.plots.vcontoursf1D(x1=-20, x2=20, nx=100, levels=20, color="C0", figsize=(10, 3));

Number of elements, Number of equations: 7 , 22
.
.
.
.
.
.
.

solution complete
../_images/66aefc56b33480b7d305b9b898f469d79bbaafea17bbb51dda4d8a30242b84dd.png 
ml = tml.ModelXsection(naq=2)
tml.XsectionMaq(
    ml,
    x1=-np.inf,
    x2=-50,
    kaq=[1, 2],
    z=[4, 3, 2, 1, 0],
    c=[1000, 1000],
    npor=0.3,
    topboundary="semi",
    hstar=15,
)
tml.XsectionMaq(
    ml,
    x1=-50,
    x2=50,
    kaq=[1, 2],
    z=[4, 3, 2, 1, 0],
    c=[1000, 1000],
    npor=0.3,
    topboundary="semi",
    hstar=13,
)
tml.XsectionMaq(
    ml,
    x1=50,
    x2=np.inf,
    kaq=[1, 2],
    z=[4, 3, 2, 1, 0],
    c=[1000, 1000],
    npor=0.3,
    topboundary="semi",
    hstar=11,
)
ml.solve()

Number of elements, Number of equations: 7 , 8
.
.
.
.
.
.
.

solution complete

from timml.linesink1d import FluxDiffLineSink1D, HeadDiffLineSink1D

ml = tml.ModelMaq(kaq=[10], z=[0, -10], topboundary="conf")
ls = tml.HeadLineSink1D(ml, xls=0, hls=1, wh="H", layers=0)
ls = tml.HeadLineSink1D(ml, xls=200, hls=0, wh="H", layers=0)
hd = HeadDiffLineSink1D(ml, xls=100)
fd = FluxDiffLineSink1D(ml, xls=100)
ml.solve()

x = np.linspace(0, 200, 101)
h = ml.headalongline(x, np.zeros_like(x))
plt.plot(x, h[0], label="layer 0")
# plt.plot(x, h[1], label="layer 1")
plt.legend(loc="best")

Number of elements, Number of equations: 2 , 2
.
.

solution complete

<matplotlib.legend.Legend at 0x70abf03af290>
../_images/d584cd399662b528e26498f9a1489bd53126db9fd6fdbea38bb935ab8c180b2a.png