2. Synthetic Pumping Test - 2 aquifers

import matplotlib.pyplot as plt
import numpy as np

import ttim

Head data is generated for a pumping test in a two-aquifer model. The well starts pumping at time \(t=0\) with a discharge \(Q=800\) m\(^3\)/d. The head is measured in an observation well 10 m from the pumping well. The thickness of the aquifer is 20 m. Questions:

1. Determine the optimal values of the hydraulic conductivity and specific storage coefficient of the aquifer when the aquifer is approximated as confined. Use a least squares approach and make use of the fmin function of scipy.optimize to find the optimal values. Plot the data with dots and the best-fit model in one graph. Print the optimal values of \(k\) and \(S_s\) to the screen as well as the root mean squared error of the residuals.

2. Repeat Question 1 but now approximate the aquifer as semi-confined. Plot the data with dots and the best-fit model in one graph. Print to the screen the optimal values of \(k\), \(S_s\) and \(c\) to the screen as well as the root mean squared error of the residuals. Is the semi-cofined model a better fit than the confined model?

def generate_data():
    # 2 layer model with some random error
    ml = ttim.ModelMaq(
        kaq=[10, 20],
        z=[0, -20, -22, -42],
        c=[1000],
        Saq=[0.0002, 0.0001],
        tmin=0.001,
        tmax=100,
    )
    ttim.Well(ml, 0, 0, rw=0.3, tsandQ=[(0, 800)])
    ml.solve()
    t = np.logspace(-2, 1, 100)
    h = ml.head(10, 0, t)
    plt.figure()
    r = 0.01 * rnd.random(100)
    n = np.zeros_like(r)
    # alpha = 0.8
    for i in range(1, len(n)):
        n[i] = 0.8 * n[i - 1] + r[i]
    ho = h[0] + n
    plt.plot(t, ho, ".")
    data = np.zeros((len(ho), 2))
    data[:, 0] = t
    data[:, 1] = ho
    # np.savetxt('pumpingtestdata.txt', data, fmt='%2.3f', header='time (d), head (m)')
    return data

rnd = np.random.default_rng(11)
data = generate_data()
to = data[:, 0]
ho = data[:, 1]

self.neq  1
solution complete

def func(p, to=to, ho=ho, returnmodel=False):
    k = p[0]
    S = p[1]
    ml = ttim.ModelMaq(kaq=k, z=[0, -20], Saq=S, tmin=0.001, tmax=100)
    ttim.Well(ml, 0, 0, rw=0.3, tsandQ=[(0, 800)])
    ml.solve(silent=True)
    if returnmodel:
        return ml
    h = ml.head(10, 0, to)
    return np.sum((h[0] - ho) ** 2)

from scipy.optimize import fmin

lsopt = fmin(func, [10, 1e-4])
print("optimal parameters:", lsopt)
print("rmse:", np.sqrt(func(lsopt) / len(ho)))

Optimization terminated successfully.
         Current function value: 0.251681
         Iterations: 41
         Function evaluations: 85
optimal parameters: [1.16063293e+01 1.27878469e-04]
rmse: 0.05016784816354969

ml = func(lsopt, returnmodel=True)
plt.figure()
plt.plot(data[:, 0], data[:, 1], ".", label="observed")
hm = ml.head(10, 0, to)
plt.plot(to, hm[0], "r", label="modeled")
plt.legend()
plt.xlabel("time (d)")
plt.ylabel("head (m)")

Text(0, 0.5, 'head (m)')

cal = ttim.Calibrate(ml)
cal.set_parameter(name="kaq0", layers=0, initial=10, pmin=0.1, pmax=1000)
cal.set_parameter(name="Saq0", layers=0, initial=1e-4, pmin=1e-5, pmax=1e-3)
cal.series(name="obs1", x=10, y=0, layer=0, t=to, h=ho)
cal.fit(report=False)
print("rmse:", cal.rmse())

......................
Fit succeeded.
rmse: 0.050167861160009285

cal.parameters

	layers	optimal	std	perc_std	pmin	pmax	initial	inhoms	parray
kaq0_0_0	0	11.606187	0.118774	1.023369	0.10000	1000.000	10.0000	None	[[11.60618735610534]]
Saq0_0_0	0	0.000128	0.000007	5.559354	0.00001	0.001	0.0001	None	[[0.0001278868669011676]]

Model as semi-confined

def func2(p, to=to, ho=ho, returnmodel=False):
    k = p[0]
    S = p[1]
    c = p[2]
    ml = ttim.ModelMaq(
        kaq=k, z=[2, 0, -20], Saq=S, c=c, topboundary="semi", tmin=0.001, tmax=100
    )
    ttim.Well(ml, 0, 0, rw=0.3, tsandQ=[(0, 800)])
    ml.solve(silent=True)
    if returnmodel:
        return ml
    h = ml.head(10, 0, to)
    return np.sum((h[0] - ho) ** 2)

lsopt2 = fmin(func2, [10, 1e-4, 1000])
print("optimal parameters:", lsopt2)
print("rmse:", np.sqrt(func2(lsopt2) / len(ho)))

Optimization terminated successfully.
         Current function value: 0.004826
         Iterations: 119
         Function evaluations: 236
optimal parameters: [1.02193607e+01 2.01296811e-04 1.56256292e+03]
rmse: 0.006947175980032902

ml = func2(lsopt2, returnmodel=True)
plt.figure()
plt.plot(data[:, 0], data[:, 1], ".", label="observed")
hm = ml.head(10, 0, to)
plt.plot(to, hm[0], "r", label="modeled")
plt.legend()
plt.xlabel("time (d)")
plt.ylabel("head (m)")

Text(0, 0.5, 'head (m)')

ml = ttim.ModelMaq(
    kaq=10, z=[2, 0, -20], Saq=1e-4, c=1000, topboundary="semi", tmin=0.001, tmax=100
)
w = ttim.Well(ml, 0, 0, rw=0.3, tsandQ=[(0, 800)])
ml.solve(silent=True)

cal = ttim.Calibrate(ml)
cal.set_parameter(name="kaq0", layers=0, initial=10)
cal.set_parameter(name="Saq0", layers=0, initial=1e-4)
cal.set_parameter(name="c0", layers=0, initial=1000)
cal.series(name="obs1", x=10, y=0, layer=0, t=to, h=ho)
cal.fit(report=False)
cal.parameters

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

	layers	optimal	std	perc_std	pmin	pmax	initial	inhoms	parray
kaq0_0_0	0	10.219379	0.023065	0.225699	-inf	inf	10.0000	None	[[10.219379490020255]]
Saq0_0_0	0	0.000201	0.000002	0.880785	-inf	inf	0.0001	None	[[0.000201297027266808]]
c0_0_0	0	1562.587540	37.886728	2.424615	-inf	inf	1000.0000	None	[[1562.5875401213193]]

cal.rmse(), ml.aq.kaq

(np.float64(0.006947177191584833), array([10.21937949]))

plt.figure()
plt.plot(data[:, 0], data[:, 1], ".", label="observed")
hm = ml.head(10, 0, to)
plt.plot(to, hm[0], "r", label="modeled")
plt.legend()
plt.xlabel("time (d)")
plt.ylabel("head (m)")

Text(0, 0.5, 'head (m)')