<jupyter_start><jupyter_text>2.Test for an anisotropic water-table aquifer - Moench Example**This test is taken from examples presented in MLU tutorial.**  Introduction and Conceptual ModelThis test is based on a synthetic example reported by Barlow and Moench (1999), utilizing an analytical solution developed by Moench and Allen (1997) for the transient flow of partially-penetrating wells in unconfined aquifers. The data reported has been used in MLU (Carlson and Randall, 2012) to check the model performance.We will reproduce the work of Yang (2020) to explore the performance of anisotropic water table aquifer modelling with TTim and compare the results with the published values and the MLU solution.The conceptual model of the test is of an aquifer, partially saturated with water (10 m water table). A pumping well is screened from 5 to 10 m depth. The well and the well-casing radius is 0.1 m. Drawdown is recorded at the pumping well and four piezometers located at two different distances and two different depths. Two piezometers, PS1 and PS2, are located at one-meter depth below the water table and 3.16 and 31.6 m distance, respectively. Another two (PD1 and PD2) piezometers are at 7.5 m depth below the water table and the same distances,  directly below the previous piezometers. The figure below shows the location of the well and the piezometersPumping starts at time t = 0 at a constant rate of 172.8 m3/d. Drawdown is recorded until t = 3 days.<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,2), width = 70, height = 5, fc = 'b', zorder=0, alpha=0.1)
ax.add_patch(sky)

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



well = plt.Rectangle((-0.5,-10), width = 1, height = 12, fc = np.array([200,200,200])/255, zorder=1)
ax.add_patch(well)

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

#Screen for the well:
screen = plt.Rectangle((-0.5,-10), width = 1, height = 5, 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 = 1.1,y = 3, dx = 3, dy = 0, color = "#00035b")
ax.add_patch(pumping_arrow)
ax.text(x = 4.2, y = 3, s = r'$ Q = 172.8 \frac{m^3}{d}$', fontsize = 'large' )

#Piezometers
piez1 =plt.Rectangle((31.4,-10), width = 0.5, height = 12, fc = np.array([200,200,200])/255, zorder=1)
screen_piez_1 = plt.Rectangle((31.4,-8), width = 0.5, height = 1, fc = np.array([200,200,200])/255, alpha=1, zorder = 2, ec = "k", ls = '--')
screen_piez_1.set_linewidth(2)
screen_piez_2 = plt.Rectangle((31.4,-1.5), width = 0.5, height = 1, fc = np.array([200,200,200])/255, alpha=1, zorder = 2, ec = "k", ls = '--')
screen_piez_2.set_linewidth(2)

ax.add_patch(piez1)
ax.add_patch(screen_piez_1)
ax.add_patch(screen_piez_2)

piez2 =plt.Rectangle((2.75,-10), width = 0.5, height = 12, fc = np.array([200,200,200])/255, zorder=1)
screen_piez_3 = plt.Rectangle((2.75,-8), width = 0.5, height = 1, fc = np.array([200,200,200])/255, alpha=1, zorder = 2, ec = "k", ls = '--')
screen_piez_3.set_linewidth(2)
screen_piez_4 = plt.Rectangle((2.75,-1.5), width = 0.5, height = 1, fc = np.array([200,200,200])/255, alpha=1, zorder = 2, ec = "k", ls = '--')
screen_piez_4.set_linewidth(2)

ax.add_patch(piez2)
ax.add_patch(screen_piez_3)
ax.add_patch(screen_piez_4)

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

#Water table
line2 = plt.Line2D(xdata = [-200,1200], ydata = [0,0], color = 'b')
ax.add_line(line2)

ax.text(-9,0.5, s = 'Water Table', fontsize = 'large', color = 'b',bbox = {'fc' : 'w'})
ax.text(4.5, -1, s = 'PS1', bbox = {'fc' : 'w'})
ax.text(34, -1, s = 'PS2', bbox = {'fc' : 'w'})
ax.text(4.5, -7.5, s = 'PD1', bbox = {'fc' : 'w'})
ax.text(34, -7.5, s = 'PD2', bbox = {'fc' : 'w'})

ax.set_xlim([-10,50])
ax.set_ylim([-10,4])
ax.set_xlabel('Distance [m]')
ax.set_ylabel('Relative height [m]')
ax.set_title('Conceptual Model - Vennebulten Example');<jupyter_output><empty_output><jupyter_text>Step 1. Import required libraries<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from ttim import *
import pandas as pd<jupyter_output><empty_output><jupyter_text>Step 2. Set basic parameters<jupyter_code>b = 10 #aquifer thickness in m
Q = 172.8 #constant discharge rate in m^3/d
rw = 0.1 #well radius in m
rc = 0.1 #casing radius in m
r1 = 3.16 # distance of closer wells in m 
r2 = 31.6 # distance of wells more far away in m<jupyter_output><empty_output><jupyter_text>Step 3. Load datasets of observation wellsThe dataset for each well consists of a column with the time data in seconds and drawdown in meters. We are loading it and converting it to days and meters.<jupyter_code>data0 = np.loadtxt('data/moench_pumped.txt', skiprows=1)
t0 = data0[:, 0] / 60 / 60 / 24 #convert time from seconds to days
h0 = -data0[:, 1] #converting drawdown to heads
data1 = np.loadtxt('data/moench_ps1.txt', skiprows=1)
t1 = data1[:, 0] / 60 / 60 / 24 #convert time from seconds to days
h1 = -data1[:, 1]
data2 = np.loadtxt('data/moench_pd1.txt', skiprows=1)
t2 = data2[:, 0] / 60 / 60 / 24 #convert time from seconds to days
h2 = -data2[:, 1]
data3 = np.loadtxt('data/moench_ps2.txt', skiprows=1)
t3 = data3[:, 0] / 60 / 60 / 24 #convert time from seconds to days
h3 = -data3[:, 1]
data4 = np.loadtxt('data/moench_pd2.txt', skiprows=1)
t4 = data4[:, 0] / 60 / 60 / 24 #convert time from seconds to days
h4 = -data4[:, 1]<jupyter_output><empty_output><jupyter_text>Step 4. Creating a TTim modelWe will create an initial model divided in the same way as in the MLU documentation. A first layer with 0.1 m thick and phreatic storage, followed by a 2 m thick layer where the shallow piezometers are located, a 3 m layer and finally, a 5 m layer where the pump is placed and the last piezometers are screened. Additionally, we will set the model parameters with the given ones in Barlow and Moench (1999) and compare the results with the heads given in the paper.<jupyter_code>#Set kaq, Saq, Sy and kzoverkh as given in Moench (1997)
kaq = 1e-4 * 60 * 60 * 24 #convert from m/s to m/d
Sy = 0.2
Saq = 2e-5
zh = 0.5 #kzoverkh
ml1 = Model3D(kaq=kaq, z=[0, -0.1, -2.1, -5.1, -10.1], Saq=[Sy, Saq, Saq, Saq], \
              kzoverkh=zh, tmin=1e-5, tmax=3)
w1 = Well(ml1, xw=0, yw=0, rw=rw, rc=rc, tsandQ=[(0, Q)], layers=3)
ml1.solve()<jupyter_output>self.neq  1
solution complete<jupyter_text>Step 5. Check the TTim model with values obtained by Moench<jupyter_code>hm1 = ml1.head(r1, 0, t1, layers=1)[0]
hm2 = ml1.head(r1, 0, t2, layers=3)[0]
hm3 = ml1.head(r2, 0, t3, layers=1)[0]
hm4 = ml1.head(r2, 0, t4, layers=3)[0]
hm0 = ml1.head(0, 0, t0, layers=3)[0]
plt.figure(figsize=(8, 5))
plt.semilogx(t0, h0, '.', label='pumped')
plt.semilogx(t0, hm0, label='ttim pumped')
plt.semilogx(t1, h1, '.', label='PS1')
plt.semilogx(t1, hm1, label='ttim PS1')
plt.semilogx(t2, h2, '.', label='PD1')
plt.semilogx(t2, hm2, label='ttim PD1')
plt.semilogx(t3, h3, '.', label='PS2')
plt.semilogx(t3, hm3, label='ttim PS2')
plt.semilogx(t4, h4, '.', label='PD2')
plt.semilogx(t4, hm4, label='ttim PD2')
plt.xlabel('time [d]')
plt.ylabel('head [m]')
plt.title('Model Results - Simulation 1')
plt.legend();<jupyter_output><empty_output><jupyter_text>Visually, TTim's solution is in good agreement with the model, but we can see that the simulated values in the pumping well and PD1 are a little below the values obtained in the analytical solution.  Step 5.1. Calculate RMSE of TTim model<jupyter_code>rmse = np.sqrt(np.sum(np.sum((h1-hm1)**2)+np.sum((h2-hm2)**2)+np.sum((h3-hm3)**2)+np.sum((h4-hm4)**2)+np.sum((h0-hm0)**2))/\
    (len(h1)+len(h2)+len(h3)+len(h4)+len(h0)))

print('RMSE:', rmse)<jupyter_output>RMSE: 0.06131792435676353<jupyter_text>The model has obtained a very close result with a good approximation to the analytical solution.  Step 6. Model CalibrationNow, instead of using the given values, we will try to find them through the TTim optimization framework. One can learn more about the calibration procedure and how to set the parameters in the following notebook: [Unconfined Test - Vennebulten](unconfined1_vennebulten.ipynb)<jupyter_code>ml2 = Model3D(kaq=1, z=[0, -0.1, -2.1, -5.1, -10.1], Saq=[0.1, 1e-4, 1e-4, 1e-4], \
              kzoverkh=1, tmin=1e-5, tmax=3)
w2 = Well(ml2, xw=0, yw=0, rw=rw, rc=rc, tsandQ=[(0, Q)], layers=3)
ml2.solve()
ca2 = Calibrate(ml2)
ca2.set_parameter(name='kaq0_3', initial=1)
ca2.set_parameter(name='Saq0', initial=0.2)
ca2.set_parameter(name='Saq1_3', initial=1e-4, pmin=0)
ca2.set_parameter_by_reference(name='kzoverkh', parameter=ml2.aq.kzoverkh, \
                               initial=0.1, pmin=0)
ca2.series(name='pumped', x=0, y=0, t=t0, h=h0, layer=3)
ca2.series(name='PS1', x=r1, y=0, t=t1, h=h1, layer=1)
ca2.series(name='PD1', x=r1, y=0, t=t2, h=h2, layer=3)
ca2.series(name='PS2', x=r2, y=0, t=t3, h=h3, layer=1)
ca2.series(name='PD2', x=r2, y=0, t=t4, h=h4, layer=3)
ca2.fit()
display(ca2.parameters)
print('RMSE:', ca2.rmse())<jupyter_output><empty_output><jupyter_text>The values are close to the values in Moench.<jupyter_code>hm0_2 = ml2.head(0, 0, t0, layers=3)[0]
hm1_2 = ml2.head(r1, 0, t1, layers=1)[0]
hm2_2 = ml2.head(r1, 0, t2, layers=3)[0]
hm3_2 = ml2.head(r2, 0, t3, layers=1)[0]
hm4_2 = ml2.head(r2, 0, t4, layers=3)[0]
plt.figure(figsize=(8, 5))
plt.semilogx(t0, h0, '.', label='pumped')
plt.semilogx(t0, hm0_2, label='ttim pumped')
plt.semilogx(t1, h1, '.', label='PS1')
plt.semilogx(t1, hm1_2, label='ttim PS1')
plt.semilogx(t2, h2, '.', label='PD1')
plt.semilogx(t2, hm2_2, label='ttim PD1')
plt.semilogx(t3, h3, ',', label='PS2')
plt.semilogx(t3, hm3_2, label='ttim PS2')
plt.semilogx(t4, h4, '.', label='PD2')
plt.semilogx(t4, hm4_2, label='ttim PD2')
plt.xlabel('time [d]')
plt.ylabel('head [m]')
plt.title('Model Results - Calibration 1')
plt.legend();<jupyter_output><empty_output><jupyter_text>Step 7. Creating and calibrating a model with more layers The calibration resulted in a lower error than the simulated model, which could mean that we could improve the conceptualization of the model in TTim. In this next step, we will increase the number of layers to better account for the vertical flow component. We will create an 18-layers model with 0.5 m thick layers, except at the layers in the piezometers where a 1 m thickness has been established to make sure the piezometers lie at the centre of the layer:<jupyter_code>## Model Configuration:

nlay = 18
zlayers = (np.concatenate((np.array([0,-0.5,-1.5]),np.linspace(-2,-7,11),np.array([-8,-8.5,-9,-9.5,-10])))) #elevation of each layer
Saq_n = 1e-4 * np.ones(nlay)
Saq_n[0] = 0.1 # Setting the first storage as specific yield
print(zlayers)
print(Saq)
ml3 = Model3D(kaq=1, z=zlayers, Saq=Saq_n, \
              kzoverkh=1, tmin=1e-5, tmax=3)
w3 = Well(ml3, xw=0, yw=0, rw=rw, rc=rc, tsandQ=[(0, Q)], layers=range(9,18))
ml3.solve()
ca3 = Calibrate(ml3)
ca3.set_parameter(name='kaq0_19', initial=1, pmin=0)
ca3.set_parameter(name='Saq0', initial=0.2, pmin=0)
ca3.set_parameter(name='Saq1_19', initial=1e-4, pmin=0)
ca3.set_parameter_by_reference(name='kzoverkh', parameter=ml3.aq.kzoverkh, \
                               initial=0.1, pmin=0)
ca3.series(name='pumped', x=0, y=0, t=t0, h=h0, layer=15)
ca3.series(name='PS1', x=r1, y=0, t=t1, h=h1, layer=1)
ca3.series(name='PD1', x=r1, y=0, t=t2, h=h2, layer=13)
ca3.series(name='PS2', x=r2, y=0, t=t3, h=h3, layer=1)
ca3.series(name='PD2', x=r2, y=0, t=t4, h=h4, layer=13)
ca3.fit()
display(ca3.parameters)
print('RMSE:', ca3.rmse())
hm0_3 = ml3.head(0, 0, t0, layers=15)[0]
hm1_3 = ml3.head(r1, 0, t1, layers=1)[0]
hm2_3 = ml3.head(r1, 0, t2, layers=13)[0]
hm3_3 = ml3.head(r2, 0, t3, layers=1)[0]
hm4_3 = ml3.head(r2, 0, t4, layers=13)[0]
plt.figure(figsize=(8, 5))
plt.semilogx(t0, h0, '.', label='pumped')
plt.semilogx(t0, hm0_3, label='ttim pumped')
plt.semilogx(t1, h1, '.', label='PS1')
plt.semilogx(t1, hm1_3, label='ttim PS1')
plt.semilogx(t2, h2, '.', label='PD1')
plt.semilogx(t2, hm2_3, label='ttim PD1')
plt.semilogx(t3, h3, ',', label='PS2')
plt.semilogx(t3, hm3_3, label='ttim PS2')
plt.semilogx(t4, h4, '.', label='PD2')
plt.semilogx(t4, hm4_3, label='ttim PD2');<jupyter_output><empty_output><jupyter_text>Increasing the number of layers has not improved the model calibration performance. Let's see how the model behaves with the same parameters from the analytical solution:  Step 7.1. Simulating multi-layered model with Moench's parameters<jupyter_code>ml4 = Model3D(kaq=kaq, z=zlayers, Saq=[Sy]+list(np.repeat(Saq, nlay-1)), \
              kzoverkh=zh, tmin=1e-5, tmax=3)
w4 = Well(ml4, xw=0, yw=0, rw=rw, rc=rc, tsandQ=[(0, Q)], layers=range(9,18))
ml4.solve()
hm0_4 = ml4.head(0, 0, t0, layers=16)[0]
hm1_4 = ml4.head(r1, 0, t1, layers=1)[0]
hm2_4 = ml4.head(r1, 0, t2, layers=13)[0]
hm3_4 = ml4.head(r2, 0, t3, layers=1)[0]
hm4_4 = ml4.head(r2, 0, t4, layers=13)[0]
plt.figure(figsize=(10, 7))
plt.semilogx(t0, h0, '.', label='pumped')
plt.semilogx(t0, hm0_4, label='ttim pumped')
plt.semilogx(t1, h1, '.', label='PS1')
plt.semilogx(t1, hm1_4, label='ttim PS1')
plt.semilogx(t2, h2, '.', label='PD1')
plt.semilogx(t2, hm2_4, label='ttim PD1')
plt.semilogx(t3, h3, ',', label='PS2')
plt.semilogx(t3, hm3_4, label='ttim PS2')
plt.semilogx(t4, h4, '.', label='PD2')
plt.semilogx(t4, hm4_4, label='ttim PD2');<jupyter_output><empty_output><jupyter_text>Step 7.1.1. Checking RMSE performance improvement<jupyter_code>rmse_n = np.sqrt(np.sum(np.sum((h1-hm1_4)**2)+np.sum((h2-hm2_4)**2)+np.sum((h3-hm3_4)**2)+np.sum((h4-hm4_4)**2)+np.sum((h0-hm0_4)**2))/\
    (len(h1)+len(h2)+len(h3)+len(h4)+len(h0)))

print('RMSE:', rmse_n)<jupyter_output>RMSE: 0.013118477167630852<jupyter_text>Step 8. Summary of and analysis of calibrated values and model errors<jupyter_code>ta = pd.DataFrame(columns=['Moench','Moench - more layers' ,'TTim', 'TTim - more layers'],\
                 index=['k[m/d]', 'Sy[-]', 'Ss[1/m]', 'kz/kh'])
ta.loc[:, 'TTim - more layers'] = ca3.parameters['optimal'].values
ta.loc[:, 'TTim'] = ca2.parameters['optimal'].values
ta.loc[:, 'Moench'] = [kaq, Sy, Saq, zh]
ta.loc[:, 'Moench - more layers'] = [kaq,Sy, Saq, zh]
ta.loc['RMSE'] = [rmse,rmse_n, ca2.rmse(), ca3.rmse()]
ta<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17457044673539518, 'fraction_numerical': 0.07711340206185567, 'mean_word_length': 4.965969659696597, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 50, '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': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '217057', 'license': 'MIT', 'n_tokens_mistral': 6636, 'n_tokens_neox': 6147, 'n_words': 1574, 'path': 'pumpingtest_benchmarks/unconfined2_moench.ipynb', 'repo_name': 'vcantarella/ttim', 'size': '14550'}