4. Confined Aquifer Test - Schroth#

This test is taken from examples presented in MLU tutorial.

Step 1. Import Libraries#

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

import ttim

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

This test example was a pumping test conducted near San Francisco, US, and reported by Brian et al. (1997). The system consists of a confined system of two aquifers separated by an aquitard layer. The upper aquifer layer is located from 46 to 49 m depth, followed by an aquitard layer from 49 to 52 m depth and the second aquifer at 52 to 55 m depth.

The lower aquifer is pumped by a well, named EW-712 that fully penetrates the formation. An observation well is placed 46 m away from the well, in the lower aquifer formation, and it is named MW-616. The last observation well is placed in the upper aquifer right on top of the pumping well. However, data was not available for this well. The radius of all wells was 0.05 m.

All wells before pumping had water levels around 20 m depth, which means that the system can be characterized as confined.

The time-drawdown data for the observation well MW-616 and pumping well EW-712 was obtained from MLU documentation (Duffield, 2007).

In this example, we are reproducing the results obtained by Yang (2020). We will use TTim to test two hypotheses: the first is that the lower aquifer is, by confinement, disconnected to the upper aquifer, and the second is there is enough leakage from the upper aquifer to consider a leaky aquifer relation between them.

A simplified cross-section of the model area can be seen below:

Hide code cell source

 
import matplotlib.pyplot as plt

##Now printing the conceptual model figure:

fig = plt.figure(figsize=(14, 9))
ax = fig.add_subplot(1, 1, 1)
# sky
sky = plt.Rectangle((-20, 0), width=100, height=3, fc="b", zorder=0, alpha=0.1)
ax.add_patch(sky)

# Aquifer 1:
ground = plt.Rectangle(
    (-20, -55),
    width=100,
    height=3,
    fc=np.array([209, 179, 127]) / 255,
    zorder=0,
    alpha=0.9,
)
ax.add_patch(ground)

# Aquifer 2:

ground = plt.Rectangle(
    (-20, -49),
    width=100,
    height=3,
    fc=np.array([209, 179, 127]) / 255,
    zorder=0,
    alpha=0.9,
)
ax.add_patch(ground)

# Confining bed:
confining_unit = plt.Rectangle(
    (-20, -52),
    width=100,
    height=3,
    fc=np.array([100, 100, 100]) / 255,
    zorder=0,
    alpha=0.7,
)
ax.add_patch(confining_unit)

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

# Wellhead
wellhead = plt.Rectangle(
    (-2.5, 0), width=5, height=1.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, -55),
    width=3,
    height=3,
    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=5, y=0.75, dx=3, dy=0, color="#00035b")
ax.add_patch(pumping_arrow)
ax.text(x=9, y=0.75, s=r"$ Q = 220 \frac{gal}{min}$", fontsize="large")
# Piezometers
piez1 = plt.Rectangle(
    (45, -55), width=2, height=55, fc=np.array([200, 200, 200]) / 255, zorder=1
)
screen_piez_1 = plt.Rectangle(
    (45, -55),
    width=2,
    height=3,
    fc=np.array([200, 200, 200]) / 255,
    alpha=1,
    zorder=2,
    ec="k",
    ls="--",
)
screen_piez_1.set_linewidth(2)

ax.add_patch(piez1)

ax.add_patch(screen_piez_1)

# last line
line = plt.Line2D(xdata=[-20, 100], ydata=[0, 0], color="k")
ax.add_line(line)
ax.text(x=45, y=0.75, s="Obs Well 1", fontsize="large")

ax.text(
    x=-17,
    y=-40,
    s="Upper Formations",
    fontsize="large",
    bbox={"facecolor": "w", "alpha": 0.9},
)
ax.text(x=-17, y=-48, s="Upper Aquifer", fontsize="large")
ax.text(x=-17, y=-51, s="Aquitard", fontsize="large")
ax.text(x=-17, y=-54, s="Lower Aquifer", fontsize="large")
# Complete the figure:

upper_formations = plt.Rectangle(
    (-20, -46), width=100, height=46, fc="white", hatch="-/", zorder=0, alpha=0.9
)
ax.add_patch(upper_formations)

ax.set_xlim([-20, 80])
ax.set_ylim([-55, 3])
ax.set_xlabel("Distance [m]")
ax.set_ylabel("Relative height [m]")
ax.set_title("Conceptual Model - Schroth Example");
../_images/d2f755f643cff31ba329e6b388305e5caeddd6d5d556842348bb2c9fadb69ec3.png

Step 2. Set basic parameters#

Q = 82.08  # constant discharge in m^3/d
zt0 = -46  # top boundary of upper aquifer in m
zb0 = -49  # bottom boundary of upper aquifer in m
zt1 = -52  # top boundary of lower aquifer in m
zb1 = -55  # bottom boundary of lower aquifer in m
rw = 0.05  # well radius in m

Step 3. Load data of the pumping and observation well#

The preferred method of loading data into TTim is to use numpy arrays.

The data is in a text file where the first column is the time data in days and the second column is the drawdown in meters

The observation well is referred to as Well 1 and the pumping well as Well 3.

For each piezometer, we will load the data as a numpy array. We further split the data into two different 1d arrays, one for time and another for drawdown. Finally, we convert the time data from minutes to days

# Loading data for the pumping well
data1 = np.loadtxt("data/schroth_obs1.txt", skiprows=1)
t1 = data1[:, 0]
h1 = data1[:, 1]
r1 = 0  # Pumping well is at distance 0 to pumping

# Loading data for the observation well
data2 = np.loadtxt("data/schroth_obs2.txt", skiprows=1)
t2 = data2[:, 0]
h2 = data2[:, 1]
r2 = 46  # distance between observation well2 and pumping well

Step 4. Create TTim model#

We begin by considering the underlying aquifer as a single confined aquifer (overlying aquifer and aquitard are excluded).

In this example, we are using the ModelMaq model to conceptualize our aquifer. ModelMaq defines the aquifer system as a stacked vertical sequence of aquifers and leaky layers (aquifer-leaky layer, aquifer-leaky layer, etc). A thorough explanation of the ModelMaq and TTim one-layer modelling conceptualization is given in the notebook: [Confined 1 - Oude Korendijk](confined1_oude_korendijk

ml_0 = ttim.ModelMaq(z=[zt1, zb1], kaq=10, Saq=1e-4, tmin=1e-4, tmax=1)
w_0 = ttim.Well(ml_0, xw=0, yw=0, rw=rw, tsandQ=[(0, Q), (1e08, 0)])
ml_0.solve()

self.neq  1
solution complete

Step 5. Calibration#

The calibration workflow has been described in detail in the notebook: Confined 1 - Oude Korendijk

ca_0 = ttim.Calibrate(ml_0)
ca_0.set_parameter(name="kaq0", initial=10)
ca_0.set_parameter(name="Saq0", initial=1e-4)
ca_0.series(name="well", x=r1, y=0, t=t1, h=h1, layer=0)
ca_0.series(name="obs_well", x=r2, y=0, t=t2, h=h2, layer=0)
ca_0.fit(report=True)

...................................
......................................
.....................................
......
Fit succeeded.
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 113
    # data points      = 40
    # variables        = 2
    chi-square         = 111.249415
    reduced chi-square = 2.92761620
    Akaike info crit   = 44.9158085
    Bayesian info crit = 48.2935674
[[Variables]]
    kaq0_0_0:  1.03182460 +/- 0.10473289 (10.15%) (init = 10)
    Saq0_0_0:  0.04018095 +/- 0.02031336 (50.55%) (init = 0.0001)
[[Correlations]] (unreported correlations are < 0.100)
    C(kaq0_0_0, Saq0_0_0) = -0.9309

/tmp/ipykernel_1192/1435324636.py:2: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'kaq0' to silence this warning. The parameter name can still include layer info, but this will be ignored in a future version of TTim.
  ca_0.set_parameter(name="kaq0", initial=10)
/tmp/ipykernel_1192/1435324636.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_0.set_parameter(name="Saq0", initial=1e-4)

display(ca_0.parameters)
print("RMSE:", ca_0.rmse())
layers optimal std perc_std pmin pmax initial inhoms parray
kaq0_0_0 None 1.031825 0.104733 10.150261 -inf inf 10.0000 None [[1.03182460080969]]
Saq0_0_0 None 0.040181 0.020313 50.554706 -inf inf 0.0001 None [[0.04018094856463222]] 
RMSE: 1.6677036263520377

Now let’s plot the model with our observation data:

hm1_0 = ml_0.head(r1, 0, t1)
hm2_0 = ml_0.head(r2, 0, t2)
plt.figure(figsize=(8, 5))
plt.semilogx(t1, h1, ".", label="well")
plt.semilogx(t2, h2, ".", label="obs_well")
plt.semilogx(t1, hm1_0[-1], label="ttim model - well")
plt.semilogx(t2, hm2_0[-1], label="ttim model - obs_well")
plt.xlabel("time [d]")
plt.ylabel("head [m]")
plt.legend();
../_images/f3d25adbcdc69fdd2f8e46b09dd356053ff5fa6631a43fcca7fcaa87e3bd47c7.png

The figure shows a poor fit specially for the well. Probably well effects might be relevant in this case, so we will try to calibrate them next.

Step 5.2. Model Calibration with skin resistance and wellbore storage#

To improve the fit, we will try to calibrate the model adding wellbore storage and skin resistance to the pumping well.

For this, we create a new model and add two extra parameters to the Well object:

  • The radius of the caisson rc, which we use to simulate wellbore storage. In this case, we use a value in meters (float);

  • The skin resistance res, a float value with the unit in days.

ml_1 = ttim.ModelMaq(z=[zt1, zb1], kaq=10, Saq=1e-4, tmin=1e-4, tmax=1)
w_1 = ttim.Well(ml_1, xw=0, yw=0, rw=rw, rc=0, res=5, tsandQ=[(0, Q), (1e08, 0)])
ml_1.solve()

self.neq  1
solution complete

Here we use the method .set_parameter_by_reference to calibrate the rc and res parameters in our well.

.set_parameter_by_reference takes the following arguments:

  • name: a string of the parameter name

  • parameter: numpy-array with the parameter to be optimized. It should be specified as a reference, for example, in our case: w1.rc[0:] referencing to the parameter rc in object w1.

  • initial: float with the initial guess for the parameter value.

  • pmin and pmax: floats with the minimum and maximum values allowed. If not specified, these will be -np.inf and np.inf.

ca_1 = ttim.Calibrate(ml_1)
ca_1.set_parameter(name="kaq0", initial=10)
ca_1.set_parameter(name="Saq0", initial=1e-4)
ca_1.set_parameter_by_reference(name="rc", parameter=w_1.rc[:], initial=0.2)
ca_1.set_parameter_by_reference(name="res", parameter=w_1.res[:], initial=3)
ca_1.series(name="well", x=r1, y=0, t=t1, h=h1, layer=0)
ca_1.series(name="obs_well", x=r2, y=0, t=t2, h=h2, layer=0)
ca_1.fit(report=True)

.................................
....
.................................
.....
................................
......
...............................
.......
.............................
.........
.............................
.........
............................
.........
.............................
.........
............................
.........
............................
........
.............................
........
.............................
........
................
Fit succeeded.
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 462
    # data points      = 40
    # variables        = 4
    chi-square         = 13.6497878
    reduced chi-square = 0.37916077
    Akaike info crit   = -35.0062192
    Bayesian info crit = -28.2507013
[[Variables]]
    kaq0_0_0:  1.95210264 +/- 0.05268909 (2.70%) (init = 10)
    Saq0_0_0:  1.1463e-04 +/- 3.3125e-05 (28.90%) (init = 0.0001)
    rc:        0.00289266 +/- 0.02556637 (883.84%) (init = 0.2)
    res:       31.8297934 +/- 564.722578 (1774.19%) (init = 3)
[[Correlations]] (unreported correlations are < 0.100)
    C(rc, res)            = -1.0000
    C(kaq0_0_0, Saq0_0_0) = -0.8643
    C(Saq0_0_0, rc)       = -0.1061
    C(Saq0_0_0, res)      = +0.1048

/tmp/ipykernel_1192/1453838537.py:2: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'kaq0' 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", initial=10)
/tmp/ipykernel_1192/1453838537.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_1.set_parameter(name="Saq0", initial=1e-4)

display(ca_1.parameters)
print("RMSE:", ca_1.rmse())
layers optimal std perc_std pmin pmax initial inhoms parray
kaq0_0_0 None 1.952103 0.052689 2.699094 -inf inf 10.0000 None [[1.9521026384548064]]
Saq0_0_0 None 0.000115 0.000033 28.897788 -inf inf 0.0001 None [[0.00011462817288296093]]
rc NaN 0.002893 0.025566 883.837522 -inf inf 0.2000 NaN [[0.002892655531761506]]
res NaN 31.829793 564.722578 1774.194921 -inf inf 3.0000 NaN [[31.82979339345705]] 
RMSE: 0.5841615314746983

hm1_1 = ml_1.head(r1, 0, t1)
hm2_1 = ml_1.head(r2, 0, t2)
plt.figure(figsize=(10, 7))
plt.semilogx(t1, h1, ".", label="well")
plt.semilogx(t2, h2, ".", label="obs_well")
plt.semilogx(t1, hm1_1[-1], label="ttim model - well")
plt.semilogx(t2, hm2_1[-1], label="ttim model - obs_well")
plt.xlabel("time [d]")
plt.ylabel("drawdown [m]")
plt.legend()
plt.title("Model Results - One Aquifer + Well Parameters");
../_images/7899ca6416a9489988127c1b865c45304a0120f317f81bd9c251c22244094514.png

The model has improved significantly. However, we can visually see that the fit is not very good for late time data. Let us investigate if we can do better with a three-layer model such as the one defined in our conceptual model.

Step 6. Create a two-aquifer conceptual model#

Until so far, we have only considered horizontal flow in our TTim models. The assumption is sufficient in fully-penetrated wells in confined aquifers. However, it might not represent the situation so well in a more complex scenario, where the vertical flow is a relevant component of the groundwater flow. In TTim, this is done by assigning more layers to the model.

An advantage of TTim is the ability to create a model with more than one aquifer layer. We will explore this feature under the ModelMaq model as in Step4.

When we use ModelMaq, the leaky layers are located in between the aquifers. These leaky layers only have vertical flow and are characterized by the parameters resistance to vertical flow (c) and storage (Sll). The specific flux is computed as (Bakker, 2013):

\[q_n = \frac{h_n-h_{n-1}}{c_n}\]

where \(q_n\) is the vertical flux from layer \(n\) to layer \(n-1\), \(h_n\) is the head in layer \(n\) and \(c_n\) is the vertical resistance to flow. \(c_n\) is computed as: \(H_n/k_n\) where, \(H_n\) is the leaky-layer thickness and \(k_n\) the vertical hydraulic conductance. \(c_n\) is the inverse of the parameter Leakance (\(L_n = 1/c_n\)), that is used in MODFLOW (Harbaugh, 2005) or analytical solutions of leaky-layers, such as in Hantush (1955).

Alternatively, we can also model the interface of two aquifers that are not separated by a leaky layer. In that case, the leaky layer is set to 0 m thickness, and the model calculates the resistance to vertical flow by a finite differences scheme (Bakker, 2013).

In the model construction, we have to set the parameters for each layer (which consists of an aquifer layer and an aquitard layer).

For our two-layer 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].

  • 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 aquitard layers are defined by the space between the aquifer layers. The top of the first aquitard is the bottom of the first aquifer, the bottom of the first aquitard is the top of the second aquifer, and so on. But they can also have 0 thickness. In case the parameter topboundary is set to semi, where we assume we have semi-confined conditions, we have to add one more element to the beginning of the list, which is the top of the aquitard layer that is above of the first aquifer.

  • The specific storage: Saq. It is a list/array with a float element for every aquifer, for example: [Saq0, Saq1].

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

  • And the maximum time: tmax, float.

  • topboundary is the parameter that sets whether the upper layer is a confined unit or a semi-confined unit. We can assign this parameter as:

    • 'conf': This means that the upper layer is sealed from above in a confined condition; *'semi': This means that the upper layer receives leakage from phreatic storage (phreatictop = True) or from a fixed head above the upper leaky-layer (phreatictop = False).

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

  • The resistance to vertical flow: c of the aquitard layers, which is a list with length n-1 where n is the number of aquifer layers, setting the resistance of each aquitard layer. In the case of semi-confined conditions (topboundary = " semi"), we also add a first element representing the resistance of the aquitard above the first aquifer.

  • The storage Sll of the aquitard layers: float/list/array with the specific storage values for each aquitard layer. In case a float is defined, the same storage is assumed for every layer. In case topboundary = semi and phreatictop = True, the first element of Sll is the specific yield (see example Unconfined - 1 - Vennebulten).

  • phreatictop: Is a boolean (True/False). If True, the first element in Saq is considered phreatic storage (Specific Yield), and it is not multiplied by the layer thickness. The default value is False in ModelMaq. This parameter is normally True only in unconfined aquifers.

ml_2 = ttim.ModelMaq(
    kaq=[17.28, 2],
    z=[zt0, zb0, zt1, zb1],
    c=200,
    Saq=[1.2e-4, 1e-5],
    Sll=3e-5,
    topboundary="conf",
    tmin=1e-4,
    tmax=0.5,
)
w_2 = ttim.Well(ml_2, xw=0, yw=0, rw=rw, tsandQ=[(0, Q), (1e08, 0)], layers=1)
ml_2.solve()

self.neq  1
solution complete

Step 7. Calibrate Multi-Aquifer Model#

Now we follow the procedures in step 5 again, but calibrating the parameters of both aquifers.

Step 7.1. Calibrate model without wellstorage and skin resistance#

Complementing the prior calibration, we add the hydraulic parameters of both layers 0 and 1 and the parameters of the aquitard layer in between. The procedure is the same as explained in Step 5.

ca_2 = ttim.Calibrate(ml_2)
ca_2.set_parameter(name="kaq0", initial=20, pmin=0, pmax=200)
ca_2.set_parameter(name="kaq1", initial=1, pmin=0)
ca_2.set_parameter(name="Saq0", initial=1e-4, pmin=0)
ca_2.set_parameter(name="Saq1", initial=1e-5, pmin=0)
ca_2.set_parameter_by_reference(
    name="Sll", parameter=ml_2.aq.Sll[:], initial=1e-4, pmin=0
)
ca_2.set_parameter(name="c1", initial=100, pmin=0)
ca_2.series(name="well", x=r1, y=0, t=t1, h=h1, layer=1)
ca_2.series(name="obs_well", x=r2, y=0, t=t2, h=h2, layer=1)
ca_2.fit(report=True)

............................
...
.............................
...
............................
...
............................
...
............................
...
............................
...
............................
...
............................
...
.............................
..
.............................
..
..............................
.
..............................
.
..............................
.
..............................
.
..............................
.
..............................
.
..............................
.
..............................
.
..............................
.
..............................
.
..............................
.
..............................
.
..............................
.
..............................
.
..............................
.
..............................
.
..............................
.
..............................
.
..............................
.
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
................................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...............................
...................

/tmp/ipykernel_1192/3079571136.py:2: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'kaq0' 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", initial=20, pmin=0, pmax=200)
/tmp/ipykernel_1192/3079571136.py:3: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'kaq1' 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="kaq1", initial=1, pmin=0)
/tmp/ipykernel_1192/3079571136.py:4: 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=1e-4, pmin=0)
/tmp/ipykernel_1192/3079571136.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-5, pmin=0)
/tmp/ipykernel_1192/3079571136.py:9: DeprecationWarning: Setting layers in the parameter name is deprecated. Set the layers= keyword argument for parameter 'c1' 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="c1", initial=100, pmin=0)

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[14], line 12
     10 ca_2.series(name="well", x=r1, y=0, t=t1, h=h1, layer=1)
     11 ca_2.series(name="obs_well", x=r2, y=0, t=t2, h=h2, layer=1)
---> 12 ca_2.fit(report=True)

File ~/checkouts/readthedocs.org/user_builds/ttim/envs/stable/lib/python3.11/site-packages/ttim/fit.py:378, in Calibrate.fit(self, report, printdot, **kwargs)
    375 def fit(self, report=False, printdot=True, **kwargs):
    376     # current default fitting routine is lmfit
    377     # return self.fit_least_squares(report) # does not support bounds by default
--> 378     return self.fit_lmfit(report, printdot, **kwargs)

File ~/checkouts/readthedocs.org/user_builds/ttim/envs/stable/lib/python3.11/site-packages/ttim/fit.py:319, in Calibrate.fit_lmfit(self, report, printdot, **kwargs)
    317     self.lmfitparams.add(name, value=p["initial"], min=p["pmin"], max=p["pmax"])
    318 fit_kws = {"epsfcn": 1e-4}  # this is essential to specify step for the Jacobian
--> 319 self.fitresult = lmfit.minimize(
    320     self.residuals_lmfit,
    321     self.lmfitparams,
    322     method="leastsq",
    323     kws={"printdot": printdot},
    324     **fit_kws,
    325     **kwargs,
    326 )
    327 print("", flush=True)
    328 print(self.fitresult.message)

File ~/checkouts/readthedocs.org/user_builds/ttim/envs/stable/lib/python3.11/site-packages/lmfit/minimizer.py:2616, in minimize(fcn, params, method, args, kws, iter_cb, scale_covar, nan_policy, reduce_fcn, calc_covar, max_nfev, **fit_kws)
   2476 """Perform the minimization of the objective function.
   2477 
   2478 The minimize function takes an objective function to be minimized,
   (...)   2610 
   2611 """
   2612 fitter = Minimizer(fcn, params, fcn_args=args, fcn_kws=kws,
   2613                    iter_cb=iter_cb, scale_covar=scale_covar,
   2614                    nan_policy=nan_policy, reduce_fcn=reduce_fcn,
   2615                    calc_covar=calc_covar, max_nfev=max_nfev, **fit_kws)
-> 2616 return fitter.minimize(method=method)

File ~/checkouts/readthedocs.org/user_builds/ttim/envs/stable/lib/python3.11/site-packages/lmfit/minimizer.py:2355, in Minimizer.minimize(self, method, params, **kws)
   2352         if (key.lower().startswith(user_method) or
   2353                 val.lower().startswith(user_method)):
   2354             kwargs['method'] = val
-> 2355 return function(**kwargs)

File ~/checkouts/readthedocs.org/user_builds/ttim/envs/stable/lib/python3.11/site-packages/lmfit/minimizer.py:1674, in Minimizer.leastsq(self, params, max_nfev, **kws)
   1672 result.call_kws = lskws
   1673 try:
-> 1674     lsout = scipy_leastsq(self.__residual, variables, **lskws)
   1675 except AbortFitException:
   1676     pass

File ~/checkouts/readthedocs.org/user_builds/ttim/envs/stable/lib/python3.11/site-packages/scipy/optimize/_minpack_py.py:439, in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    437     if maxfev == 0:
    438         maxfev = 200*(n + 1)
--> 439     retval = _minpack._lmdif(func, x0, args, full_output, ftol, xtol,
    440                              gtol, maxfev, epsfcn, factor, diag)
    441 else:
    442     if col_deriv:

File ~/checkouts/readthedocs.org/user_builds/ttim/envs/stable/lib/python3.11/site-packages/lmfit/minimizer.py:540, in Minimizer.__residual(self, fvars, apply_bounds_transformation)
    537     self.result.success = False
    538     raise AbortFitException(f"fit aborted: too many function evaluations {self.max_nfev}")
--> 540 out = self.userfcn(params, *self.userargs, **self.userkws)
    542 if callable(self.iter_cb):
    543     abort = self.iter_cb(params, self.result.nfev, out,
    544                          *self.userargs, **self.userkws)

File ~/checkouts/readthedocs.org/user_builds/ttim/envs/stable/lib/python3.11/site-packages/ttim/fit.py:278, in Calibrate.residuals_lmfit(self, lmfitparams, printdot)
    276 vals = lmfitparams.valuesdict()
    277 p = np.array([vals[k] for k in self.parameters.index])
--> 278 return self.residuals(p, printdot)

File ~/checkouts/readthedocs.org/user_builds/ttim/envs/stable/lib/python3.11/site-packages/ttim/fit.py:254, in Calibrate.residuals(self, p, printdot, weighted, layers, series)
    251     layers = range(self.model.aq.naq)
    253 for i, k in enumerate(self.parameters.index):
--> 254     parraylist = self.parameters.loc[k, "parray"]
    255     for parray in parraylist:
    256         # [:] needed to do set value in array
    257         parray[:] = p[i]

File ~/checkouts/readthedocs.org/user_builds/ttim/envs/stable/lib/python3.11/site-packages/pandas/core/indexing.py:1194, in _LocationIndexer.__getitem__(self, key)
   1192 @final
   1193 def __getitem__(self, key):
-> 1194     check_dict_or_set_indexers(key)
   1195     if type(key) is tuple:
   1196         key = (list(x) if is_iterator(x) else x for x in key)

File ~/checkouts/readthedocs.org/user_builds/ttim/envs/stable/lib/python3.11/site-packages/pandas/core/indexing.py:2785, in check_dict_or_set_indexers(key)
   2780 def check_dict_or_set_indexers(key) -> None:
   2781     """
   2782     Check if the indexer is or contains a dict or set, which is no longer allowed.
   2783     """
   2784     if isinstance(key, set) or (
-> 2785         isinstance(key, tuple) and any(isinstance(x, set) for x in key)
   2786     ):
   2787         raise TypeError(
   2788             "Passing a set as an indexer is not supported. Use a list instead."
   2789         )
   2791     if isinstance(key, dict) or (
   2792         isinstance(key, tuple) and any(isinstance(x, dict) for x in key)
   2793     ):

File ~/checkouts/readthedocs.org/user_builds/ttim/envs/stable/lib/python3.11/site-packages/pandas/core/indexing.py:2785, in <genexpr>(.0)
   2780 def check_dict_or_set_indexers(key) -> None:
   2781     """
   2782     Check if the indexer is or contains a dict or set, which is no longer allowed.
   2783     """
   2784     if isinstance(key, set) or (
-> 2785         isinstance(key, tuple) and any(isinstance(x, set) for x in key)
   2786     ):
   2787         raise TypeError(
   2788             "Passing a set as an indexer is not supported. Use a list instead."
   2789         )
   2791     if isinstance(key, dict) or (
   2792         isinstance(key, tuple) and any(isinstance(x, dict) for x in key)
   2793     ):

KeyboardInterrupt: 

display(ca_2.parameters)
print("RMSE:", ca_2.rmse())
optimal std perc_std pmin pmax initial parray
kaq0 199.161550 13053.628430 6554.291438 0 200.0 20.00000 [199.1615501616947]
kaq1 0.097464 0.323649 332.068998 0 inf 1.00000 [0.09746438855592698]
Saq0 5.610347 2477.773317 44164.354256 0 inf 0.00010 [5.610346530366302]
Saq1 1.737358 5.461490 314.355977 0 inf 0.00001 [1.737358478821435]
Sll 0.063543 15.739979 24770.624116 0 inf 0.00010 [0.06354292456986621, 0.06354292456986621]
c1 0.001425 0.006811 477.970244 0 inf 100.00000 [0.0014250704284297644] 
RMSE: 0.744305083413874

The more complex model has improved the fit significantly from the first model. However, it showed worse RMSE, AIC and BIC values than the second model, with wellbore storage and skin resistance. We also have to note the unrealistic values found for the hydraulic conductivity of the upper layer and the storage parameters.

Next, we plot this model results.

hm1_2 = ml_2.head(r1, 0, t1)
hm2_2 = ml_2.head(r2, 0, t2)
plt.figure(figsize=(8, 5))
plt.semilogx(t1, h1, ".", label="well")
plt.semilogx(t2, h2, ".", label="obs_well")
plt.semilogx(t1, hm1_2[-1], label="ttim model - well")
plt.semilogx(t2, hm2_2[-1], label="ttim model - obs_well")
plt.xlabel("time [d]")
plt.ylabel("drawdown [m]")
plt.legend()
plt.title("Model Results - Two Aquifer Model");

Text(0.5, 0.98, 'Model Results - Two Aquifer Model')
../_images/8441354485f47594b6f337c989ccf2b31dfecdd277a654d4d6442ee7208d71f1.png

Step 7.2. Calibrate the two-layer model with wellparameters#

We will add the skinresistance and wellbore storage parameters to our well and resume calibration.

ml_3 = ttim.ModelMaq(
    kaq=[19, 2],
    z=[zt0, zb0, zt1, zb1],
    c=200,
    Saq=[4e-4, 1e-5],
    Sll=1e-4,
    topboundary="conf",
    tmin=1e-4,
    tmax=0.5,
)
w_3 = ttim.Well(
    ml_3, xw=0, yw=0, rw=rw, rc=0.2, res=0, tsandQ=[(0, Q), (1e08, 0)], layers=1
)
ml_3.solve()

ca_3 = ttim.Calibrate(ml_3)
ca_3.set_parameter(name="kaq0", initial=20, pmin=0)
ca_3.set_parameter(name="kaq1", initial=1, pmin=0)
ca_3.set_parameter(name="Saq0", initial=1e-4, pmin=1e-7)
ca_3.set_parameter(name="Saq1", initial=1e-5, pmin=1e-7)
ca_3.set_parameter_by_reference(
    name="Sll", parameter=ml_3.aq.Sll[:], initial=1e-4, pmin=1e-7
)
ca_3.set_parameter(name="c1", initial=100, pmin=0)
ca_3.set_parameter_by_reference(name="res", parameter=w_3.res[:], initial=0, pmin=0)
ca_3.set_parameter_by_reference(name="rc", parameter=w_3.rc[:], initial=0.2, pmin=0)
ca_3.series(name="obs1", x=r1, y=0, t=t1, h=h1, layer=1)
ca_3.series(name="obs2", x=r2, y=0, t=t2, h=h2, layer=1)
ca_3.fit(report=True)

display(ca_3.parameters)
print("RMSE:", ca_3.rmse())

Now we see a much better fit and AIC BIC values from the previous simulation. One thing to consider is the tiny storage values of the first aquifer and the aquitard. It might mean that it is troubling to fit all these parameters together. The small result of the skin resistance and its very high standard deviation might also mean it is irrelevant.

hm1_3 = ml_3.head(r1, 0, t1)
hm2_3 = ml_3.head(r2, 0, t2)
plt.figure(figsize=(8, 5))
plt.semilogx(t1, h1, ".", label="well")
plt.semilogx(t2, h2, ".", label="obs_well")
plt.semilogx(t1, hm1_3[-1], label="ttim model - well")
plt.semilogx(t2, hm2_3[-1], label="ttim model - observation well")
plt.xlabel("time [d]")
plt.ylabel("drawdown [m]")
plt.title("Model Results - Two Aquifers + Well Parameters")
plt.legend();

The plot showed the significant fit improvement with the new model.

Step 7.3. Calibrate Model with some fixed parameters#

In this new model, we will fix some of the parameters. Thus, we expect to have better confidence in estimating the others.

We will fix the hydraulic parameters of the first aquifer layer, the storage coefficient of the leaky layer and the skin resistance of the well. The fixed parameters were defined from the estimates of Brian et al. (1997).

ml_4 = ttim.ModelMaq(
    kaq=[17.28, 2],
    z=[zt0, zb0, zt1, zb1],
    c=200,
    Saq=[1.2e-4, 1e-5],
    Sll=3e-5,
    topboundary="conf",
    tmin=1e-4,
    tmax=0.5,
)
w_4 = ttim.Well(
    ml_4, xw=0, yw=0, rw=rw, rc=None, res=0, tsandQ=[(0, Q), (1e08, 0)], layers=1
)
ml_4.solve()

self.neq  1
solution complete

ca_4 = ttim.Calibrate(ml_4)
ca_4.set_parameter(name="kaq1", initial=1, pmin=0)
ca_4.set_parameter(name="Saq1", initial=1e-5, pmin=0)
ca_4.set_parameter(name="c1", initial=100, pmin=0)
ca_4.set_parameter_by_reference(name="rc", parameter=w_4.rc[:], initial=0.2, pmin=0)
ca_4.series(name="well", x=r1, y=0, t=t1, h=h1, layer=1)
ca_4.series(name="obs_well", x=r2, y=0, t=t2, h=h2, layer=1)
ca_4.fit(report=True)

...............................................................
Fit succeeded.
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 60
    # data points      = 40
    # variables        = 4
    chi-square         = 0.56419452
    reduced chi-square = 0.01567207
    Akaike info crit   = -162.449426
    Bayesian info crit = -155.693908
[[Variables]]
    kaq1:  1.93430644 +/- 0.01232844 (0.64%) (init = 1)
    Saq1:  1.3182e-05 +/- 2.3131e-06 (17.55%) (init = 1e-05)
    c1:    239.727821 +/- 20.5021300 (8.55%) (init = 100)
    rc:    0.05268527 +/- 4.0799e-04 (0.77%) (init = 0.2)
[[Correlations]] (unreported correlations are < 0.100)
    C(kaq1, c1)   = +0.8394
    C(Saq1, rc)   = -0.5818
    C(kaq1, Saq1) = -0.5034
    C(Saq1, c1)   = -0.2139
    C(c1, rc)     = -0.1116

display(ca_4.parameters)
print("RMSE:", ca_4.rmse())
optimal std perc_std pmin pmax initial parray
kaq1 1.934306 0.012328 0.637357 0 inf 1.00000 [1.93430643512287]
Saq1 0.000013 0.000002 17.547945 0 inf 0.00001 [1.3181571984599572e-05]
c1 239.727821 20.502130 8.552253 0 inf 100.00000 [239.72782149888832]
rc 0.052685 0.000408 0.774393 0 inf 0.20000 [0.05268527152690705] 
RMSE: 0.11876389590424886

The new estimates had much better confidence intervals when we used part of the parameters fixed. The RMSE, BIC and AIC indicators also showed better results.

hm1_4 = ml_4.head(r1, 0, t1)
hm2_4 = ml_4.head(r2, 0, t2)
plt.figure()
plt.semilogx(t1, h1, ".", label="well")
plt.semilogx(t2, h2, ".", label="obs_well")
plt.semilogx(t1, hm1_4[-1], label="ttim model - well")
plt.semilogx(t2, hm2_4[-1], label="ttim model - observation well")
plt.xlabel("time [d]")
plt.ylabel("drawdown [m]")
plt.title("Model Results - Two Aquifer Model + Well Parameters")
plt.legend(bbox_to_anchor=(1.05, 1))
plt.grid()
../_images/de47abd4445f26aa6e7ec69609dffef9fb1c6afc84fc83afdad55575b7537fe0.png

Step 8. Comparison of Results#

The results simulated with the two aquifer models are present below. Furthermore, Yang (2020) compared TTim results with the results obtained from MLU (Carlson & Randall, 2012) and the published results obtained by Brian et al. (1997). The authors of the original paper applied a finite-difference model of radial flow with a trial and error approach to find the optimal parameters.

t = pd.DataFrame(
    columns=[
        "k0[m/d]",
        "k1[m/d]",
        "Ss0[1/m]",
        "Ss1[1/m]",
        "Sll[1/m]",
        "c[d]",
        "res",
        "rc",
    ],
    index=["MLU", "Brian et al.", "ttim", "ttim-rc", "ttim-fixed upper"],
)
t.loc["ttim-rc"] = ca_3.parameters["optimal"].values
t.iloc[2, 0:6] = ca_2.parameters["optimal"].values
t.iloc[4, 5] = ca_4.parameters["optimal"].values[2]
t.iloc[4, 7] = ca_4.parameters["optimal"].values[3]
t.iloc[4, 0] = 17.28
t.iloc[4, 1] = ca_4.parameters["optimal"].values[0]
t.iloc[4, 2] = 1.2e-4
t.iloc[4, 3] = ca_4.parameters["optimal"].values[1]
t.iloc[4, 4] = 3e-5
t.iloc[0, 0:6] = [17.424, 6.027e-05, 1.747, 6.473e-06, 3.997e-05, 216]
t.iloc[1, 0:6] = [17.28, 3.456, 1.2e-04, 1.5e-5, 3e-05, np.nan]
t["RMSE"] = [0.023452, np.nan, ca_2.rmse(), ca_3.rmse(), ca_4.rmse()]
t
k0[m/d] k1[m/d] Ss0[1/m] Ss1[1/m] Sll[1/m] c[d] res rc RMSE
MLU 17.424 0.00006 1.747 0.000006 0.00004 216 NaN NaN 0.023452
Brian et al. 17.28 3.456 0.00012 0.000015 0.00003 NaN NaN NaN NaN
ttim 199.16155 0.097464 5.610347 1.737358 0.063543 0.001425 NaN NaN 0.744305
ttim-rc 5.492296 1.251603 0.0 0.000045 0.000003 0.745506 0.000018 0.055127 0.178735
ttim-fixed upper 17.28 1.934315 0.00012 0.000013 0.00003 239.719879 NaN 0.052688 0.118764

The first TTim model showed unrealistic results. The second model, although improving the fit, had wide standard deviations (and confidence intervals). That indicates that this solution, with more parameters, in multi-aquifer configuration, is non-unique.

When we kept some of the data fixed with the values from Brian et al., we improved the fit and the confidence interval of the estimates of the model.

MLU results have very low RMSE. However, the specific storage of the upper aquifer (Ss0) is probably too high.

Step 8.1. Comparison of estimates and model error#

# Preparing the DataFrame:
t1 = pd.DataFrame(
    columns=["kaq - opt -l1", "kaq - 95% -l1", "kaq - opt -l2", "kaq - 95% -l2"],
    index=["MLU", "Brian", "TTim - rc", "TTim - fixed upper"],
)
simulation = ["MLU", "Brian", "TTim - rc", "TTim - fixed upper"]
t1.loc["MLU"] = [17.424, 124.160 * 1e-2 * 17.424, 1.747, 4.521 * 1e-2 * 1.747]
t1.loc["Brian"] = [17.28, np.nan, 3.456, np.nan]
t1.loc["TTim - fixed upper"] = [
    17.28,
    np.nan,
    ca_4.parameters.loc["kaq1", "optimal"],
    2 * ca_4.parameters.loc["kaq1", "std"],
]
t1.loc["TTim - rc"] = [
    ca_3.parameters.loc["kaq0", "optimal"],
    2 * ca_3.parameters.loc["kaq0", "std"],
    ca_3.parameters.loc["kaq1", "optimal"],
    2 * ca_3.parameters.loc["kaq1", "std"],
]

# Plotting

plt.figure(figsize=(10, 7))

plt.errorbar(
    x=t1.index,
    y=t1["kaq - opt -l1"],
    yerr=[t1["kaq - 95% -l1"], t1["kaq - 95% -l1"]],
    marker="o",
    linestyle="",
    markersize=12,
    label="upper aquifer",
)
plt.errorbar(
    x=t1.index,
    y=t1["kaq - opt -l2"],
    yerr=[t1["kaq - 95% -l2"], t1["kaq - 95% -l2"]],
    marker="o",
    linestyle="",
    markersize=12,
    label="lower aquifer",
)

plt.legend()
plt.title("Error bar plot of calibrated hydraulic conductivities")
plt.ylabel("K [m/d]")
# plt.ylim([278,289])
plt.xlabel("Model")
../_images/4fd75eb425845551ed1e963cca4f2db3f3c091a952861fe49e02ddbd49515b96.png

The MLU model has a larger confidence interval for the first aquifer layer, while TTim got lower values for the generic model (TTim - rc) and a relatively small range. The model with some fixed values has a small confidence interval. Confidence intervals were not reported in Brian et al. (1997).

References#

  • Bakker, M. Semi-analytic modeling of transient multi-layer flow with TTim. Hydrogeol J 21, 935–943 (2013). https://doi.org/10.1007/s10040-013-0975-2

  • Brian, Schroth, T., N., Narasimhan, 1997. Application of a numerical model in the interpretation of a leaky aquifer test. Ground Water .

  • 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.

  • Newville, M.,Stensitzki, T., Allen, D.B., Ingargiola, A. (2014) LMFIT: Non Linear Least-Squares Minimization and Curve Fitting for Python.https://dx.doi.org/10.5281/zenodo.11813. https://lmfit.github.io/lmfit-py/intro.html (last access: August,2021).

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