Application of Numerical Groundwater Model to Determine Spatial Configuration of Confining Unit Breaches near a Municipal Well Field in Memphis, Tennessee

The aquifer system of interest is part of the Mississippi embayment and adjoins the Gulf Coastal aquifer system to its south. The Mississippi embayment is a geologic syncline that descends toward the Gulf of Mexico and whose axis approximates the pathway of the present-day Mississippi River (Graham and Parks 1986). It is filled with more than 1,000 m of sediments of Cretaceous, Tertiary, and Quaternary age with thinning of the geologic units toward updip boundaries to the north, east, and west (Fig. 2). Many of the shallower aquifer systems in the embayment provide fresh water, whereas deeper saturated units are saline (Waldron et al. 2011).

Centrally located in the Mississippi embayment, the City of Memphis, Shelby County, Tennessee, obtains all of its potable water from groundwater. The Memphis aquifer is the primary aquifer used for drinking water and industrial use, with some supplemental supply from a deeper aquifer, the Fort Pillow aquifer. Overlying the Memphis aquifer is an aquitard, an unconfined aquifer, and a surficial unit (Fig. 2). The surficial unit is wind-blown silt or loess, which blankets Shelby County in thicknesses ranging from 25 m near the bluffs of the Mississippi River along the western boundary to 3 m approaching the eastern boundary (Gentry et al. 2006b; Pell et al. 2005). Beneath the loess is an unconfined aquifer called the shallow aquifer. The shallow aquifer is composed mostly of sand and gravel (Van Arsdale et al. 2007; Parks 1990); its horizontal hydraulic conductivity ranges from 1.5 to $45\mathrm{m}/\mathrm{day}$ (Graham and Parks 1986).

Below the shallow aquifer is an aquitard, called the upper Claiborne confining unit (UCCU), that is mostly clay with silt, fine sand, and lignite. Its thickness ranges from 0 to 110 m, where zero thickness relates to its subcrop along the eastern border of Shelby County (toward the recharge zone) and breaches (Graham and Parks 1986). The vertical hydraulic conductivity of the UCCU ranges from $1.5\times {10}^{-6}$ to $3.0\times {10}^{-4}\mathrm{m}/\mathrm{day}$ in the absence of breaches (Robinson et al. 1997). The hydraulic conductivity of breaches in UCCU has not been thoroughly determined; however, a falling-head slug test conducted by Gentry et al. (2006a) at Well MS-12 in Shelby Farms in Shelby County, where the confining unit was found absent, measured a hydraulic conductivity of approximately $0.1524\mathrm{m}/\mathrm{day}$ ($0.5\mathrm{ft}/\mathrm{day}$). The UCCU confines the underlying Memphis aquifer beneath most of Shelby County, providing a layer of protection from contamination.

The Memphis aquifer is a prolific freshwater aquifer composed mostly of fine to coarse sand and ranging in thickness from 152 to 271 m. Its horizontal hydraulic conductivity ranges from 8.5 to $47\mathrm{m}/\mathrm{day}$ (Brahana and Broshears 2001; Gentry et al. 2006a; Parks and Carmichael 1990; Waldron et al. 2011). Separating the Memphis and deeper Fort Pillow aquifers is the Flour Island Formation, composed mostly of clay, with a thickness range of 49–94 m in southwest Tennessee. The Fort Pillow aquifer is also composed of sand but is thinner than the Memphis aquifer, with a thickness range of 38–93 m (Graham and Parks 1986).

Prior to the countywide groundwater model developed by Villalpando-Vizcaíno (2019), many numerical models were developed to simulate hydrologic budgets and evaluate groundwater availability the Mississippi embayment aquifer system (Arthur and Taylor 1990, 1998; Brahana and Broshears 2001; Brahana 1982; Clark et al. 2013; Clark and Hart 2009); however, these models were unable to accurately represent leakage through local breaches due to the absence or misrepresentation of the shallow aquifer (Villalpando-Vizcaíno 2019). Clark and Hart (2009) and Clark et al. (2013) did not simulate the shallow aquifer in western Tennessee. Brahana (1982) and Brahana and Broshears (2001) represented the shallow aquifer using constant-head values over its entire extent. Arthur and Taylor (1990, 1998) represented the shallow aquifer as river cells with a degree of hydraulic connection with the Memphis aquifer based on streambed conductance. Recognizing these limitations, the simulated hydrologic budgets determined by Arthur and Taylor (1990) estimate overall downward leakage through the UCCU of about 35% of the water withdrawn from the Memphis aquifer. Leakage through the UCCU is representative of the UCCU as a whole and not exclusive to breaches. Similarly, Brahana and Broshears (2001) estimated that leakage through the UCCU to the Memphis aquifer accounted for as much as 54% of its total inflow, identifying it as a major component of the hydrologic budget in the Memphis area that is not uniformly distributed. Other sources for the water withdrawn from the Memphis aquifer are recharge in the outcrop area, leakage from the Fort Pillow aquifer, and water released from storage (Brahana and Broshears 2001).

A multilayered three-dimensional (3D) groundwater model of the upper Mississippi embayment aquifer system beneath Shelby County was developed by Villalpando-Vizcaíno (2019) to simulate groundwater in the shallow, Memphis and Fort Pillow aquifers, and the exchange through the intervening aquitards. Unlike its predecessors, Villalpando-Vizcaíno’s (2019) model treated hydrologic conditions in the shallow aquifer as variable head under transient conditions, and it also included areas of enhanced leakage through the UCCU based on eight spatial configurations for breaches inferred from Parks (1990), one of which is located in the area of interest of this study and will be analyzed further (Fig. 3). Specific to leakage between the shallow aquifer and the Memphis aquifer, Villalpando-Vizcaíno (2019) determined that overall leakage of $\mathrm{87,840}{\mathrm{m}}^3/\mathrm{day}$ occurred from the UCCU into the Memphis aquifer. This accounted for 10% of the total recharge to the upper 60 m of the aquifer, with as much as 20% in some localized areas. In comparison to previous leakage estimates, these were based on the 2016 simulated flow budgets for subregions within Shelby County, of smaller area than those previously discussed. As such, leakage from the UCCU that could flow laterally from adjacent subregions is neglected, leading to smaller leakage averages than those achieved if all subregions were analyzed as one. However, because no previous studies have thoroughly determined the hydraulic conductivity through the breaches, the hydraulic conductivity of the UCCU was treated as a calibration parameter to achieve the best fit between simulated and observed hydraulic heads (Villalpando-Vizcaíno 2019). The Villalpando-Vizcaíno (2019) numerical model was used in the present investigation of the Sheahan well field in conjunction with particle tracking and groundwater age dating to advance the understanding of groundwater exchange through these breaches and to evaluate previously proposed spatial configurations.

The Sheahan well field (Fig. 3), in operation since 1932, is one of MLGW’s 10 well fields and consists of 24 production wells screened in the Memphis aquifer along with a water treatment plant with a capacity of $0.13\mathrm{million}{\mathrm{m}}^3/\mathrm{day}$. The production wells have screens that are generally 24–30.5 m in length at depths ranging from 91 to 236 m below the ground surface (Ivey et al. 2008; Larsen et al. 2013). Downward leakage from the shallow aquifer into the Memphis aquifer near the Sheahan well field has been suspected since 1965. According to Nyman (1965), a segment of Nonconnah Creek south of the Sheahan well field behaved as a losing stream during the fall or dry season when groundwater levels are typically at their lowest and perennial stream baseflows are sustained through groundwater discharge (Ogletree 2016; Parks 1990). This segment of Nonconnah Creek was observed to be completely dry, while the channel had some flow in the upstream and downstream reaches (Nyman 1965). A later investigation by Graham and Parks (1986) measured tritium concentrations in observation wells Sh:K-073 and Sh:K-074 (Fig. 3), located in the southern part of the Sheahan well field and screened in the Memphis aquifer, that indicate a component of modern water had entered the aquifer. Enhanced leakage between the shallow and Memphis aquifer due to the compromised nature of the aquitard was also substantiated by a depression in the water table and a distorted temperature gradient in the Memphis aquifer in the southern part of the well field (Graham and Parks 1986). In the Sheahan well field, Parks (1990) inferred the location of a potential breach through the identification of a water table depression and thinning of the UCCU west of the Sheahan well field using an isopach map (Fig. 3).

In more recent studies, Larsen et al. (2003) determined 6% to 32% of the water being extracted by production wells in the Sheahan area is from the shallow aquifer according to major solute chemistry and tritium and tritium/helium-3 data. In contrast with leakage estimates from models previously discussed, this approach is based on the actual composition of water samples extracted from the production wells at Sheahan during fall 2000, with the location and depth interval of the production well screens playing a significant role in the fraction of downward leakage through the UCCU that is extracted, with variability at a much finer scale than that captured by the models. Additionally, leakage estimates based on tritium/helium-3 data only account for downward leakage with detectable tracer amounts (waters less than 60 years old), making these expectedly lower. S-wave reflection surveys were conducted in the southern part of the well field, finding an erosional feature or paleochannel in the confining unit (Pell et al. 2005). The rotosonic drilling at MLGW-99s in the central part of the well field (Fig. 3) indicated the absence of clay-rich deposits up to 45 m below the ground surface where the top of the UCCU was expected to be found (Gentry et al. 2006b). Ivey et al. (2008) used inverse age-distribution modeling with environmental tracer data to determine the most probable location of a leakage source impacting the Sheahan production wells. The location identified by Ivey et al. (2008) is in partial agreement with the breach inferred from Parks (1990), though differing in size, and includes the location drilled by Gentry et al. (2006b). Finally, a hydrologic study of Nonconnah Creek as a probable recharge source to the Sheahan area was conducted by Larsen et al. (2013). They interpreted a leakage flow path from the creek into the shallow aquifer that flows along a paleochannel toward the well field.

The regional 3D groundwater flow model by Villalpando-Vizcaíno (2019) consists of a transient finite-difference model discretizing the study area into a grid of 388 rows by 391 columns and eight layers, with uniform square 250-m cells. The first layer represents the shallow aquifer, the second layer represents the UCCU, Layers 3 to 6 represent the Memphis aquifer, and Layers 7 and 8 represent the Flour Island confining unit and Fort Pillow aquifer, respectively. The water table is the upper boundary of the model, while its base, representing the Old Breastworks confining unit underlying the Fort Pillow aquifer, was treated as a no-flow boundary. The shallow aquifer is bounded along the west by the Mississippi River as a transient boundary, and bounded along the east by the outcrop of the unconfined Memphis aquifer. The eastern and southwestern boundaries of the Memphis aquifer were defined as constant-head boundaries consistent with historic water level data (Villalpando-Vizcaíno 2019). The eastern boundary of the Fort Pillow aquifer was set as constant head to accommodate forced patterns adjacent to the boundary, while all the remaining boundaries in the model were defined as no-flow boundaries because they are parallel to hydrologic gradients (Villalpando-Vizcaíno 2019).

Variable pumping rates and river stages were defined for the simulation period of January 2005 to December 2016. The model was built using the USGS MODFLOW-NWT program (Niswonger et al. 2011). MODFLOW-NWT was chosen because of its ability to handle dry cells and rewetting (Villalpando-Vizcaíno 2019), both of which were common occurrences when simulating the shallow aquifer. MODFLOW-NWT is executed within the pre- and postprocessing software GMS, v.10.3 developed by Aquaveo. Calibration was carried out using observations from 74 monitoring wells distributed among the shallow, Memphis, and Fort Pillow aquifers (15, 46, and 13 observation wells, respectively), and additional monitoring points were extracted from available water table maps (Kingsbury 2018; Konduru 2007; Ogletree 2016; Schrader 2008). The model was parametrized with the use of pilot points of varying density and placement to allow for heterogeneity, leading to a total of 938 calibrated parameters including the hydraulic conductivity, specific storage, riverbed conductance, and recharge. For hydraulic conductivity alone, more than 100 pilot points were used for every hydrogeologic unit except for the Flour Island confining unit, which used 60 pilot points. The values of the pilot points were calibrated using automated parameter estimation (PEST) (Doherty and Hunt 2010), and targeting the reduction of root-mean-square error (RMSE) residuals between calculated and observed heads below $\pm 2.5\mathrm{m}$ for the observation wells in the shallow, Memphis, and Fort Pillow aquifers because a further reduction of the error was found to yield unrealistic parameters (Villalpando-Vizcaíno 2019). Hydraulic conductivities resulting from calibration are shown in Fig. 4, while plots of calculated versus observed heads over time for two observation wells in the Sheahan area are presented in Fig. 5. Additionally, measures of error for the calibrated model are shown in Table 1.

To simulate flow paths of modern water into the Memphis aquifer and more specifically the Sheahan well field, the regional model was extended back to January 1960. This date was chosen because groundwater age dating (discussed subsequently) identifies waters $<60\mathrm{years}$ old. The monthly pumping rates for all the production wells in MLGW well fields for years prior to 2005 were calculated from the cumulative monthly pumping rates for each well field in MLGW’s summary of operations reports, where available, and evenly distributing the total rate among all the active production wells by well field. The production wells were considered to be active for all the months following the month in which they were constructed according to their construction logs. For the years when monthly pumping rates per well field were not available, rates were estimated based on the monthly pumping rates assigned to the previous year, and the changes in the average daily pumpage of all MLGW production wells and the number of active wells with respect to the previous year.

Using MODFLOW-NWT, Villalpando-Vizcaíno (2019) observed that the numerical model remained very sensitive to cells exhibiting low saturation or desaturation/rewetting. Such conditions occur proximal to breaches where water loss from the shallow aquifer to the Memphis aquifer partially to completely drains the shallow aquifer. Because breaches would locally drain the shallow aquifer, modeling the transient water table as far back as 1960 was identified as a challenge and beyond the scope of the Villalpando-Vizcaíno (2019) study due to data constraints, mainly because the oldest water table map for the area was for fall 1987 (Parks 1990). Due to this limitation, the water levels in the shallow aquifer were changed to constant head using Ogletree’s (2016) values for 2005, based on measurements made in fall of 2005 by Konduru (2007), although it is recognized that these conditions are not representative of the extended simulation period. The 2005 water table surface is considered to be the oldest reliable water table data for the Memphis area because the 1987 water table surface (Parks 1990) used measurements from 1944 to 1987, assuming all water levels were unaffected during the 43-year period by increased pumping from the Memphis aquifer (Konduru 2007).

Villalpando-Vizcaíno (2019) simulated exchange between the shallow aquifer and Nonconnah Creek. River stages vary and needed to be extended back in time as well. The river stages defined for the period of 2005–2017 were repeated back to January 1960 at 12-year intervals. The influence of variable river stage is minimized due to the constant head in the shallow aquifer because significant interaction with the Memphis aquifer only occurred at the outcrop of the Memphis aquifer, more than 11 km east of the well field, and then with a negligible impact on the study area.

The error of the extended model was assessed through the mean absolute residual of the simulated versus observed heads in the Memphis aquifer around the Sheahan well field, based on the available water level data for Wells Sh:P-061, Sh:P-076, Sh:K-021, Sh:K-066, Sh:K-110, and Sh:K-122 in the USGS National Water Information System (USGS 2019). Using these six observation points, the total mean absolute residual for the extended model (January 1960 to January 2017) was found to be 3.91 m, whereas the mean absolute residual for the unmodified model by Villalpando-Vizcaíno (2019) was 3.09 m. When conducting the same analysis for two separate time periods, prior to and after January 1, 2005, the mean absolute residual of the extended model was found to be 5.97 m for the 1960–2005 period and 1.85 m for the 2005–2017 period. Thus, a smaller error was found around Sheahan for this period when compared to the original model (1.85 against 3.09 m). The higher error in the earlier period is attributed to the data limitations, such as the unknown pumping rates for some of the simulated years that had to be estimated and the constant-head assumption for the shallow aquifer.

The extended model was used to determine likely spatial configurations of breaches in the Sheahan well field by comparing parameter estimates obtained through the analysis of flow paths to published groundwater data for the production wells at Sheahan. Five breach spatial configurations were analyzed (Figs. 6 and 7): (1) Parks’s (1990) breach to the west of the well field (referred to as PB); (2) breach inferred from Ivey et al. (2008) in the central part of the well field (referred to as IB); (3) a combination of the Parks (1990) and Ivey et al. (2008) breaches that includes additional areas where stratigraphic control is scarce (termed large breach and referred to as LB); (4) a paleochannel-like configuration (termed paleochannel and referred to as PC) interpreted from the interpolated surface of the top elevation of the UCCU in the vicinity of the well field; and (5) a wider interpretation of the paleochannel from the previous configuration (termed large paleochannel and referred to as LPC) extending further from the well field. These two paleochannel configurations (Fig. 7) do not represent breach configurations proposed by previous studies; however, they are based on the presence of a paleochannel within the UCCU suggested by Pell et al. (2005) and their locations are partially in agreement with the paleovalley identified by Larsen et al. (2013). The plan view areas of the breaches PB, IB, and LB range from 2.19 to ${6.06\mathrm{km}}^2$ (Fig. 6), while the plan view areas of the paleochannels PC and LPC are 5.19 and ${11.31\mathrm{km}}^2$, respectively (Fig. 7). The areas of the spatial configurations for the breach LB, PC, and LPC are significantly larger than the areas of PB and IB, and thus represent greater vertical leakage in the extended model.

To represent the breach spatial configurations in the finite-difference grid of the extended model, all the cells contained or intersected by the configurations were converted to breach cells, sharing a different hydraulic conductivity value than the rest of the cells in the second layer of the model representing the UCCU. Because the cells have an uniform size of $250\times 250\mathrm{m}$, the actual representation of linear features such as the paleochannels in the model is akin to channels that are at least 250 m wide, whereas the paleochannel identified by Pell et al. (2005) is approximately 140–350 m wide.

Physical measures of hydraulic conductivities of breaches in Shelby County are almost nonexistent. Gentry et al. (2006a) investigated a breach proximal to the former Shelby County landfill approximately 8 km from the Sheahan well field. From a single falling-head slug test, Gentry et al. (2006a) measured the horizontal hydraulic conductivity of the breach to be approximately $0.1524\mathrm{m}/\mathrm{day}$ ($0.5\mathrm{ft}/\mathrm{day}$). Villalpando-Vizcaíno (2019) attempted to use model calibration to ascertain breach properties, such as hydraulic conductivity, based on sparse measured values. From his calibration, Villalpando-Vizcaíno (2019) determined an average horizontal hydraulic conductivity of $1.4\times {10}^{-3}\mathrm{m}/\mathrm{day}$ for all the breaches inferred from Parks (1990), and an average value of $1.8\times {10}^{-3}\mathrm{m}/\mathrm{day}$ for the pilot points representing the breach in the Sheahan area, specifically.

For this investigation, the breaches in the five scenarios were simulated with a horizonal hydraulic conductivity matching that of Gentry et al. (2006a) that was then varied an order of magnitude above and below in order to capture unknown variability and account for its uncertainty. Given this, a total of 15 models were analyzed representing five different breach spatial configurations under three different hydraulic conductivity values. The regional model by Villalpando-Vizcaíno (2019) is considered to be optimally calibrated because its time frame was defined considering data availability for the shallow aquifer. Therefore, recalibration of the extended model was discarded due to the added uncertainty from the data shortcomings previously discussed. However, particle tracking analysis revealed that the hydraulic conductivity values determined for the breach at Sheahan in the calibration were unable to provide pathways of modern water to the well field. Consequently, the value determined by Gentry et al. (2006a), approximately two orders of magnitude greater than the average calibrated value for the breaches by Villalpando-Vizcaíno (2019) of $1.4\times {10}^{-3}\mathrm{m}/\mathrm{day}$ [Fig. 4(b)], was favored to represent the breach because it allowed for greater downward leakage and a significant presence of modern water in the Memphis aquifer more in line with what is supported by the age-dating data (Gentry et al. 2006b; Larsen et al. 2003, 2016). The vertical hydraulic conductivity of the breaches was assumed to be 1/10th of its horizontal value (Freeze and Cherry 1979). Because predominantly vertical flow through a confining unit breach is governed by the vertical hydraulic conductivity, the focus will be shifted to this parameter instead of its horizontal component.

Flow paths were produced using USGS MODPATH v7.2.01 (Pollock 2016) in GMS. For this analysis, 20 particles were placed within each cell containing a well and traced backward in time from January 2017 to January 1960. In the case of production wells MLGW-88, MLGW-87A, MLGW-63A, MLGW-98, MLGW-93, and MLGW-96, when their screen length penetrated two vertical grid cells, particles were placed in both cells. Flow paths were analyzed in order to calculate estimates for the modern water percentage and the apparent age of the modern water extracted by the production wells in the Sheahan well field. These two estimates were calculated following different approaches.

Estimates for the apparent age of the modern water extracted by the production wells were obtained as the average travel time of the particles (i.e., advection) (Pollock 2016) that would have originated from the shallow aquifer, passed through a breach, and then reached a production well. This study recognizes the possibility of a particle going through a cell in the Memphis aquifer in which the presence of modern water could have been suggested by other flow paths, such as a particle coming from the shallow aquifer that previously reached that cell. However, the apparent age estimates represent average travel times of particles with modern water fractions coming from the shallow aquifer. The apparent age of modern water fractions involving the pathway of more than one particle and not directly traceable to the shallow aquifer, such as those in overlapping well capture zones, were not determined because they were considered to introduce greater uncertainty.

To obtain an estimate of modern water percentage, a mixed model using number of particles and normalized cell volumes was developed in Microsoft Excel. A grid mirroring the finite-difference grid proximal at Sheahan was recreated in Excel comprising 26 rows by 26 columns and six layers. A custom VBA script processed the path line report from MODPATH of the groundwater models (MODPATH output) and tracked the particles as they moved forward in time, calculating new modern water percentages for each cell. Initial conditions were split between the following scenarios: (1) modern water (100%) is assigned only to shallow aquifer cells while the remaining cells are 0% (lower bound); and (2) modern water (100%) was assigned to the shallow aquifer and breaches while the remaining cells are 0% (upper bound). Volumetric calculations were based on cell size (i.e., constant $250\times 250\mathrm{m}$) and saturated thickness, and then were normalized (Table 2). The modern water percentage for all cells was calculated using Eq. (1)where ${\mathrm{MW}\mathrm{P}}_{\mathrm{out}}$ = modern water percentage of the cell the particle is leaving; ${\mathrm{MWP}}_{\mathrm{in}-1}$ = modern water percentage of the cell the particle is moving into, before mixing; ${\mathrm{MWP}}_{\mathrm{in}}$ = modern water percentage of the cell the particle is moving into, after mixing; ${\mathrm{NCV}}_{\mathrm{out}}$ = normalized cell volume of the cell the particle is leaving; and ${\mathrm{N}\mathrm{C}\mathrm{V}}_{\mathrm{in}}$ = normalized cell volume of the cell the particle is moving into. For the wells with screens split across two cells, the upper- and lower-bound values were calculated as weighted averages using the fraction of the screen length in each cell as the weight (Table 3).

The potential spatial configurations for the breach around the Sheahan well field were ranked against one another and scored according to the mean residuals and mean absolute residuals when comparing (1) the simulated water levels to historic water level measurements around Sheahan available in the USGS National Water Information System (USGS 2019) for Wells Sh:P-061, Sh:P-076, Sh:K-021, Sh:K-066, Sh:K-110, and Sh:K-122; (2) simulated versus published apparent ages of the modern water in Sheahan production wells (Gentry et al. 2006b; Larsen et al. 2003, 2016); and (3) simulated modern water percentages against published tracer-based mixing percentages (Larsen et al. 2003, 2016). The model with the smallest residual in a given criterion, of the three previously described, was ranked as first for that criterion and the score of the model was calculated as the sum of the ranks obtained for the three criteria, with a lower score being more favorable.

Of the 24 production wells in the Sheahan well field, only Wells MLGW-78B, MLGW-55B, MLGW-80A, MLGW-99, MLGW-87A, MLGW-88, and MLGW-86R have published data on both their modern water percentages (Larsen et al. 2003, 2016) and apparent ages (Gentry et al. 2006b; Larsen et al. 2003, 2016), thus restricting the assessment to the calculated residuals for these wells (Table 4). The absence of modern water in any of these target wells was considered as additional error for a model, adding five points to its score. The scores of models sharing the same spatial configuration for the breach were summed to determine the spatial configuration that would be more likely to be present at the study area based on having smaller residuals and scores.

The parameter estimates (Table 5) were not within the published ranges for both the modern water percentages and apparent ages in any of the target wells (Gentry et al. 2006b; Larsen et al. 2003, 2016). However, in some target wells, one of these parameter estimates was able to match (i.e., fall within) the published ranges (indicated with boldface font in Table 5). Overall, the LPC configuration had the most matches, two with a conductivity of $0.01524\mathrm{m}/\mathrm{day}$ and four with a conductivity of $0.1524\mathrm{m}/\mathrm{day}$. The PC configuration had five matches, while the PB configuration achieved two matches and the remaining configurations tied at one each.

The models with PB and IB configurations were unable to generate pathways for modern water toward the seven target wells under any of the conductivity values tested. Target well MLGW-86R in the southeast section of the well field did not receive modern water in any of the models simulating the PB configuration. This was also the case for target well MLGW-80 in the northern section of the well field, which did not pull water through the breach in any of the models simulating the IB configuration. The LB, PC, and LPC configurations allowed for the generation of modern water pathways toward the seven target production wells under one of the hydraulic conductivity values tested. Due to limited pumping rates for individual wells over time, production wells were assigned equal rates. From discussions with MLGW, not all wells were active month to month; hence, the presence of modern water in the target wells may be impacted by these data limitations.

The mean and mean absolute residuals based on the difference between the simulated water levels and the historic water levels around Sheahan were used as the first criterion of the scoring scheme to compare the model results. The simulated water levels in the Memphis aquifer were consistently higher than the historical water level measurements available for the area. The mean and mean absolute residuals (Tables 6–11) were between 2.95 and 13.75 m, increasing as the extension and the hydraulic conductivity of the breach became greater because it allowed for more downward leakage from the constant-head shallow aquifer. This indicates an excess of water in the Memphis aquifer caused by the shallow aquifer acting as an unlimited source of water, although it was observed to be dry across much of the well field during the early 2000s (Larsen et al. 2003, 2013).

The mean and mean absolute residuals when comparing the mean advective travel times to published tritium/helium-3 apparent ages were significantly impacted by the number of wells that could be analyzed because some of the target wells did not pull water through the breach in many models. With a vertical hydraulic conductivity for the breach of $0.001524\mathrm{m}/\mathrm{day}$ (Tables 6 and 7), the residuals could only be calculated for Well MLGW-87A, which was found to pull water through the breach in the PC and LPC spatial configurations with an average travel time of approximately 43 years. Although Tables 6 and 7 indicate that more than one well was found to have modern water in these models, only the particles for Well MLGW-87A could be used to calculate the apparent age residuals because no other particles placed at the cells containing the target well screens indicated a pathway through the breach. However, this was not the only well found to pull water through the breach under the lower hydraulic conductivity value. Linkage between breaches in all spatial configurations and MLGW-54B were observed. Other linkages were observed between MLGW-63A and IB and LB, MLGW-76A and LB, MLGW-88 and PC, and MLGW-98 and LPC.

Using a hydraulic conductivity of $0.01524\mathrm{m}/\mathrm{day}$ (Tables 8 and 9), the residuals and number of wells for which the analysis was feasible tend to favor the LPC spatial configuration, which extends more widely across the well field, allowing for the presence of modern water in a greater number of wells while also achieving smaller residuals. The IB model was found to have the highest residuals because the enhanced leakage through a narrower area led to the presence of modern water in a smaller number of wells while also reaching these few wells at a quicker pace than what is suggested by the published data. For example, simulating IB and focusing on MLGW-99, particles had an average travel time of 16.4 years, which is, respectively, 23.5 and 35.2 years lower than the minimum and maximum published values (Table 3).

Although it is inferred from the water level residuals that the LPC configuration allows for greater downward leakage into the Memphis aquifer followed by the PC configuration, this was found not to lead to faster travel times into the well field. As observed when using the largest simulated vertical hydraulic conductivity value of $0.1524\mathrm{m}/\mathrm{day}$, the configurations that contribute larger areas of leakage (LB, PB, and IB configurations) had smaller average travel times and less favorable mean absolute residuals (Tables 10 and 11). As seen in the models with a conductivity value of $0.01524\mathrm{m}/\mathrm{day}$, concentrated areas of leakage resulted in faster pathways than those created by leakage features that spread more widely across the well field (i.e., PC and LPC). Though travel times represent advective flow, the impact of probable heterogeneities in the breaches and Memphis aquifer are assumed to be captured by the orders of magnitude variation in hydraulic conductivity from Gentry et al. (2006a). Additionally, calculated tritium/helium-3 ages are expected to be within 10% of the true ages (Larsen et al. 2016).

With the lowest vertical hydraulic conductivity value tested of $0.001524\mathrm{m}/\mathrm{day}$, modern water percentage residuals could only be calculated for the LB, PC, and LPC models (Tables 6 and 7). However, none of these models supported the presence of modern water in the seven target wells. This, in conjunction with the negative sign of the mean residuals, suggests that vertical hydraulic conductivity values for the breach below or equal to $0.001524\mathrm{m}/\mathrm{day}$ are unable to match the extent and magnitude of the presence of modern water around the well field supported by previous studies (Gentry et al. 2006b; Larsen et al. 2003, 2016).

With vertical hydraulic conductivity values of 0.01524 and $0.1524\mathrm{m}/\mathrm{day}$, the mean and mean absolute residuals for the modern water percentage estimates (Tables 8–11) were more favorable for the PB and LPC configurations. However, these residuals were higher than 15%, indicating a consistent overestimation of the calculated modern water percentages. This is attributed to the unrestrained downward leakage of the shallow aquifer due to it being modeled as a constant-head aquifer. A transient simulation of the shallow aquifer is expected to result in smaller residuals because thinning of its saturated thickness would decrease the volume of modern water leaking into the Memphis aquifer. However, due to the limited capabilities of MODPATH and MODFLOW-NWT when dealing with partially or completely dry cells (Villalpando-Vizcaíno 2019) and the limited water table data, this was not possible for this study.

To determine the most likely spatial configuration, composite scores were calculated as the sum of the scores of mean and mean absolute residuals for each modeled breach and for each value of breach hydraulic conductivity. These were then summed into a final score (Table 12). Most of the composite scores of the models decrease as the vertical hydraulic conductivity of the breach increases due to modern water reaching more of the target wells and not losing five points for every target well that did not capture modern water. The exception to this was the LB configuration; its score increased when raising the vertical hydraulic conductivity of the breach from 0.01524 to $0.1524\mathrm{m}/\mathrm{day}$. This was due to both the large breach models with hydraulic conductivity values of 0.01524 and $0.1524\mathrm{m}/\mathrm{day}$ already having modern water at the seven target wells, and the ranks becoming less favorable when increasing the vertical hydraulic conductivity to $0.1524\mathrm{m}/\mathrm{day}$. The LPC model had the most favorable final score, followed by the PC model. As previously mentioned, better results were found for those configurations that are more distributed throughout the well field, which allows for a greater diversity of alternative flow paths for shallow aquifer leakage into the Memphis aquifer.