In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from ttim import *
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:
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.
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?
In [2]:
def generate_data():
# 2 layer model with some random error
ml = ModelMaq(kaq=[10, 20], z=[0, -20, -22, -42], c=[1000],
Saq=[0.0002, 0.0001], tmin=0.001, tmax=100)
w = 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 * np.random.randn(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
In [3]:
np.random.seed(11)
data = generate_data()
to = data[:, 0]
ho = data[:, 1]
In [4]:
def func(p, to=to, ho=ho, returnmodel=False):
k = p[0]
S = p[1]
ml = ModelMaq(kaq=k, z=[0, -20], Saq=S, tmin=0.001, tmax=100)
w = 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)
In [5]:
from scipy.optimize import fmin
lsopt = fmin(func, [10, 1e-4])
print('optimal parameters:', lsopt)
print('rmse:', np.sqrt(func(lsopt) / len(ho)))
In [6]:
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)');
In [7]:
cal = Calibrate(ml)
cal.set_parameter(name='kaq0', initial=10, pmin=0.1, pmax=1000)
cal.set_parameter(name='Saq0', 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())
In [8]:
cal.parameters.style.set_precision(3)
Out[8]:
In [9]:
def func2(p, to=to, ho=ho, returnmodel=False):
k = p[0]
S = p[1]
c = p[2]
ml = ModelMaq(kaq=k, z=[2, 0, -20], Saq=S, c=c, topboundary='semi',
tmin=0.001, tmax=100)
w = 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)
In [10]:
lsopt2 = fmin(func2, [10, 1e-4, 1000])
print('optimal parameters:', lsopt2)
print('rmse:', np.sqrt(func2(lsopt2) / len(ho)))
In [11]:
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)');
In [12]:
ml = ModelMaq(kaq=10, z=[2, 0, -20], Saq=1e-4, c=1000, topboundary='semi', tmin=0.001, tmax=100)
w = Well(ml, 0, 0, rw=0.3, tsandQ=[(0, 800)])
ml.solve(silent=True)
In [15]:
cal = Calibrate(ml)
cal.set_parameter(name='kaq0', initial=10)
cal.set_parameter(name='Saq0', initial=1e-4)
cal.set_parameter(name='c0', initial=1000)
cal.series(name='obs1', x=10, y=0, layer=0, t=to, h=ho)
cal.fit(report=False)
cal.parameters.style.set_precision(5)
Out[15]:
In [16]:
cal.rmse(), ml.aq.kaq
Out[16]:
In [17]:
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)');