<jupyter_start><jupyter_text>1. Slug Test - Pratt County**This test is taken from AQTESOLV examples.**  Introduction and Conceptual Model In this notebook, we reproduce the work of Yang (2020) to check the TTim performance in analysing slug tests. We later compare the solution in TTim with the KGS analytical model (Hyder et al. 1994) implemented in AQTESOLV (Duffield, 2007).This slug test was conducted in Pratt County Monitoring Site, US, and reported by Butler (1998). A partially penetrating well is screened in unconsolidated alluvial deposits, consisting of sand and gravel interbedded by clay. The total thickness of the aquifer is 47.87 m. The screen is located at 16.77 m depth and has a screen length 1.52 m the well radius is 0.125 m, and the casing radius 0.064 m.The slug displacement is 0.671 m. Head change has been recorded at the slug well.The conceptual model can be seen in the figure below.<jupyter_code>import matplotlib.pyplot as plt
import numpy as np

##Now printing the conceptual model figure:

fig = plt.figure(figsize=(14, 9))
ax = fig.add_subplot(1,1,1)
#sky
sky = plt.Rectangle((-20,0), width = 50, height = 20, fc = 'b', zorder=0, alpha=0.1)
ax.add_patch(sky)

#Aquifer:
ground = plt.Rectangle((-20,-47.87), width = 50, height = 47.87, fc = np.array([209,179,127])/255, zorder=0, alpha=0.9)
ax.add_patch(ground)

well = plt.Rectangle((-1,-(16.77+1.52)), width = 2, height = (16.77+1.52), fc = np.array([200,200,200])/255, zorder=1)
ax.add_patch(well)

#Wellhead
wellhead = plt.Rectangle((-1.25,0),width = 2.5, height = 10, fc = np.array([200,200,200])/255, zorder=2, ec='k')
ax.add_patch(wellhead)

#Screen for the well:
screen = plt.Rectangle((-1,-(16.77+1.52)), width = 2, height = 1.52, fc = np.array([200,200,200])/255, alpha=1, zorder = 2, ec = "k", ls = '--')
screen.set_linewidth(2)
ax.add_patch(screen)
pumping_arrow = plt.Arrow(x = 0,y = 10, dx = 0, dy = 3, color = "#00035b")
ax.add_patch(pumping_arrow)
ax.text(x = 0.5, y = 11, s = r'$ D = 0.671 m $', fontsize = 'large' )

#last line
line = plt.Line2D(xdata= [-200,1200], ydata = [0,0], color = "k")
ax.add_line(line)


ax.set_xlim([-20,30])
ax.set_ylim([-47,20])
ax.set_xlabel('Distance [m]')
ax.set_ylabel('Relative height [m]')
ax.set_title('Conceptual Model - Pratt County Example');<jupyter_output><empty_output><jupyter_text>Step 1. Load required libraries<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from ttim import *<jupyter_output><empty_output><jupyter_text>Step 2. Set basic parameters<jupyter_code>rw = 0.125 # well radius
rc = 0.064 # well casing radius
L = 1.52 # screen length
b = -47.87 # aquifer thickness
zt = -16.77 # depth to top of screen
H0 = 0.671 # initial displacement in the well
zb = zt - L # bottom of screen<jupyter_output><empty_output><jupyter_text>Step 3. Calculate the added volumeAs we will see later, the input for the slug test in TTim is the added or removed volume. Therefore we must first convert our measured displacement into volume.<jupyter_code>Q = np.pi * rc ** 2 * H0
print('slug:', round(Q, 5), 'm^3')<jupyter_output>slug: 0.00863 m^3<jupyter_text>Step 4. Load data from well Data from the slug well is loaded from a text file, where the first column is the time in seconds and the second column is the head.<jupyter_code>data = np.loadtxt('data/slug.txt', skiprows = 1)
t = data[:, 0] / 60 / 60 / 24 #convert time to days
h = data[:, 1]<jupyter_output><empty_output><jupyter_text>Step 5. Create a conceptual model in TTimWe conceptualize the aquifer as a three-layer model, one layer above the screen, one layer at the screen top and bottom and another layer just below it.We use ```Model3D``` method to build this model. Details on how to set the model can be seen in notebook: [Unconfined - 1 - Vennebulten](unconfined1_vennebulten.ipynb).The setting of the slug well is slightly different from the pumping well. We detail the differences below:* the ```tsandQ``` argument in the ```Well``` object has a different meaning. Instead of meaning the time of start or shutdown and the pumping rate of the pumping well, it means the time a volume is instantaneously added or removed from the well. In our case, we defined it as ```[(0, -Q)]``` where ```0``` is the moment in time when we added the slug and ```Q``` is the added volume. A negative sign means a volume is added. Otherwise, it would mean an extracted volume.* the ```wbstype``` argument is set to ```' slug'```, so TTim knows the ```tsandQ``` argument means time and instant volumes, instead of pumping rates.<jupyter_code>ml = Model3D(kaq=10, z=[0, zt, zb, b], Saq=1e-4, kzoverkh=1, tmin=1e-6, tmax=0.01)
w = Well(ml, xw=0, yw=0, rw=rw,rc=rc, tsandQ=[(0, -Q)], layers=1, wbstype='slug')
ml.solve()<jupyter_output>self.neq  1
solution complete<jupyter_text>Step 6. Model calibrationWe calibrate both hydraulic conductivity and specific storage, considering uniform parameters in all layers.<jupyter_code>ca = Calibrate(ml)
ca.set_parameter(name='kaq0_2', initial=10)
ca.set_parameter(name='Saq0_2', initial=1e-4)
ca.series(name='obs', x=0, y=0, layer=1, t=t, h=h)
ca.fit(report=True)
display(ca.parameters)
print('RMSE:', ca.rmse())<jupyter_output><empty_output><jupyter_text>Here we plot the results in heads over initial heads, as done for the graphical solution.<jupyter_code>hm = ml.head(0, 0, t, layers=1)
plt.figure(figsize = (10, 7))
plt.semilogx(t, h/H0, '.', label='obs')
plt.semilogx(t, hm[-1]/H0, 'r', label='ttim')
plt.ylim([0, 1])
plt.xlabel('time [d]')
plt.ylabel('h/H0')
plt.title('Model Results - Three layers')
plt.legend();<jupyter_output><empty_output><jupyter_text>Step 7. A conceptual model with more layersNow we check the model performance in a model with more layers.We use the same division that we set before, but now we divide the first layers into 18 layers, the second layer into 3 layers and the third layer into 29 layers. So each layer is roughly 1 m thick.<jupyter_code>n1 = 18
n2 = 3
n3 = 29
nlay = n1 + n2 + n3 #number of layers
zlay1 = np.linspace(0, zt, n1 + 1)
zlay2 = np.linspace(zt, zb, n2 + 1)
zlay3 = np.linspace(zb, b, n3 + 1)
layers = np.append(zlay1[:-1], zlay2[:-1])
layers = np.append(layers, zlay3) #elevation of each layer
Saq = 1e-4 * np.ones(nlay)
Saq[0] = 0.1
M_nlay = Model3D(kaq=10, z=layers, Saq=Saq, kzoverkh=1, phreatictop=True,\
                 tmin=1e-6, tmax=0.01)
W_nlay = Well(M_nlay, xw=0, yw=0, rw=rw, tsandQ=[(0, -Q)], layers=[18,19,20], \
              rc=rc, wbstype='slug')
M_nlay.solve()<jupyter_output>self.neq  3
solution complete<jupyter_text>Step 8. Calibration of multi-layer model<jupyter_code>cM = Calibrate(M_nlay)
cM.set_parameter(name='kaq0_49', initial=10)
cM.set_parameter(name='Saq0_49', initial=1e-4, pmin=1e-7)
cM.series(name='obs', x=0, y=0, layer=[18,19,20], t=t, h=h)
cM.fit()
display(cM.parameters)
print('RMSE:', cM.rmse())
hM = M_nlay.head(0, 0, t, layers = 20)
plt.figure(figsize = (10, 7))
plt.semilogx(t, h / H0, '.', label = 'obs')
plt.semilogx(t, hM[0] / H0, label = 'ttim')
plt.ylim([0, 1])
plt.xlabel('time(d)')
plt.ylabel('h/H0')
plt.title("Model Results - more layers")
plt.legend();<jupyter_output><empty_output><jupyter_text>Parameters varied significantly, but the AIC value has dropped sharply.  Step 9. Final Model calibration with well skin resistance Now we test if the skin resistance of the well has an impact on model calibration. For this, we add the ```res``` parameter in the calibration settings. We use the same multi-layer model.<jupyter_code>cR = Calibrate(M_nlay)
cR.set_parameter(name='kaq0_49', initial=10)
cR.set_parameter(name='Saq0_49', initial=1e-4, pmin=1e-7)
cR.set_parameter_by_reference(name='res', parameter = W_nlay.res, initial = 0.2, pmin = 0)
cR.series(name='obs', x=0, y=0, layer=[18,19,20], t=t, h=h)
cR.fit()
display(cR.parameters)
print('RMSE:', cR.rmse())<jupyter_output><empty_output><jupyter_text>The RMSE has improved slightly, and AIC has also decreased, which means adding skin resistance improved the model. However, the improvement is minimal, as the ```res``` calibrated value is tiny. Therefore, one can justify a calibration without ```res```.  Step 10. Analysis and comparison of simulated valuesWe now compare the values in TTim and add the results of the AQTESOLV modelling reported by Yang (2020).<jupyter_code>ta = pd.DataFrame(columns=['k [m/d]', 'Ss [1/m]', 'res [1/d]'],\
                 index=['AQTESOLV', 'ttim-three', 'ttim-multi', 'ttim-res'])
ta.loc['ttim-three'] = np.concatenate((ca.parameters['optimal'].values,[np.nan]))
ta.loc['ttim-multi'] = np.concatenate((cM.parameters['optimal'].values,[np.nan]))
ta.loc['ttim-res'] = cR.parameters['optimal'].values
ta.loc['AQTESOLV'] = [4.034, 3.834E-04]+[np.nan]
ta['RMSE'] = [0.002976, ca.rmse(), cM.rmse(), cR.rmse()]
ta.style.set_caption('Comparison of parameter values and error under different models')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14322682268226822, 'fraction_numerical': 0.048492349234923494, 'mean_word_length': 5.046938775510204, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 46, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '84001', 'license': 'MIT', 'n_tokens_mistral': 3503, 'n_tokens_neox': 3216, 'n_words': 1034, 'path': 'pumpingtest_benchmarks/slug1_pratt_county.ipynb', 'repo_name': 'vcantarella/ttim', 'size': '8888'}