<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Slug Test for Confined Aquifer
This test is taken from examples of AQTESOLV.
Step1:  Set background parameters
Step2:  Slug
Step3:  Load data
Step4:  Create single layer conceptual model
Step5:  Calibrate with two datasets simultaneously
Step6:  Try adding well skin resistance res
Step7:  Adding well screen resistance does not improve the performance obviously. While the AIC value increases. Thus, res should be removed from the model.
Try multilayer conceptual model
Step8:  Calibrate with two datasets simultaneously
Step9:  Summary of values presented by AQTESOLV & MLU
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
from ttim import *
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
Explanation: Slug Test for Confined Aquifer
This test is taken from examples of AQTESOLV.
End of explanation
H0 = 2.798 #initial displacement in m
b = -6.1 #aquifer thickness
rw1 = 0.102 #well radius of Ln-2 Well
rw2 = 0.071 #well radius of observation Ln-3 Well
rc1 = 0.051 #casing radius of Ln-2 Well
rc2 = 0.025 #casing radius of Ln-3 Well
r = 6.45 #distance from observation well to test well
Explanation: Set background parameters:
End of explanation
Q = np.pi * rc1 ** 2 * H0
print('Slug:', round(Q, 5), 'm^3')
Explanation: Slug:
End of explanation
data1 = np.loadtxt('data/ln-2.txt')
t1 = data1[:, 0] / 60 / 60 / 24 #convert time from seconds to days
h1 = data1[:, 1]
data2 = np.loadtxt('data/ln-3.txt')
t2 = data2[:, 0] / 60 / 60 / 24
h2 = data2[:, 1]
Explanation: Load data:
End of explanation
ml_0 = ModelMaq(kaq=10, z=[0, b], Saq=1e-4,  \
           tmin=1e-5, tmax=0.01)
w_0 = Well(ml_0, xw=0, yw=0, rw=rw1, rc=rc1, tsandQ=[(0, -Q)], layers=0, wbstype='slug')
ml_0.solve()
Explanation: Create single layer conceptual model:
End of explanation
#unknown parameters: kaq, Saq
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='Ln-2', x=0, y=0, layer=0, t=t1, h=h1)
ca_0.series(name='Ln-3', x=r, y=0, layer=0, t=t2, h=h2)
ca_0.fit(report=True)
display(ca_0.parameters)
print('RMSE:', ca_0.rmse())
hm1_0 = ml_0.head(0, 0, t1, layers=0)
hm2_0 = ml_0.head(r, 0, t2, layers=0)
plt.figure(figsize=(8, 5))
plt.semilogx(t1, h1/H0, '.', label='obs ln-2')
plt.semilogx(t1, hm1_0[0]/H0, label='ttim ln-2')
plt.semilogx(t2, h2/H0, '.', label='obs ln-3')
plt.semilogx(t2, hm2_0[0]/H0, label='ttim ln-3')
plt.xlabel('time(d)')
plt.ylabel('h/H0')
plt.legend();
Explanation: Calibrate with two datasets simultaneously:
End of explanation
ml_1 = ModelMaq(kaq=10, z=[0, b], Saq=1e-4,  \
           tmin=1e-5, tmax=0.01)
w_1 = Well(ml_1, xw=0, yw=0, rw=rw1, res=0, rc=rc1, tsandQ=[(0, -Q)], layers=0, wbstype='slug')
ml_1.solve()
#unknown parameters: kaq, Saq, res
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='res', parameter=w_1.res, initial=0)
ca_1.series(name='Ln-2', x=0, y=0, layer=0, t=t1, h=h1)
ca_1.series(name='Ln-3', x=r, y=0, layer=0, t=t2, h=h2)
ca_1.fit(report=True)
display(ca_1.parameters)
print('RMSE:', ca_1.rmse())
hm1_1 = ml_1.head(0, 0, t1, layers=0)
hm2_1 = ml_1.head(r, 0, t2, layers=0)
plt.figure(figsize=(8, 5))
plt.semilogx(t1, h1/H0, '.', label='obs ln-2')
plt.semilogx(t1, hm1_1[0]/H0, label='ttim ln-2')
plt.semilogx(t2, h2/H0, '.', label='obs ln-3')
plt.semilogx(t2, hm2_1[0]/H0, label='ttim ln-3')
plt.xlabel('time(d)')
plt.ylabel('h/H0')
plt.legend();
Explanation: Try adding well skin resistance res:
End of explanation
#Determine elevations of each layer.
#Thickness of each layer is set to be 0.5 m.
z = np.arange(0, b, -0.5)
zlay = np.append(z, b)
nlay = len(zlay) - 1
Saq_2 = 1e-4 * np.ones(nlay)
n = np.arange(0, 13,1)
ml_2 = Model3D(kaq=10, z=zlay, Saq=Saq_2, kzoverkh=1, tmin=1e-5, tmax=0.01, \
              phreatictop=True)
w_2 = Well(ml_2, xw=0, yw=0, rw=rw1, tsandQ=[(0, -Q)], layers=n, rc=rc1, \
          wbstype='slug')
ml_2.solve()
Explanation: Adding well screen resistance does not improve the performance obviously. While the AIC value increases. Thus, res should be removed from the model.
Try multilayer conceptual model:
End of explanation
ca_2 = Calibrate(ml_2)
ca_2.set_parameter(name='kaq0_12', initial=10)
ca_2.set_parameter(name='Saq0_12', initial=1e-4, pmin=0)
ca_2.series(name='Ln-2', x=0, y=0, layer=n, t=t1, h=h1)
ca_2.series(name='Ln-3', x=r, y=0, layer=n, t=t2, h=h2)
ca_2.fit(report=True)
display(ca_2.parameters)
print('RMSE:', ca_2.rmse())
hm1_2 = ml_2.head(0, 0, t1, layers=n)
hm2_2 = ml_2.head(r, 0, t2, layers=n)
plt.figure(figsize=(8, 5))
plt.semilogx(t1, h1/H0, '.', label='obs ln-2')
plt.semilogx(t1, hm1_2[0]/H0, label='ttim ln-2')
plt.semilogx(t2, h2/H0, '.', label='obs ln-3')
plt.semilogx(t2, hm2_2[0]/H0, label='ttim ln-3')
plt.xlabel('time(d)')
plt.ylabel('h/H0')
plt.legend();
Explanation: Calibrate with two datasets simultaneously:
End of explanation
t = pd.DataFrame(columns=['k [m/d]', 'Ss [1/m]'], \
                          index=['MLU', 'AQTESOLV', 'ttim-single', 'ttim-multi'])
t.loc['AQTESOLV'] = [1.166, 9.368E-06]
t.loc['MLU'] = [1.311, 8.197E-06]
t.loc['ttim-single'] = ca_0.parameters['optimal'].values
t.loc['ttim-multi'] = ca_2.parameters['optimal'].values
t['RMSE'] = [0.010373, 0.009151, ca_0.rmse(), ca_1.rmse()]
t
Explanation: Summary of values presented by AQTESOLV & MLU
End of explanation
<END_TASK>