<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 test is taken from examples presented in MLU tutorial.
Step1:  The test is condected at a fully confined two-aquifer system. Both the pumping well and the observation piezometer are screened at the second aquifer.
Set basic parameters
Step2:  Load data of two observation wells
Step3:  Create single layer model (overlying aquifer and aquitard are excluded)
Step4:  To improve model's performance, rc & res are adding
Step5:  Create three-layer conceptual model
Step6:  Try adding res & rc
Step7:  Calibrate with fitted characters for upper aquifer
Step8:  The optimized value of res is very close to the minimum limitation, thus res has little effect on the performance of the model. res is removed in this calibration.
Step9:  Summary of values simulated by MLU
Results of calibrations done with three-layer model of ttim are presented below.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from ttim import *
Explanation: Confined Aquifer Test
This test is taken from examples presented in MLU tutorial.
End of explanation
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
Explanation: The test is condected at a fully confined two-aquifer system. Both the pumping well and the observation piezometer are screened at the second aquifer.
Set basic parameters:
End of explanation
data1 = np.loadtxt('data/schroth_obs1.txt', skiprows = 1)
t1 = data1[:, 0]
h1 = data1[:, 1]
r1 = 0 
data2 = np.loadtxt('data/schroth_obs2.txt', skiprows = 1)
t2 = data2[:, 0]
h2 = data2[:, 1]
r2 = 46 #distance between observation well2 and pumping well
Explanation: Load data of two observation wells:
End of explanation
ml_0 = ModelMaq(z=[zt1, zb1], kaq=10, Saq=1e-4, tmin=1e-4, tmax=1)
w_0 = Well(ml_0, xw=0, yw=0, rw=rw, tsandQ = [(0, Q), (1e+08, 0)])
ml_0.solve()
ca_0 = Calibrate(ml_0)
ca_0.set_parameter(name='kaq0', initial=10)
ca_0.set_parameter(name='Saq0', initial=1e-4)
ca_0.series(name='obs1', x=r1, y=0, t=t1, h=h1, layer=0)
ca_0.series(name='obs2', x=r2, y=0, t=t2, h=h2, layer=0)
ca_0.fit(report=True)
display(ca_0.parameters)
print('RMSE:', ca_0.rmse())
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='obs1')
plt.semilogx(t2, h2, '.', label='obs2')
plt.semilogx(t1, hm1_0[-1], label='ttim1')
plt.semilogx(t2, hm2_0[-1], label='ttim2')
plt.xlabel('time(d)')
plt.ylabel('head(m)')
plt.legend()
plt.savefig('C:/Users/DELL/Python Notebook/MT BE/Fig/schroth_one1.eps');
Explanation: Create single layer model (overlying aquifer and aquitard are excluded):
End of explanation
ml_1 = ModelMaq(z=[zt1, zb1], kaq=10, Saq=1e-4, tmin=1e-4, tmax=1)
w_1 = Well(ml_1, xw=0, yw=0, rw=rw, rc=0, res=5, tsandQ = [(0, Q), (1e+08, 0)])
ml_1.solve()
ca_1 = 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='obs1', x=r1, y=0, t=t1, h=h1, layer=0)
ca_1.series(name='obs2', x=r2, y=0, t=t2, h=h2, layer=0)
ca_1.fit(report=True)
display(ca_1.parameters)
print('RMSE:', ca_1.rmse())
hm1_1 = ml_1.head(r1, 0, t1)
hm2_1 = ml_1.head(r2, 0, t2)
plt.figure(figsize = (8, 5))
plt.semilogx(t1, h1, '.', label='obs1')
plt.semilogx(t2, h2, '.', label='obs2')
plt.semilogx(t1, hm1_1[-1], label='ttim1')
plt.semilogx(t2, hm2_1[-1], label='ttim2')
plt.xlabel('time(d)')
plt.ylabel('head(m)')
plt.legend()
plt.savefig('C:/Users/DELL/Python Notebook/MT BE/Fig/schroth_one2.eps');
Explanation: To improve model's performance, rc & res are adding:
End of explanation
ml_2 = 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 = Well(ml_2, xw=0, yw=0, rw=rw, tsandQ = [(0, Q), (1e+08, 0)], layers=1)
ml_2.solve()
ca_2 = Calibrate(ml_2)
ca_2.set_parameter(name= 'kaq0', initial=20, pmin=0)
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='obs1', x=r1, y=0, t=t1, h=h1, layer=1)
ca_2.series(name='obs2', x=r2, y=0, t=t2, h=h2, layer=1)
ca_2.fit(report=True)
display(ca_2.parameters)
print('RMSE:',ca_2.rmse())
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='obs1')
plt.semilogx(t2, h2, '.', label='obs2')
plt.semilogx(t1, hm1_2[-1], label='ttim1')
plt.semilogx(t2, hm2_2[-1], label='ttim2')
plt.xlabel('time(d)')
plt.ylabel('head(m)')
plt.legend()
plt.savefig('C:/Users/DELL/Python Notebook/MT BE/Fig/schroth_three1.eps');
Explanation: Create three-layer conceptual model:
End of explanation
ml_3 = 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 = Well(ml_3, xw=0, yw=0, rw=rw, rc=None, res=0, tsandQ = [(0, Q), (1e+08, 0)], \
           layers=1)
ml_3.solve()
ca_3 = 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=0)
ca_3.set_parameter(name='Saq1', initial=1e-5, pmin=0)
ca_3.set_parameter_by_reference(name='Sll', parameter=ml_3.aq.Sll[:],\
                                initial=1e-4, pmin=0)
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())
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='obs1')
plt.semilogx(t2, h2, '.', label='obs2')
plt.semilogx(t1, hm1_3[-1], label='ttim1')
plt.semilogx(t2, hm2_3[-1], label='ttim2')
plt.xlabel('time(d)')
plt.ylabel('head(m)')
plt.legend()
plt.savefig('C:/Users/DELL/Python Notebook/MT BE/Fig/schroth_three2.eps');
Explanation: Try adding res & rc:
End of explanation
ml_4 = 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 = Well(ml_4, xw=0, yw=0, rw=rw, rc=None, res=0, tsandQ = [(0, Q), (1e+08, 0)], \
           layers=1)
ml_4.solve()
Explanation: Calibrate with fitted characters for upper aquifer:
End of explanation
ca_4 = 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='obs1', x=r1, y=0, t=t1, h=h1, layer=1)
ca_4.series(name='obs2', x=r2, y=0, t=t2, h=h2, layer=1)
ca_4.fit(report=True)
display(ca_4.parameters)
print('RMSE:', ca_4.rmse())
hm1_4 = ml_4.head(r1, 0, t1)
hm2_4 = ml_4.head(r2, 0, t2)
plt.figure(figsize = (8, 5))
plt.semilogx(t1, h1, '.', label='obs1')
plt.semilogx(t2, h2, '.', label='obs2')
plt.semilogx(t1, hm1_4[-1], label='ttim1')
plt.semilogx(t2, hm2_4[-1], label='ttim2')
plt.xlabel('time(d)')
plt.ylabel('head(m)')
plt.legend()
plt.savefig('C:/Users/DELL/Python Notebook/MT BE/Fig/schroth_three3.eps');
Explanation: The optimized value of res is very close to the minimum limitation, thus res has little effect on the performance of the model. res is removed in this calibration.
End of explanation
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', 'MLU-fixed k1','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] = [2.020e-04, 9.110e-04, 3.456, 6.214e-05, 7.286e-05, 453.5]
t['RMSE'] = [0.023452, 0.162596, ca_2.rmse(), ca_3.rmse(), ca_4.rmse()]
t
Explanation: Summary of values simulated by MLU
Results of calibrations done with three-layer model of ttim are presented below.
End of explanation
<END_TASK>