Due 4/25/17
5 pts
Please submit your assignment as an html export, and for written responses, please type them in a cell that is of type Markdown. The final part of the exercise involves drawing a flow net by hand (actually, you could tackle this part of the assignment first). For that, you may turn in a hard copy of your answer.
In [1]:
# Import numerical tools
import numpy as np
#Import pandas for reading in and managing data
import pandas as pd
# Import pyplot for plotting
import matplotlib.pyplot as plt
#Import seaborn (useful for plotting)
import seaborn as sns
# Magic function to make matplotlib inline; other style specs must come AFTER
%matplotlib inline
%config InlineBackend.figure_formats = {'svg',}
This investigation of groundwater flow and drawdown in wells is based on the lawsuit described in the book and movie A Civil Action (a true story). For the background behind this story, read the Wikipedia page. Then look at the map and animation of the study site here. Wells G and H shown on the map are the municipal wells from which Woburn’s water supply was withdrawn.
In this assignment you will be investigating some of the very problems that the hydrologists hired to the case worked on. As usual, you may work with another person, but each person must turn in separate assignments and specify whom they worked with. Below I provide additional information specific to this assignment.
One of the important questions raised during the trial centered on the drawdown of the water table when wells G and H (see map linked above) were operating and the effects of the Riley tannery well (also included in the model). If, when wells G and/or H operated, the water table was drawn down on the east side of the Aberjona River, then contaminants from west side of the river valley (Beatrice 15 acres and Olympia Trucking) could flow under the river to wells G and H. If the water table was not drawn down on the east side of the river when wells G, H, and the Riley Tannery well operated, then it is unlikely that contaminants on the west side of the river valley could flow to wells G and H. Thus, the interaction of these three wells with the Aberjona River and the groundwater flow system was critical to determining which of the contaminated properties contributed TCE and PCE to wells G and H.
Proving the connection between the contaminants migrating from the defendants’ properties to wells G and H during the first phase of the trial was a requirement to proceed to the next phase. To make arguments regarding the potential sources of TCE to wells G and H, the plaintiffs' expert had to define the subsurface geology, the surface and subsurface, and show that TCE was used at both the Riley/Beatrice and W.R. Grace properties. Counsel for the plaintiffs was able to identify Grace employees involved with on-site TCE disposal activities. However, the plaintiffs' counsel was not able to directly link disposal of wastes on the Riley 15-acre property with activities at the tannery itself. The ATSDR report, which was prepared three years after the trial, is a good illustration of how difficult it is to conclusively assess risk from a specific contaminant source at the Woburn site.
The U.S. Geological Survey first performed a seven-hour pumping test to characterize aquifer properties, and then, in December 1985-January 1986, a 30-day pumping test to evaluate the effects of pumping on the distribution of hydraulic head. This was the first time wells G and H operated since their closure in May 1979. It was also the first time synoptic sets of water levels had ever been measured in the network of monitoring wells surrounding wells G and H. During the 30-day pumping test, wells G and H operated at their average rates of 700 gpm (well G) and 400 gpm (well H). However, pumping records kept by the City of Woburn showed that between 1964 and 1979 the wells rarely operated together and that the two municipal wells were frequently not in use for months at a time. As a result of changes in the pumping rates of wells G and H and their periodic (discontinuous) use, the water table in the buried valley aquifer was dynamic, rising and falling as wells G and H were turned on and off and pumped heavily or lightly. These dynamic changes in water levels result in the groundwater flow system being transient in character – as opposed to being steady state in character.
Under steady-state conditions, water levels at any location in the flow system do not change with time, which results in hydraulic gradients and flow velocities at any location not changing with time, and therein no net change in the amount of groundwater in storage. Under transient conditions, water levels at any location in the flow system change with time, which causes hydraulic gradients, flow velocities, and the amount of groundwater in storage (water levels) to change with time. In this exercise we assume that the wells G and H have been pumping together for a long enough time that the water table is stationary and flow is steady state. This is an assumption that we recognize is not true over periods of time greater than a few months. This simplistic assumption lets us learn about drawdown created by a single pumping well, the formation of a cone of depression, superposition of drawdowns from one or more pumping wells, and the creation of groundwater divides between pumping wells, all of which are important to understanding the responses of the groundwater flow system to the pumping history of wells G and H.
The operation of the Riley well was estimated by de Lima and Olimpio (USGS, 1989) to operate at an average pumping rate of 200 gallons per minute, continuously. This was based on the wastewater discharge emanating from the tannery; actual pumping rates were not available. It is believed that pumping volumes fluctuated depending on the number of shifts working, the volume of leather orders, actual type of leather being processed, etc.
Helpful commands for linear regression: np.polyfit(), np.polyval().
See an example of their usage in the Feb 2 tutorial, found under Files/Materials on bCouses.
In [2]:
# The data
time_minutes = np.arange(170,570,10) #Times (minutes) at which drawdown was observed after the start of pumping
s = np.array([0.110256896, 0.122567503, 0.180480293, 0.214489433, 0.356304352, 0.554603882, 0.49240574, 0.524209467, 0.562727644, 0.754849464, 0.718713002, 0.752910019, 0.73903658, 0.89009263, 0.967712464, 0.910807162, 0.986238396, 1.042178187, 1.081114186, 1.080825045, 1.196565491, 1.264971986, 1.430805272, 1.377858223, 1.376787182, 1.340970634, 1.466832052, 1.528405086, 1.44136224, 1.610936499, 1.503519725, 1.55398487, 1.628028833, 1.675649643, 1.672772239, 1.730410501, 1.730935188, 1.756850444, 1.731143013, 1.818924173])#Drawdown in observation well, ft
With knowledge of the hydraulic conductivity, you can now solve the Thiem equation for the equilibrium potentiometric surface under different combinations of well pumping rates. For this part of the assignment, you will be working with a module that solves the Thiem equation over a two-dimensional spatial domain.
The module requires the user to input several parameters, as described in the red text. When you first run it, use the hydraulic conductivity that you just solved for. Assume that the Riley well (QR) pumps at an average rate of 200 gallons per minute (gpm), Well G (QG) pumps at an average rate of 700 gpm, and Well H (QH) pumps at an average rate of 400 gpm. When you run the module, you will see two plots. The first plot is a filled contour plot that depicts the drawdown relative to non-pumping conditions. The value shown is the composite amount of drawdown related to pumping from each well (drawdowns are additive). Note how the values increase in the cells immediately around the wells and decrease toward the margins of the graph.
The second plot is a cross section that shows the drawdown associated with the north-south transect that passes through wells G and H under steady state conditions.
2. Take a look at the model below. What is a key assumption that is being made in the implementation of the Theim equation, beyond the assumptions we talked about in class? Hint: compare the formulas under the Calculate Drawdown for Each Well section to Eq. 5.49 in the Fetter handout. [1/3 pt.]
In [31]:
def PlotWoburnDD(K,b,QG,QH,QR, returnval=0):
"""
This routine uses the Thiem equation for unCONFINED aquifers to generate a plot of drawdown
contours around the Aberjona River in Woburn, Massachusetts. The Riley well is the source of
the contamination. Wells G and H are wells for the town's municipal water supply. Users can
evaluate how the potentiometric surface changes with changes in hydraulic conductivity (K),
saturated thickness (b), and well pumping rates (QG for well G, QH for well H, and QR for the
Riley well). Note that K is in units of ft/day, b is in units of ft, and QG, QH, and QR are
in units of gallons per minute. returnval is an optional argument. If you set it equal to 1,
the function will return the value of drawdown at well G. If unspecified, the function returns
just the plots.
"""
#CONVERT GALLONS PER MINUTE TO CUBIC FEET PER DAY
QG = QG*192.5
QH = QH*192.5
QR = QR*192.5
#OTHER CONSTANTS
delx = 100 #Cell size in feet
r0 = 2700 #You figure this out! (It is in feet)
#SET UP THE EXPERIMENTAL GRID
nrows = 25 #number of rows
ncols = 27 #number of columns
rowvect = np.arange(nrows) #A vector of row coordinates
colvect = np.arange(ncols) #A vector of column coordinates
colcoords, rowcoords = np.meshgrid(colvect,rowvect) #Creates two matrices. In colcoords, each
#entry is the column index of the cell (i.e., point in space). In rowcoords, each entry is the
#row index of the cell (i.e., point in space).
#SPECIFY WELL LOCATIONS AND RIVER COORDINATES WITHIN THE EXPERIMENTAL GRID
loc_R = np.array([20,8]) #Index locations for the Riley well (row, column)
loc_H = np.array([10,13]) #Index locations for H well
loc_G = np.array([17,13]) #Index locations for G well
river_row = np.array([0, 1, 2, 3, 4, 5, 6, 6.4, 7, 8, 9, 9.9, 10, 11, 11.5, 12, 13, 14, 15, 16, 16.1, 17, 17.5, 18, 18.4, 19, 19.6, 20, 20.3, 21, 21.7, 22, 22.4, 23, 23.1, 23.9, 24, 24.8, 25])
river_col = np.array([8.6, 9, 9.1, 9.3, 9.5, 9.6, 9.9, 10, 10.3, 11, 11.7, 12, 12.1, 12.7, 12.9, 12.8, 12.6, 12, 11.6, 11.1, 11, 10.8, 10.7, 10.8, 11, 11.4, 12, 12.6, 13, 14, 15, 15.5, 16, 16.9, 17, 18, 18.1, 19, 19.2])
#CALCULATE DISTANCES BETWEEN EACH WELL AND EVERY OTHER CELL
d_R = np.sqrt((rowcoords-loc_R[0])**2+(colcoords-loc_R[1])**2)*delx #Solve for distance (feet) using the Pythagorean theorem
d_H = np.sqrt((rowcoords-loc_H[0])**2+(colcoords-loc_H[1])**2)*delx #Solve for distance (feet) using the Pythagorean theorem
d_G = np.sqrt((rowcoords-loc_G[0])**2+(colcoords-loc_G[1])**2)*delx #Solve for distance (feet) using the Pythagorean theorem
#SET DISTANCES TO WELL IN CELLS WITH WELL TO 25 FT (TO AVOID SINGULARITY AT 0)
d_R[20,8] = 25 #feet
d_H[10,13] = 25 #feet
d_G[17,13] = 25 #feet
#CALCULATE DRAWDOWN FROM EACH WELL
s_R = b-np.sqrt(b**2-QR*np.log(r0/d_R)/(np.pi*K)) #Drawdown from the Riley well
s_H = b-np.sqrt(b**2-QH*np.log(r0/d_H)/(np.pi*K)) #Drawdown from well H
s_G = b-np.sqrt(b**2-QG*np.log(r0/d_G)/(np.pi*K)) #Drawadown from well G
s_total = s_R+s_H+s_G #Combined drawdown from all of the wells, feet
#NOW GENERATE PLOTS
plt.figure(figsize=(10,10))
s_plot = plt.contourf(np.transpose((colcoords+0.5)*delx), np.transpose((rowcoords+0.5)*delx), np.transpose(s_total), cmap=plt.cm.plasma)
plt.gca().invert_yaxis() #This puts the origin of the plot on the upper left
cb = plt.colorbar(s_plot, orientation='horizontal')
cb.set_label('Drawdown, ft')
riv_plot = plt.plot(river_col*delx, river_row*delx, 'b-',linewidth=3.0)
wH = plt.plot((loc_H[1]+0.5)*delx, (loc_H[0]+0.5)*delx, 'ko')
ax = plt.gca()
ax.annotate('H', xy = ((loc_H[1]+0.5)*delx, (loc_H[0]+0.5)*delx))
wG = plt.plot((loc_G[1]+0.5)*delx, (loc_G[0]+0.5)*delx, 'ko')
wR = plt.plot((loc_R[1]+0.5)*delx, (loc_R[0]+0.5)*delx, 'ko')
ax.annotate('G', xy = ((loc_G[1]+0.5)*delx, (loc_G[0]+0.5)*delx))
ax.annotate('Riley', xy = ((loc_R[1]+0.5)*delx, (loc_R[0]+0.5)*delx))
#plt.axis('equal') #Uncomment this to make the x- and y- axis display on the same scale.
plt.title('Map of Aberjona River and Drawdown Contours')
plt.show()
plt.figure()
plt.fill_between(rowvect*delx, max(np.ceil(s_total[:,13])), s_total[:,13])
plt.plot([loc_H[0]*delx, loc_H[0]*delx], [min(s_total[:,13]), max(s_total[:,13])], 'r')
plt.plot([loc_G[0]*delx, loc_G[0]*delx], [min(s_total[:,13]), max(s_total[:,13])], 'r')
plt.gca().invert_yaxis()
plt.xlabel('y distance, feet')
plt.ylabel('Drawdown, feet')
plt.title('Cross-section across column with wells G and H')
if returnval==1:
return s_total[17,13]
3. Now run the model! (You might want to check back to the Jan 31 class tutorial for a refresher on how to run modules/functions.) Using the pumping rates and aquifer thickness specified above, and the hydraulic conductivity that you solved for (in feet per day), how do pumping from the Riley well and Well H affect the shape of the drawdown contours? Why does the contour shape become uniformly round as you move away from the wells? [1/3 pt.]
4. Set the pumping rates for wells G and H to 0. How does this change the shape of drawdowns on the map and the contours? How would this pumping condition affect the direction and rate of movement of a contaminant plume? [1/3 pt.]
5. Now reset the pumping rate for well G to 700 gpm and keep the Riley well value at 200 gpm.
a. How does this pumping configuration change the shape of drawdowns on the cross section and contour map? [0.25 pts.]
b. What is the net change in drawdown in well G when the value used in question 3 (0 GPM) is compared to that used for question 4 (700 gpm)? Note if you want the figures from both runs of the module to show up together, you need to insert the line plt.figure() between your first call to PlotWoburnDD and your second call. [0.5 pts.]
c. Could the new pumping conditions (well G = 700 gpm) pull TCE-contaminated water from the Riley property to well G? Explain why or why not. [0.25 pts.]
6. How sensitive are your results to hydraulic conductivity? To aquifer thickness? In other words, how does the potentiometric surface change when you use higher or lower values of these parameters? Write your answer in paragraph format. Hint: If you generate several plots in a row, you might want to separate them with print statements so that you will know which is which. [0.5 pts.]
7. Although this Thiem module can be an effective tool for examining hypothetical pumping rates and groundwater conditions, there are some severe limitations to using this simplistic treatment in a trial. What are these limitations? Hint: review the assumptions of the Thiem equation and assess the geology/geography of the area, carefully considering the role of the river and the selection of the model’s domain. [0.25 pts.]
8. Numerical groundwater flow models commonly are used today to address many of the limitations assumptions of Thiem Equation as applied to the flow system at Woburn (see MODFLOW information). Why would a numerical groundwater flow model, such as the USGS model, MODFLOW, which was in development during the trial, have been advantageous to the experts? [0.25 pts]
9. Let's say you are not allowed to access any data about well pumping rates. However, you are able to make measurements in piezometers around the well of interest, which generates a very nice contour plot that looks like this:
Assume that K = 100 ft/day, and the saturated thickness of the aquifer away from the cone of depression is 80 ft. Drawing a flow net by hand, estimate the volumetric pumping rate of the well. (It may be difficult to do this precisely, so you also might want to describe what your intent is.) The person whose estimate is closest gets extra credit points! [1 pt, + 1 pt E.C.]
In [ ]: