Application of Numerical Tools to Investigate a Leaky Aquitard beneath Urban Well Fields

Abstract

Memphis aquifer is the primary drinking water source in Shelby County (Tennessee, USA), and it supplies industrial, commercial, and residential water. Memphis aquifer is separated from the Shallow aquifer by a clayey layer known as the Upper Claiborne confining unit (UCCU). All of the production wells in the Memphis area are screened in the Memphis aquifer, or even deeper in the Fort Pillow aquifer. Traditionally, it was assumed that the UCCU could fully protect the Memphis aquifer from the contaminated Shallow aquifer groundwater. However, recent studies show that at some locations, the UCCU is thin or absent, which possibly leads to the contribution of Shallow aquifer to the Memphis aquifer. Accurately locating the breaches demands expensive and difficult geological or geochemical investigations, especially within an urban area. Hence, a pre-field investigation to identify the locations where the presence of breaches is likely can significantly reduce the cost of field investigations and improve their results. In this study, to identify the locations where the presence of breaches in the UCCU is likely, we develop a reliable MODFLOW-based numerical model, and use three different analyses: (1) pilot-point calibration (PPC), (2) velocity and flow budget (VFB), and (3) particle tracking (PT), to post-process the developed groundwater model results. These pre-field numerical investigations provide relevant and defensible explanations for groundwater flow anomalies in an aquifer system for informed decision-making and future field investigations. In this study, we identify five specific zones within the broad study area which are reasonable candidates for the future field investigations. Finally, we test the results of each analysis against other evidence for breaches, to demonstrate that the results of the numerical analyses are reliable and supported by previous studies.

1. Introduction

The City of Memphis (Shelby County, Tennessee, USA) has undergone significant urban development from 1950 to 2000, partially due to the availability of inexpensive, high-quality water from Memphis aquifer, a regional groundwater resource. Local utilities and industries in Shelby County, Tennessee extract more than 680 thousand cubic meters of groundwater per day that requires minimal treatment for municipal use [1]. The groundwater quality of Memphis is not only crucial for general consumption, but also for the economy of the region [2]. After decades of continuous extraction, water quality degradation and contamination of the Memphis aquifer has become a challenging social, economic, and environmental concern. Hence, local utilities and citizen groups need to take more proactive approaches to address such serious water resources issues, rather than the traditional reactive approaches.

Memphis is within the Mississippi embayment aquifer system, and underlain by approximately one kilometer of unconsolidated sediments (Figure 1). The subsurface is divided into a series of alternating sand and clay layers, grouped into hydrostratigraphic units that include several water-bearing units. As shown in Figure 1, four sand aquifers beneath Shelby County, Tennessee include the Shallow aquifer, the Memphis aquifer, the Fort Pillow aquifer, and the Cretaceous aquifer, each of which are separated from adjacent aquifers by aquitard units, known as Upper Claiborne confining units (UCCU), Flour Island, Tertiary–Cretaceous confining unit, and other Cretaceous units [3,4]. The Shallow aquifer consists of the Mississippi River Valley alluvial (MRVA) aquifer, and the Gulf Coastal Plain (GCP) fluvial and alluvial aquifer. The MRVA aquifer serves as a water source for industrial, domestic, and irrigation wells in the Mississippi alluvial valley, and it is typically not utilized in urban settings. Water from the Shallow aquifer is of poorer quality than that in the underlying Memphis and Fort Pillow aquifers, due to the absence of an overlying aquitard; thus, it is locally contaminated in places. The Memphis aquifer serves as the primary source of drinking water in Shelby County, Tennessee, with minor withdrawals coming from the deeper Fort Pillow aquifer. However, the production well screened in the Fort Pillow aquifer is located in other cities (Cordova, and Millington) far from the City of Memphis. The Cretaceous aquifer is not utilized in Shelby County, due to the depth and poorer water quality.

The UCCU layer separates the Shallow and the Memphis aquifers throughout most of the Shelby County, acting as an upper protective layer to the Memphis aquifer. However, hydrogeologic breaches, also termed windows, within the UCCU locally create hydraulic connections between these aquifers [5]. The quantity of water passing through these hydraulic connections increases with greater groundwater withdrawal and groundwater head depression within the Memphis aquifer [6]. Hence, through these leaky hydrogeologic breaches where the hydraulic conductance is relatively high compared to the rest of the UCCU, the Shallow aquifer recharges the Memphis aquifer. Figure 2 demonstrates a schematic of a breach in the aquitard layer that acts as a preferential flow path from the Shallow aquifer to a screened production well in the Memphis aquifer.

Breaches in the UCCU aquitard threaten the water quality of the City of Memphis; therefore, identifying locations of the breaches are essential for source-water protection, and the development and implementation of wellhead protection plans. However, the seemingly straightforward question is—Where are the breaches? So far, two different approaches have been employed to answer this question. The first approach is the “geochemical investigation” [7,8,9,10,11,12,13] where groundwater samples are analyzed to identify geochemical signatures that are attributed to surface water within the groundwater. Within the Memphis area, age-dating of groundwater by using tritium-helium-3 and sulfur hexafluoride support geochemical data, and reflects the influence of modern water (<50 years) in the Memphis aquifer. Subsequent studies have shown that modern water contributes up to 40% of groundwater withdrawn from production wells—in one case, production water contains as much as 75% modern water [9,10,11]. The second approach is the “geological investigation”, in which techniques, such as surveying stratigraphic thickness and water table elevation, electrical resistivity surveys, shear-wave seismic reflection, and other geophysical analyses are used to identify gaps in impermeable layers [5,14,15,16,17]. Geological investigations show how hydrogeologic units are connected in the subsurface, and where preferential leakage pathways may exist. Parks [5] created an isopach map of the UCCU throughout Shelby County, and identified several areas where clay layers within the UCCU are thin or absent. Parks [5] also prepared a water-table elevation map for Shelby County, and showed that several areas with anomalous water-table depressions also had thin or absent confining unit clay, consistent with vertical leakage through a breach in the UCCU. Waldron et al. [16] used an S-wave seismic reflection survey to map the size and the orientation of localized breaches in the UCCU near an unlined landfill.

The abovementioned approaches have several limitations. First, mapping breaches are complicated, costly, and time-consuming owing to their depth, unknown location, and randomness of origin causality. Second, mapping studies with geophysical methods are technically limited to a small study area, and their measuring accuracy dissipates as the study area extends. Third, approaches involving extensive drilling and geophysical studies are inapplicable or prohibitive in urban areas, due to infrastructure, and compact residential and commercial areas. Therefore, when the presence of a breach in the aquitard layer is likely, numerical investigations conducted prior to field studies (pre-field investigation) can be extremely helpful to identify locations where a breach is likely to exist. Pre-field numerical investigations provide relevant and defensible explanations for groundwater flow anomalies in an aquifer system for informed decision-making.

Several numerical tools are developed to simulate the groundwater flow within the system. Groundwater models provide invaluable information about what cannot be adequately measured or observed, in both space and time domains (e.g., geological conceptualization, flow rates and flow pathways). Such numerical tools apply mathematical equations to calculate aquifer system responses to forcing conditions such as pumping, stream-stage variation, and rainfall events [18]. Besides groundwater models, several other numerical tools and techniques are developed to post-process the flow equations and the results of numerical models, to improve our understanding of the studied groundwater system [19,20,21,22,23,24,25,26,27].

In this work, we show how a combination of forward and backward numerical modeling, and post-processing numerical tools can be used to produce a more satisfactory representation of the groundwater system to make informed field investigation decisions. Herein, three numerical analyses are utilized: (1) pilot-point calibration analysis (PPC); (2) velocity and flow budget analyses (VFB); and (3) particle tracking analysis (PT), to reasonably identify the locations where the presence of breaches is likely. Finally, we test the results of each analysis against other evidence for breaches, to demonstrate that the results of the numerical analyses are reliable and supported by previous studies.

2. Methodology

This study incorporates multiple Memphis Light, Gas and Water (MLGW) well fields, including the Allen, Davis, and Palmer well fields. These three well fields are located east of the Mississippi River bluff line in southwestern Shelby County (Figure 3), including 49 production wells (26 wells in Allen; 19 wells in Davis, and four wells in Palmer) and they collectively produce more than 50 million cubic meter per year over 2005 to 2016 [28]. Herein, we consider only the upper three geological units: (1) Shallow aquifer, (2) UCCU, and (3) Memphis aquifer (Figure 1). This approach is justified because all the production wells within the studying well fields are screened in the upper and middle Memphis aquifer, and little evidence exists for the upward flow of water from deeper hydrogeologic units [3]. Therefore, we neglect the interaction between the lower Memphis aquifer and the Flour Island aquitard, by assuming a no-flow boundary condition at the base of the Memphis aquifer (Figure 1).

2.1. Extent of the Study Area

In this study, we first assemble and analyze the existing relevant hydrogeological field data to define a proper extent for the study area, and we build an advisory conceptual model. In groundwater studies, the delineation of the extent of the study area is quite critical and challenging. To evaluate a proper extent for this study, we use the available potentiometric data of the Memphis aquifer in 2005, 2007, 2010 and 2015 [29,30]. Figure 4 shows the types and locations of boundary conditions that encompass Allen, Davis, and Palmer well fields. The purple, blue, green, and orange lines are associated with the years of 2005, 2007, 2010, and 2015, respectively. Dashed lines symbolize temporally varying boundary conditions in the Memphis aquifer (i.e., the groundwater heads along the dashed lines are spatially constant while their magnitudes vary with time). The solid line symbolizes the no-flow boundary condition at the groundwater divides between Allen, Sheahan, and Mallory well fields (Figure 3 and Figure 4). As shown in Figure 4, locations of both no-flow and temporarily varying boundary conditions remain nearly consistent from 2005 to 2015. Hence, we estimate the study area to encompass the average extent of the boundary conditions in different years, which is shown in black in Figure 4. The same extent is delineated for the Shallow aquifer and UCCU; however, with relevant spatially and temporally varying boundary conditions.

2.2. Conceptual Model

Figure 5 shows the conceptual model for the site-specific hydrogeological setting. In the conceptual model, we simplify the system conditions and capture the essential features of the site. Figure 5a shows a plan view of the conceptual model. A bluff line of 15 to 45 m high separates the alluvial and fluvial Shallow aquifers. Consequently, the hydraulic characteristics of the Shallow aquifer vary significantly on either side of the bluff line. Figure 5b shows a side view of the conceptual model which is made up of three layers: Shallow aquifer, UCCU, and the Memphis aquifer. Stratigraphic unit elevations were derived from [4]. Production well screen intervals are illustrated in Figure 5b, taking note that they vary in depth. This aquifer system conceptualization was very helpful in developing the numerical model.

2.3. Groundwater model

To describe the transient groundwater flow in the study area, we develop a high-resolution three-dimensional numerical model using [31]:

$$S\frac{\partial h}{\partial t}=\frac{\partial }{\partial x}\left(K_x\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_y\frac{\partial h}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_z\frac{\partial h}{\partial z}\right)+w,$$

(1)

where h (L) is groundwater head, Kx, Ky, and Kz (L/T) are hydraulic conductivity in x, y, and z directions, respectively; t is time (T); S (1/L) is the storage coefficient of the aquifer and w (1/T) is a source and sink parameter (volumetric flux per unit volume). We used MODFLOW-NWT to solve Equation (1) [32]. The numerical model has six layers with a uniform $250\times 250$ m² cells, totaling 20,773 active cells. The first and second numerical layers represent the Shallow aquifer and UCCU, respectively. The other four numerical layers uniformly divide the Memphis aquifer. The time domain of the model begins from 2005 to 2015, a span for which we have relatively sufficient data for Shallow and Memphis aquifer groundwater head within the study area. This duration is discretized into 121 monthly stress periods, each subdivided into five smaller time steps.

2.4. Initial Conditions and Boundary Conditions

In this study, the initial head conditions of the Shallow aquifer were set to the water-table surface (Shallow aquifer) as measured in 2005. There was no information available regarding the groundwater head in UCCU. Therefore, we assumed the same initial conditions within the UCCU in 2005 [33]. The initial head conditions for the Memphis aquifer were derived from the 2005 Memphis aquifer potentiometric surface [30]. Boundary conditions of the Shallow aquifer, UCCU, and Memphis aquifer were estimated using available potentiometric surfaces from 2005 to 2015. The boundary condition of the Shallow aquifer was assumed to be spatially and temporally varying during 2005–2015, except beneath the Mississippi River, where it is controlled by the river stage variations. Therefore, the boundary condition under the Mississippi River is estimated according to the monthly measurements of river stage at the nearby Beale Street gage during the period of 2005 to 2015 [34]. We neglected the horizontal inflow and outflow at the UCCU boundary condition; therefore, the boundary condition of the UCCU was assumed to be a no-flow boundary condition. Two boundary conditions were applied to the Memphis aquifer: no-flow (which is perpendicular to equipotential lines) and temporarily varying boundary conditions (Figure 3). The Memphis aquifer head variation at the varying boundary condition, at many points, remained within a relatively small range of fluctuation during 2005 to 2015. Therefore, the temporally varying boundary condition of the Memphis aquifer was estimated by interpolating available 2005 to 2015 Memphis aquifer heads [30].

2.5. Internal Forcing Conditions

Internal forcing conditions within the study area included McKellar Lake, Nonconnah Creek, production wells, and surface recharge. Since McKellar Lake has an open connection to the Mississippi River, the lake stage follows Mississippi River stages. The Nonconnah Creek stage variation was modeled using monthly mean stages measured at upstream and downstream stations, Farris View Drive stage and Riverport stage, during the period 2011 to 2017, when the stage data were recorded (Rivergages.com, 2018). The exact withdrawal rates of the production wells are unavailable. Therefore, monthly production for the entire well fields are uniformly distributed to each active well during that month [34]. The pumping rates of each well field and the number of active wells in each month are shown in Figure 6a–c. Precipitation recharge was applied to the upper layer, representing the Shallow aquifer surface recharge. Unfortunately, recharge rates to the shallow aquifer are unknown, and they must be inferred from available precipitation records. However, since precipitation events occur almost in all months of each year in the study area [35] we assumed a temporarily constant recharge rate over the modelling time domain and a spatially varying recharge rate over the study area (See Appendix C). Finally, the recharge rates will be estimated during model calibration processes.

3. Numerical Analysis Results

In this study, we located zones where the presence of leaky aquitard (UCCU) breaches was likely, by applying three numerical analyses: (1) PPC performed against horizontal and vertical hydraulic conductivities, (2) VFB, and (3) PT. As it pertains the identifying plausible aquitard breach locations, the following observances were assumed: (1) hydraulic conductivities at breaches were expected to be locally high, (2) convergence of flow to a specific area was a sign of the presence of a breach, and (3) vertical flow at breaches was expectedly high, compared with other parts of the study area. These conditions will be derived from the PPC and VFB. In contrast, PT should show that groundwater particles at the breaches vertically transport further than in the surrounding area.

3.1. Pilot Point Calibration Analyses (PPC)

Although the boundary conditions and withdrawal rates were estimated, several parameters required for groundwater modeling, including horizontal and vertical hydraulic conductivities, storage coefficients, surface recharge, and the conductance of the hyporheic zone at Nonconnah Creek bed were remain unknown. To estimate the unknowns, we used an advanced nonlinear method of spatial parameter characterization to develop an inverse groundwater model [36]. An inverse groundwater model includes a calibration process to achieve a satisfactory fit between the model predictions and field measurements. However, most of the time, the solution to the groundwater problem is inherently non-unique [36,37,38]. Based on the groundwater literature, one would address the non-uniqueness issue by introducing additional “simplifying” information such as the use of zonal calibration, in which the model is subdivided into a limited number of zones with constant parameter values [39]. The values of parameters in each zone can be restricted to a range that is estimated from the actual field data. Furthermore, additional information can be introduced to the model through advanced mathematical regularization in the parameter calibration process, like truncated singular value decomposition (SVD) [36,37,38,40]. Finally, like all predicting models, the model predictions are required to fit the collected observation data, and to produce reasonable results, based on previous knowledge of the system. In this study, we used all of the abovementioned techniques to address the challenging non-uniqueness issue.

We first used uniform properties for each layer to estimate proper starting calibration values within a meaningful range of parameters, as reported in previous studies conducted in the Mississippi embayment area [3,39,41,42,43,44,45,46]. In the groundwater literature, this approach is known as zonal calibration [47,48]. We considered two zones at the east and west of the bluff for Shallow aquifer. Similarly, we considered two zones at either side of the bluff for UCCU. Finally, we considered one zone for each numerical layer (four layers) of the Memphis aquifer. The range of parameters estimated in previous studies within and outside of the study area are summarized in

Table 1. Minimum and maximum of hydraulic parameters reported in previous studies within and nearby the study area [3,39,41,42,43,44,45,46].

	Minimum				Maximum
	Horizontal hydraulic conductivity (m/d)	Vertical hydraulic conductivity (m/d)	Storage coefficients (1/m)	Specific yield (-)	Horizontal hydraulic conductivity (m/d)	Vertical hydraulic conductivity (m/d)	Storage coefficients (1/m)	Specific yield (-)
Shallow aquifer - East of the bluff	1	2.60 × 10−5	-	0.22	46	7.00 × 10−1	-	0.48
Shallow aquifer - West of the bluff	36	-	-	0.28 *	119	-	-	0.28 *
UCCU	6.10 × 10−6	1.50 × 10−6	9.84 × 10−4 *	-	9.14 × 10−1	6.10 × 10−3	9.84 × 10−4 *	-
Memphis aquifer	8	-	3.28 × 10−4	-	48	-	9.84 × 10−3	-

* Only the estimation is available.

. However, it should be noted that inherent uncertainty exists in the estimated data, due to the limited number of estimations, and due to the different test protocols [49,50]. This type of uncertainty increases with depth, especially in the Memphis aquifer, which needs to be considered during modeling and calibration.

Parameters used for this pre-calibration process were horizontal and vertical hydraulic conductivities, storage coefficients, the Nonconnah Creek bed conductance, and surface recharge. To describe the heterogeneity in hydraulic conductivity and surface recharge within each zone, we used “pilot points” accompanied by advanced singular value decomposition (SVD) regularization functionality available through PEST [36,37]. This model included 69 pilot points distributed in all layers. The inverse distance weighted method is used to spatially distribute pilot points values to each cell of the model. We used PEST to estimate a new set of the model parameters (i.e., horizontal and vertical hydraulic conductivities, storage coefficients, the Nonconnah Creek bed conductance, and surface recharge) using a robust mathematical technique by comparing simulation results to a set of observations [36,37]. Observation points were derived from the available Shallow and Memphis aquifers potentiometric surfaces maps [30,33]. To evaluate the calibration results, we compared the simulated transient groundwater head with the observed data at 23 observation points within the Shallow and Memphis aquifers. Shallow aquifer observation data included water levels collected in 2005 and 2015. The Memphis aquifer observation data included groundwater heads collected during 2005 to 2015. Groundwater head fluctuation in each observation point was not high; therefore, the monthly varying head at each observation point was estimated by linear interpolation between the years for when the observation data were available. To evaluate the goodness of calibration, the observed and simulated groundwater heads were compared with each other. The observed versus simulated groundwater heads are shown in Figure 7. The goodness-of-fit coefficients between the simulated and calibrated results of the model were computed for all observation points within the Shallow and the Memphis aquifers [51]. $R^2$ values of Shallow and Memphis aquifers were estimated to be 0.98 and 0.84, respectively. Root Mean Square Error (RMSE) values of Shallow and Memphis aquifers were estimated to be 1.6 m and 1.0 m, respectively. The estimated $R^2$ and RMSE indicated a satisfactory goodness of fit.

The calibrated horizontal and vertical hydraulic conductivities of the UCCU are shown in Figure 8a,b, respectively. Figure 8 shows that the patterns of the distribution of horizontal and vertical hydraulic conductivities were quite similar, however, with different values (see Appendixe A, Appendixe B and Appendixe C which summarize all six numerical layers’ calibrated hydraulic conductivities, storage coefficients, and surface recharge, respectively). As shown in Figure 8, PPC revealed three major zones (dashed lines) in which the vertical and horizontal hydraulic conductivities of the UCCU were relatively high. Since we expected locally high hydraulic conductivities at the breaches, we concluded that the presence of the breaches at these three zones were more probable than other parts of the study area.

3.2. Validation of PPC Results

The increase of hydraulic conductivities represented within the three identified zones (Figure 8) likely exist due to lithologic variation or erosion of the UCCU, with an increase of coarser grained sediments [52,53]. The first identified area by PPC is located in the Memphis Depot Industrial Park (MDIP) (Figure 8). Recent investigations have already indicated that the Shallow groundwater level in this area is locally lower than that in the surrounding areas [33,54]. Figure 9 shows the approximate flow directions of Shallow aquifer in 2005 and 2015, which converge from the surrounding areas toward the MDIP. There are no pumping wells in this area, which otherwise would have resulted in a local cone of depression.

The second zone is located on President’s Island (Figure 8). Previous studies and field investigation in this area have indicated that the thickness of the UCCU beneath President’s Island significantly thins [41]. Calibration results in the President’s Island area affirm the suspicion of the thinning or absence of the UCCU. Moreover, in the area immediately south of President’s Island, which includes the Tennessee Valley Authority (TVA) combined cycle plant, a hydraulic connection of the Shallow and Memphis aquifer is identified [54] (Figure 8). This recent study supports the likely presence of a breach in that area.

The third identified zone is located at the southern part of the study area adjacent to the bluff line. This area is an undiscovered zone, where field investigations are lacking. Therefore, this area is considered as an area of interest, where further investigations are required to validate the PPC results.

3.3. Velocity and Flow Budget Analyses (VFB)

High hydraulic conductivity alone at a breach does not necessarily lead to a vertical flow, and it needs sufficient vertical gradient to create flow. Determining the velocity and flow directions of the Shallow aquifer can provide valuable information regarding the likely presence of a breach. Figure 10a–c shows the horizontal velocity vector of the Shallow aquifer in the study area in 2005, 2010, and 2015 derived from the developed numerical model. The magnitude of the velocity is proportional to the arrows lengths. A longer arrow refers to a higher horizontal velocity. As shown in Figure 10a–c, the VFB approach identified three zones where the velocity vectors converged to specific areas. The horizontal velocity of groundwater decreases inside the convergence zones, suggesting downward flow (assuming upward flow is negligible due to the regional downward vertical gradient). As shown in Figure 11a–c, VFB analyses also show vertical downward flow budgets at the bottom of the Shallow aquifer in 2005, 2010, and 2015, respectively. Negative values refer to an upward flow, and positive values refer to a downward flow. As shown in Figure 11a–c, the vertical flows in all three identified zones by VFB are locally high. Figure 11a–c shows a weak correlation between the extent and location of the identified zones, and the extent and location of the production wells (note that all wells are screened within the Memphis aquifer). It might be because the UCCU weakens the influence of wells on the vertical flow of the Shallow aquifer. Therefore, it can be concluded that the local increase in the vertical flow of the shallow aquifer is not associated with the production wells. The main difference between the results of the PPC and VFB analyses is that VFB did not show any local flow convergence and an increase in vertical flow budget at the TVA area. Using this approach, we can conclude that the presence of breaches at the three identified zones is more likely than other parts of the study area.

3.4. Validation of VFB Results

As shown in Figure 10, the velocity vectors converge and focalize to three specific zones where the vertical flow of the shallow aquifer is locally high, as shown in Figure 11. These zones are almost the same as the three zones identified by PPC analysis, except in the TVA area; therefore, VFB analysis only verifies the previous findings in the MDIP and President’s Island. However, it should be noted that VFB identifies smaller breach areas, compared with PPC. The third identified zone at the southern part of the study area is an undiscovered zone, where further local and detailed investigations are required for verification. It should be noted that the convergence of the flow near the Nonconnah Creek is not necessarily, because of the presence of a breach in this area. Flow budget analysis shows that the convergence of the flow in this area is mainly because of the drainage of the Shallow aquifer to the creek.

3.5. Particle Tracking Analyses (PT)

In this approach, we used the numerical tool, MODPATH [20], to conduct forward particle tracking. Figure 12a shows the cross-section A-A’ of the Shallow aquifer and UCCU, an east-to-west transect of the study area that is shown in context in Figure 12b. Forward particle tracking represents the location of imaginary particles at the Shallow aquifer surface, and at the top of the UCCU from 2005 (green marks) to 2015 (red marks). Figure 12a shows the vertical downward flow transport from the Shallow aquifer and UCCU to the Memphis aquifer. A sudden change was observed at the bluff line, due to differences in alluvial (west of the bluff) and fluvial (east of the bluff) deposits characteristics. The particles in the MDIP area also locally penetrated further from 2005 to 2015. As shown, particles west of the bluff within the President’s Island moved downward significantly faster than the east of the bluff, consistent with the higher hydraulic conductivity of the MRVA and fluvial deposits at either side of the bluff line.

Figure 13a shows placement of north–south cross-section B–B’ of the Shallow aquifer and UCCU, which are shown in context in Figure 13b. Cross-section B–B’ intersects the bluff at three points. Figure 13a shows a sudden change in particle migration proximal to the bluff, due to hydraulic characteristics differences in alluvial and fluvial deposits at either side of the bluff line. Particles at the southern part of the cross-section B–B’ transport further at depth from 2005 to 2015, compared with particles at the northern part of the cross-section.

Figure 14a shows the cross-section C–C’ of the Shallow aquifer and UCCU, an east-to-west transect of the study area that is shown in context in Figure 14b. The cross-section intersects the bluff at one point. Similarly, we place imaginary particles on the Shallow aquifer and top of UCCU in the year of 2005 (green marks). Forward particle tracking represents the location of each particle in 2015 (red marks). Same as in cross-section A–A’, an abrupt change is observed at the bluff; however, PPT results show a relatively high interaction between the layers adjacent the bluff line. As shown in Figure 14a, PPT shows no sign of a breach at the TVA area. However, there is a broad area east of the bluff where particles transport paths are relatively long. Therefore, there is a high possibility for the presence of a breach in this area.

3.6. Validation of PT Results

PT results confirmed all the zones identified by PPC and VFB analyses. The results show that the particles in the President’s Island and MDIP transport locally further, which is consistent with previous studies [14,34,54] in these areas. PT shows that particles in the TVA area transport slowly downward, except adjacent to the bluff line. This finding is consistent with a recent study [54]. However, similar to PPC analysis, there is an undiscovered area in the southern part of the study area where PT identifies it as the locations where the presence of breaches is likely.

4. Summary and Conclusions

The Memphis aquifer is the primary drinking water source in the Shelby County (Tennessee, USA), and it supplies industrial, commercial, and residential water. Memphis aquifer is one of several deep-sand aquifers separated from the Shallow aquifer by a clayey confining layer, UCCU. All of the production wells in the Memphis area are screened in the Memphis aquifer or even deeper in Fort Pillow aquifer, to avoid Shallow aquifer contamination. Traditionally, the public water supply authority in the Memphis area assumed that the UCCU could fully protect the water in the Memphis aquifer from downward leakage of the contaminated or poor-quality groundwater from the overlying Shallow aquifer. However, recent studies show that at some locations, the UCCU is thin or absent, which leads to the contribution of water from the Shallow aquifer to the Memphis aquifer [9,10,11,14]. Hence, sustainable water management in the Memphis area requires the locations of breaches in the UCCU to ensure the protection of water pumped from the production wells. Locating breaches demands expensive geological or geochemical investigations. These studies are even more costly and complicated in densely developed urban areas. Therefore, other techniques are desired that can help to determine the locations of breaches, so that complicated and costly field investigations can be focused in limited areas.

In this study, to identify the locations where the presence of breaches in the UCCU is likely, we use three different analyses: (1) pilot point calibration (PPC), (2) velocity and flow budget (VFB), and (3) particle tracking (PT). As shown in Figure 15, these approaches identify five zones where the presence of breaches is expected. As summarized in

Table 2. Summary of the PPC, VFB, and PT analyses identified zones.

	Type of Analysis			Is it Consistent with Previous Investigations?
	PPC	VFB	PT
Zone-1	☑	☑	☑	Yes
Zone-2	☑	☑	☑	Yes
Zone-3	☑	☑	☑	N/A (no previous study)
Zone-4	☑	□	☑	Yes
Zone-5	☑	□	☑	N/A (no previous study)

, all three approaches distinguish Zones 1, 2, and 3 from other parts of the study area. However, VFB does not identify Zones 4 and 5 as the locations of possible breaches. Previous investigations results are concurrent with this study, in that the presence of breaches in Zone 1, 2, and 4 is likely. These findings validate the reliability of the numerical model. However, the presence of breaches in Zones 3 and 5 have never been explored, and they need further investigation. Five zones identified in this study are reasonable candidates for future field investigations, based on the current understanding of the study area. The results of this study can reduce the cost and increase the efficiency of investigations all over the broad study area, by directing future studies.

In this study, we show how a numerical pre-field investigation can tackle complex hydrogeological problems. A similar approach and analyses that have been developed in this work can be applied to study the interaction between the geological layers of an aquifer, provided that enough information about the boundary conditions, hydraulic properties, and geologic characteristics are available. As in all numerical models, the consistency, stability, and accuracy of the model depends on the quantity and quality of the available data. Gathering hydraulic and geologic data is costly and time-consuming, and it usually needs several years to fulfill. This study is an example of reliable models that have been built on information collected over several years.