<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Confined Aquifer Test
This example is taken from Kruseman and de Ridder (1970)
Step1:  Set basic parameters for the model
Step2:  Create conceptual model
Step3:  Load data of two observation wells
Step4:  Calibrate using only the data from observation well 1
Step5:  Calibrate using only the data from observation well 2
Step6:  Calibrate model with two datasets simultaneously
Step7:  Investigate whether adding well bore storage improves the fit
Step8:  Calibrate using only the data from observation well 1
Step9:  Calibrate using only the data from observation well 2
Step10:  Calibrate model with two datasets simultaneously
Step11:  Summary of values presented in Kruseman and de Ridder (1970)
Compare effect of rc
Step12:  Adding wellbore storage improve the performance when use drawdown data of two observation wells respectively. However, when calibrate model with two datasets simultaneously, rc approaches minimum value. Adding rc does not improve the performance much.
Compare ttim to results of K&dR, AQTEOLV and MLU
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from ttim import *
import pandas as pd
Explanation: Confined Aquifer Test
This example is taken from Kruseman and de Ridder (1970)
End of explanation
H = 7 #aquifer thickness
zt = -18 #top boundary of aquifer
zb = zt - H #bottom boundary of aquifer
Q = 788 #constant discharge
Explanation: Set basic parameters for the model:
End of explanation
#unkonwn parameters: kaq, Saq
ml = ModelMaq(kaq=60, z=[zt, zb], Saq=1e-4, tmin=1e-5, tmax=1)
w = Well(ml, xw=0, yw=0, rw=0.2, tsandQ=[(0, Q)], layers=0)
ml.solve(silent='True')
Explanation: Create conceptual model:
End of explanation
#time and drawdown of piezometer 30m away from pumping well
data1 = np.loadtxt('data/piezometer_h30.txt', skiprows = 1)
t1 = data1[:, 0] / 60 / 24 #convert min to days
h1 = data1[:, 1]
r1 = 30
#time and drawdown of piezometer 90m away from pumping well
data2 = np.loadtxt('data/piezometer_h90.txt', skiprows = 1)
t2 = data2[:, 0] / 60 / 24 #convert min to days
h2 = data2[:, 1]
r2 = 90
Explanation: Load data of two observation wells:
End of explanation
ca1 = Calibrate(ml)
ca1.set_parameter(name='kaq0', initial=10)
ca1.set_parameter(name='Saq0', initial=1e-4)
ca1.series(name='obs1', x=r1, y=0, t=t1, h=h1, layer=0)
ca1.fit(report=True)
display(ca1.parameters)
print('rmse:', ca1.rmse())
hm1 = ml.head(r1, 0, t1)
plt.figure(figsize=(8, 5))
plt.semilogx(t1, h1, '.', label='obs at 30 m')
plt.semilogx(t1, hm1[0], label='ttim at 30 m')
plt.xlabel('time(d)')
plt.ylabel('drawdown(m)')
plt.title('ttim analysis for Oude Korendijk')
plt.legend();
Explanation: Calibrate using only the data from observation well 1:
End of explanation
ca2 = Calibrate(ml)
ca2.set_parameter(name='kaq0', initial=10)
ca2.set_parameter(name='Saq0', initial=1e-4)
ca2.series(name='obs2', x=r2, y=0, t=t2, h=h2, layer=0)
ca2.fit(report=True)
display(ca2.parameters)
print('rmse:', ca2.rmse())
hm2 = ml.head(r2, 0, t2)
plt.figure(figsize=(8, 5))
plt.semilogx(t2, h2, '.', label='obs at 90 m')
plt.semilogx(t2, hm2[0], label='ttim at 90 m')
plt.xlabel('time(d)')
plt.ylabel('drawdown(m)')
plt.title('ttim analysis for Oude Korendijk')
plt.legend();
Explanation: Calibrate using only the data from observation well 2:
End of explanation
ca = Calibrate(ml)
ca.set_parameter(name='kaq0', initial=10)
ca.set_parameter(name='Saq0', initial=1e-4)
ca.series(name='obs1', x=r1, y=0, t=t1, h=h1, layer=0)
ca.series(name='obs2', x=r2, y=0, t=t2, h=h2, layer=0)
ca.fit(report=True)
display(ca.parameters)
print('rmse:', ca.rmse())
hs1 = ml.head(r1, 0, t1)
hs2 = ml.head(r2, 0 ,t2)
plt.figure(figsize = (8, 5))
plt.semilogx(t1, h1, '.', label='obs at 30m')
plt.semilogx(t1, hs1[0], label='ttim at 30 m')
plt.semilogx(t2, h2, '.', label='obs at 90m')
plt.semilogx(t2, hs2[0], label = 'ttim at 90m')
plt.xlabel('time(d)')
plt.ylabel('drawdown(m)')
plt.title('ttim analysis for Oude Korendijk')
plt.legend();
Explanation: Calibrate model with two datasets simultaneously:
End of explanation
#unknown parameters: kaq, Saq and rc
ml1 = ModelMaq(kaq=60, z=[zt, zb], Saq=1e-4, tmin=1e-5, tmax=1)
w1 = Well(ml1, xw=0, yw=0, rw=0.2, rc=0.2, tsandQ=[(0, Q)], layers=0)
ml1.solve(silent='True')
Explanation: Investigate whether adding well bore storage improves the fit:
Create conceptual model:
End of explanation
ca3 = Calibrate(ml1)
ca3.set_parameter(name='kaq0', initial=10)
ca3.set_parameter(name='Saq0', initial=1e-4)
ca3.set_parameter_by_reference(name='rc', parameter=w1.rc[0:], initial=0.2, pmin=0.01)
ca3.series(name='obs1', x=r1, y=0, t=t1, h=h1, layer=0)
ca3.fit(report=True)
display(ca3.parameters)
print('rmse:', ca3.rmse())
hm3 = ml1.head(r1, 0, t1)
plt.figure(figsize=(8, 5))
plt.semilogx(t1, h1, '.', label='obs at 30 m')
plt.semilogx(t1, hm3[0], label='ttim at 30 m')
plt.xlabel('time(d)')
plt.ylabel('drawdown(m)')
plt.title('ttim analysis for Oude Korendijk')
plt.legend();
Explanation: Calibrate using only the data from observation well 1:
End of explanation
ca4 = Calibrate(ml1)
ca4.set_parameter(name='kaq0', initial=10)
ca4.set_parameter(name='Saq0', initial=1e-4)
ca4.set_parameter_by_reference(name='rc', parameter=w1.rc[0:], initial=0.2, pmin=0.01)
ca4.series(name='obs2', x=r2, y=0, t=t2, h=h2, layer=0)
ca4.fit(report=True)
display(ca4.parameters)
print('rmse:', ca4.rmse())
hm4 = ml1.head(r2, 0, t2)
plt.figure(figsize=(8, 5))
plt.semilogx(t2, h2, '.', label='obs at 90 m')
plt.semilogx(t2, hm4[0], label='ttim at 90 m')
plt.xlabel('time(d)')
plt.ylabel('drawdown(m)')
plt.title('ttim analysis for Oude Korendijk')
plt.legend();
Explanation: Calibrate using only the data from observation well 2:
End of explanation
ca0 = Calibrate(ml1)
ca0.set_parameter(name='kaq0', initial=10)
ca0.set_parameter(name='Saq0', initial=1e-4)
ca0.set_parameter_by_reference(name='rc', parameter=w1.rc[0:], initial=0.2, pmin=0.01)
ca0.series(name='obs1', x=r1, y=0, t=t1, h=h1, layer=0)
ca0.series(name='obs2', x=r2, y=0, t=t2, h=h2, layer=0)
ca0.fit(report=True)
display(ca0.parameters)
print('rmse:', ca0.rmse())
hs1 = ml1.head(r1, 0, t1)
hs2 = ml1.head(r2, 0 ,t2)
plt.figure(figsize = (8, 5))
plt.semilogx(t1, h1, '.', label='obs at 30m')
plt.semilogx(t1, hs1[0], label='ttim at 30 m')
plt.semilogx(t2, h2, '.', label='obs at 90m')
plt.semilogx(t2, hs2[0], label = 'ttim at 90m')
plt.xlabel('time(d)')
plt.ylabel('drawdown(m)')
plt.title('ttim analysis for Oude Korendijk')
plt.legend();
Explanation: Calibrate model with two datasets simultaneously:
End of explanation
t0 = pd.DataFrame(columns=['obs 30 m', 'obs 90 m', 'obs simultaneously'], index=['without rc', 'with rc'])
t0.loc['without rc', 'obs 30 m'] = ca1.rmse()
t0.loc['without rc', 'obs 90 m'] = ca2.rmse()
t0.loc['without rc', 'obs simultaneously'] = ca.rmse()
t0.loc['with rc', 'obs 30 m'] = ca3.rmse()
t0.loc['with rc', 'obs 90 m'] = ca4.rmse()
t0.loc['with rc', 'obs simultaneously'] = ca0.rmse()
print('RMSE of two conceptual models:')
t0
Explanation: Summary of values presented in Kruseman and de Ridder (1970)
Compare effect of rc:
End of explanation
t = pd.DataFrame(columns=['k [m/d]', 'Ss [1/m]', 'RMSE'], \
                 index=['K&dR', 'ttim', 'AQTESOLV', 'MLU'])
t.loc['ttim'] = np.append(ca.parameters['optimal'].values, ca.rmse())
t.loc['AQTESOLV'] = [66.086, 2.541e-05, 0.05006]
t.loc['MLU'] = [66.850, 2.400e-05, 0.05083]
t.loc['K&dR'] = [55.71429, 1.7E-4, '-']
t
Explanation: Adding wellbore storage improve the performance when use drawdown data of two observation wells respectively. However, when calibrate model with two datasets simultaneously, rc approaches minimum value. Adding rc does not improve the performance much.
Compare ttim to results of K&dR, AQTEOLV and MLU:
End of explanation
<END_TASK>