1. Unconfined Aquifer Test - Vennebulten#

This example is taken from Kruseman et al. (1970).

Import packages#

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import ttim as ttm

plt.rcParams["figure.figsize"] = [5, 3]

Introduction and Conceptual Model#

In aquifer tests in unconfined aquifers, there is also the vertical component to flow to the well. The drawdown data shows the delayed water table response, a distinguishable S-shape in the log-log plot. In the early times of the drawdown, the drawdown behaves as a confined aquifer: when the aquifer releases the elastic storage. However, as pumping continues, the water table storage begins to be released, generating further drawdown and the S-shape.

This test conducted in Vennebulten, the Netherlands, is reported in Kruseman et al. (1970). The cross-section consists of a first layer up to 6 m depth of very fine and loamy sands, followed by coarse sands until 21 m deep.

In this example, we will reproduce the work of Xinzhu (2020) that compared different conceptualizations in TTim to various solutions presented in the original report (Kruseman et al., 1970) and in other software, MLU (Carson & Randall, 2012) and AQTESOLV (Duffield, 2007).

The screen of the pumping well is placed between 10 and 21 meters depth, and pumping has taken place for 25 hours at a rate of 873 m3/d. The available drawdown data comes from two piezometers, a shallow one, screened at 3 m depth, and a deeper one, screened in the depths between 12 to 19 m. Both wells are located 90 m from the pumping well.

The conceptual model of the aquifer is displayed below:

Hide code cell source

 
import matplotlib.pyplot as plt

##Now printing the conceptual model figure:

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
# sky
sky = plt.Rectangle((-20, 2), width=150, height=5, fc="b", zorder=0, alpha=0.1)
ax.add_patch(sky)

# Aquifer:
ground = plt.Rectangle(
    (-20, -6),
    width=150,
    height=8,
    fc=np.array([209, 179, 127]) / 255,
    zorder=0,
    alpha=0.9,
    hatch=".",
)
ax.add_patch(ground)

# Aquifer 2:
ground2 = plt.Rectangle(
    (-20, -21),
    width=150,
    height=15,
    fc=np.array([209, 179, 127]) / 255,
    zorder=0,
    alpha=0.9,
    hatch="o",
)
ax.add_patch(ground2)


well = plt.Rectangle(
    (-1.5, -21), width=3, height=23, fc=np.array([200, 200, 200]) / 255, zorder=1
)
ax.add_patch(well)

# Wellhead
wellhead = plt.Rectangle(
    (-2, 2), width=4, height=2.5, fc=np.array([200, 200, 200]) / 255, zorder=2, ec="k"
)
ax.add_patch(wellhead)

# Screen for the well:
screen = plt.Rectangle(
    (-1.5, -21),
    width=3,
    height=11,
    fc=np.array([200, 200, 200]) / 255,
    alpha=1,
    zorder=2,
    ec="k",
    ls="--",
)
screen.set_linewidth(2)
ax.add_patch(screen)
pumping_arrow = plt.Arrow(x=2, y=3.5, dx=5, dy=0, color="#00035b")
ax.add_patch(pumping_arrow)
ax.text(x=7, y=3.5, s=r"$ Q = 873$ m$^3$/d", fontsize="large")

# Piezometers
piez1 = plt.Rectangle(
    (89, -21), width=2, height=23, fc=np.array([200, 200, 200]) / 255, zorder=1
)
screen_piez_1 = plt.Rectangle(
    (89, -19),
    width=2,
    height=7,
    fc=np.array([200, 200, 200]) / 255,
    alpha=1,
    zorder=2,
    ec="k",
    ls="--",
)
screen_piez_1.set_linewidth(2)
screen_piez_2 = plt.Rectangle(
    (89, -3),
    width=2,
    height=0.5,
    fc=np.array([200, 200, 200]) / 255,
    alpha=1,
    zorder=2,
    ec="k",
    ls="--",
)
screen_piez_2.set_linewidth(2)

ax.add_patch(piez1)
ax.add_patch(screen_piez_1)
ax.add_patch(screen_piez_2)


# last line
line = plt.Line2D(xdata=[-200, 1200], ydata=[2, 2], color="k")
ax.add_line(line)

# Water table
line2 = plt.Line2D(xdata=[-200, 1200], ydata=[0, 0], color="b")
ax.add_line(line2)

ax.text(-18, 0.5, s="Water Table", fontsize="large", color="b", bbox={"fc": "w"})
ax.text(93, -3, s="Shallow piezometer", bbox={"fc": "w"})
ax.text(93, -16, s="Deeper piezometer", bbox={"fc": "w"})

ax.set_xlim([-20, 130])
ax.set_ylim([-21, 7])
ax.set_xlabel("Distance [m]")
ax.set_ylabel("Relative height [m]")
ax.set_title("Conceptual Model - Vennebulten Example");
../_images/4485d28ceac5444a372ea5219b9a29dafb286727c2da0fc2e15ed9d86f278108.png

Set basic parameters#

b = -21  # aquifer thickness, m
r = 90  # distance from observation wells to pumping well, m
Q = 873  # constant discharge, m^3/d

Load data of the two piezometers#

data1 = np.loadtxt("data/venne_shallow.txt", skiprows=1)
ts = data1[:, 0] / 60 / 24  # convert min to days
hs = data1[:, 1]

data2 = np.loadtxt("data/venne_deep.txt", skiprows=1)
td = data2[:, 0] / 60 / 24  # convert min to days
hd = data2[:, 1]

Create a conceptual one-layer model#

Both Kruseman et al. (1970) and AQTESOLV solutions that use the Neuman method (Neumann, 1969) assume a one layer unconfined model. To compare TTim with both, we begin by modelling a one-layer aquifer.

For the unconfined test, the preferred method for modelling is to use the Model3D class. Model3D assumes the system is a vertical stacking of aquifer layers. Vertical flow is computed between layers by calculating the vertical resistance between layers. Vertical resistance between the aquifer layers is determined as the resistance from the middle of one layer to the middle of the next layer. The vertical anisotropy can be specified for each layer.

Model construction is similar to the ModelMaq class. We detail it below:

For our Model3D model, we have to set:

  • The hydraulic conductivity: kaq. It is a list/array with a float element for every aquifer, for example: [kaq0,kaq1]. We can also set a float value. In this case, the same kaq is assumed for every layer.

  • The top and bottom of each aquifer: z defined by a list/array [zt0,zb0,zt1,zb1,...], where the inputs are a sequence of top and bottoms of the aquifer layers.

  • The specific storage: Saq. It is a list/array with a float element for every aquifer, for example: [Saq0, Saq1]. We can also set a float value. In this case, the same Saq is assumed for every layer.

  • The minimum time for which TTim solve the groundwater flow: tmin, a float.

  • And the maximum time: tmax, float.

  • TTim automatically assumes the topboundary is confined. In this case, we also assume the topboundary is confined, so we do not need to set this parameter. In the code example, the parameter is set for clarity.

  • The vertical anisotropy, defined by the parameter: kzoverkh, which means the vertical hydraulic conductivity divided by the horizontal conductivity. This parameter is a list/array with a float element for every aquifer, for example: [kzoverkh0,kzoverkh1]. We can also set a float value. In this case, the same kzoverkh is assumed for every layer. If one does not set this parameter, a isotropic model is considered: kzoverkh = 1

  • phreatictop: Is a boolean (True/False). If True, the first element in Saq is considered phreatic storage (Specific Yield) and is not multiplied by the layer thickness. The default value is True. This parameter is relevant for the unconfined aquifer test, and we will show underneath how to set it.

To reproduce the one-layer aquifer model in Kruseman et al. (1970), we will build a two-layer Model3D model. The first layer is a very thin (0.1 m thick) layer with phreatic storage, followed by the 21 m thick aquifer layer. This thin layer is how TTim accounts for the water table storage in the unconfined situation. The first conceptual model is represented in the image below.

# Model Figure - One - layer model
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)


# Aquifer:
ground = plt.Rectangle(
    (-20, -21), width=150, height=23, fc="w", zorder=0, alpha=0.9, hatch="."
)
ax.add_patch(ground)

well = plt.Rectangle((-1.5, -21), width=3, height=23, fc="w", zorder=1, ec="k")
ax.add_patch(well)

# Wellhead
wellhead = plt.Rectangle((-2, 2), width=4, height=2.5, fc="w", zorder=2, ec="k")
ax.add_patch(wellhead)

# Screen for the well:
screen = plt.Rectangle(
    (-1.5, -21), width=3, height=11, fc="w", alpha=1, zorder=2, ec="k", ls="--"
)
screen.set_linewidth(2)
ax.add_patch(screen)

# Piezometers
piez1 = plt.Rectangle((89, -21), width=2, height=23, fc="w", zorder=1, ec="k")
screen_piez_1 = plt.Rectangle(
    (89, -19), width=2, height=7, fc="w", alpha=1, zorder=2, ec="k", ls="--"
)
screen_piez_1.set_linewidth(2)
screen_piez_2 = plt.Rectangle(
    (89, -3), width=2, height=0.5, fc="w", alpha=1, zorder=2, ec="k", ls="--"
)
screen_piez_2.set_linewidth(2)

ax.add_patch(piez1)
ax.add_patch(screen_piez_1)
ax.add_patch(screen_piez_2)


# last line
line = plt.Line2D(xdata=[-200, 1200], ydata=[2, 2], color="k")
ax.add_line(line)

# Water table
line2 = plt.Line2D(xdata=[-200, 1200], ydata=[0, 0], color="b")
ax.add_line(line2)

ax.text(-18, 0.5, s="Water Table", fontsize="large", color="b", bbox={"fc": "w"})
ax.text(93, -3, s="Shallow piezometer", bbox={"fc": "w"})
ax.text(93, -16, s="Deeper piezometer", bbox={"fc": "w"})

ax.set_xlim([-20, 130])
ax.set_ylim([-21, 7])
ax.set_xlabel("Distance [m]")
ax.set_ylabel("Relative height [m]")
ax.set_title("Conceptual Model 1 - Vennebulten Example");
../_images/d2c180f46a2452bccdcee123f46b9141019f1c14d9fa73f7ab05ae4e5af8256c.png 
ml_1 = ttm.Model3D(
    kaq=10,
    z=[0, -0.1, b],
    Saq=[0.1, 1e-4],
    tmin=1e-4,
    tmax=1.1,
    kzoverkh=1,
    phreatictop=True,
)
w_1 = ttm.Well(ml_1, xw=0, yw=0, rw=0.1, tsandQ=[(0, Q)])
ml_1.solve()

self.neq  1
solution complete

Calibrate the one layer model with the shallow piezometer#

We begin the initial model by adding the shallow observation well as the observation for the residuals calibration. And we calibrate hydraulic conductivity, specific yield and specific storage of our one layer unconfined aquifer:

# calibrate with data of shallow piezometer
# unknown parameters: kaq, Saq
ca_1 = ttm.Calibrate(ml_1)
ca_1.set_parameter(name="kaq0_1", initial=10)
ca_1.set_parameter(name="Saq0", initial=0.2)
ca_1.set_parameter(name="Saq1", initial=1e-4, pmin=0)
ca_1.set_parameter(name="kzoverkh0_1", initial=1, pmin=1e-5)
ca_1.series(name="obs", x=r, y=0, t=ts, h=hs, layer=0)  # shallow piezometer
ca_1.fit()

........................
...........................
...........................
..........................
...............
Fit succeeded.

/tmp/ipykernel_1462/3076412918.py:4: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'kaq0_1' to silence this warning. The parameter name can still include layer info, but this will be ignored in a future version of TTim.
  ca_1.set_parameter(name="kaq0_1", initial=10)
/tmp/ipykernel_1462/3076412918.py:5: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'Saq0' to silence this warning. The parameter name can still include layer info, but this will be ignored in a future version of TTim.
  ca_1.set_parameter(name="Saq0", initial=0.2)
/tmp/ipykernel_1462/3076412918.py:6: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'Saq1' to silence this warning. The parameter name can still include layer info, but this will be ignored in a future version of TTim.
  ca_1.set_parameter(name="Saq1", initial=1e-4, pmin=0)
/tmp/ipykernel_1462/3076412918.py:7: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'kzoverkh0_1' to silence this warning. The parameter name can still include layer info, but this will be ignored in a future version of TTim.
  ca_1.set_parameter(name="kzoverkh0_1", initial=1, pmin=1e-5)

display(ca_1.parameters)
print(f"RMSE: {ca_1.rmse():.4f} m")
layers optimal std perc_std pmin pmax initial inhoms parray
kaq0_1_0_1 None 137.281122 8.094764 5.896487 -inf inf 10.0000 None [[137.28112218891238, 137.28112218891238]]
Saq0_0_0 None 0.006004 0.144598 2408.481102 -inf inf 0.2000 None [[0.006003688838702366]]
Saq1_1_1 None 0.000505 0.006880 1361.917152 0.00000 inf 0.0001 None [[0.0005052051840837013]]
kzoverkh0_1_0_1 None 43.607887 1438.054461 3297.693524 0.00001 inf 1.0000 None [[43.6078868512209, 43.6078868512209]] 
RMSE: 0.0035 m

hs_1 = ml_1.head(r, 0, ts)
hd_1 = ml_1.head(r, 0, td)
plt.semilogx(ts, hs, ".", label="shallow obs")
plt.semilogx(ts, hs_1[0], label="shallow ttim")
plt.semilogx(td, hd, ".", label="deep obs")
plt.semilogx(td, hd_1[0], "--", label="deep ttim")
plt.xlabel("time [d]")
plt.ylabel("drawdown [m]")
plt.title("TTim Unconfined Model Results - Shallow Piezometer")
plt.legend()
plt.grid()
../_images/0b4b07ef982acf22450e563bf075f8a910f04b1750e228fbb12acb1f60b0e69a.png

Calibrate the one layer model with the deeper piezometer#

In this second approach, we adjust the model to the deeper piezometer, as done by Kruseman and de Ridder (1970).

# calibrate with data of deeper piezometer
# unknown parameters: kaq, Saq, kzoverkh
ca_2 = ttm.Calibrate(ml_1)
ca_2.set_parameter(name="kaq0_1", initial=10, pmin=1e-8)
ca_2.set_parameter(name="Saq1", initial=1e-4, pmin=1e-5)
ca_2.set_parameter(name="Saq0", initial=0.2, pmin=1e-8)
ca_2.set_parameter(name="kzoverkh0_1", initial=1, pmin=1e-5)
ca_2.series(name="obs", x=r, y=0, t=td, h=hd, layer=0)  # deep piezometer
ca_2.fit()

.
.........................
.
.........................
.
..........................
..........................
..........................
.........
Fit succeeded.

/tmp/ipykernel_1462/994403479.py:4: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'kaq0_1' to silence this warning. The parameter name can still include layer info, but this will be ignored in a future version of TTim.
  ca_2.set_parameter(name="kaq0_1", initial=10, pmin=1e-8)
/tmp/ipykernel_1462/994403479.py:5: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'Saq1' to silence this warning. The parameter name can still include layer info, but this will be ignored in a future version of TTim.
  ca_2.set_parameter(name="Saq1", initial=1e-4, pmin=1e-5)
/tmp/ipykernel_1462/994403479.py:6: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'Saq0' to silence this warning. The parameter name can still include layer info, but this will be ignored in a future version of TTim.
  ca_2.set_parameter(name="Saq0", initial=0.2, pmin=1e-8)
/tmp/ipykernel_1462/994403479.py:7: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'kzoverkh0_1' to silence this warning. The parameter name can still include layer info, but this will be ignored in a future version of TTim.
  ca_2.set_parameter(name="kzoverkh0_1", initial=1, pmin=1e-5)

display(ca_2.parameters)
print("RMSE:", ca_2.rmse())
layers optimal std perc_std pmin pmax initial inhoms parray
kaq0_1_0_1 None 115.607587 5.249664 4.540934 1.000000e-08 inf 10.0000 None [[115.60758713535701, 115.60758713535701]]
Saq1_1_1 None 0.000010 0.000008 78.943286 1.000000e-05 inf 0.0001 None [[1.0108519128304572e-05]]
Saq0_0_0 None 0.000146 0.000156 107.158126 1.000000e-08 inf 0.2000 None [[0.0001460064170965003]]
kzoverkh0_1_0_1 None 13.678577 838.702739 6131.505601 1.000000e-05 inf 1.0000 None [[13.678577394965371, 13.678577394965371]] 
RMSE: 0.010081870854008614

hd_2 = ml_1.head(r, 0, td)
hs_2 = ml_1.head(r, 0, ts)
plt.semilogx(td, hd, ".", label="deep obs")
plt.semilogx(td, hd_2[0], label="deep ttim")
plt.semilogx(ts, hs, ".", label="shallow obs")
plt.semilogx(ts, hs_2[0], "--", label="shallow ttim")
plt.xlabel("time [d]")
plt.ylabel("drawdown [m]")
plt.title("TTim Unconfined Model Results - Deeper Piezometer")
plt.legend();
../_images/a04a5f16866ba55c87078cc1ce09d6d0f3fa45f04e1be76405f5d698b6e81726.png

Create a conceptual model with n-layers#

As we can see in the examples of step 5, the single-layer simplification does not represent the system well as we have a vertical component to flow, shown in the head difference between both piezometers.

We now explore the feature of TTim to create a multi-layer model to represent better the unconfined system and simulate the vertical flow component. We will discretize the aquifer in a 21 layer model, with 1 m thick each.

# Model Figure - One - layer model
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)


# Aquifer:
for i in range(21, -2, -1):
    ground = plt.Rectangle(
        (-20, -i),
        width=150,
        height=1,
        fc="w",
        zorder=0,
        alpha=0.9,
        hatch="..",
        ec="k",
        ls="--",
    )
    ax.add_patch(ground)


well = plt.Rectangle((-1.5, -21), width=3, height=23, fc="w", zorder=1, ec="k")
ax.add_patch(well)

# Wellhead
wellhead = plt.Rectangle((-2, 2), width=4, height=2.5, fc="w", zorder=2, ec="k")
ax.add_patch(wellhead)

# Screen for the well:
screen = plt.Rectangle(
    (-1.5, -21),
    width=3,
    height=11,
    fc="w",
    alpha=1,
    zorder=2,
    ec="k",
    ls="--",
    hatch="-",
)
screen.set_linewidth(2)
ax.add_patch(screen)

# Piezometers
piez1 = plt.Rectangle((89, -21), width=2, height=23, fc="w", zorder=1, ec="k")
screen_piez_1 = plt.Rectangle(
    (89, -19), width=2, height=7, fc="w", alpha=1, zorder=2, ec="k", ls="--", hatch="-"
)
screen_piez_1.set_linewidth(2)
screen_piez_2 = plt.Rectangle(
    (89, -3), width=2, height=0.5, fc="w", alpha=1, zorder=2, ec="k", ls="--", hatch="-"
)
screen_piez_2.set_linewidth(2)

ax.add_patch(piez1)
ax.add_patch(screen_piez_1)
ax.add_patch(screen_piez_2)


# last line
line = plt.Line2D(xdata=[-200, 1200], ydata=[2, 2], color="k")
ax.add_line(line)

# Water table
line2 = plt.Line2D(xdata=[-200, 1200], ydata=[0, 0], color="b")
ax.add_line(line2)

ax.text(-18, 0.5, s="Water Table", fontsize="large", color="b", bbox={"fc": "w"})
ax.text(93, -3, s="Shallow piezometer", bbox={"fc": "w"})
ax.text(93, -16, s="Deeper piezometer", bbox={"fc": "w"})

ax.set_xlim([-20, 130])
ax.set_ylim([-21, 7])
ax.set_xlabel("Distance [m]")
ax.set_ylabel("Relative height [m]")
ax.set_title("Conceptual Model 2 - Multi-layer Model 3D - Vennebulten Example");
../_images/912fd9489982b43a52baa2690ab6c4f75c7c6def61f14bee768979242fd7f28c.png 
nlay = 21  # number of layers
zlayers = np.linspace(0, b, nlay + 1)  # elevation of each layer
Saq = 1e-4 * np.ones(nlay)
Saq[0] = 0.1  # Setting the first storage as specific yield

The model is created just as in the previous step, however with the new parameters defined above:

ml_2 = ttm.Model3D(
    kaq=10, z=zlayers, Saq=Saq, kzoverkh=0.1, phreatictop=True, tmin=1e-4, tmax=1.1
)
w_2 = ttm.Well(ml_2, xw=0, yw=0, rw=0.1, tsandQ=[(0, Q)], layers=range(nlay))
ml_2.solve()

self.neq  21
solution complete

Calibrate multi-layer model with the two piezometers simultaneously#

In the TTim multi-layer model, we can fit the parameters using data from both piezometers simultaneously. For this initial assumption, we assume the aquifer has one hydraulic conductivity and storage parameter.

The unknown parameters are kaq, Saq, kzoverkh.

Now, on the series method, we have to remember to set a different layer for each piezometer, corresponding to the depth of the screen.

ca_3 = ttm.Calibrate(ml_2)
ca_3.set_parameter(name="kaq0_20", initial=10)
ca_3.set_parameter(name="Saq0", initial=0.2)
ca_3.set_parameter(name="Saq1_20", initial=1e-4)
ca_3.set_parameter(name="kzoverkh0_20", initial=0.1, pmin=1e-5, pmax=0.5)
ca_3.series(name="obs1", x=r, y=0, layer=1, t=ts, h=hs)
ca_3.series(name="obs2", x=r, y=0, layer=15, t=td, h=hd)
ca_3.fit(report=True)

......
.
......
.
......
.
......
.
......
.
......
.
......
.
......
.
......
.
......
.
......
.
......
.
......
.
......
.
Fit succeeded.
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 95
    # data points      = 48
    # variables        = 4
    chi-square         = 0.00242958
    reduced chi-square = 5.5218e-05
    Akaike info crit   = -466.779339
    Bayesian info crit = -459.294535
[[Variables]]
    kaq0_20_0_20:       56.6087015 +/- 3.65165543 (6.45%) (init = 10)
    Saq0_0_0:           0.03170812 +/- 0.00331430 (10.45%) (init = 0.2)
    Saq1_20_1_20:       3.4568e-05 +/- 2.0607e-06 (5.96%) (init = 0.0001)
    kzoverkh0_20_0_20:  0.00183709 +/- 4.5738e-04 (24.90%) (init = 0.1)
[[Correlations]] (unreported correlations are < 0.100)
    C(kaq0_20_0_20, kzoverkh0_20_0_20) = -0.9753
    C(kaq0_20_0_20, Saq0_0_0)          = -0.7921
    C(Saq0_0_0, kzoverkh0_20_0_20)     = +0.7509
    C(kaq0_20_0_20, Saq1_20_1_20)      = -0.4841
    C(Saq1_20_1_20, kzoverkh0_20_0_20) = +0.4175
    C(Saq0_0_0, Saq1_20_1_20)          = +0.3713

/tmp/ipykernel_1462/4244544214.py:2: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'kaq0_20' to silence this warning. The parameter name can still include layer info, but this will be ignored in a future version of TTim.
  ca_3.set_parameter(name="kaq0_20", initial=10)
/tmp/ipykernel_1462/4244544214.py:3: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'Saq0' to silence this warning. The parameter name can still include layer info, but this will be ignored in a future version of TTim.
  ca_3.set_parameter(name="Saq0", initial=0.2)
/tmp/ipykernel_1462/4244544214.py:4: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'Saq1_20' to silence this warning. The parameter name can still include layer info, but this will be ignored in a future version of TTim.
  ca_3.set_parameter(name="Saq1_20", initial=1e-4)
/tmp/ipykernel_1462/4244544214.py:5: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'kzoverkh0_20' to silence this warning. The parameter name can still include layer info, but this will be ignored in a future version of TTim.
  ca_3.set_parameter(name="kzoverkh0_20", initial=0.1, pmin=1e-5, pmax=0.5)

display(ca_3.parameters)
print("RMSE:", ca_3.rmse())
layers optimal std perc_std pmin pmax initial inhoms parray
kaq0_20_0_20 None 56.608702 3.651655 6.450696 -inf inf 10.0000 None [[56.60870151083327, 56.60870151083327, 56.608...
Saq0_0_0 None 0.031708 0.003314 10.452526 -inf inf 0.2000 None [[0.031708118804235]]
Saq1_20_1_20 None 0.000035 0.000002 5.961170 -inf inf 0.0001 None [[3.456801812921833e-05, 3.456801812921833e-05...
kzoverkh0_20_0_20 None 0.001837 0.000457 24.897285 0.00001 0.5 0.1000 None [[0.0018370874318220005, 0.0018370874318220005... 
RMSE: 0.007114515717856613

hs_3 = ml_2.head(x=r, y=0, t=ts, layers=1)
hd_3 = ml_2.head(x=r, y=0, t=td, layers=15)
plt.semilogx(ts, hs, ".", label="shallow obs")
plt.semilogx(td, hd, ".", label="deep obs")
plt.semilogx(ts, hs_3[0], label="shallow ttim")
plt.semilogx(td, hd_3[0], label="deep ttim")
plt.xlabel("time [d]")
plt.ylabel("drawdown [m]")
plt.title("TTim Multi - Layer Unconfined Model Results")
plt.legend()
plt.grid()
../_images/0ce6afee5dcdc797bf5d69e4b7561469edbe01823f02bc3b92b9a8230ec34883.png

We already see significant improvement from the previous single-layer model. The fit is better and the AIC and BIC indicators have also significantly improved.

What if we take into account the described stratification of the aquifer? In that case, we could try to stratify our model into two: The first 6 m and the deeper layers.

Calibration of the Stratified Model#

In this final example, we will assume the storage is distributed according to the sediment stratification in the aquifer. We will adjust two different Saqvalues, one for the first 6 m of the aquifer and another for the deeper layers. We assume the hydraulic conductivity and the anisotropy is constant.

# Model Figure - One - layer model
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)


# Aquifer 1:
for i in range(21, 6, -1):
    ground = plt.Rectangle(
        (-20, -i),
        width=150,
        height=1,
        fc="w",
        zorder=0,
        alpha=0.9,
        hatch="oo",
        ec="k",
        ls="--",
    )
    ax.add_patch(ground)
# Aquifer 2:
for i in range(6, -2, -1):
    ground = plt.Rectangle(
        (-20, -i),
        width=150,
        height=1,
        fc="w",
        zorder=0,
        alpha=0.9,
        hatch="..",
        ec="k",
        ls="--",
    )
    ax.add_patch(ground)


well = plt.Rectangle((-1.5, -21), width=3, height=23, fc="w", zorder=1, ec="k")
ax.add_patch(well)

# Wellhead
wellhead = plt.Rectangle((-2, 2), width=4, height=2.5, fc="w", zorder=2, ec="k")
ax.add_patch(wellhead)

# Screen for the well:
screen = plt.Rectangle(
    (-1.5, -21),
    width=3,
    height=11,
    fc="w",
    alpha=1,
    zorder=2,
    ec="k",
    ls="--",
    hatch="-",
)
screen.set_linewidth(2)
ax.add_patch(screen)

# Piezometers
piez1 = plt.Rectangle((89, -21), width=2, height=23, fc="w", zorder=1, ec="k")
screen_piez_1 = plt.Rectangle(
    (89, -19), width=2, height=7, fc="w", alpha=1, zorder=2, ec="k", ls="--", hatch="-"
)
screen_piez_1.set_linewidth(2)
screen_piez_2 = plt.Rectangle(
    (89, -3), width=2, height=0.5, fc="w", alpha=1, zorder=2, ec="k", ls="--", hatch="-"
)
screen_piez_2.set_linewidth(2)

ax.add_patch(piez1)
ax.add_patch(screen_piez_1)
ax.add_patch(screen_piez_2)


# last line
line = plt.Line2D(xdata=[-200, 1200], ydata=[2, 2], color="k")
ax.add_line(line)

# Water table
line2 = plt.Line2D(xdata=[-200, 1200], ydata=[0, 0], color="b")
ax.add_line(line2)

ax.text(-18, 0.5, s="Water Table", fontsize="large", color="b", bbox={"fc": "w"})
ax.text(93, -3, s="Shallow piezometer", bbox={"fc": "w"})
ax.text(93, -16, s="Deeper piezometer", bbox={"fc": "w"})

ax.set_xlim([-20, 130])
ax.set_ylim([-21, 7])
ax.set_xlabel("Distance [m]")
ax.set_ylabel("Relative height [m]")
ax.set_title("Conceptual Model 3 - Multi-layer Model 3D - Vennebulten Example");
../_images/d9942cab215ef97b4b0e5ce0fc21e73089f717d02bf541ec74dbd93f305cbb7b.png 
ml_3 = ttm.Model3D(
    kaq=10, z=zlayers, Saq=Saq, kzoverkh=0.1, phreatictop=True, tmin=1e-4, tmax=1.1
)
w_3 = ttm.Well(ml_3, xw=0, yw=0, rw=0.1, tsandQ=[(0, Q)], layers=range(nlay))
ml_3.solve()

self.neq  21
solution complete

ca_4 = ttm.Calibrate(ml_3)
ca_4.set_parameter(name="kaq0_20", initial=50)
ca_4.set_parameter(name="Saq0", initial=0.1)
ca_4.set_parameter(name="Saq1_7", initial=1e-4, pmin=0)
ca_4.set_parameter(name="Saq7_20", initial=1e-4, pmin=0)
ca_4.set_parameter(name="kzoverkh0_20", initial=0.1, pmin=1e-5, pmax=0.5)
ca_4.series(name="obs1", x=r, y=0, layer=1, t=ts, h=hs)
ca_4.series(name="obs2", x=r, y=0, layer=15, t=td, h=hd)
ca_4.fit()

......
.
......
.
......
.
......
.
.....
.
.....
.
.....
.
.....
.
.....
.
.....
.
......
.
......
.
......
.
......
.
......
.
......
.
.....
..
.....
..
.....
..
.....
..
....
Fit succeeded.

/tmp/ipykernel_1462/4291606139.py:2: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'kaq0_20' to silence this warning. The parameter name can still include layer info, but this will be ignored in a future version of TTim.
  ca_4.set_parameter(name="kaq0_20", initial=50)
/tmp/ipykernel_1462/4291606139.py:3: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'Saq0' to silence this warning. The parameter name can still include layer info, but this will be ignored in a future version of TTim.
  ca_4.set_parameter(name="Saq0", initial=0.1)
/tmp/ipykernel_1462/4291606139.py:4: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'Saq1_7' to silence this warning. The parameter name can still include layer info, but this will be ignored in a future version of TTim.
  ca_4.set_parameter(name="Saq1_7", initial=1e-4, pmin=0)
/tmp/ipykernel_1462/4291606139.py:5: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'Saq7_20' to silence this warning. The parameter name can still include layer info, but this will be ignored in a future version of TTim.
  ca_4.set_parameter(name="Saq7_20", initial=1e-4, pmin=0)
/tmp/ipykernel_1462/4291606139.py:6: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'kzoverkh0_20' to silence this warning. The parameter name can still include layer info, but this will be ignored in a future version of TTim.
  ca_4.set_parameter(name="kzoverkh0_20", initial=0.1, pmin=1e-5, pmax=0.5)

display(ca_4.parameters)
print("RMSE:", ca_4.rmse())
layers optimal std perc_std pmin pmax initial inhoms parray
kaq0_20_0_20 None 73.422023 2.928747 3.988921 -inf inf 50.0000 None [[73.42202304879444, 73.42202304879444, 73.422...
Saq0_0_0 None 0.021280 0.002075 9.749326 -inf inf 0.1000 None [[0.021279910580349958]]
Saq1_7_1_7 None 0.000420 0.000060 14.405552 0.00000 inf 0.0001 None [[0.0004195699808910991, 0.0004195699808910991...
Saq7_20_7_20 None 0.000024 0.000001 4.946894 0.00000 inf 0.0001 None [[2.365855114194737e-05, 2.365855114194737e-05...
kzoverkh0_20_0_20 None 0.000433 0.000122 28.189746 0.00001 0.5 0.1000 None [[0.0004334769044372613, 0.0004334769044372613... 
RMSE: 0.003437332288521483

hs_4 = ml_3.head(x=r, y=0, t=ts, layers=1)
hd_4 = ml_3.head(x=r, y=0, t=td, layers=15)
plt.semilogx(ts, hs, ".", label="shallow obs")
plt.semilogx(td, hd, ".", label="deep obs")
plt.semilogx(ts, hs_4[0], label="shallow ttim")
plt.semilogx(td, hd_4[0], label="deep ttim")
plt.xlabel("time [d]")
plt.ylabel("drawdown [m]")
plt.title("TTim Stratified Unconfined Model Results")
plt.legend()
plt.grid()
../_images/09357bb31fc8e06e96de04c461eeb507bce755e9b12d79d64b6251245fcdab79.png

Here we see that the model fit has significantly improved. AIC and BIC indicators are lower than the previous multi-layer model.

Comparison of results#

t2 = pd.DataFrame(
    columns=["k [m/d]", "Sy [-]", "Ss [1/m]", "kzoverkh"],
    index=["MLU", "ttim-multilayer", "ttim-stratified Ss"],
)
t2.loc["MLU"] = [62.657, 0.0012, 2.790e-05, 0.002595]
t2.loc["ttim-multilayer"] = ca_3.parameters["optimal"].values
t2.iloc[2, 0:2] = ca_4.parameters["optimal"].values[0:2]
t2.iloc[2, 2:4] = ca_4.parameters["optimal"].values[3:5]
t2.loc[:, "RMSE"] = pd.Series([0.013540, ca_3.rmse(), ca_4.rmse()], index=t2.index)
t2
k [m/d] Sy [-] Ss [1/m] kzoverkh RMSE
MLU 62.657 0.0012 0.000028 0.002595 0.013540
ttim-multilayer 56.608702 0.031708 0.000035 0.001837 0.007115
ttim-stratified Ss 73.422023 0.02128 0.000024 0.000433 0.003437

The multi-layer approach allowed us to fit both piezometers and better represent the vertical component of flow. However, the parameters were sensitive to the conceptualization applied. The stratified model had much larger hydraulic conductivity in comparison to the multi-layer model. In the stratified approach, the fit has significantly improved.

References#

  • Carlson F, Randall J (2012) MLU: a Windows application for the analysis of aquifer tests and the design of well fields in layered systems. Ground Water 50(4):504–510

  • Duffield, G.M., 2007. AQTESOLV for Windows Version 4.5 User’s Guide, HydroSOLVE, Inc., Reston, VA.

  • Kruseman, G.P., De Ridder, N.A., Verweij, J.M., 1970. Analysis and evaluationof pumping test data. volume 11. International institute for land reclamation and improvement The Netherlands.

  • Neuman, S.P., Witherspoon, P.A., 1969. Applicability of current theories of flow in leaky aquifers. Water Resources Research 5, 817–829.

  • Yang, Xinzhu (2020) Application and comparison of different methodsfor aquifer test analysis using TTim. Master Thesis, Delft University of Technology (TUDelft), Delft, The Netherlands.