Assignment steps
Separate steps to successfully carry out the assignments.
Steps to take
The assignment can be broken up into several pieces from which to build the result.
1 A well near a river 1 Well with varying monthly extraction 1 Infiltration from the river
A well near a river
This means that the river keeps the head at its banks constant. This situation can be modeled using the well and its mirror well. So you need the real well and its mirror well on the other side and at the same distance to the riverbank.
Well with varying monthly extraction
As always this is modeled by means of superpostion in time. At each new month a new well at the same location is started with the difference of flow compared to the previous one and the results are added to the effect of all the previous well for the times after the moment that the new starts. The best way is to use a vector of times, say t, at which you want to calculate the result. This may be one point per day of 13 points per hour our whatever you find enough to get a nice graph with sufficient detail. Then you need a vector of times at which the extaction rate changes, that is, at which a new well starts. This vector may be called change-times. Change times have nothing to do with the times you calculate the results, they are independent.
Infiltration from the river
The flow towards any well in a confined or unconfined aquifer at a radius $r$ from the well equals the extraction rate divided by the surface are of the aquifer with a radius $r$: $qr = Q_r / (2 \pi r D)$. Then you need to compute the component perpendicular to the riverbank to get the infiltration rate. Don't forget to add the contribution of the mirror well.
Advice: first do it for steady-state situation, so that $Q_r = Q_0$. Then replace the $Q_r$ by the time-dependent value, which can be found in the syllabus.
Compute the total inflow from the river (by integrating along the river bank)
The easiest way to use coordinate here is, of course, to take $x=0$ at the river bank and let $x$ be the distance perpendicular to the river bank. The x-axis is perpendicular to the river bank. Then $y$ is along the river bank. Take $y=0$ at the point where the x-axis intersects the riverbank. Then divide the y-axis between say -a and +a into a fairly large number of pieces dy. Compute the infiltration velocity at these points. Multiply by the length dy and sum over all dy to get the total inflow.
Show this graphically.
The extraction from the well
Is seems most convenient to set up an Excel spreadsheet with the extraction in every month over e number of years. This will be table that has columns, date and Q, and possibly other columns which we may need to set up the spreadsheet but don't need here in Python.
Then we import the data using Pandas
In [6]:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.special import exp1
In [7]:
def newfig(title='forgot title?', xlabel='forgot xalbel?', ylabel='forgot ylabel?', xlim=None, ylim=None,
xscale='linear', yscale='linear', size_inches=(12, 6)):
fig, ax = plt.subplots()
fig.set_size_inches(size_inches)
ax.set_title(title)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_xscale(xscale)
ax.set_yscale(yscale)
if not xlim is None: ax.set_xlim(xlim)
if not ylim is None: ax.set_ylim(ylim)
ax.grid()
return ax
#newfig()
In [8]:
ayear = 365.25 # Convenence parameter, to convert between days and years
# Aquifer propertise
kD = 600 # m2/d
S = 0.2
In [9]:
# Demo of how to read data from an Excel workbook
# A workbook is convenient for preparing data and then reading
# it into Python (as a Pandas DataFrame sort of table) at once
welldata = pd.read_excel('PumpingStationQ.xlsx', sheet_name='Sheet1', parse_dates=True, dayfirst=True, index_col=0)
welldata.columns
Out[9]:
In [165]:
subttl = f'\nkD={kD:.0f} m2/d, S={S:4g}'
# Generate two axes for our visualization:
ax0 = newfig(title='Extraction from well' + subttl, xlabel='time [y]', ylabel='Q [m3/d]')
ax1 = newfig(title='head due to extraction of one well',
xlabel='time [y]', ylabel='s [m]')
xw, yw = -300., 0. # The well coordinates
xp, yp = 0., 0. # The coordinates of the observation point
r = np.sqrt((xw - xp) ** 2 + (yw - yp) ** 2) # Distance between well and observation point
s = np.zeros_like(time) # Initialize the head to all zeros
Q = np.zeros_like(time) # Same for the discharge
# Loop over all wells in the well data DataFrame
for iw in range(len(welldata)): # iw is the line number in the DataFrame
dQ, tch = welldata.iloc[iw][['dQ', 'tch']] # From this line select flow change dQ and change time tch
u = r ** 2 * S / (4 * kD * (time[time > tch] - tch))
ds = dQ/(4 * np.pi * kD) * exp1(u) # head change due to this single dQ
Q[time > tch] += dQ # Add this dQ to all times > tch
ax1.plot(time[time > tch] / ayear, ds) # plot ds versus time for t > tch
s[time > tch] += ds # Add this contribution to the total hean change
ax0.plot(time / ayear, Q, 'r', lw=5, alpha=0.3, label='total') # plot total Discharge
ax1.plot(time / ayear, s, 'r', lw=5, alpha=0.3, label='total') # Plot overall head change
Out[165]:
In [194]:
subttl = f'\nkD={kD:.0f} m2/d, S={S:4g}'
ax0 = newfig(title='Extraction from well' + subttl, xlabel='time [y]',
ylabel='Q [m3/d] infiltration -->')
ax1 = newfig(title='head due to extraction of one well',
xlabel='time [y]', ylabel='s [m]')
x1, y1 = -200, 0
x2, y2 = 200, 0
xp, yp = 100, 100
r1 = np.sqrt((x1 - xp) ** 2 + (y1 - yp) ** 2)
r2 = np.sqrt((x2 - xp) ** 2 + (y2 - yp) ** 2)
s = np.zeros_like(time)
Q = np.zeros_like(time)
for iw in range(len(welldata)):
dQ, tch = welldata.iloc[iw][['dQ', 'tch']]
n = np.sum(time > tch)
u1 = r1 ** 2 * S / (4 * kD * (time[time > tch] - tch))
u2 = r2 ** 2 * S / (4 * kD * (time[time > tch] - tch))
Q[time > tch] += dQ
ds = -dQ/(4 * np.pi * kD) * (exp1(u1) - exp1(u2))
#ax1.plot(time[time > tch] / ayear, ds)
s[time > tch] += ds
ax0.plot(time / ayear, Q, 'r', lw=5, alpha=0.3, label='total')
ax1.plot(time / ayear, s, 'r', lw=5, alpha=0.3, label='total')
ax0.legend()
ax1.legend()
Out[194]:
In [205]:
subttl = f'\nkD={kD:.0f} m2/d, S={S:4g}, xw = {xw:.0f} m'
ax0 = newfig(title='Extraction from well' + subttl, xlabel='time [y]',
ylabel='Q [m3/d] infiltration -->')
ax1 = newfig(title='Infiltration [m2/d] at y=0' + subttl, xlabel='time [y]',
ylabel='Qin [m2/d] --> in')
x1, y1 = -300, 0
x2, y2 = -x1, 0
xp, yp = 0, 0
r1 = np.sqrt((x1 - xp) ** 2 + (y1 - yp) ** 2)
r2 = np.sqrt((x2 - xp) ** 2 + (y2 - yp) ** 2)
Qin = np.zeros_like(time)
Q = np.zeros_like(time)
for iw in range(len(welldata)):
dQ, tch = welldata.iloc[iw][['dQ', 'tch']]
n = np.sum(time > tch)
u1 = r1 ** 2 * S / (4 * kD * (time[time > tch] - tch))
u2 = r2 ** 2 * S / (4 * kD * (time[time > tch] - tch))
Q[time > tch] += dQ
dQin = -dQ * (np.exp(-u1)/ (2 * np.pi * r1) + np.exp(-u2) / (2 * np.pi * r2))
#ax1.plot(time[time > tch] / ayear, dQin)
Qin[time > tch] += dQin
ax0.plot(time / ayear, Q, 'r', lw=5, alpha=0.3, label='total')
ax1.plot(time / ayear, Qin, 'r', lw=5, alpha=0.3, label='total')
ax2 = ax1.twiny()
ax2.plot(time / ayear, Q/np.abs(Q), 'g', lw=5, alpha=0.3, label='season')
ax3 = ax0.twiny()
ax3.plot(time / ayear, Q/np.abs(Q), 'g', lw=5, alpha=0.3, label='season')
ax0.legend()
ax1.legend()
Out[205]:
In [201]:
secaxyaxis
Out[201]:
In [ ]:
ax1.twiny(