Development of a Numerical Multi-Layered Groundwater Model to Simulate Inter-Aquifer Water Exchange in Shelby County, Tennessee

Abstract

Inter-aquifer water exchange between the shallow and Memphis aquifers in Shelby County, Tennessee may pose a contamination threat due to the downward migration of younger, poor quality groundwater into deeper, more pristine aquifer. Discontinuities (breaches) in the upper Claiborne confining unit (UCCU) allow for leakage into the Memphis aquifer, a sand-dominated aquifer that provides about 95% of the groundwater used in the Memphis area. This study created a multi-layered 3D groundwater model for Shelby County using the United States Geological Survey’s MODFLOW-NWT program to evaluate water exchange for a simulation period from January 2005 to December 2016. Results indicate an overall leakage through the UCCU of 61 m³/min into the Memphis aquifer in Shelby County, accounting for 10% of its water budget inflow, with localized areas experiencing as much as 20% water exchange. As young water tends to stay in the upper part of the Memphis aquifer, water budget assessment for the upper 60 m of the Memphis aquifer revealed leakage representing 29% of the zone inflow, and as much as 53% in certain areas. More localized studies must be conducted to understand the location, characteristics, and orientation of the confining unit breaches, as well as the inter-aquifer water exchange.

1. Introduction

Water exchange in groundwater systems comprised of alternating aquifers and leaky confining units can influence the water characteristics and may pose a contamination threat due to downward leakage of low-quality water, especially in urban areas where groundwater degradation is common [1,2,3], as is the case in Shelby County, Tennessee [4,5,6,7,8]. Due to inter-aquifer water exchange, the Memphis aquifer is experiencing localized water quality degradation and contamination [9,10,11,12,13,14].

Numerical modeling of inter-aquifer water exchange is commonly difficult because of the varying methods for characterizing the heterogeneity leading to the exchange [15,16], the disconnections between surface water and alluvial aquifer water from deeper aquifer systems [17], and changes in inter-aquifer water exchange with pumping [18,19].

Analytical methods that evaluate inter-aquifer water exchange can be applied to simple hydrogeologic systems [20]; however, more complex settings with time variable stressors, such as pumping, require numerical simulations [21].

Local variability in inter-aquifer exchange can be evaluated numerically [22], but the impacts on water budget under transient conditions in a complex hydrogeologic and pumping setting are not well established. Previous numerical modeling studies of the Memphis aquifer in Shelby County [23,24,25,26,27,28] have not adequately addressed inter-aquifer water exchange to the degree needed to be predictive of current hydrogeologic conditions and useful as a sub-regional water management tool.

Shelby County lies in the northern half of the Mississippi embayment in the southwestern corner of Tennessee (Figure 1a–c). The embayment is a trough-shaped basin filled with more than a thousand meters of Cretaceous, Tertiary and Quaternary age sediments forming a series of aquifers and confining units [29,30,31,32,33,34]. The surficial Pliocene, Pleistocene and Holocene fluvial terrace and alluvial deposits in the County form a shallow (unconfined) aquifer. Withdrawal from the shallow aquifer is minimal and primarily used for irrigation and domestic use [12]. These deposits are overlain by Pleistocene loess and underlain by silt, sand, and clay of the Eocene upper Claiborne Group (i.e., upper Claiborne confining unit, UCCU). The UCCU provides partial to complete confinement for the underlying Memphis aquifer, a sand-dominated aquifer that provides about 95% of the groundwater used in the Memphis area [35,36].

The UCCU includes an extensive thickness of interstratified silt and clay throughout the Memphis area [4,6], but local discontinuities within the silt and clay, termed breaches, provide avenues for inter-aquifer water exchange [4,32]. The presence of these breaches, even at a local scale, can have a region-wide costly impact on water quality when younger, potentially contaminating water, and causing leaks into deeper, more pristine Memphis aquifer reservoirs [12] with residence times between 1000–3000 years [32]. Age dating of local production wells have persistently shown presence of modern water, with ages between 20 and 45 years in certain areas, and with tracer-based mixing percentages indicating as much as 30% modern water in the upper Memphis aquifer [37].

Total groundwater withdrawals in Shelby County increased consistently from pre-development conditions (1887), reaching a maximum of 830,000 m³/day in 2005 [38], but show a decreasing trend in recent years with a total pumping of 696,000 m³/day in 2015 [39]. Groundwater level observations for the Memphis aquifer indicate that a broad regional cone of depression developed in the potentiometric surface during the mid-1920s [35] and continues today (Figure 2a). Though the Memphis aquifer initially experienced an upward gradient through the UCCU when flowing artesian conditions existed, continual withdrawals since 1887 have resulted in a gradient reversal, whereby the shallow aquifer now contributes water by downward leakage through the UCCU, exacerbated where breaches exist [4,6,23,24,25,26,37,40]. In addition to stratigraphic interpretations of breaches, water-table maps of the shallow aquifer [4,32,41,42,43,44] have identified anomalous depressions (Figure 2b) that indicate probable locations where the confining unit is also thin or absent and where less restricted downward leakage occurs, being that no high-capacity production wells for public supply pump from the shallow aquifer.

Groundwater modeling of water resources in the Mississippi embayment has been used to assist with water management and to understand flow through the system under pre-development and modern conditions, mainly focusing on the Mississippi River Valley Alluvial (MRVA) aquifer in Arkansas (e.g., [45,46,47]). Most recently, a 13-layer, 1.6-km cell, three-dimensional model known as the Mississippi Embayment Regional Aquifer Study (MERAS) was developed by Clark and Hart [27], receiving updates by Clark et al. [28]. It is a regional-scale tool, mainly focused on agriculture, intended for quantifying groundwater availability within the Mississippi embayment.

Although models have been used extensively, only a handful have specifically evaluated inter-aquifer water exchange and flow in Shelby County. Regionally, Arthur and Taylor [24] created a five-layer, 8-km cell, steady-state model that evaluated 1980 conditions. Simulation of water budgets showed that about 35% of the water pumped in the Memphis area comes from downward leakage through the UCCU, moving at a rate of as much as 10 cm/year. About 0.25 cm/year flow moves upward from the Flour Island confining unit into the Memphis aquifer. Arthur and Taylor [25] created a modified version of the model to evaluate the Mississippi embayment in a piece-wise manner for 1987 conditions, stating that for the northern section of the embayment (i.e., above the 35th parallel, 54,000 km²), net vertical downward flow contributed 180 m³/min to the Memphis aquifer and 130 m³/min to the Fort Pillow aquifer.

Focusing on Shelby County, Brahana [23] created a one-layer model for the Memphis aquifer, with cells ranging from 30-km at the margins of the model to 1-km in the Memphis area. During calibration, zones of high leakage in the UCCU along rivers were essential in simulating observed water levels. The water budget for the period of 1966–1970 showed a vertical leakage into the Memphis of 160 m³/min, or 37% of the total budget. An updated version of the same model, including the Memphis and Fort Pillow aquifers [26], stated that water exchange between the aquifers and confining units is a major component of the hydrologic budget in the Memphis area, and is not distributed uniformly. Leakage from the shallow aquifer contributed 266 m³/min at a rate of 0.35 cm/year to the 1980 water budget of the Memphis aquifer, which represents 54% of the total inflow. Leakage of 3.6 m³/min coming from the Fort Pillow represented 1% of the Memphis aquifer inflow.

Even though previous models have attempted to address inter-aquifer water exchange, each has shortcomings that present challenges yet to be resolved. The models do not simulate the shallow aquifer in Memphis [27,28], or they only add it as constant head nodes for inflow to the system [23,26], or as river cells with a head dependent flux component [24,25]. In previous models of the study area [23,24,25,26,27,28], local breaches were not addressed due to coarse grid sizes or lack of data, and pumping was based on water-use reports and distributed among approximated well locations. Model layers were treated as confined under the assumption that saturated thickness is constant in time and equal to the entire layer. Except for the MERAS model [27,28], the other models [23,24,25,26] are considered quasi-3D models because vertical flow through confining units is not modelled explicitly, instead using a leakance term that divides the vertical hydraulic conductivity by the thickness of the unit; an assumption that neglects storage and could introduce significant error [50].

Because of these assumptions, inconsistencies among the models, and the fact that the last “local” model that addressed water exchange was developed for 1980 conditions [26], the purpose of this study is to develop a detailed, multilayer, three-dimensional numerical groundwater-flow model that more accurately simulates the actual hydrogeologic conditions of water exchange between external stressors to the shallow, Memphis, and Fort Pillow aquifers and their respective intervening confining units in Shelby County.

Hydrogeologic Units

Following prior investigations, the first 500 m of the subsurface in Shelby County is divided into the following six hydrogeologic units in descending order: shallow aquifer (Quaternary units), UCCU, Memphis aquifer, Flour Island confining unit, Fort Pillow aquifer, and the Old Breastworks confining unit (

Table 1. Geologic and hydrostratigraphic units underlying Shelby County, Tennessee, with minimum and maximum values for hydraulic conductivity and model layer designation. (Adapted from Graham and Parks [32], and modified from Larsen et al. [37].

System and Series		Group	Stratigraphic Unit	Hydro- Stratigraphic Unit	Thickness (m)	Lithology	Model Layer	Hydraulic Conductivity (m/day)
								Min	Max
Quaternary	Holocene and Pleistocene		Alluvium	shallow (alluvial) aquifer	0–50	Sand, gravel, silt, and clay. Underlies the Mississippi River alluvial plain and alluvial plains of tributary streams in western Tennessee. Thickest beneath the Mississippi River alluvial plain, where commonly between 30 and 45 m thick; generally, less than 15 m thick elsewhere.	1	2.5 × 10−5 f	45 f
	Pleistocene		Loess	Leaky confining unit	0–20	Silt, silty clay, and minor sand. Principal unit at the surface in upland areas of western Tennessee. Thickest on the bluffs that border the Mississippi alluvial plain; thinner eastward from the bluffs.
Quaternary and Tertiary	Pleistocene and Pliocene (?)		Fluvial terrace deposits	shallow (Fluvial) aquifer	0–30	Sand, gravel, minor clay, and ferruginous sandstone. Generally, underlies the loess in upland areas, but locally absent. Thickness varies greatly because of erosional surfaces at top and base.
Tertiary	Eocene	Claiborne	Jackson Formation	upper Claiborne confining unit	0–110	Clay, silt, sand, and lignite. Because of similarities in lithology, the Jackson Formation and upper part of the Claiborne Group cannot be reliably subdivided based on available information. Most of the preserved sequence is the Cockfield and Cook Mountain Formations undivided.	2	1.5 × 10−6 f	6.0 × 10−3 f
			Cockfield and Cook Mountain formations
			Memphis Sand	Memphis aquifer	150–270	Sand, clay, and minor lignite. Thick body of sand with lenses of clay at various stratigraphic horizons and minor lignite. Thickest in the southwestern part of the Memphis area; thinnest in the northeastern part.	3–6	8.5 d	47 g
	Eocene(?)	Wilcox	Flour Island Formation	Flour Island confining unit	50–95	Clay, silt, sand, and lignite. Consists primarily of silty clays and sandy silts with lenses and interbedded fine sand and lignite.	7	3.5 × 10−3 b
	Paleocene		Fort Pillow Sand	Fort Pillow aquifer	40–95	Sand with minor clay and lignite. Sand is fine to medium grained. Thickest in the southwestern part of the Memphis area; thinnest in the northern and northeastern parts.	8	4 e	28 a,c
		Midway	Old Breastworks Formation	Old Breastworks confining unit	55–110	Clay, silt, sand, and lignite. Consists primarily of silty clays and clayey silts with lenses and interbedded fine sand and lignite.	No-Flow Boundary

^a Schneider and Cushing [54], ^b Graham and Parks [32], ^c Parks and Carmichael [55], ^d Parks and Carmichael [56], ^e Carmichael et al. [33], ^f Robinson et al. [5], and ^g Gentry et al. [57].

and Figure 3).

West of the Mississippi River bluff, the water table corresponds to the MRVA aquifer [49]. Quaternary units east of the bluff are included in the shallow aquifer. The water table is lower than the land surface (except at springs and seeps), but it generally conforms to the topography [4]. Commonly, and during most of the year, the shallow aquifer discharges as baseflow to the three main tributaries of the Mississippi River: Loosahatchie River, Wolf River, and Nonconnah Creek [4,9,23,32], although there are locations where the rivers lose water into the aquifer [9,10,23,51,52,53].

Even though only limited and sporadic water level records are available for the shallow aquifer, observed seasonal fluctuation is less than 3 m [26]. According to Ogletree [44], the average overall change in head between 2005 and 2015 was only 15 cm, even though there are areas that change and behave independently. Due to the complex heterogeneity of the shallow aquifer, its hydraulic properties span over several orders of magnitude (see Table 1).

The UCCU is composed of the Eocene Cockfield and Cook Mountain formations. It underlies the shallow aquifer and overlays much of the Memphis aquifer with clay, silt, and sand, ranging in thickness from 0 to 60 m but thins toward the eastern part of Shelby County [4] (See Figure 3). Several investigations have shown the presence of hydrogeologic breaches where the shallow aquifer is in direct hydraulic connection with the Memphis aquifer owing to thin or absent clay intervals [4,6,7,9,12,13,22,32,33,40,59,60,61].

The Memphis aquifer and its hydrogeologic equivalents in adjacent states are a predominantly fine- to coarse-grained, sand-dominated regional aquifer [56] that serves as the main water source in the Memphis area and western Tennessee. It ranges in thickness from 150 to 270 m, tapering to zero thickness along the margins of the Mississippi embayment. It has an average hydraulic conductivity of 15 m/day [12] with common storage coefficients ranging from 0.003 to 0.01 m⁻¹ [56]. Even though it becomes unconfined in east-southeast Shelby County, it is generally overlain by a thin (<5 m) veneer of fluvial gravel and Pleistocene loess [62]. The potentiometric surface in the Memphis aquifer is strongly influenced by pumping, causing a complex flow pattern (see Figure 2a) with documented head drawdowns of more than 45 m [25,26,35,48].

The Flour Island Formation is a relatively thick and widespread confining unit in most of western Tennessee, forming the confining unit between the Memphis and Fort Pillow aquifers [56]. The Fort Pillow aquifer is composed of fine- to medium-grained sand with minor clay lenses and is commonly about 60 m thick beneath Shelby County [55]. It is used to supplement extraction from the Memphis aquifer in Shelby County and serves as the major source of water for the City of West Memphis, Arkansas [26,35]. The Old Breastworks confining unit underlies the Fort Pillow aquifer (See Figure 3) as a thick sequence of silty clay beds and interbedded fine sand that hydraulically separate the Tertiary from the underlying Cretaceous aquifers [49].

2. Materials and Methods

2.1. Conceptual Model

To simulate the complex groundwater flow beneath Shelby County, a conceptual model was assembled based on geological, geophysical, hydrological, hydrogeochemical, and other ancillary information (e.g., [50]). Boundary conditions conform with observed historic water levels around Shelby County. Additionally, extensive effort was taken to accurately represent the dynamic pumping regime employed by the five utilities in Shelby County and those in Mississippi, as well as surface water–groundwater exchange stressors, and selected water-demanding industrial and recreational users.

The United States Geological Survey (USGS) MODFLOW-NWT program [63] was chosen to simulate groundwater conditions beneath Shelby County by solving the partial differential equation that describes the three-dimensional movement of groundwater through porous material as defined by Harbaugh [64].

$$S_s\frac{\partial h}{\partial t}=\frac{\partial }{\partial x}\left(K_{xx}\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yy}\frac{\partial h}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_{zz}\frac{\partial h}{\partial z}\right)+W$$

(1)

where Kxx, Kyy and Kzz correspond to hydraulic conductivity values (L/T) along the x, y and z coordinate axes, h is the potentiometric head (L), Ss is the specific storage of the material (L⁻¹), t is time, and W is a volumetric flux per unit volume that represents sources and/or sinks (T⁻¹).

This Newton formulation variant of MODFLOW-2005 [64] has an improved algorithm for handling cell desaturation, a common problem that would create instability and cell deactivation in models with water-table aquifers where saturated thicknesses thin or, in our case, groundwater drains through confining unit breaches. The transformation of the conceptual model into a MODFLOW-NWT model [63] was performed using the pre- and post-processing software GMS (Groundwater Modeling System, v.10.3), developed by Aquaveo [65].

Hydrogeologic layers (see Table 1 and Figure 3) were created from digital surfaces developed for the MERAS model using a cell spacing of 1.6 km [58]. The top of the model is overlain by a resampled 500-m digital elevation model (DEM), extracted from a composite 1/3 arc-second (10-m) DEM obtained from the USGS 3D Elevation Program (3DEP) [66,67].

The model area, as shown in Figure 1c and Figure S3 in the Supplementary Material, was chosen to capture all of the well fields in Shelby County as well as pumping centers in nearby neighboring municipalities that may influence flow conditions. The western side of the model was truncated to maintain no-flow conditions in the northwestern corner of the Memphis aquifer as it is perpendicular to the general flow, and to match a groundwater divide [48,68] with a constant head in the southwestern corner, both considered invariant hydraulic boundaries.

The model discretizes the study area into a finite-difference grid (Figure 4) with 388 rows by 391 columns and eight layers, rotated to conform to the principal west-northwest groundwater flow direction in the Memphis aquifer [69]. The cells have uniform horizontal grid dimensions of 250 m, which approximates the average distance between Memphis well field production wells; therefore, no more than one well falls inside any single cell in most cases, avoiding over-lumping pumping.

The model is transient to simulate temporal variations in pumping and the river stage. The simulation period extends from January 2005 to December of 2016. This timeframe was chosen to match the period with the best availability of consistent water level information in all three aquifers. Stressors to the system were discretized into monthly periods, resulting in a total of 144 stress periods. The model structure has just under 600,000 active cells. The Quaternary units were grouped into a single layer due to their complex lithology, varying thicknesses, and poorly defined transitions. The Memphis aquifer was split into four layers to better represent the vertical distribution of heads and variations of hydraulic conductivity, as well as to better approximate the actual screen positions of the production wells.

The average total depth for public supply wells in the model, including screen length, is 140 m in the Memphis aquifer, and 420 m in the Fort Pillow aquifer (

Table 2. Average depth and screen lengths for wells, listed by well field. The well fields of Millington, Shaw, and DeSoto County municipalities have wells screened in both the Memphis and Fort Pillow aquifers.

Location	Wells 1	Aquifer	Average Depth (m)	Screen Length (m)
Bartlett	11	M	150	30
Collierville	12	M	100	20
Germantown	20	M	120	20
Millington	2	M	105	20
	4	F.P	460	20
Allen	27	M	140	30
Davis	19	M	140	25
Lichterman	23	M	145	30
LNG	2	M	115	30
Mallory	24	M	190	30
McCord	27	M	165	30
Morton	17	M	165	35
Palmer	4	M	120	25
Shaw	14	M	195	35
	3	F.P	360	35
Sheahan	24	M	170	30
DeSoto	84	M	105	20
	9	F.P	445	30

¹ Industrial and recreational use wells (e.g., golf courses) are not included. M = Memphis aquifer. F.P = Fort Pillow aquifer.

). More than 80% of the wells are screened in the upper two-thirds of the Memphis aquifer. Additionally, about 55% are screened in the upper third, where intervals of coarse-grained, transmissive sand are commonly found [4,35].

MODFLOW considers a well to be screened over the entire cell thickness, regardless of the actual screen length, and when a screen is split among layers, the total pumping will be apportioned using the percentage of screen located in that layer. To minimize the errors generated by these conditions, an analysis was performed to compute the optimum number of layers needed in the Memphis aquifer for which most of the wells will have the majority of their screen in only one layer (

Table 3. Percentage of wells included in the model, for which 80% or more of their screen will fall within a single layer. Thicknesses are approximated from the hydrogeologic surfaces developed by Hart et al. [58].

Vertical Discretization	Average Layer Thickness (m)	Wells > 80% Screen
3 layers	77	70%
4 layers	58	72%
5 layers	46	63%
6 layers	39	50%
7 layers	33	57%

).

Using the principle that a higher percentage of a well screen is optimum within a single layer, an individual threshold of 80% was set. Using four layers in the model with an average thickness around 60 m, results in 72% of the wells in the model had at least 80% of their screen within a single layer; this reduced the need for pumping apportionment and closed the gap between layer thickness and average screen length (30 m). Even though the use of more layers could reduce this gap (e.g., seven layers of ~30 m), each layer will mean an additional 150,000 cells; hence, increasing the computation time significantly.

2.1.1. Boundary Conditions

The upper boundary of the model corresponds to the water table, whereas its base (i.e., Old Breastworks confining unit) is treated as a no-flow boundary due to assumed negligible flow exchange between the Old Breastworks confining unit and older units [24,25,26,27,70]. Physically, the shallow aquifer is bounded along the west by the MRVA (not modelled west of the Mississippi River) and the Mississippi River, represented as a transient boundary with the MODFLOW river (RIV) package [64]. Specified heads were implemented for the shallow aquifer in the eastern part of Shelby County (see Figure 1c and Figure 4), following the outcrop of the unconfined Memphis aquifer where saturated thickness thins to zero in order to aid in model convergence.

The Memphis aquifer was set to a constant-head along the eastern boundary of the model area (Figure 4) to conform with the predominant flow direction and to introduce water into the system based on regional data presented in Schrader [48] which is consistent with historic water level readings between 2005 and 2016 at the nearby monitoring well Fa:R-002 (Figure 5), indicating an average annual deviation of only 20 cm. The southwestern boundary is also set as a constant head as supported by observations in the nearby wells Ar:C-001 and Ar:H-002 (see Figure 5), and based on a persistent groundwater divide mapped by Schrader [48,68]. The north, south, west, and northwestern boundaries correspond to no-flow conditions as model lines are perpendicular to the regional equipotential lines shown by Schrader [48,68] (See Figure 2a) and they are positioned far enough from the main aquifer stressors as to add or remove significant flow across the system.

The Fort Pillow aquifer was originally modelled with a constant-head boundary at the east, as supported by information from well Fa:R-001. Due to unnatural, forced patterns noted during simulation and lack of additional monitoring data, it was changed to a no-flow condition in subsequent models, without affecting simulation results. All the remaining boundaries in the model, including those of confining units, are set as no-flow because hydrologic gradients are parallel to the model boundaries.

2.1.2. Initial Conditions

Initial conditions for the shallow aquifer were obtained from the water table map developed by Narsimha [41] and modified by Ogletree [44] for conditions in Fall 2005 (see Figure 2b). Low water levels observed at certain breach locations were slightly increased to avoid initiating the model with a high probability of cells going dry in the first iteration, since in certain areas the shallow aquifer is so low that head levels are either dry or near the base of the formation [6] creating issues with model convergence as MODFLOW is unable to cope with a simulation that has dry cells in the first time step. Due to lack of data, head in the shallow aquifer outside Shelby County was estimated by subtracting 10 m from the DEM elevation as the average depth to the water table [41,44]. The initial head values in the Memphis aquifer were specified using data in Schrader [48] (see Figure 2a). Initial heads for the Fort Pillow aquifer were derived from observation well data from Fall 2005 [71].

2.1.3. Hydraulic Properties and Recharge

Using the information in Table 1, starting values for the parameters of hydraulic conductivity and specific storage for each hydrogeologic unit were set to averages. A specific yield was set to 0.25 and the anisotropy ratio was assumed to be 0.1, so that the vertical conductivity equals 10% of the horizontal [11,72], although higher anisotropy has been reported in the system [27]. Precipitation recharge was applied at the uppermost layer of the model, comprised of the shallow aquifer and the unconfined area of the Memphis aquifer [32]. An initial recharge rate of 1 cm/year was used following the average regional recharge extracted from the MERAS model of Clark and Hart [27].

2.1.4. Rivers and Streams

Exchange between groundwater and surface water was simulated using the MODFLOW (RIV) package [64]. The four rivers modeled are the Mississippi River and three of its tributaries: Wolf River, Loosahatchie River, and Nonconnah Creek (Figure 6). Model cells representing rivers were intersected using the stream location obtained from the USGS National Hydrography Dataset (NHD) [73] (see Figure 6).

To include stages at cells containing river locations, readings from ten river gauges were processed from the NWIS database [71]. As shown in Figure 7, month-to-month height deviations of less than 2 m are common for all rivers except the Mississippi, with changes between 1.5 and 3.5 m. Seasonal variations are marked by higher stages in winter and spring, with maxima usually between April and May and minima around September and October. Downstream gauges of the Wolf River and Nonconnah Creek (i.e., gauges 7031740 and 7032252) are affected by backwater conditions along the Mississippi River.

The rivers in Shelby County and surrounding area flow through unconsolidated sediments; hence, the riverbed is comprised of a mixture of sand, gravel, and clay. As there are no measures of riverbed conductivity nor a good determination of its thickness, a conceptual sensitivity analysis was performed to estimate the conductance per unit length required by the MODFLOW (RIV) package [64]. The groundwater model was run 16 times with varying conductance values and the observed heads in 15 shallow monitoring wells (see green markers in Figure 6) were compared to model heads, as well as extracting the flow budget computed by MODFLOW for the shallow aquifer for an average stage month (i.e., February, 2016) for an area inside Shelby County, delimited at the west by the bluff, and at the east by the outcrop of the Memphis aquifer under unconfined conditions.

The analysis was performed to determine a probable range of values for conductance per unit length whose range could accommodate both losing and gaining conditions along rivers in the model (Figure 8). According to calculated root-mean-square error (RMSE) and maximum and minimum leakage, conductance values between 0.5 and 10 m/day likely accommodate the possible ranges for the model. The RMSE does not increase significantly after 2.5 m/day, whereas cells adjacent to the rivers flood at conductance values lower than 0.1 m/day. It seems that values around 2.5 m/day will result in the lowest RMSE, and in the best combination of gaining–losing conditions, allowing the river to gain water from the aquifers during the dry season, and to lose it during certain times in the wet season as is expected in Shelby County.

2.1.5. Wells and Groundwater Pumpage

There are 336 high-capacity production wells inside the model area (Figure 1c and Figure 4) with a majority of them located within the ten Memphis Light, Gas, and Water (MLGW) well fields. Well depths range between 30 and 270 m in the Memphis aquifer, and 350 and 470 m in the Fort Pillow aquifer. Screen lengths range from less than 10 to more than 30 m. Industrial and recreational (e.g., golf courses) well depths and screen lengths were assumed to be the same, on average, as the MLGW production wells inside Shelby County as no construction data exist for them. Efforts to access well and pumping information for the other areas surrounding Shelby County were unsuccessful with the exception of DeSoto County, Mississippi.

Monthly pumping rates were obtained from the various well owners and state agencies. Other models used in the region typically distribute withdrawal rates over an area representative of the well field boundary using approximate locations, and most rely on 5-year water-use reports for yearly apportioning, subsequently limiting the model capabilities and temporal resolution [23,24,25,26,27,28]. This study simulates pumping at individual wells in the Memphis area, which affords greater accuracy in describing well activity (i.e., active, inactive, abandoned) within a well field.

As pumping rates for individual wells are not known, the monthly pumping rate for a well field was distributed among only the active production wells during that month. Inconsistent record keeping, data collection schemes, and database disparity created challenges in compiling pumping data for so many different wells and owners. The records were compiled into a single format, and missing data were populated using average values and values from trend analysis, mostly determined using custom scripts created for that task.

Monthly values were transformed to daily distributions for MODFLOW; this was caried out parsimoniously assuming a constant, daily pumping volume during each month. Withdrawal rates per well vary from less than 1000 m³/day to more than 5000 m³/day. Average monthly withdrawals over the study period reveal that the highest water demand occurs during summer, usually peaking in July–August (Figure 9). A slight decreasing trend in withdrawals over the course of the years modeled was also observed.

2.2. Calibration

Calibration of the model was achieved during the period January 2005 to December 2016. Water levels from monitoring wells were obtained through the NWIS database [71]. A total of 74 monitoring wells were used in the calibration process (Figure 10), distributed among the shallow, Memphis, and Fort Pillow aquifers by count: 15, 46, and 13, respectively. Additional monitoring points were created from water level maps [41,44,48,74], especially in the shallow aquifer, to provide some control in those areas lacking continuous information as recommended by Hill and Tiedeman [75].

The model was calibrated using automated Parameter ESTimation (PEST) [76] with a pre-process, manual, trial-and-error analysis to ensure convergence. The model was parametrized at discrete locations with the use of pilot points following guidelines by Fienen et al. [77] and Doherty and Hunt [78], providing an intermediate approach for characterizing heterogeneity [79]. Due to the ill-posed inverse problem (i.e., more parameters than observations), differences in aquifer characteristics, and a high number of pilot points, singular value decomposition (SVD-Assist) [78,80] and Tikhonov regularization [81] were used to ensure the most stable and geologically reasonable approach [82,83]. The calibration operation was expedited by using Parallel-PEST [76], which distributes model runs between computer cores, thereby reducing iteration time.

Each layer in the model incorporated a different number of pilot points depending on available monitoring data, distributing them first within a regularly spaced grid, to avoid “deserts”, then creating additional points to increase density in areas of interest and where monitoring density is high [79,84]. A total of 938 calibration points were used. The number of pilot points used for hydraulic conductivity alone is as follows: 119 in the shallow aquifer, 160 in the UCCU, 110 in the Memphis aquifer, 60 in the Flour Island confining unit, and 110 in the Fort Pillow aquifer. Additional parameters were riverbed conductance, specific storage, and recharge.

Following Hill and Tiedeman [75], composite scale sensitivities for individual pilot points were reviewed before re-running the process, holding them constant or rearranging them if their sensitivities represented less than 1% of the largest value. The best fit was obtained when the RMSE residuals between calculated and observed heads were reduced below a pre-set tolerance of ±2.5 m defined during preliminary calibration for the monitoring wells in all three aquifers. Further reduction in the error resulted in the loss of “hydrosense” [85], yielding geologically unrealistic parameters [78] or bulls-eye patterns, due to the inappropriate pursuit of a zero misfit [77].

3. Results and Discussion

The model was run enumerable times before using PEST to ensure convergence and a smooth simulation. Aquifer thickness interpolation errors were corrected in the grid and among layers, as well as in pumping and river stage assignment. Calibration using PEST was carried out a few dozen times, varying the density of pilot points, their placement, and the parameters calibrated per model run in an attempt to reach the best approach. Initially, less than 20 iterations taking between three to four days of continuous processing were required to evaluate the resulting parameters and make decisions for subsequent runs. After the parameters approached their optimal ranges, between five and six iterations were needed to yield the lowest error with the best parameter values.

3.1. Hydraulic Parameters

Resulting calibrated hydraulic conductivities are shown in Figure 11a,b with the majority falling within the range of previous studies (see Table 1) and generally accepted values [72]. Modelled values for the shallow aquifer have a slightly higher average than those in the Memphis aquifer (i.e., 20.9 m/day in contrast to 18.4 m/day), in accordance to the presence of gravel in the fluvial terrace and alluvial deposits. On the other hand, the Fort Pillow aquifer has a comparatively lower average hydraulic conductivity with 9.8 m/day—an expected result because of a higher proportion of finer sediments in this geologic unit [35].

Conductivity values for confining units (Figure 11b) are consistent with values for aquitards composed mostly of clay [72]. Both the UCCU and Flour Island confining unit have median hydraulic conductivities in the order of hundred thousandths m/day, although the Flour Island is less leaky and has a smaller range of conductivities. Higher values in the UCCU, mostly shown as outliers in the plot, are in accordance with those calibrated inside the breaches delimited by Parks [4] (e.g., average of 1.4 × 10⁻³ m/day). These results suggest the presence of areas with high conductance not previously identified as a breach location, which may warrant further investigation to determine their plausible role in inter-aquifer water exchange. In such cases, confining unit hydraulic conductivity values as high as 0.015 m/day suggest a greater presence of silt and sand than clay.

Calibrated, average specific storage values are 0.0014 m⁻¹ for the UCCU, 0.0033 m⁻¹ for the Memphis aquifer, 0.0018 m⁻¹ for the Flour Island confining unit, and 0.0019 m⁻¹ for the Fort Pillow aquifer, which all fall within the ranges described by Criner et al. [59] and Parks and Carmichael [55,56]. Final calibrated values of riverbed conductance per unit length ranged between 1 and 11.6 m/day, with an average of 3 m/day, all of which are within the expected values obtained in the pre-analysis (see Figure 8).

Modeled recharge values range from less than 1 cm/yr to 8 cm/yr with an average value of 1.8 cm/yr. Water balance studies on the unconfined area of the Memphis aquifer have reported recharge rates of less than 1.6 cm/yr [86] to more than 30 cm/yr [87,88,89]. Previously modeled values for the study area range from less than 1 cm/yr to more than 10 cm/yr [11,25,27,57,90,91].

Although the hydraulic values are within expected ranges, the modelled spatial heterogeneity of the pilot points does not necessarily reflect the real distribution of the parameters (e.g., hydraulic conductivity, storage, recharge). More physical measurements and aquifer testing are needed to reduce parameter uncertainty, especially where only head data are used for calibration [75]. Additionally, uncertainty exists regarding the validity of measured hydraulic conductivity values in the region, based on questionable adherence to aquifer test protocols [12].

3.2. Model Evaluation

The resulting measures of error of the calibration (i.e., mean error, mean absolute error and root-mean-square error [RMSE]) are shown in

Table 4. Average measure of error for each aquifer in this study and previous investigations.

Measure of Error (m)		Entire Model	Shallow	Memphis	Fort Pillow
Monitoring wells	This model	74	15	46	13
Mean error		0.0	0.1	0.4	−0.8
Mean absolute error		1.6	1.8	1.4	1.8
RMSE		2.0	2.2	1.8	2.3
		Reference	Shallow	Memphis	Fort Pillow
Previous models RMSE		h	NM	±1.5 *	NM
		i	NM	~16	~11
		j	NM	11.5	10.3
		k	NM	4.3	3.0
		l	NM †	10.9	NR
		m	NM †	9.8	NR

h Brahana [23], i Arthur and Taylor [24], j Arthur and Taylor [25], k Brahana and Broshears [26], l Clark and Hart [27] and m Clark et al. [28]. * For 75% of the control points. ^† Modelled the MRVA aquifer, RMSE ~ 4.5 m. NM = Not Modelled. NR = Not Reported.

, including the RMSE for previous models reviewed earlier.

The RMSE for this model represents an equivalent of only between 3% and 6% of the total change in head, when observed over a range of hydraulic head of 49 m for the shallow, 54 m for the Memphis, and 34 m for the Fort Pillow aquifers. Although the errors obtained are the lowest of all the models available for the Shelby County area, it is important to acknowledge that the nature of most of the models [24,25,27,28] is of regional scale, compromising any robust local comparison.

Comparison of time-series of simulated versus observed hydraulic heads through hydrographs allows for the visualization on the goodness-of-fit at those locations with monitoring wells. Figure 12 compares the monitoring wells previously circled in Figure 10. If the model is appropriately calibrated and successful at recreating the flow conditions of an area, the simulated head values will be within the ±2.5 m of the target shown as green bars (e.g., wells Sh:J-126 and Sh:M-040 screened in the Memphis aquifer, and Sh:U-001 screened in the Fort Pillow aquifer). Figure 12a, and to a lesser degree Figure 12b, converge towards the mean of the observations and demonstrate the model’s ability to simulate the pumping regime and the water level trends as the aquifer responds to stressors. Figure 12c does follow the observation heads but is unable to respond to the seasonal pumping as the model does not include the deep production wells in the area of West Memphis, where extensive withdrawal from the Fort Pillow aquifer takes place [26,55].

Nearly all the monitoring wells in the model follow the same pattern of a good fit seen in Figure 12a–c with a few exceptions. Well Sh:P-099 screened in the shallow aquifer (Figure 12d) is a good example when simulation followed a continuously diverging trend as compared with the observations, resulting in simulation errors between ±2.5 m and ±5 m (i.e., yellow bar), and above ±5 m (i.e., red bar), thereby suggesting a poor recreation of flow conditions. The well is located in a forested area in Overton Park (Memphis) right beside the Memphis Zoo and its numerous, aging (some circa 1906) water features, therefore allowing the shallow aquifer to receive water at a continual rate. However, in the model, these local water features are not included (because their characteristics are unknown) and recharge is applied at a constant rate in time. Therefore, the surface area proximal to this well contributes water so the model is weak in this area. Potentially, other areas with similar conditions to well Sh:P-099 might be poorly simulated due to the complex interactions that small water bodies, tributaries, and seasonal recharge create with the shallow aquifer.

3.3. Flow Budget and Inter-Aquifer Water Exchange

The simulated multi-layer aquifer system in Shelby County can be examined through the flow budget computed by MODFLOW, where the model calculates the volume of water that is added or removed from each cell. Inflows include recharge by precipitation, seepage from water bodies, lateral/vertical flow, constant heads, and water released from storage. Outflows consist of groundwater discharge to surface water bodies, pumping, lateral/vertical flow, constant heads or water entering storage. Lateral flow is defined as the total amount of water that is being released or added to a cell/area as a sum of the flow going through all the faces except for the upper and lower (Equation (2)).

$$\pm Lateral Flow=Right face+Front face+Left face+Back face$$

(2)

Change in storage occurs when inflow is not balanced by outflow. Flow into and out of storage is also considered part of the overall budget inasmuch as storage release effectively adds water to the flow system and it is considered an inflow, while accumulation in storage effectively removes water from the flow system and it is considered an outflow, even though neither process, unto itself, involves the transfer of water into or out of the groundwater system [92].

Even though inflows and outflows vary spatially and temporally, the model indicates a strong, seasonal pattern owing to monthly, simulated variations in river stage and pumping—year-to-year flow budgets remain nearly the same. To illustrate the complexity of inflows/outflows throughout the groundwater system, the month of August when pumping withdrawals typically peak, in 2016, was chosen to analyze the average daily flow budget in Shelby County (Figure 13).

Each layer and its flow components are examined in Figure 13 by showing the percentage that they contribute to the overall inflow (+) or outflow (−) of the hydrogeologic unit. Analysis of the shaded area in the shallow aquifer west of the bluff corresponds to the MRVA. This was not included because simulation showed that the flow in that area is controlled by the Mississippi River, resulting in a flow that is four times larger than the total flow of the area east of the bluff, which masks the shallow aquifer components, as also noted by Brahana and Broshears [26].

3.4. Multi-Layer Average Daily Flow Budget

Figure 13 shows that the main inflow component of the shallow aquifer is water released from storage (57.6%) followed by recharge and constant head which contribute 101 m³/min or 30%, combined. Upward leakage coming from the UCCU is less than 1% or 1.5 m³/min, an expected value since present head differences between the shallow and Memphis aquifers favor the downward movement of water. Downward leakage of approximately 35.9 m³/min represents 10.5% of the total outflows.

The upper left corner of Figure 13 shows a plot of the interaction between the shallow aquifer and the rivers during 2016. During the months of higher stage (i.e., April and May), the rivers lose water into the groundwater system, around 40 m³/min, a condition also noted by Graham and Parks [32] and Brahana and Broshears [26]. These results corroborate observations of areas and times when the rivers are under losing conditions [9,10,23,43,51,52,53]. Groundwater discharge from the shallow aquifer to the rivers (i.e., gaining conditions) accounts for over 52% of the overall outflows with a flow of 177 m³/min occurring during the highlighted month of August. This is close to the yearly peak of 229 m³/min in September and is bracketed by months when rivers are most commonly under baseflow conditions (see also Figure 7); this is corroborated by the opposing gaining conditions mentioned in previous studies during the rest of the year [4,9,23,32].

Downward movement of 36 m³/min from the shallow aquifer contributes a little over 57% of the total inflows to the UCCU; the rest of the inflows are supplied by the water released from storage at a rate of 26.7 m³/min, the composite outflow of 1.5 m³/min moving upward into the shallow aquifer, and 61.1 m³/min discharged and exchanged with the Memphis aquifer. The latter represents 97.5% of the total outflows from the UCCU. This reduction in storage water and flow through confining units has been documented by Czarnecki et al. [93] and Masterson et al. [94], where confining unit depletions due to storage release are of concern in certain areas of the U.S. because the water removed from clayey sediments cannot be replenished as these units gradually compress [94].

The Memphis aquifer inflow comes mainly from lateral contributions at a rate of 331 m³/min (56.4%). In eastern Shelby County (under unconfined conditions), rivers lose water to the Memphis aquifer at a rate of 15.5 m³/min. Downward vertical exchange from the UCCU represents an overall contribution of a little over 10% or 61.1 m³/min of the total inflow to the Memphis aquifer in Shelby County. According to Larsen et al. [37], the rate of movement from the UCCU to the Memphis aquifer is an estimated 1.6 cm/yr over an area of 2033 km². Larsen et al. [37] also found that young water tends to stay in the upper part of the Memphis aquifer resulting in higher mixing percentages. Isolating calculations to the upper 60 m of the Memphis aquifer (i.e., model layer 4) resulted in UCCU leakage (i.e., Memphis aquifer inflow) of almost 30%.

The main outflow and stressor to the Memphis aquifer is generated by pumping which accounts for 85.3% of the total extractions or 500 m³/min. Layer apportioning shows that the upper and middle areas are more stressed (i.e., 117 and 256 m³/min, respectively), and accounts for approximately 70% of the extractions. Upward exchange from the Flour Island confining unit to the Memphis aquifer is negligible (<0.1%), and downward leakage to the confining unit is 4.5 m³/min. Exchanges between the four layers simulated in the Memphis aquifer represent around 590 m³/min of internal water redistribution.

The Flour Island confining unit receives almost 95% of its inflow from the Memphis aquifer (4.5 m³/min), exchanging 1 m³/min of the previous inflow as leakage to the underlying Fort Pillow aquifer. Opposite to the conditions of the shallow aquifer and UCCU, and because of the difference between inflow and outflow, around 3.7 m³/min of groundwater is being removed from the flow system and going into storage accumulation to replenish the confining unit (referred to as outflow in Figure 13).

Regarding the Fort Pillow aquifer, downward leakage from the Flour Island confining unit, of about 1 m³/min, accounts for less than 2% of the inflow to the Fort Pillow aquifer, which receives most of its inflow via lateral flow. Outflow is mostly driven by lateral flow, and a discharge of 16.5 m³/min or 28% due to pumping. Upward leakage to the Flour Island confining unit is negligible.

3.5. Water Exchanged between the UCCU and the Memphis Aquifer

Figure 14 includes only the average daily water exchanged between the UCCU and the Memphis aquifer for August 2016, as the water coming from and into the UCCU in relation with the aquifer. Shelby County was segmented into quadrants (A-X) and the breaches Parks [4] enumerated (1–8) to locally compute the leakage in each of them.

The analysis of water exchanges, as shown in Figure 14, was compiled to quantify the interaction that the UCCU has with the Memphis aquifer at individual areas, selecting those cells within the zone of interest and using MODFLOW’s flow budget to compute flow values. Segmentation of the county into quadrants indicates that present-day gradients favor downward movement of water into the Memphis aquifer over much of Shelby County, with greater interactions occurring in quadrants with breaches or within areas of higher conductivity.

Quadrants near the Mississippi River, and where the confining unit is mostly absent (i.e., quadrants E, R, S, T, and X), show contributions of less than 1 m³/min, which indicates that leakage might not be significant. However, it is recognized that quadrants E, R, S, and X do not have much coverage by the UCCU which would be a contributing factor to a lower contribution. Areas where incoming leakage accounts for more than 5% of the total inflow budget are in quadrants A, B, K, L, O, P, Q, U, and V, concentrating more than 70% of the total leakage with most of these containing a breach mapped by Parks [4] or an anomalous water depression (Figure 2), except for quadrants A, B, and K. If only the upper 60 m of the Memphis aquifer are considered, the percent of exchanges can fluctuate between less than 1%, to more than 35% in quadrant V with its 8.6 m³/min moving at a rate of 4.5 cm/yr. An area of interest that does not contain any mapped breaches by Parks [4], but has an anomalous water depression, is quadrant O showing a 7% water exchange between the UCCU and the Memphis aquifer, thus reaffirming the plausibility of a breach inferred by Bradshaw [42].

The areas delimited as breaches by Parks [4] are analyzed on the right bottom side of Figure 14. Water contributions from the UCCU to the Memphis aquifer can range from less than 0.3% (<0.05 m³/min) to more than 20.8% of the total inflows in breaches 8 and 4, respectively. It is also notable how the flow through breach 4 exceeds the combined flow from the rest of the breaches, with a water exchange of 3.2 m³/min. This high flow rate is attributed mostly to the larger size of the breach as compared to the others, realizing that breach 5, though also large, resides to the west of the bluff where contribution by the overlying shallow aquifer does not occur. Total flow from the UCCU through Parks [4] breaches (i.e., 4.5 m³/min) accounts for only 7.4% of the total water exchange occurring in Shelby County (i.e., 61.1 m³/min) with the Memphis aquifer, suggesting the presence of more breaches within the UCCU.

4. Conclusions and Implications

Shelby County relies almost entirely on the Memphis aquifer to satisfy its water demand, and although the UCCU protects and limits the exchange of water between the Memphis aquifer and the overlying shallow aquifer, localized discontinuities (i.e., breaches) in the confining unit provide preferential avenues for increased inter-aquifer water exchange. As there are numerous potential surficial contaminant sources proximal to the breaches [37], this water exchange adds concern about water quality degradation. A multi-layered three-dimensional model was developed to simulate inter-aquifer water exchange in a complex hydrogeologic system, with unconfined and confined conditions within individual hydrogeologic units, large seasonal detailed pumping withdrawals, surface water interactions, and recharge.

Even though the simulation shows that the vast majority of the rivers are gaining water from the shallow aquifer during most of the year, losing conditions still locally exist during specific times of the year, especially the wetter months of April and May. With the rivers also traversing the unconfined portion of the Memphis aquifer, similar concerns for water quality degradation exist with the general gradient of the Memphis aquifer being east to west toward the large depression beneath much of downtown Memphis.

Downward leakage from the UCCU into the Memphis aquifer across Shelby County accounts for more than 10% of all the inflow to the aquifer, totaling 61 m³/min moving at a rate of 1.6 cm/yr. Evaluation of the upper 60 m of the Memphis aquifer indicates that leakage through the UCCU results in mixing of as much as 29% young water in the upper part of the aquifer. Water exchange between the Memphis aquifer and the Flour Island shows a downward leakage of 4.5 m³/min from the aquifer into the confining unit, with a negligible flow contribution going upwards into the aquifer.

Spatially discrete water exchange analysis using quadrants shows that leakage is mostly concentrated in areas with previously mapped breaches in the UCCU and water-table anomalies in the shallow aquifer. Incoming flow from the UCCU into the Memphis aquifer ranges from less than 0.05 m³/min to as much as 8.7 m³/min. Considering the individual breaches mapped by Parks [4], flow contribution through the UCCU into the Memphis aquifer is generally less than 0.5 m³/min, with a maximum of 3.2 m³/min or 20.8% of the total breach area inflow. This value can be as much as 53% of the inflow when looking at the upper 60 m of the Memphis aquifer in a particular breach area.

The simulated flows coming from the UCCU and the Flour Island confining unit into the Memphis aquifer are lower than the those estimated in previous models [23,24,25,26], yet are in accordance with estimates based on age-dating studies and tracer-based mixed percentages obtained for the Memphis aquifer [6,37,95]. The present results argue that simplifications in previous models, such as (1) the use of a coarse grid, (2) regionalization, (3) indirect simulation of the shallow aquifer, rivers and confining units, (4) use of manual calibration techniques, and (5) inaccurate (or overly generalized) pumping rates and well locations, resulted in overestimated leakage to the Memphis aquifer in Shelby County.

Although this study identified areas of inter-aquifer water exchange proportions that are within the range of published values, the rates and percentages should be treated with caution. Most of the model uncertainties can be related back to confining unit insensitivity during the calibration process, absence of information on riverbed characteristics, missing data on pumping and well design characteristics (especially in DeSoto County), and lack of monitoring data outside Shelby County and generally in the shallow aquifer. Regarding the monitoring data, long periods of record and observation proximal to breaches are lacking.

The correspondence of localized leakage in parts of Shelby County with previously identified breach areas and anomalous water-table depressions suggests that determinative methods of assessing heterogeneity in confining units may provide useful model constraints for modeling inter-aquifer water exchange. Although other methods also have proven useful for identifying locations and extent of inter-aquifer exchange through the UCCU in Shelby County [22,61,96], determinative methods have the advantage of potentially providing more predictive behavior.

Although locally unconfined conditions in the Memphis aquifer exist within Shelby County [4,97], the resulting partial saturation of the aquifer, even beneath local streams, appears to have limited influence on leakage. This assertion is exemplified by the greatest amount of leakage being simulated adjacent to the Lichterman well field, where the Memphis aquifer is unconfined [97]. Similarly, annual variations in pumping stress and the slight decrease in total annual pumping do not appear to affect overall leakage rates and percentages, although the limited monitoring data constrain the time period of calibration and application of the results.

This newly developed model is a great improvement upon previous efforts [23,24,25,26,27,28], as it is capable of explicitly modeling breaches in the UCCU and does not use a leakance term like previous quasi-3D models [23,24,25,26]. At the same time, this model was not developed with a regional intent [24,25,27,28], which increases its guidance capabilities for more localized studies and can provide useful information from a well field, to a city level.

This is the first model in the Mississippi embayment area that completely simulates the shallow aquifer within Shelby County as past models used constant cells for the entire area [23,26], used river cells to simulate the head as a dependent flux [24,25], or it was simply not modeled [27,28]. Since the shallow aquifer is connected to the Memphis aquifer in certain areas, it is of great importance to be able to explicitly simulate its flow and the inter-aquifer interactions that occur with the other units, the rivers, and the multiple recharge processes.

This model takes advantage of using PEST and pilot points during the calibration process, using 74 monitoring wells with historical data and resulting in the lowest error of all the models covering Shelby County [23,24,25,26,27,28], with hydraulic parameters well in accordance to those published for the area. Furthermore, this investigation provides an update for inter-aquifer water exchange quantification, as the most recent flow analysis was developed for 1980 conditions [26] and does not reflect current water uses, nor does it include the level of detail needed to characterize the hydrogeologic system at present.

The implication of the modeled inter-aquifer water exchange is that more monitoring needs to occur in areas where exchange may lead to leakage of poor-quality water from the shallow aquifer into the Memphis aquifer. This should be coordinated with recurrent investigations, focusing on the state of the groundwater system (e.g., water levels, contaminant potential, declines in saturated thickness, etc.), and better characterization of and/or information on the external stressors, such as the rivers, pumping rates, and recharge.

The benefit of numerical models, such as that developed in this study, is that they are dynamic and can be updated when new information becomes available. Potential additional applications of this model include pumping rate projections or configurations, delineation of capture zones and wellhead protection areas, and coupling with a solute-transport model, among other uses. The successful application of models in the Memphis area suggests that similar detailed modeling endeavors may be useful in assessing inter-aquifer water exchange in other layered aquifer systems with similar stressors.