Assessing Vulnerability of Regional-Scale Aquifer-Aquitard Systems in East Gulf Coastal Plain of Alabama by Developing Groundwater Flow and Transport Models

: Groundwater vulnerability assessment helps subsurface water resources management by providing scientiﬁc information for decision-makers. Rigorous, quantitative assessment of groundwater vulnerability usually requires process-based approaches such as groundwater ﬂow and transport modeling, which have seldom been used for large aquifer-aquitard systems due to limited data and high model uncertainty. To quantify the vulnerability of regional-scale aquifer-aquitard systems in the East Gulf Coastal Plain of Alabama, a three-dimensional (3D) steady-state groundwater ﬂow model was developed using MODFLOW, after applying detailed hydrogeologic information to characterize seven main aquifers bounded by aquitards. The velocity ﬁeld calibrated by observed groundwater heads was then applied to calculate groundwater age and residence time for this 3D aquifer-aquitard system via backward/forward particle tracking. Radioactive isotope data ( 14 C and 36 Cl) were used to calibrate the backward particle tracking model. Results showed that shallow groundwater (<300 ft below the groundwater table) in southern Alabama is mainly the Anthropocene age (25–75 years) and hence susceptible to surface contamination, while the deep aquifer-aquitard systems (700 ft or deeper below the groundwater table) contain “fossil” waters and may be safe from modern contamination if there is no artiﬁcial recharge/discharge. Variable horizontal and vertical vulnerability maps for southern Alabama aquifer-aquitard systems reﬂect hydrologic conditions and intermediate-scale aquifer-aquitard architectures in the regional-scale models. These large-scale ﬂow/transport models with coarse resolutions reasonably characterize the broad distribution and vertical ﬂuctuation of groundwater ages, probably due to aquifer-aquitard structures being captured reliably in the geology model. Parameter sensitivity analysis, vadose zone percolation time, wavelet analysis, and a preliminary extension to transient ﬂow were also discussed to support the aquifer vulnerability assessment indexed by groundwater ages for southern Alabama.


Introduction
Groundwater vulnerability (GWV) to contamination, which refers to the likelihood that an aquifer will become contaminated as a function of pollutant and medium properties, has been a research topic in water resources management for decades [1][2][3][4][5].This is in part because groundwater serves as the primary source of drinking water for 1.5 billion people worldwide [6], and this valuable resource quality is degrading [7].Extensive literature reviews [8,9] have shown that GWV has been assessed by qualitative, statistical, and process-based approaches.Despite the tremendous efforts mentioned above, these methods for assessing aquifer vulnerability are still in their infancy [10].Particularly, the process-based approach is known to provide a more detailed quantification of both sustainability of Alabama's coastal and state-wide groundwater resources management.No GWV has been determined using rigorous, costly, process-based approaches for Alabama.
To reach these goals, this study is organized as follows.Section 2 describes the area of study and process-based numerical modeling approach.Section 3 presents a GWV assessment for aquifer-aquitard systems in southern Alabama using the process-based approach.Section 4 evaluates the GWV assessment using groundwater age isotope data and discusses the groundwater age distribution.Section 5 presents the main conclusions.Appendix A shows the local-scale result for vulnerability analysis (critical for the local coastal aquifer).Appendix B contains a wavelet analysis that correlates precipitation, river stage, and groundwater levels, to explore the response time of subsurface water to surface inputs and supplement the GWV analysis in Section 4. Appendix C lists the vadose zone age and model parameters used in the physical model shown in the main text.

Methodology and Model Development 2.1. Study Site: Geology and Hydrogeology Conditions
The study site is the East Gulf Coastal Plain (EGCP) of Alabama (Figure 1a).The EGCP of Alabama encompasses 40 of Alabama's 67 counties (whose population has grown by 2.6% since 2010) that have landscapes varying from flat areas to rounded/eroded hills.The total groundwater withdrawal in Alabama was 496 million gallons per day in 2015, and more than half of the groundwater withdrawal was for public consumption [36].Groundwater levels have been decreasing for two decades in some areas in Alabama, due to extraction and discharge exceeding recharge [37,38].Land subsidence and groundwater quality deterioration have also been detected in parts of the state, affecting long-term groundwater sustainability [39].
Water 2023, 15, x FOR PEER REVIEW 3 of 25 society, environment, and economy; however, this hidden resource has not received enough attention in Alabama for decades.Large data gaps and incomplete quantifications (of groundwater resources) at almost all spatiotemporal levels are a challenge for the sustainability of Alabama's coastal and state-wide groundwater resources management.No GWV has been determined using rigorous, costly, process-based approaches for Alabama.
To reach these goals, this study is organized as follows.Section 2 describes the area of study and process-based numerical modeling approach.Section 3 presents a GWV assessment for aquifer-aquitard systems in southern Alabama using the process-based approach.Section 4 evaluates the GWV assessment using groundwater age isotope data and discusses the groundwater age distribution.Section 5 presents the main conclusions.Appendix A shows the local-scale result for vulnerability analysis (critical for the local coastal aquifer).Appendix B contains a wavelet analysis that correlates precipitation, river stage, and groundwater levels, to explore the response time of subsurface water to surface inputs and supplement the GWV analysis in Section 4. Appendix C lists the vadose zone age and model parameters used in the physical model shown in the main text.

Study Site: Geology and Hydrogeology Conditions
The study site is the East Gulf Coastal Plain (EGCP) of Alabama (Figure 1a).The EGCP of Alabama encompasses 40 of Alabama's 67 counties (whose population has grown by 2.6% since 2010) that have landscapes varying from flat areas to rounded/eroded hills.The total groundwater withdrawal in Alabama was 496 million gallons per day in 2015, and more than half of the groundwater withdrawal was for public consumption [36].Groundwater levels have been decreasing for two decades in some areas in Alabama, due to extraction and discharge exceeding recharge [37,38].Land subsidence and groundwater quality deterioration have also been detected in parts of the state, affecting long-term groundwater sustainability [39].The Coastal Plain of Alabama is composed mostly of sediment (sand, gravel, and silt) and sedimentary rocks (e.g., chalk, limestone, and sandstone), and these sedimentary deposits range in age from Late Cretaceous to the Middle Eocene [40].Aquifers in southern Alabama and coastal regions can be vulnerable to point/non-point source pollutants due to seawater intrusion, groundwater mineralization, and agricultural activities [41][42][43].The coastal area, such as Baldwin County, Alabama, has a relatively flat topography with low rolling hills, which minimizes surface runoff and enhances infiltration of water and pollutants, for example, soluble nitrate, ability to infiltrate through permeable soils [44].Shallow, unconfined aquifers receive mainly surface recharge (notably, the Alabama climate is temperate, with annual average precipitation of 142 cm evenly distributed throughout the year and concentrated to the south), while deep and layered confined aquifers distributed in EGCP of Alabama may receive recharge from both precipitation (via leakage, with a potential delay) and upgradient inflow (due to the regional groundwater flow direction in coastal aquifers), resulting in complex flow paths presenting challenges for any GWV assessment.Agricultural and animal wastes are significant contaminant sources of surface water and groundwater contamination in this area [42].Pioneering work has been conducted recently to evaluate the shallow groundwater recharge potential using an improved statistical index method "DRASTIC" [45], while rigorous quantification of GWV for deep aquifers requires a process-based model described below.Deep aquifer vulnerability can differ significantly from that for shallow aquifers in aquifer-aquitard complexes because the travel time from water table to well takes can be orders of magnitudes longer than that in the vadose zone [11].

Process-Based Numerical Models to Assess GWV
The process-based modeling approach to assess GWV used in this study has three steps.
Step 1 develops the hydrogeologic aquifer model using the software suite Groundwater Modeling System (GMS) (Figure 1b).Step 2 calculates groundwater heads by developing a steady-state groundwater flow model for the EGCP using the U.S. Geological Survey's modular 3D finite-difference groundwater flow model (MODFLOW) [46].Step 3 evaluates GWV for these aquifers by constructing both forward and backward contaminant transport models using the particle-tracking post-processing model "MODPATH" [47].These three steps are introduced in the following three sub-sections.

Geology Model for Southern Alabama Aquifers
There are seven main aquifers in the EGCP.Their dominant lithology is listed in Table 1.These aquifers mainly consist of fine-medium sand and interlayered clay layers.The irregular 3D spatial distributions of these aquifers are captured by the hydrogeologic model depicted in Figure 1b.Specific characteristics and properties are discussed below: (1) Coker aquifer.It is the lower unit of the Tuscaloosa Group, whose thickness varies from 250 to 500 feet near the fall line and reaches 900 feet downdip [40].It provides a significant source of groundwater in northwest to east-central Alabama [48].(2) Gordo aquifer.The upper part of the Tuscaloosa Group is composed of the Gordo Formation, which is an important source of groundwater from northwest to eastcentral Alabama [49].The Gordo Formation has an average thickness of 300 feet, but it thickens to 500 feet downdip [40].A nonmarine clay layer in the upper section of the Gordo Formation serves as the confining unit above it [40].(3) Eutaw aquifer.Provides a significant amount of water for Alabama.Outcrops range in a thickness from 100 to150 ft in the eastern part of the state and 350 to 400 ft in the western part [49].The Mooreville and Demopolis Chalks form the confining unit above this aquifer.The chalks extend until around Bullock County where they transition to sands and clays of the Blufftown Formation [40].(4) Providence-Ripley aquifer.The Ripley Formation has a thickness varying from 150 to 250 feet, and the Providence Sand exhibits a thickness increasing from less than 50 feet in Lowndes County to around 300 feet at the eastern boundary of Alabama [40].It provides a significant groundwater source in south-central and east-central Alabama.
There are 117 screened wells in the Ripley aquifer, and their depths vary from 18 ft (below land surface, or bls) in the outcrop area to 1045 ft bls downdip.The confining unit above this aquifer is made up of the Prairie Bluff Chalk, Clayton Formation, and Porters Creek Formation.Meanwhile, in the eastern part of the state, a marine clay layer located in the lower section of the Clayton Formation serves as the confining unit above [40].(5) Nanafalia-Clayton aquifer.The thickness of this aquifer varies from 250 ft in southcentral and southwestern Alabama to 75 ft in southeast Alabama [50].To the southwest of the state, the confining unit above consists of the Yazoo Clay.In other regions, the confining unit is made up of silt, clay, and clayey sand located near the middle of the Tuscahoma Formation [40].(6) Lisbon aquifer.Provides the significant public, domestic, agricultural, and industrial water source for Alabama's EGCP [48].The Lisbon Formation is 75~165 ft thick from east to west [50].The confining unit above is located near the middle of the Tuscahoma Formation.It consists of silt, clay, and clayey sand [40].(7) Gulf-coastal lowland aquifer.This aquifer can be found in southern Mobile and Baldwin Counties and consists of clastic sediments in the Miocene undifferentiated, where a complete Miocene section exists stratigraphic interval of the Miocene section is progressively abbreviated farther north due to erosion [49].Where present, the entire Miocene thickness ranges from less than 50 to approximately 2500 ft [40].The confining unit above this aquifer in the southwestern region of the state is the Yazoo Clay.However, in the south-central and southeastern parts, the Yazoo Clay transitions into the Ocala Limestone toward the east.In these regions, the confining unit is a gray clay that is dense and has a soft texture [40].
Most of the aquifers' information mentioned above, such as the spatial location and dimension for each aquifer, are incorporated into the geology model built in this study.

Aquifer Lithology
Coker aquifer Cross-bedded sand, light-colored micaceous, very fine to medium sand, and varicolored micaceous clay [49].Deposited during a time of marine transgression, the Coker Formation was deposited [51].

Gordo aquifer
Cross-bedded sand, gravelly sand, and lenticular beds of locally carbonaceous clay that are partially mottled moderate-red and pale-red-purple.Pale-yellowish-orange, poorly sorted, cross-bedded gravelly fine to very coarse quartz sand, containing irregular beds of moderate-reddish-brown to pale-red-purple sandy clay [47].The boundary between marine sediments of the Coker Formation (which consists of massive, marine clay with thin beds of fine-grained sand) and nonmarine sediments of the Gordo Formation formed during a period of significant sea level regression [51].

Eutaw aquifer
The western part is described as light-greenish-gray well-sorted micaceous cross-bedded fine to medium sand.The eastern part is described as light-greenish-gray to yellowish-gray, well-sorted, micaceous, partly fossiliferous fine to medium quartz sand interbedded with dark-gray carbonaceous clay, greenish-gray micaceous sandy clay, and thin beds of glauconitic, fossiliferous sandstone [49].
The Eutaw Formation was primarily formed in a marginal marine setting connected to a barrier island and deltaic environment [51].

Providence-Ripley aquifer
The eastern part consists of light gray to pale-olive massive, micaceous, glauconitic, fossiliferous fine sand, sandy calcareous clay, and thin indurated beds of fossiliferous sandstone.The western part of this formation contains micaceous fine to medium quartz sand, cross-bedded in the upper part, and sandy calcareous clay [49].Gray, fossiliferous, silty Demopolis chalk in central and western Alabama.
The chalk overlies the Mooreville Chalk and grades into the Blufftown and Ripley Formations in the east of the region [51].

Nanafalia-Clayton aquifer
Clayton Formation in eastern Alabama comprises fine sand, medium-gray silty, calcareous clay, sandy fossiliferous limestone, and gravelly, medium to coarse sand containing clay pebbles.Glauconitic sand, massive clay, and fossiliferous sandy clay form the Nanafalia Formation [49,50] The aquifer includes the basal sand of Tuscahoma Formation, and the whole of the Nanafalia and equivalent Baker Hill (in eastern Alabama), Naheola, Porters Creek, and Clayton Formations.However, one or more of these formations is absent at any one geographical location.The aquifer consists mostly of unconsolidated sand and clay beds, but locally includes carbonate rocks [52].

Lisbon aquifer
Sand, limestone, and sandy limestone, highly fossiliferous, glauconitic, quartz sand, and lenses of greenish-gray clay [53].According to Toulmin and LaMoreaux (1963) [54], the Lisbon Formation in southeast Alabama consists primarily of sand but also contains significant amounts of limestone and sandy limestone.The Gosport Sand is only mapped in the west and central Alabama, between the Alabama River and the Alabama-Mississippi state line, and is comprises highly fossiliferous, glauconitic, quartz sand and lenses of greenish-gray clay [53], with an outcrop thickness ranging from 17 to 30 feet [50].

Groundwater Flow Model Built by MODFLOW
The hydrogeologic model for south Alabama developed in the MODFLOW package with built-in GMS contains detailed hydrogeologic and geologic information, including six cross-sections (after stratigraphic analysis) that used aquifer-formation data from 48 boreholes across the state.The borehole data were collected by Davis [40] and the Geological Survey of Alabama (GSA) [48] (Figure 2).Davis [40] used electric logs from water wells and oil test holes to describe the Coastal Plain's overall structure and geology.The spontaneous potential (SP) and resistivity curves from well logs have been widely employed to confirm the similarity of curve shapes on logs from other wells.Drillers' and sample logs were used in locations where electric logs were not available.The hydrogeologic model developed for this study consists of seven aquifers listed in Table 1, with a grid array of 260 × 260 × 15 along x/y/x directions, and with a total grid number of ~1 million.
Boundary conditions of the steady-state groundwater flow model are defined as follows.The top layer is an active recharge boundary.The bottom layer is bedrock [40], and a no-flow boundary is defined for the entire layer at the bottom.The other boundaries are constant head boundaries and variable head boundaries, conveniently defined by cell types in MODFLOW.The variable heads are defined with the reference heads obtained from 30 observation wells in the unconfined aquifers using November 2020 for the steady-state flow model.The constant heads are defined using the reference river stages, with data from November 2020.The pumping rates were obtained from the GSA November 2020 data.
Thirteen rivers and their riverbed conductance were also incorporated into the groundwater model [48,55].From Clark and Hart [55], the final streambed conductance values were calculated from PEST.Faults were incorporated using the Horizontal Flow Barrier (HFB) package, either by simulating enhanced flow or by acting as a flow barrier.The recharge potential-intrinsic (RPI) layer in the 30-year range (1989-2019) from Guthrie et al. [45] was used to define the recharge area, and the 1981-2010 annual average precipitation (AAP) from the U.S. Department of Agriculture (USDA, https://datagateway.nrcs.usda.gov/,accessed on 1 September 2022), the precipitation rate from November 2020 (NOAA, https://hdsc.nws.noaa.gov/hdsc/pfds/,accessed on 1 September 2022), and GSA data [48] were used to calculate the potential recharge rate.Hydraulic conductivity for each sedimentary rock is approximated initially using representative values from Domenico and Schwartz [56], GSA [48], Faye and Smith [57], Martin and Whiteman [58], Sun and Johnston [59], and Mallory [60].Further parameter calibration in GMS is shown in Section 3.1.
Notably, the accurate modeling of large aquifer systems is very challenging given the paucity of data.Hence, it is impossible to obtain detailed hydrologic dynamics by such models, except for the overall pattern of groundwater flow and age distribution focused by this study.An enhanced sub-area model can help to identify information missed by the large-scale coarse resolution model, which will be the focus of the next study.Thirteen rivers and their riverbed conductance were also incorporated into the groundwater model [48,55].From Clark and Hart [55], the final streambed conductance values were calculated from PEST.Faults were incorporated using the Horizontal Flow Barrier (HFB) package, either by simulating enhanced flow or by acting as a flow barrier.The recharge potential-intrinsic (RPI) layer in the 30-year range (1989-2019) from Guthrie et al. [45] was used to define the recharge area, and the 1981-2010 annual average precipitation (AAP) from the U.S. Department of Agriculture (USDA, https://datagateway.nrcs.usda.gov/,accessed on 1 September 2022), the precipitation rate from November 2020 (NOAA, https://hdsc.nws.noaa.gov/hdsc/pfds/,accessed on 1 September 2022), and GSA data [48] were used to calculate the potential recharge rate.Hydraulic conductivity for each sedimentary rock is approximated initially using representative values from Domenico and Schwartz [56], GSA [48], Faye and Smith [57], Martin and Whiteman [58], Sun and Johnston [59], and Mallory [60].Further parameter calibration in GMS is shown in Section 3.1.Notably, the accurate modeling of large aquifer systems is very challenging given the paucity of data.Hence, it is impossible to obtain detailed hydrologic dynamics by such models, except for the overall pattern of groundwater flow and age distribution focused by this study.An enhanced sub-area model can help to identify information missed by the large-scale coarse resolution model, which will be the focus of the next study.Note that the stratigraphic columns using lithologic correlation (with also uniformities) shown in (b,c) end at the bottom of each borehole, which does not mean that the aquifer terminates abruptly (GMS matches the rock types below boreholes, which are not shown here since the plots (b,c) only illustrate the cross-sections derived by boreholes).

Contaminant Transport Model Using MODPATH
The velocity calculated by MODFLOW is then used by MODPATH to track 3D advective trajectories either forward or backward [47].In MODPATH, the calculated travel path and history of random-walking particles that represent water parcels in the EGCP aquifers provide critical information for GWV assessment.On the one hand, the forward particle tracking in MODPATH simulates a water parcel's "life expectancy" and residence time in the aquifer by moving particles along streamlines.It tracks their forward position from a specified location in the aquifer, such as the well screen, before they are sampled by one of the outflows such as pumping or groundwater discharge to streams or oceans [61].The backward particle tracking calculates the "age" of the water parcel, which represents the time elapsed since entering the aquifer and can be used to characterize the aquifer's vulnerability to non-point source contamination from land surface or water table.The sum of groundwater "life expectancy" and "age" is the water parcel's total travel time [61], which defines the aquifer renewal time (note it does not account for the impact of future climate change or the change of pumping on transport if these factors are not modeled).Therefore, the forward and backward particle tracking schemes embedded in MODPATH can lead to useful hydrogeologic information for assessing GWV and the aquifer renewal time framework.

Model Calibration
Groundwater heads for all 41 monitoring wells (in either unconfined or confined aquifers) observed in November 2020 were fitted by the steady-state flow model.Three parameters-hydraulic conductivity, recharge rates, and riverbed conductance-were selected as calibration parameters, since there were no measurements in this area for these parameters, and they may significantly affect groundwater head distributions.To achieve the best calibration, we combined a trial-and-error approach (for a preliminary calibration) and then automated parameter estimation using Pilot PEST.The simulated groundwater heads are similar to the measured ones (Figure 3), where the root mean square error (RMSE) for all 41 monitoring wells is 9.38 ft.The model error is unbiased and is relatively low compared to the overall range of the observed groundwater heads (~200 ft).
ble.The sum of groundwater "life expectancy" and "age" is the water parcel's total travel time [61], which defines the aquifer renewal time (note it does not account for the impact of future climate change or the change of pumping on transport if these factors are not modeled).Therefore, the forward and backward particle tracking schemes embedded in MODPATH can lead to useful hydrogeologic information for assessing GWV and the aquifer renewal time framework.

Model Calibration
Groundwater heads for all 41 monitoring wells (in either unconfined or confined aquifers) observed in November 2020 were fitted by the steady-state flow model.Three parameters-hydraulic conductivity, recharge rates, and riverbed conductance-were selected as calibration parameters, since there were no measurements in this area for these parameters, and they may significantly affect groundwater head distributions.To achieve the best calibration, we combined a trial-and-error approach (for a preliminary calibration) and then automated parameter estimation using Pilot PEST.The simulated groundwater heads are similar to the measured ones (Figure 3), where the root mean square error (RMSE) for all 41 monitoring wells is 9.38 ft.The model error is unbiased and is relatively low compared to the overall range of the observed groundwater heads (~200 ft).The simulated groundwater table contour map (plotted in Figure 4b) captures the overall pattern, including the locations of high and low heads, of the groundwater table contour map interpolated by the 30 observation wells (whose well screens are located in the unconfined aquifer), using the ordinary kriging method (shown in Figure 4a), although these two contour maps exhibit a subtle difference in the local groundwater head distribution.Notably, the interpolated contour shown in Figure 4a may be too smooth to The simulated groundwater table contour map (plotted in Figure 4b) captures the overall pattern, including the locations of high and low heads, of the groundwater table contour map interpolated by the 30 observation wells (whose well screens are located in the unconfined aquifer), using the ordinary kriging method (shown in Figure 4a), although these two contour maps exhibit a subtle difference in the local groundwater head distribution.Notably, the interpolated contour shown in Figure 4a may be too smooth to capture the local variation of the real groundwater head distribution, due to the relatively small number of observation points in a large aquifer/aquitard system.

Parameter Sensitive Analysis
PEST calibrates 158 model parameters (in three types, as mentioned above) that may affect groundwater flow magnitude and/or direction (Figure 5).The PEST simulation also provides the parameter sensitivity calculation.Sensitivity analyses showed that the recharge rate is the most sensitive parameter, followed by horizontal hydraulic conductivity and riverbed conductance, respectively (Figure 5).PEST keeps track of each parameter's composite and relative composite sensitivity.At the end of parameter calibration, PEST performs sensitivity analysis [62].To calculate parameter sensitivity, Hill and Tiedeman [63] utilized the following equation to obtain the (dimensionless) composite scaled sensitivity (CSS) of parameter i: where ND is the total number of observations, y i is a simulated value for the i-th observation, b j is the j-th estimated parameter, and w 1/2 is the square root of the weight matrix determined.This data analysis helps to determine which parameters have the greatest impact on the model output and which ones have the least impact.
Water 2023, 15, x FOR PEER REVIEW 9 of 25 capture the local variation of the real groundwater head distribution, due to the relatively small number of observation points in a large aquifer/aquitard system.

Parameter Sensitive Analysis
PEST calibrates 158 model parameters (in three types, as mentioned above) that may affect groundwater flow magnitude and/or direction (Figure 5).The PEST simulation also provides the parameter sensitivity calculation.Sensitivity analyses showed that the recharge rate is the most sensitive parameter, followed by horizontal hydraulic conductivity and riverbed conductance, respectively (Figure 5).PEST keeps track of each parameter's composite and relative composite sensitivity.At the end of parameter calibration, PEST performs sensitivity analysis [62].To calculate parameter sensitivity, Hill and Tiedeman [63] utilized the following equation to obtain the (dimensionless) composite scaled sensitivity (CSS) of parameter i: where ND is the total number of observations,  is a simulated value for the i-th observation,  is the j-th estimated parameter, and  / is the square root of the weight matrix determined.This data analysis helps to determine which parameters have the greatest impact on the model output and which ones have the least impact.A3).A3).The groundwater age sampled from 24 wells in Coastal Plain Alabama ranged from less than 100 years to 80,000 years [64].Penny et al. [65] used 36 Cl/Cl ratios to calculate the groundwater age differences and flow velocities in the Eutaw and Tuscaloosa aquifers in EGCP.The 36 Cl/Cl ratio method is particularly effective for dating groundwater that is significantly older than that which can be dated with 14 C [65]. 14 C possesses a half-life of approximately 5730 years, whereas 36 Cl has a half-life of nearly 301,000 years.This means that 14 C is more useful for dating relatively young fossil samples (up to around 50,000 years), while 36 Cl can be more reliable for dating groundwater whose age is in the range of 60,000 to 1 million years.

Groundwater Age and
Groundwater in Moundville (well E1, marked in Figure 6a) and Greensboro (well E4, marked in Figure 6a) have a 36 Cl age difference of approximately 110,000 years.Groundwater in South Macon's well T1 and Troy's well T8 has a 36 Cl age difference of almost 459,000 years (Figure 6a).Penny et al. [65] found that the groundwater age in the confined aquifers of the Coastal Plain was significantly different from the groundwater age in the unconfined aquifers.The large age differences observed between the confined and unconfined aquifers and the distance from the recharge area to the sampling location are the result of the different recharge rates and flow processes that occur in each aquifer.In addition, significant mixing of groundwater can reduce the 36 Cl/Cl ratio, which may result in an underestimation of the groundwater recharge rates when interpreting groundwater ages based on field data.
Water 2023, 15, x FOR PEER REVIEW 11 of 25 50,000 years), while 36 Cl can be more reliable for dating groundwater whose age is in the range of 60,000 to 1 million years.
Groundwater in Moundville (well E1, marked in Figure 6a) and Greensboro (well E4, marked in Figure 6a) have a 36 Cl age difference of approximately 110,000 years.Groundwater in South Macon's well T1 and Troy's well T8 has a 36 Cl age difference of almost 459,000 years (Figure 6a).Penny et al. [65] found that the groundwater age in the confined aquifers of the Coastal Plain was significantly different from the groundwater age in the unconfined aquifers.The large age differences observed between the confined and unconfined aquifers and the distance from the recharge area to the sampling location are the result of the different recharge rates and flow processes that occur in each aquifer.In addition, significant mixing of groundwater can reduce the 36 Cl/Cl ratio, which may result in an underestimation of the groundwater recharge rates when interpreting groundwater ages based on field data.Figure 6.These images show 24 wells of 14 C samples from Solder [60] (circles) and 4 wells of 36 Cl samples from Penny and Lee [61] (squares) (a).Comparison between the isotopic-dated ages (using 14 C and 36 Cl) and the backward particle tracking ages using MODPATH in a log-log plot (where the unit of age is year) (b).
Figure 6b illustrates that the backward particle tracking ages calculated using the calibrated MODPATH model generally match the measured 14 C and 36 Cl ages for groundwater.It is noteworthy that the particle tracking method by MODPATH calculates the advective time only (without dispersion), implying that regional-scale advection may play a more important role than well bore mixing and local dispersion in defining the mean groundwater age in a state-wide flow/transport model.Figure 7 shows the modeled Figure 6.These images show 24 wells of 14 C samples from Solder [60] (circles) and 4 wells of 36 Cl samples from Penny and Lee [61] (squares) (a).Comparison between the isotopic-dated ages (using 14 C and 36 Cl) and the backward particle tracking ages using MODPATH in a log-log plot (where the unit of age is year) (b).
Figure 6b illustrates that the backward particle tracking ages calculated using the calibrated MODPATH model generally match the measured 14 C and 36 Cl ages for groundwater.It is noteworthy that the particle tracking method by MODPATH calculates the advective time only (without dispersion), implying that regional-scale advection may play a more important role than well bore mixing and local dispersion in defining the mean groundwater age in a state-wide flow/transport model.Figure 7 shows the modeled groundwater age using MODPATH for the regional scale aquifers-aquitards at different depths.This hypothesis needs to be validated by further tests.It is also noteworthy that the observed hydraulic heads were used to calibrate the groundwater flow model in Section 3.1, and the calibrated flow model is used here to run particle tracking and fit the isotope ages independently.This second stage only needs to calibrate the effective porosity, since it is the only new parameter required for particle tracking in MODPATH.In addition, the shallow aquifer is more vulnerable and remains the major concern for water usage, while the deep aquifers are added here to explore the 3D fluctuation of groundwater vulnerability which can be affected by shallow-deep aquifer mass exchange.(where the unit of age is year).

Forward Particle Tracking Calculates Residence Time
MODPATH, with the parameter effective porosity calibrated by the isotopic-dated ages, was then applied to predict "life expectancy" or "residence time" of groundwater using forward-in-time particle tracking.Here each particle was tracked forward along streamlines, from the monitoring well to its outlet (i.e., the drain sink cell) (assuming that the average flow field does not change significantly in the future; otherwise, a multi-million-year-long transient and predictive flow model is needed, which is however not easy, if not impossible, to build for many sites).Figure 8 shows the forward particle tracking results for water located at three different starting depths in the aquifer, which are 300, 700, and 1000 ft below the groundwater table, respectively.In this study, the residence time (calculated by forward particle tracking) is variable along the vertical direction with the times ranging from younger than 100 years to older than 90,000 years for aquifersaquitards at the depth of 300, 700, and 1000 ft below the groundwater table.

Forward Particle Tracking Calculates Residence Time
MODPATH, with the parameter effective porosity calibrated by the isotopic-dated ages, was then applied to predict "life expectancy" or "residence time" of groundwater using forward-in-time particle tracking.Here each particle was tracked forward along streamlines, from the monitoring well to its outlet (i.e., the drain sink cell) (assuming that the average flow field does not change significantly in the future; otherwise, a multi-millionyear-long transient and predictive flow model is needed, which is however not easy, if not impossible, to build for many sites).Figure 8 shows the forward particle tracking results for water located at three different starting depths in the aquifer, which are 300, 700, and 1000 ft below the groundwater table, respectively.In this study, the residence time (calculated by forward particle tracking) is variable along the vertical direction with the times ranging from younger than 100 years to older than 90,000 years for aquifers-aquitards at the depth of 300, 700, and 1000 ft below the groundwater table.
lion-year-long transient and predictive flow model is needed, which is however not easy, if not impossible, to build for many sites).Figure 8 shows the forward particle tracking results for water located at three different starting depths in the aquifer, which are 300, 700, and 1000 ft below the groundwater table, respectively.In this study, the residence time (calculated by forward particle tracking) is variable along the vertical direction with the times ranging from younger than 100 years to older than 90,000 years for aquifersaquitards at the depth of 300, 700, and 1000 ft below the groundwater table.

Total Travel Time
The total travel time for southern Alabama groundwater calculated by MODPATH is interpreted as the point-source age of groundwater arriving at a well screen plus the subsequent time taken by the water parcel to exit the EGCP.Groundwater's total ages typically increase with depth and range from less than 100 to more than 100,000 years in the EGCP aquifers (Figure 9).The vadose zone ages (listed in Figure A5, Appendix C), which are usually orders of magnitudes smaller than the underlying groundwater age, can be neglected in this long, total travel time.The total travel time for southern Alabama groundwater calculated by MODPATH is interpreted as the point-source age of groundwater arriving at a well screen plus the subsequent time taken by the water parcel to exit the EGCP.Groundwater's total ages typically increase with depth and range from less than 100 to more than 100,000 years in the EGCP aquifers (Figure 9).The vadose zone ages (listed in Figure A5, Appendix C), which are usually orders of magnitudes smaller than the underlying groundwater age, can be neglected in this long, total travel time.

Vulnerability of Southern Alabama Aquifers with Broad Groundwater Ages
Mixed young and old groundwaters are found in southern Alabama aquifers.For example, the migration of water from the recharge areas in shallow unconfined aquifers downward through vertical leakage through the adjacent aquitards can lead to relatively young groundwater (according to the backward particle tracking model), whose spatial distribution can extend into the lower, deeper aquifers along preferential flow paths consisting of high-permeable, interconnected sediment or sedimentary rocks.Meanwhile, the 3D groundwater flow model reveals that groundwater flows upward in the northern part of the Eutaw and Gordo aquifer (northern part of Figure 1a), which is consistent with the

Vulnerability of Southern Alabama Aquifers with Broad Groundwater Ages
Mixed young and old groundwaters are found in southern Alabama aquifers.For example, the migration of water from the recharge areas in shallow unconfined aquifers downward through vertical leakage through the adjacent aquitards can lead to relatively young groundwater (according to the backward particle tracking model), whose spatial distribution can extend into the lower, deeper aquifers along preferential flow paths consisting of high-permeable, interconnected sediment or sedimentary rocks.Meanwhile, the 3D groundwater flow model reveals that groundwater flows upward in the northern part of the Eutaw and Gordo aquifer (northern part of Figure 1a), which is consistent with the finding by Gardner [66].This upward flow is probably due to the upward vertical leakage through confining beds to rivers [66].In addition, groundwater from the Nanafalia-Clayton aquifer is typically of high quality and appropriate for a wide range of uses.However, in Marengo and western Wilcox Counties, there are elevated levels of chloride, bicarbonate, and dissolved solids, which may be the result of groundwater moving upward through a fault from underlying aquifers [67,68].Therefore, in the northern region of EGCP, the shallow aquifer (i.e., the Eutaw and Gordo aquifer) has relatively older groundwater than the shallow aquifer located in southern EGCP.In addition, the average groundwater age increases substantially with increasing depth, especially when low-permeability clay layers separate shallow and deep aquifers and significantly retard the vertical movement of water and pollutants.Hence, the hydrologic conditions and intermediate-scale aquifer-aquitard architecture in the regional-scale model likely contribute to the mixed ages for EGCP groundwaters, resulting in a horizontally non-uniform and vertically variable vulnerability map for southern Alabama aquifers (shown by Figure 7).
Most groundwater above 300 ft (depth below the groundwater table) in southern Alabama is younger than 100 years (Figure 7a), which is Anthropocene-age and therefore, susceptible to surface contamination [69].Most deep aquifers (≥700 ft below the groundwater table in southern Alabama contain groundwater older than 20,000 years, representing "fossil" aquifers in the Pleistocene epoch and should be "safe" from modern contamination [70].The existence of both modern and "fossil" aquifers in southern Alabama is generally consistent with the reported isotopic ages of groundwater.For example, the estimated mean groundwater age by Solder [64] using 14 C for the South Atlantic and Gulf Coast is ~30,000 years, indicating "old" groundwater in the aquifer system (i.e., "fossil" water).Relatively young groundwater with a mean age of less than 2000 years was typically found in some unconfined parts of these aquifer systems [64].In this study area (i.e., along the Southeastern Coastal Plain and the Coastal Lowlands), there were a total of 24 samples analyzed for 14 C with ages ranging from less than 100 years to more than 80,000 years (Figure 6b).The results indicate that the age of groundwater in the study area varies in an aquifer with sampling depth.Particularly, we found that the shallow wells located in the uppermost layer (i.e., Layer 1) tend to contain younger groundwater compared to deep wells (located in deeper layers).
To further interpret the groundwater age distribution, the particle tracking ages calculated in Section 3 were transformed into probability density functions (PDFs).The physical heterogeneity of regional-scale aquifers is usually characterized by highly variable flow velocities and multiscale coherence lengths [71].Here the groundwater age is defined as the amount of time since the water parcel entered the aquifer (i.e., ignoring the delay of the vadose zone).Figures 10 and 11 plot the PDF distributions of groundwater ages employing two parameters, effective porosity n, and vertical hydraulic conductivity K v .
These PDFs behave as base-10 log normal distributions with slightly elongated latetime tails (implying the tempered stable density for the age distribution identified for natural heterogeneous aquifers [72,73]) after adjusting manually the values of effective porosity and vertical hydraulic conductivity.In terms of porosity, the PDF of groundwater ages shifts to the right (representing more old groundwater) for the case of 0.01 ≤ n ≤ 0.18 compared to that for 0.1 ≤ n ≤ 0.46, probably due to the change of flow paths (driven by the change of local velocities).In addition, the resultant PDF of groundwater ages for the case of 1 × 10 −6 ≤ K v ≤ 1 × 10 1 ft/day also shifts to the right (with a higher peak for deep aquifers) than that for 1 × 10 −4 ≤ K v ≤ 1 × 10 2 ft/day, because of the decreasing vertical velocity when K v is smaller.
To further interpret the groundwater age distribution, the particle tracking ages calculated in Section 3 were transformed into probability density functions (PDFs).The physical heterogeneity of regional-scale aquifers is usually characterized by highly variable flow velocities and multiscale coherence lengths [71].Here the groundwater age is defined as the amount of time since the water parcel entered the aquifer (i.e., ignoring the delay of the vadose zone).Figures 10 and 11 plot the PDF distributions of groundwater ages employing two parameters, effective porosity , and vertical hydraulic conductivity  .These PDFs behave as base-10 log normal distributions with slightly elongated latetime tails (implying the tempered stable density for the age distribution identified for natural heterogeneous aquifers [72,73]) after adjusting manually the values of effective porosity and vertical hydraulic conductivity.In terms of porosity, the PDF of groundwater ages shifts to the right (representing more old groundwater) for the case of 0.01 ≤  ≤ 0.18 compared to that for 0.1 ≤  ≤ 0.46, probably due to the change of flow paths (driven by the change of local velocities).In addition, the resultant PDF of groundwater ages for the case of 1 10  1 10 ft/day also shifts to the right (with a higher peak for deep aquifers) than that for 1 10  1 10 ft/day, because of the decreasing vertical velocity when  is smaller.

Extension to Transient Flow: Impact on GWV
One major limitation of this study is that the steady-state flow model was used to explore groundwater ages and vulnerability (since there were no observations to support a thousand-year-long transient model).An accurate assessment of the impacts of climate change on groundwater vulnerability and depletion would require a transient groundwater model.As a preliminary test, we developed and calibrated a one-year-long transient model.Results showed that this preliminary transient model captures the overall temporal pattern of groundwater levels, although further calibrations are needed to improve the model fit (Figure 12a).

Extension to Transient Flow: Impact on GWV
One major limitation of this study is that the steady-state flow model was used to explore groundwater ages and vulnerability (since there were no observations to support a thousand-year-long transient model).An accurate assessment of the impacts of climate change on groundwater vulnerability and depletion would require a transient groundwater model.As a preliminary test, we developed and calibrated a one-year-long transient model.Results showed that this preliminary transient model captures the overall temporal pattern of groundwater levels, although further calibrations are needed to improve the model fit (Figure 12a).
To further explore the impact of climate change and anthropogenic activities on southern Alabama groundwater resources, we extended the transient time frame to 10 years (Figure 12b) by changing recharge rates and pumping rates in the three scenarios listed below.
explore groundwater ages and vulnerability (since there were no observations to support a thousand-year-long transient model).An accurate assessment of the impacts of climate change on groundwater vulnerability and depletion would require a transient groundwater model.As a preliminary test, we developed and calibrated a one-year-long transient model.Results showed that this preliminary transient model captures the overall temporal pattern of groundwater levels, although further calibrations are needed to improve the model fit (Figure 12a).To further explore the impact of climate change and anthropogenic activities on southern Alabama groundwater resources, we extended the transient time frame to 10 years (Figure 12b) by changing recharge rates and pumping rates in the three scenarios listed below.Scenario I: Base case following historical recharge and pumping rates.Harper et al. [36] found that the total groundwater withdrawal in Alabama was 496 million gallons per day (MGD) in 2015, more than half of which was used for public consumption.Scenario I assumed that the groundwater recharge and withdrawal rates would remain constant for the next decade, and there would be 78 pumping wells (the same number of wells identified for the steady-state model) with a pumping rate of 0.85 Mft 3 /day/well.The corresponding evolution of the groundwater head was modeled until 2030.Results showed that the distance of the wells from the boundary condition (i.e., the constant head boundary) had a significant impact on the effect of pumping rates on the groundwater system.On the one hand, wells located far from the boundary in Layer 1 exhibited a mixed pattern of drawdown fluctuation, where the pumping rate dominates the recharge rate, leading to a decline in groundwater levels over time.On the other hand, wells located closer to the boundary and wells in Layer 2 were not significantly affected by the change in pumping rates (Figure 13), as expected.Scenario I: Base case following historical recharge and pumping rates.Harper et al. [36] found that the total groundwater withdrawal in Alabama was 496 million gallons per day (MGD) in 2015, more than half of which was used for public consumption.Scenario I assumed that the groundwater recharge and withdrawal rates would remain constant for the next decade, and there would be 78 pumping wells (the same number of wells identified for the steady-state model) with a pumping rate of 0.85 Mft 3 /day/well.The corresponding evolution of the groundwater head was modeled until 2030.Results showed that the distance of the wells from the boundary condition (i.e., the constant head boundary) had a significant impact on the effect of pumping rates on the groundwater system.On the one hand, wells located far from the boundary in Layer 1 exhibited a mixed pattern of drawdown fluctuation, where the pumping rate dominates the recharge rate, leading to a decline in groundwater levels over time.On the other hand, wells located closer to the boundary and wells in Layer 2 were not significantly affected by the change in pumping rates (Figure 13), as expected.Scenario II: Increased pumping rate for wetting years.The recharge rate to the aquifer system was kept unchanged, while the pumping rate in the aquifer system was increased from 0.85 to 66.3 Mft 3 /day, with a 7700% increase (representing an extreme increase in freshwater demand).The continuous decline in the simulated water levels for most wells Scenario II: Increased pumping rate for wetting years.The recharge rate to the aquifer system was kept unchanged, while the pumping rate in the aquifer system was increased from 0.85 to 66.3 Mft 3 /day, with a 7700% increase (representing an extreme increase in freshwater demand).The continuous decline in the simulated water levels for most wells in Layer 1 implies the depletion of Alabama groundwater resources under an increased pumping rate under the same climate.Wells in Layer 2, however, were relatively unaffected by the increasing pumping rate (their head remained relatively stable) (Figure 13), since Alabama deep aquifers are less sensitive to the change in the pumping rate than shallow aquifers.
Scenario III: Decreasing recharge rates for dry years.This scenario explored how groundwater in the ECGP responds to a drier climate in the next decade.The recharge rate was reduced by 10% from that of the steady-state flow condition, while the pumping rate remained constant at 0.85 Mft 3 /day/well until 2030.The whole study area experienced an overall decline in groundwater levels, which may be attributed to the decreasing recharge rates under a drier climate (Figure 13).This suggests that groundwater resources in the study area may be at risk of depletion if the current recharge rates continue to decline by 10%.
It is noteworthy that future work is needed to improve this preliminary transient flow model.For example, the model grid can be refined to increase the model resolution and the cell's corresponding parameters.The boundary conditions can be further adjusted to account for various changing environmental factors.Limitations of the transient flow model shown above also need to be accounted for when refining the groundwater models.For instance, accurate pumping rates are crucial for accurately forecasting groundwater levels; hence, efforts are needed to obtain as precise and comprehensive pumping data as possible.Furthermore, in the southern part of Alabama, seawater intrusion may affect local groundwater dynamics, which needs to be considered when predicting future groundwater quality in the coastal aquifers.These efforts will be pursued in our future studies.

Conclusions
This study developed the first process-based approach to assess the vulnerability of aquifer-aquitard systems in southern Alabama by building regional-scale flow and transport models using MODFLOW and MODPATH.Model construction, parameter sensitivity analysis, and statistical analysis revealed the following four main conclusions.
First, the process-based modeling approach can be applied to assess the vulnerability of state-wide groundwater systems, where the large-scale groundwater flow model was calibrated using observed heads and the subsequent transport model was calibrated using the measured radioactive isotope ages.The large-scale model with the coarse grid structure, where the immediate-scale aquifer-aquitard structures were well captured (i.e., 7 main aquifers separated by aquitards), can characterize the broad distribution and strong vertical fluctuation of groundwater ages, due for example to preferential flow paths.
Second, sensitivity analyses showed that the recharge rate, horizontal hydraulic conductivity, and riverbed conductance are the three dominant parameters for affecting groundwater ages and therefore aquifer vulnerability.In addition, based on the result from the probability density function, vertical hydraulic conductivity and effective porosity are the main factors controlling the groundwater age distributions.
Third, distributions of the total ages of EGCP groundwater followed the same spatial pattern as the groundwater age calculated by backward particle tracking, because the EGCP is predominantly composed of unconsolidated to semi-consolidated sediments, and the groundwater age can also reflect subsurface heterogeneity.This similarity implies that the age of groundwater plays a significant role in total travel time and that groundwater moves relatively faster to the discharge area (than that moving backward to its recharge area).For example, in Baldwin County and Mobile County, southern Alabama aquifers are relatively closer (than the northern areas) to the discharge zones and exhibit shorter residence times.
Fourth, transient groundwater flow models emphasized the importance of the well location and its distance to the boundary in designing groundwater management strategies.
For example, Scenario II (with more pumping) showed the importance of maintaining a sustainable balance between recharge and pumping.The increased pumping rates resulted in an overall decline in the groundwater levels, but the wells in Layer 2 were less sensitive to changes in pumping than shallower wells.Scenario III emphasized the potential risk of groundwater depletion due to the decreasing recharge rates under a drier climate in the next decade.Sustainable groundwater requires reliable management of aquifers vulnerable to changes in climate and other environmental factors.

Appendix B. Wavelet Analysis
To further check the groundwater age, we conducted wavelet analysis to link groundwater to precipitation and surface water.The results show that precipitation has a significant impact on groundwater head (for a well located 140 ft bls) for a period longer than a few years.The river stage can reduce groundwater level fluctuations across a timescale of 30 to 365 days.However, precipitation has a more significant influence than the river stage in controlling the groundwater level change.Therefore, the shallow groundwater response to surface signals in southern Alabama may be delayed by several years after precipitation input.This is consistent with the young groundwater (1~5 years shown in Figure 13) modeled above.
This appendix focuses on a monitoring well, MON-1, located in southern Alabama's Monroe County.Specifically, MON-1 is positioned north of Frisco City and has a depth of 140 feet.Water levels at this well have been continuously measured since 2011.Additionally, MON-1 is close to agricultural land and the Alabama River.
The groundwater level time series data in Monroe County show weak periodicity, with the periodic features focusing mostly on 4-16 days in winter each year, according to the results of wavelet analysis (Figure A2).The time frame of the river phases is primarily 8-128 days.In the springs of 2016, 2019, and 2020, the periodic behavior is more obvious.Precipitation indicates a continuous period of one year.
The percentage area of significant coherence (PASC) and average wavelet coherence (AWC) are used to evaluate the quantitative relationship between the predictor factors and the response variable AWC.The variable with higher AWC and PASC can account for a larger variety of response variable fluctuations.The bivariate wavelet coherence is presented in Figure A2.According to Figure A2a, neither the river stage nor the groundwater level shows greater periodicity.However, the results identified a region with frequent fluctuation in groundwater level and river stage during a 128-day period (Figure A3b).The AWC and PASC at various temporal scales between groundwater levels and river stages are shown in Table A1.At periods between 180 and 365 days, the river stage makes the most significant contribution to the groundwater level fluctuation.and the response variable AWC.The variable with higher AWC and PASC can account for a larger variety of response variable fluctuations.The bivariate wavelet coherence is presented in Figure A2.According to Figure A2a, neither the river stage nor the groundwater level shows greater periodicity.However, the results identified a region with frequent fluctuation in groundwater level and river stage during a 128-day period (Figure A3b).The AWC and PASC at various temporal scales between groundwater levels and river stages are shown in Table A1.At periods between 180 and 365 days, the river stage makes the most significant contribution to the groundwater level fluctuation.According to Figure A4, the greatest periodicity for the groundwater level and precipitation is between 365-768 days.The groundwater level illustrates apparent lag effects in the area that passed a significant test.Additionally, for periods longer than 365 days, precipitation is found to contribute most to changes in the groundwater level.A comparison of Tables A1 and A2 shows that the river stage has the advantage of reducing groundwater level fluctuations across timescales of 30 to 365 days.However, precipitation has a more significant influence than the river stage in controlling groundwater level changes across all temporal scales.According to Figure A4, the greatest periodicity for the groundwater level and precipitation is between 365-768 days.The groundwater level illustrates apparent lag effects in the area that passed a significant test.Additionally, for periods longer than 365 days, precipitation is found to contribute most to changes in the groundwater level.A comparison of Tables A1 and A2 shows that the river stage has the advantage of reducing groundwater level fluctuations across timescales of 30 to 365 days.However, precipitation has a more significant influence than the river stage in controlling groundwater level changes across all temporal scales.

Figure 1 .
Figure 1.The study site: East Gulf Coastal Plain (EGCP), Alabama, (a) regional map illustrates the location of 41 observation wells, 74 pumping wells, main rivers, and faults used in this study.(b) The three-dimensional geologic model for EGCP.The name of each aquifer is listed in the legend.

Water 2023 , 25 Figure 2 .
Figure 2. Example of cross-sections derived from 48 boreholes (a).Cross-section A-A′ (b), and crosssection F′-F (c).Note that the stratigraphic columns using lithologic correlation (with also uniformities) shown in (b,c) end at the bottom of each borehole, which does not mean that the aquifer terminates abruptly (GMS matches the rock types below boreholes, which are not shown here since the plots (b,c) only illustrate the cross-sections derived by boreholes).

Figure 2 .
Figure 2. Example of cross-sections derived from 48 boreholes (a).Cross-section A-A (b), and cross-section F -F (c).Note that the stratigraphic columns using lithologic correlation (with also uniformities) shown in (b,c) end at the bottom of each borehole, which does not mean that the aquifer terminates abruptly (GMS matches the rock types below boreholes, which are not shown here since the plots (b,c) only illustrate the cross-sections derived by boreholes).

Figure 3 .
Figure 3.Comparison between the observed (symbols) and the modeled groundwater head at 41 observation wells.The line represents the 1:1 line.

Figure 3 .
Figure 3.Comparison between the observed (symbols) and the modeled groundwater head at 41 observation wells.The line represents the 1:1 line.

Figure 4 .
Figure 4. (a) The simulated groundwater table contour map (in November 2020) using MODFLOW, and (b) a graph showing the observed head and the residual for each computed head.

Figure 4 . 25 Figure 5 .
Figure 4. (a) The simulated groundwater table contour map (in November 2020) using MODFLOW, and (b) a graph showing the observed head and the residual for each computed head.Water 2023, 15, x FOR PEER REVIEW 10 of 25

2 .
Groundwater Age and Residence Time 3.2.1.Backward Particle Tracking to Calculate Groundwater Age Isotope ages for groundwater documented in the literature for the study area are applied to calibrate the backward tracking model built by MODPATH, since the backward tracking time represents the groundwater age if the vadose zone transport time is relatively short.For example, Solder [64] collected 14 C samples from 231 public supply wells in the South Atlantic and Gulf Coast, e.g., from the principal aquifer systems of the Southeastern Coastal Plain, Mississippi embayment-Texas coastal uplands, and the Coastal Lowlands.

Water 2023 , 25 Figure 7 .
Figure 7.The modeled groundwater age using MODPATH (backward particle tracking) for aquifers-aquitards in EGCP at different depths: (a) 300, (b) 700, and (c) 1000 ft below the groundwater table (where the unit of age is year).

Figure 7 .
Figure 7.The modeled groundwater age using MODPATH (backward particle tracking) for aquifersaquitards in EGCP at different depths: (a) 300, (b) 700, and (c) 1000 ft below the groundwater table (where the unit of age is year).

Figure 8 .
Figure 8.The modeled residence time using MODPATH for aquifers-aquitards in EGCP at different depths: (a) 300, (b) 700, and (c) 1000 ft below the groundwater table (where the unit of age is year).

Figure 8 .
Figure 8.The modeled residence time using MODPATH for aquifers-aquitards in EGCP at different depths: (a) 300, (b) 700, and (c) 1000 ft below the groundwater table (where the unit of age is year).
Water 2023, 15, x FOR PEER REVIEW 13 of 25

Figure 9 .
Figure 9.The modeled total age for groundwater in aquifers-aquitards in EGCP using MODPATH at different depths: (a) 300, (b) 700, and (c) 1000 ft below the groundwater table.

Figure 9 .
Figure 9.The modeled total age for groundwater in aquifers-aquitards in EGCP using MODPATH at different depths: (a) 300, (b) 700, and (c) 1000 ft below the groundwater table.

Figure 10 .
Figure 10.The probability density function (PDF) of groundwater age changing with the effective porosity  between 0.1 ≤ n ≤ 0.46 and 0.01 ≤ n ≤ 0.18 for EGCP aquifers located 300 (a), 700 (b), and 1000 ft (c) below the groundwater table.

Figure 10 .
Figure 10.The probability density function (PDF) of groundwater age changing with the effective porosity n between 0.1 ≤ n ≤ 0.46 and 0.01 ≤ n ≤ 0.18 for EGCP aquifers located 300 (a), 700 (b), and 1000 ft (c) below the groundwater table.

Figure 12 .
Figure 12.(a) Map shows the location of BAL-5 and DLE-2.The observed and computed head in the transient flow model for wells BAL-5 and DLE-2 at Layers 1 and 2, respectively, during a 1-year period (b) and a 10-year period (c).

Figure 12 .
Figure 12.(a) Map shows the location of BAL-5 and DLE-2.The observed and computed head in the transient flow model for wells BAL-5 and DLE-2 at Layers 1 and 2, respectively, during a 1-year period (b) and a 10-year period (c).

Figure 13 .
Figure 13.Transient model: examples of the predicted groundwater table from Wells BAL-5 and DLE-2 at Layers 1 and 2, respectively in 3 scenarios from November 2020 to 2030.

Figure 13 .
Figure 13.Transient model: examples of the predicted groundwater table from Wells BAL-5 and DLE-2 at Layers 1 and 2, respectively in 3 scenarios from November 2020 to 2030.

Figure A2 .
Figure A2.The continuous wavelet power spectrum of groundwater level (a).The thick contours indicate 5% significance levels against red noise.The pale regions indicate a cone of influence of edge effects which might distort the results.The color code for power values varies from low values (dark blue) to high values (dark red).(b) The continuous wavelet power spectrum of river stage.(c) The continuous wavelet power spectrum of precipitation.

Figure A2 .
Figure A2.The continuous wavelet power spectrum of groundwater level (a).The thick contours indicate 5% significance levels against red noise.The pale regions indicate a cone of influence of edge effects which might distort the results.The color code for power values varies from low values (dark blue) to high values (dark red).(b) The continuous wavelet power spectrum of river stage.(c) The continuous wavelet power spectrum of precipitation.

Figure A3 .
Figure A3.Cross wavelet coherence power spectrum (a) and bivariate wavelet coherency power spectrum (b) between groundwater level and river stage.The arrows indicate the relative phase relationships: point right for in-phase, point left for anti-phase, and point upward for 90° lag effect.The color shows the correlation coefficient of the two-time series.The thick contours indicate 5% significance levels against red noise.The pale regions indicate a cone of influence of edge effects which might distort the results.The color code for power values varies from low values (dark blue) to high values (dark red).

Figure A3 .
Figure A3.Cross wavelet coherence power spectrum (a) and bivariate wavelet coherency power spectrum (b) between groundwater level and river stage.The arrows indicate the relative phase relationships: point right for in-phase, point left for anti-phase, and point upward for 90 • lag effect.The color shows the correlation coefficient of the two-time series.The thick contours indicate 5% significance levels against red noise.The pale regions indicate a cone of influence of edge effects which might distort the results.The color code for power values varies from low values (dark blue) to high values (dark red).

Table 1 .
Lithology for the seven main aquifers located in EGCP.

Table A1 .
Results of wavelet coherency between groundwater level and river stage.All scales mean the average wavelet coherence at different temporal scales.Scale <30 Days 30-90 Days 90-180 Days 180-365

Table A1 .
Results of wavelet coherency between groundwater level and river stage.All scales mean the average wavelet coherence at different temporal scales.

Table A2 .
Results of wavelet coherency between groundwater level and precipitation.

Table A2 .
Results of wavelet coherency between groundwater level and precipitation.

Table A3 .
Parameter sensitivity results for the steady-state groundwater flow model for southern Alabama.