### Background

Analytical solutions for groundwater flow in pumped aquifers are convenient for simple, quick calculations that do not warrant the complexity of a discretized numerical model. A variety of solutions (e.g., those of Theis, Theim, Hantush, etc.) familiar to hydrogeologists exist to represent wells in single, confined aquifers or leaky aquifers under steady-state as well as transient pumping conditions. Oftentimes, these solutions can be superimposed to account for the influences of multiple pumping wells, regional groundwater flow, and other effects.

Hemker, 1984 (Steady groundwater flow in leaky multiple-aquifer systems, Journal of Hydrology, 72, 355-374) developed an analytical solution for steady-state flow resulting from pumping in a system with an arbitrary number of confined aquifers (horizontal flow only), separated by aquitards (vertical flow only). Pumping can occur within any of the aquifer units. The solution is based upon extracting a set of eigenvectors from a conductance-leakance-term matrix and employs a Bessel function, so it can be unwieldy to use without writing some sort of code to execute its various components. This, along with the relative scarcity of practical field problems involving many stacked aquifers-aquitards typically encountered by hydrogeologists, probably explains why this particular model is not used widely (just my opinion here). Nonetheless, I recently came upon a problem that is very well suited to the Hemker model and needed to write a script to generate some solutions. I chose python for this task because of the ease of implementation (this did not involve heavy-duty computation), the readily available math tools from numpy and scipy, and some plotting capabilities in matplotlib that allow simple visualizations of the model systems and the results.

### Details

The problem involved injection of wastewater into a deep aquifer, some 2,000 feet beneath a potable water aquifer. Between the two are multiple high- and low-permeability sequences of sedimentary materials (e.g., sands/sandstones versus clays/shales), as indicated by geophysical logs associated with the wastewater injection wells. The question to be addressed was the time scale required for the injected wastewater to be driven upwards into the potable water aquifer, assuming that all aquifer-aquitard units in the sequence are laterally continuous.

The conceptual model used to answer this question consists of the following:

- A random, alternating aquifer-aquitard sequence, generated by random placement of contacts between aquifer and aquitards along the 2,000 ft-long sedimentary sequence. The lithology model is assumed to be binary; only “aquitards” and “aquifers” exist in the model.
- A uniform (horizontal) hydraulic conductivity is assigned to all aquifer materials, and a separate (vertical) hydraulic conductivity is assigned to all aquitards.
- Injection takes place within the lowest aquifer unit in the sequence (placed adjacent to and immediately below the 2,000-ft sequence).
- Output consists of calculated resultant changes in hydraulic head within each aquifer zone (i.e., a function of radial distance from injection well as well as position within vertical aquifer sequence).
- Head differences between aquifer layers – especially close to the injection well in the radial direction, where the head perturbation is greatest – determine the travel time for upward migration of water through the aquitards, assuming a porosity (e.g., 20 percent).

### Code

The python code used to implement Hemker’s (1984) model is listed below. The script includes classes to represent aquifer, aquitard, and well properties, some numpy and scipy functions to handle the math details, and both pandas and matplotlib libraries to process the model output. The code is configured to run in two modes: one to reproduce some of the results from Hemker’s (1984) paper (for verification), and the second to generate results for the problem described under the details section:

############################################################################### # # hemker.py # # a python script for modeling steady-state drawdown in response to pumping # in a multi-aquifer system based on Hemker (1984) analytical solution # # by Walt McNab # ############################################################################### from numpy import * from scipy.special import * from scipy.spatial.distance import * from pandas import * import matplotlib.pyplot as plt import matplotlib.patches as mpatches ################################ # # hydraulic properties classes # ################################ class HydraulicUnits: def __init__(self, K, b): self.K = K # hydraulic conductivity (horizontal or vertical) self.b = b # vertical thickness class Aquifers(HydraulicUnits): def __init__(self, K, b, z): HydraulicUnits.__init__(self, K, b) self.T = self.K * self.b # transmissivity self.z = z # depth midpoint of aquifer self.L = zeros(len(K)) # aquifer leakance def Leakance(self, w): self.L = sqrt(1.0/w) class Aquitards(HydraulicUnits): def __init__(self, K, b): HydraulicUnits.__init__(self, K, b) self.c = b/K # vertical hydraulic resistance class Well: def __init__(self, x, y, Q_tot, screened_list, aquifer): self.x = x # well location self.y = y self.Q_tot = Q_tot # total injection rate self.screened_list = array(screened_list) # screened flag self.Q = self.DistributeQ(aquifer) def DistributeQ(self,aquifer): # distribute Q among aquifers, based on aquifer transmissivity T_tot = (aquifer.T * self.screened_list).sum() Q = self.Q_tot * self.screened_list * aquifer.T/T_tot return Q def HeadChange(self, x, y, k, num_aqf, aquifers, V, Z): # head change in layer k resulting from pumping this well r = pdist([[self.x, self.y], [x, y]])[0] sum_term = zeros(num_aqf, float) for j in xrange(num_aqf): sum_term[j] = (V[k, :] * Z[j, :] * k0(r/aquifers.L)).sum() return (self.Q * sum_term/(2.0*pi*aquifers.T)).sum() ################################## # # Hemker model functions # ################################## def HemkerMatrix(aquifers, aquitards): # construct the Hemker model matrix a = 1.0/(aquifers.T * aquitards.c[:-1]) b = 1.0/(aquifers.T * aquitards.c[1:]) A = diag(a+b) + diag(-b[:-1], 1) + diag(-a[1:], -1) return A def Eigen(A): # return eigenvalues and eigenvectors of matrix A and its transpose w, V = linalg.eig(A) Z = transpose(linalg.inv(V)) return w, V, Z ################################## # # support functions # ################################## def GetLayers(b_aqf, b_aqt): # for user-specified aquifer-aquitard stratigraphy, process geometry section = b_aqf.sum() + b_aqt.sum() z_aqf = section - (b_aqt[:-1].cumsum() + b_aqf.cumsum() - 0.5*b_aqf) t_aqf = z_aqf + 0.5*b_aqf # aquifer top elevation t_aqt = t_aqf + b_aqt[:-1] # aquitard top elevation mixed_top = -concatenate((t_aqt, t_aqf)) tops = -sort(mixed_top) tops = append(tops, b_aqt[-1]) return z_aqf, tops def Strat(section, num_contacts, K_aqf0, K_aqt0, K_inj, b_inj): # create strata rand = random.random(num_contacts) * section contact = sort(rand) tops = [] # top of unit; used later for plotting K_aqf = [] b_aqf = [] z = [] K_aqt = [] b_aqt = [] for i in xrange(len(contact)): # stepping downward with depth if not i: # top layer --> aquitard, by definition K_aqt.append(K_aqt0) b_aqt.append(contact[i]) tops.append(section) else: if (i % 2 == 0): # even: aquitard K_aqt.append(K_aqt0) b_aqt.append(contact[i] - contact[i-1]) else: # odd: aquifer K_aqf.append(K_aqf0) b_aqf.append(contact[i] - contact[i-1]) z.append(section - (0.5*contact[i-1] + 0.5*contact[i])) tops.append(section - contact[i-1]) # append last zone in section (aquitard, by definition) K_aqt.append(K_aqt0) b_aqt.append(section - contact[-1]) tops.append(section - contact[-1]) # append injection zone, plus ficticious underlying aquitard K_aqf.append(K_inj) b_aqf.append(b_inj) z.append(-0.5*b_inj) tops.append(0.) K_aqt.append(K_aqt0) b_aqt.append(1000.) tops.append(-b_inj) tops.append(-b_inj - 0.01) # to set plot properly # create hydraulic unit objects aquifers = Aquifers(array(K_aqf), array(b_aqf), array(z)) aquitards = Aquitards(array(K_aqt), array(b_aqt)) return aquifers, aquitards, tops def PlotSection(tops): # draw vertical section showing aquifer and aquitard delineations plt.figure(0) clr = ['Peru', 'DeepSkyBlue'] # colors for aquitard, aquifer for (i, layer) in enumerate(tops): plt.bar(1., layer, color = clr[i%2]) f = plt.gca() f.axes.get_xaxis().set_visible(False) plt.ylabel('Elev. Above Datum') aquifer_patch = mpatches.Patch(color='DeepSkyBlue', label='Aquifer') aquitard_patch = mpatches.Patch(color='Peru', label='Aquitard') plt.legend(handles=[aquifer_patch, aquitard_patch]) plt.show() def ProcessVertProfile(x, y, num_aqf, aquifer_layer, aquifers, well, V, Z): # draw vertical head change profile plt.figure(1) dh = zeros(num_aqf, float) for layer in aquifer_layer: for pump in well: dh[layer] = +pump.HeadChange(x, y, layer, num_aqf, aquifers, V, Z) vertical = {'layer': aquifer_layer, 'z': aquifers.z, 'b': aquifers.b, 'dh': dh} vertical_df = DataFrame(vertical) vertical_df.to_csv('vertical.csv') plt.semilogx(vertical_df['dh'], vertical_df['z']) plt.xlabel('Head Change') plt.ylabel('Elev. Above Datum') plt.show() def PlotLayerProfiles(profiles_df, aquifer_layer, x): # draw radial drawdown profiles for all aquifer layers plt.figure(1) for layer in aquifer_layer: subset_df = profiles_df[profiles_df['layer']==layer] plt.plot(subset_df['x'], subset_df['dh'], label='aquifer ' + str(layer)) plt.legend(loc='lower right') plt.xlim(0., x.max()) plt.xlabel('Radial Distance') plt.ylabel('Head Change') plt.show() def ProcessLayerProfiles(num_aqf, aquifer_layer, aquifers, well, V, Z): # process head changes along x-transect; write to file & return data frame x = linspace(100., 5000., 50) # example transect y = 0. for (i, xp) in enumerate(x): dh = zeros(num_aqf, float) for layer in aquifer_layer: for pump in well: dh[layer] = +pump.HeadChange(xp, y, layer, num_aqf, aquifers, V, Z) profiles = {'layer': aquifer_layer, 'x': zeros(num_aqf)+xp, 'y': zeros(num_aqf)+y, 'z': aquifers.z, 'b': aquifers.b, 'dh': dh} if not i: profiles_df = DataFrame(profiles) else: new_profiles_df = DataFrame(profiles) profiles_df = concat([profiles_df, new_profiles_df], axis=0) profiles_df.to_csv('profiles.csv') return profiles_df, x ################################### # # script # ################################### def Hemker(mode): # specifcy aquifer-aquitard structure if mode==0: # example strata K_aqf = array([40., 50., 33.333, 40.]) b_aqf = array([50., 30., 15., 50.]) K_aqt = array([0.01, 0.01, 0.01, 0.01, 0.01]) b_aqt = array([10., 15., 10., 40., 200.]) z_aqf, tops = GetLayers(b_aqf, b_aqt) aquifers = Aquifers(K_aqf, b_aqf, z_aqf) aquitards = Aquitards(K_aqt, b_aqt) else: # create stochastic strata (example parameters) section = 2000. # total thickness num_contacts = 50 # must be an even number K_aqf0 = 30. K_aqt0 = 0.002 K_inj = 0.98 b_inj = 75. aquifers, aquitards, tops = Strat(section, num_contacts, K_aqf0, K_aqt0, K_inj, b_inj) PlotSection(tops) # visualize aquifer-aquitard section # well properties num_aqf = len(aquifers.K) well = [] if mode == 0: well.append(Well(0., 0., -10000., array([0, 1, 0, 0]), aquifers)) # example else: Qw = 18096. # example screen_list = zeros(num_aqf, int) screen_list[-1] = 1 well.append(Well(0., 0., Qw, screen_list, aquifers)) # set up Hemker matrix A = HemkerMatrix(aquifers, aquitards) print 'Set up matrix for Hemker solution.' # solve eigenvalue problem w, V, Z = Eigen(A) aquifers.Leakance(w) print 'Solved eigenvalue problem.' # model drawdown aquifer_layer = arange(0, num_aqf, 1) if mode == 0: profiles_df, x = ProcessLayerProfiles(num_aqf, aquifer_layer, aquifers, well, V, Z) # create output as a dataframe PlotLayerProfiles(profiles_df, aquifer_layer, x) # plot profiles over x else: # plot head, per layer, as a function of depth x = 5. # example y = 0. ProcessVertProfile(x, y, num_aqf, aquifer_layer, aquifers, well, V, Z) print 'Done.' ### run script ### # mode==0 --> run problem with specified layer properties; # generate drawdown curves for every aquifer # mode==1 --> generate many aquifer-aquitard zones stochastically; # generate section diagram and head change profile with depth for given r mode = 1 Hemker(mode)

Feel free to use this code for your own purposes, but note that there is no warranty implied whatsoever! 🙂

### Findings

A visual comparison between the output generated for a particular example problem described in Hemker’s (1984) paper and this script, subject to the same parameter values (e.g., transmissivity) and pumping scenario indicates a correct match:

For a more complex aquifer-aquitard sequence (e.g., 25 aquifers) distributed across the 2,000 ft thickness of the application problem, the induced head perturbation declines rapidly across the aquitard layers above the injection zone. Details are omitted here, but this particular scenario implies thousands, if not tens of thousands, of years of travel time would be required for wastewater to migrate through the sequences of aquitards (again, assuming lateral continuity of all the lithologic units).