The transboundary Hueco Bolson aquifer is the primary water supply for El Paso, Texas, and Ciudad Juárez, Mexico as well as U.S. Army Fort Bliss and smaller cities in Texas and Mexico. Binational groundwater data exchanges between the United States and Mexico that have focused on the aquifer have evolved over more than a century, shaped by scientific, political, and social dynamics. The history can be viewed through distinct periods. Efforts to understand the groundwater resource began with early reconnaissance surveys and have evolved through successive efforts to refine the characterization of groundwater flowpaths, residence times, surface water-groundwater interactions, and aquifer salinization. Lessons from these efforts highlight the importance of perseverance, mutual respect, and formal agreements, such as the 1999 memorandum between El Paso Water Utilities and the Junta Municipal de Agua y Saneamiento de Juárez, in sustaining long-term cooperation. The accumulated datasets chronicle the evolution of hydrogeologic conceptual and mathematical models while providing a foundation for ongoing research, sustainable water-use strategies, and the long-term stewardship of the Hueco Bolson aquifer shared by El Paso and Ciudad Juárez. The Hueco Bolson case demonstrates how long-term, cooperative data collection can improve scientific understanding and management of complex transboundary aquifer systems.

Groundwater is increasingly needed for water supplies but may have limited utility in some locations because of its salinity. Salinity, often expressed as total dissolved solid (TDS), is frequently estimated using specific conductance (SC) measurements. However, the commonly used proxy (0.65 multiplied by SC to indicate TDS, common in many handheld meters) can result in inaccurate TDS estimates. First, the TDS-SC relationship is not linear over the entire concentration range of groundwater. Furthermore, the TDS (and salinity)-SC relationships vary substantially depending on the major-ion composition. Here we develop a proxy method utilizing SC and major-ion water type to estimate TDS and salinity specifically for groundwaters. Compared to most surface waters, groundwater tends to have a wider range of salinity (fresh to highly saline) and higher concentrations of bedrock-derived solutes such as carbonate ions, silica, and many other ions. The dataset used to develop the proxies includes water chemistry data from 149,059 discrete groundwater samples. The groundwater proxies, which employ nonlinear log-log relations, utilize five water types (HCO, Cl, Ca-Mg-SO, Na-K-SO, and mixed waters), are accurate (median percent difference between TDS and salinity determined using the proxy compared to discrete measurements was <±0.8%) over a wide range of SC (up to 200 mS/cm), rapid, cost-effective, and can be measured on-site.

Reliable mapping and delineation of contaminant plumes and accurate estimation of contaminant mass discharge (CMD) are critical for groundwater risk assessment and planning of remedial actions at contaminated sites. However, traditional interpolation methods are often challenged by low-density sampling resulting in improper plume delineation. This study introduces a probabilistic censoring method that enhances geostatistical interpolation by incorporating comparably cheap, high-resolution, but semi-quantitative data collected from direct push-probes in the subsurface. The method converts halogen-specific detector signals into binary presence-absence indicators, which are interpolated using indicator kriging to generate a probability field of contaminant distribution. The probability field is then used to censor a spatial concentration field derived from traditional groundwater sampling, retaining interpolated concentration values only in areas where contamination is likely. We apply the method to a site contaminated with chlorinated solvents using two datasets with different sampling densities. Results show that, using our new method, plume fringes became more clearly defined and the total area with low concentrations (<10 μg L) increased by 41-85%. CMD estimates were reduced by 13-18%, while relative uncertainty remained largely unchanged. The method integrates seamlessly with traditional interpolation methods and our censoring workflow can be applied to other forms of direct-push data (e.g., relative permeability). As such, the framework offers a useful method for incorporating semi-quantitative field measurements into concentration interpolation and CMD estimation at contaminated sites.

Solute transport simulators aiming to accurately describe the transport of charged chemical species in porous media need to account for electrostatic coupling effects. Each ion in pore water possesses a specific electric charge and molecular diffusion coefficient, properties that determine their mobility and the overall charge balance of aqueous solutions. Depending on the charge, concentration and aqueous diffusion coefficient, the displacement of an ion in solution influences, and is in turn influenced by, other ions in solution by means of electrostatic interactions. This phenomenon has been studied with experiments and numerical simulations in diffusion-dominated regimes, as well as in advection-dominated flow-through systems, showing that electrostatic coupling effects play a relevant role in the spatiotemporal prediction of ion concentrations. However, there is limited availability of solute transport codes incorporating electrostatic coupling, limiting applications of multispecies ionic transport at different scales. This article elaborates on the topic of electrostatic coupling and presents a methodology for incorporating the effect into multispecies solute transport simulations with MODFLOW. The integration is achieved through the Application Programming Interface of the program (MODFLOW-API). This interface enables the access to concentrations and dispersion coefficients of all species during the simulation, which are necessary to calculate a dispersive correction that effectively incorporates electrostatic coupling into the model. Numerical results demonstrate the effectiveness of the coupling strategy, benchmarking the implementation with previously validated numerical simulators and with experimental data.

Land subsidence is widely present across the globe and brings catastrophic hazards. The well-acknowledged mechanism of subsidence is groundwater pumping, which leads to the reduction of hydraulic head and hence increases the effective stress, resulting in the vertical compaction of unconsolidated sediment. Here, we propose a hypothesis that subsidence in the coastal areas might be caused by osmotic effects, given the presence of seawater intrusion. The hypothesis is corroborated by simulating fluid flow, solute transport, and elastic deformation of multi-layered aquifer-aquitard systems. The simulations potentially cover a variety of natural environments by varying concentration, hydraulic head, thickness of aquitard, and hydraulic conductivity. We find that osmotic effects due to seawater intrusion play a non-negligible role in controlling subsidence in our studied cases, suggesting that future work on subsidence in areas impacted by seawater intrusion should fully incorporate osmotic effects to improve our understanding and prediction of subsidence.

We propose that a new term, aquitardifer, be added to the hydrogeologic nomenclature. Aquitardifer, a blend of the terms aquitard and aquifer, accounts for geologic materials that have properties of both as traditionally defined. Several examples of aquitardifers are provided, as is justification for and applicability of the term.

This study investigates the spatial and temporal sensitivity of aquitard hydraulic conductivity and specific storage on drawdowns in pumping tests. The objective is to understand which area of the aquitard is represented by drawdowns in different observation wells. A three-layered MODFLOW 6 model was used to simulate pumping tests on a circular Voronoi grid for three transmissivity scenarios and both confined and semiconfined top boundary conditions. A local sensitivity analysis was performed using PEST++ to determine how perturbations in hydraulic conductivity and specific storage of the aquitard affect head changes at observation wells in the pumped and overlying aquifer. Results indicate that for observation wells in the pumped aquifer, sensitivity forms an elliptical shape that is symmetrical around the observation well and the pumping well for all scenarios. The sensitivity map for the observation well in the overlying aquifer depends on the transmissivity ratio between both aquifers. It favors the area surrounding the pumping well if the transmissivity of the pumped aquifer is lower than that of the overlying aquifer. Conversely, with higher transmissivity in the pumped aquifer, sensitivity primarily lies around the observation well. Sensitivity patterns evolve over time, expanding the area of influence and shifting the sensitivity toward the observation well for a semiconfined top boundary. These findings are relevant for understanding the information regarding aquitard heterogeneity that is present in pumping test drawdowns and optimizing pumping test design.

Cross-sectional (2D) groundwater models are commonly applied to simulate complex processes that are challenging to capture using the coarse grids of 3D regional-scale models. 2D models are often extracted from 3D models for this purpose. However, translating groundwater properties from 3D to 2D models so that regional flow patterns are preserved poses several challenges. A methodology is presented here to maximize agreement between the heads of 2D and 3D groundwater models, considering MODFLOW models with rectilinear grids. This includes careful averaging of hydraulic properties and stresses from the 3D model to create commensurate properties and stresses in cross section. The approach was evaluated by examining the statistical match of transient heads within 10 cross sections extracted from a 3D model of the Limestone Coast (Australia). Concordance between 2D and 3D models was generally poor but was improved by incorporating lateral flow as inflows/outflows in 2D models. Lateral flows required inputs from the 3D model, which limits the application of 2D models as independent predictive tools. Pumping in the 3D model was redistributed to neighboring cells to reduce errors in the 2D model that arise from the limited capability to simulate pumping effects. Although pumping redistribution led to minimal improvement for the case study model, simpler modeling scenarios with more intense, localized pumping showed substantially better head matches between 2D and 3D models when pumping redistribution was applied. The methodology for creating cross-sectional models offered in this article provides relatively simple steps for creating 2D models that are consistent with 3D parent models, although further work is needed to develop a methodology for 2D models that are oblique to 3D model grids.

PEST++IES (White 2018; White et al. 2020) is widely used in the groundwater modeling community for its ability to perform computationally efficient history matching and uncertainty analysis in a highly parameterized context. One primary advantage of using an iterative ensemble smoother is that the number of model runs required per iteration depends on the number of realizations in an ensemble, not the number of parameters in each realization. However, this raises the question: what is the optimal number of realizations and iterations to use for any one model before the point of diminishing returns? Using a modified version of the Freyberg model (Freyberg 1988; Hunt et al. 2020), different parameter and observation scenarios were evaluated for four iterations and ensembles of 10, 25, 50, 100, 250, 500, 1000, and 2000 realizations. To match observations, PEST++IES altered hydraulic conductivity (k), both globally across the model and locally at three different pilot point densities, as well as global recharge (via a single multiplier), global river conductance, and individual well flow rates. Risk-based well capture zone results (Fienen et al. 2022a) and estimated hydraulic conductivity fields from each scenario were quantitatively and qualitatively compared against the "truth" model and its outputs. Across the cases examined, ensemble sizes of 100 to 250 realizations and two PEST++IES iterations were generally sufficient to achieve good results.

Borehole nuclear magnetic resonance (NMR) can be used to estimate the hydraulic conductivity (K) of unconsolidated materials. Various petrophysical models have been developed to predict K based on NMR response, with considerable efforts on optimizing site-specific constants. In this study, we assessed the utility of NMR logs to estimate K within highly heterogeneous glaciofluvial deposits by comparing them with K measurements from three types of co-located hydraulic testing methods, including permeameter, multi-level slug, and direct-push hydraulic profiling tool (HPT) logging tests. Four NMR models, including Schlumberger-Doll Research (SDR), Seevers, Sum-of-Echoes (SOE), and Kozeny-Godefroy (KGM), were applied to construct K profiles at four locations with model constants optimized using permeameter-based K. Model constants suitable for glaciofluvial deposits were provided. Results showed that NMR logging can provide reliable K estimates for interbedded layers of sand/gravel, silt, and clay. Through cross-hole comparison of NMR-derived K profiles, the trends and magnitudes of K for aquifers/aquitards were readily mapped. Quantitatively, the NMR-derived K coincided with hydraulic-testing K, with optimal model fits within one order of magnitude. We noticed that (1) Seevers performed similarly but no better than SDR in predicting permeameter and slug testing measurements; (2) SOE yielded slightly better predictions than SDR; (3) the removal of porosity in SDR did not deteriorate its prediction, and the optimized SDR constant resembled the literature-based values for glacial deposits; and (4) KGM yielded the optimal fits with slug-based K, demonstrating its reliable performance. Lastly, we made recommendations on selecting suitable petrophysical models.

Coastal lowlands are increasingly vulnerable to threats from sea-level and associated groundwater rise. This study introduces a categorical modeling framework that redefines groundwater depth estimation as a classification problem rather than a continuous prediction task. By dividing groundwater occurrence into multiple depth thresholds (0.9-2.0 m), the approach explicitly quantifies prediction uncertainty through Type I (false positive) and Type II (false negative) errors. A national-scale ensemble model developed at 100 m resolution using the Random Forest algorithm was trained on New Zealand's comprehensive depth-to-water database. Thirty-seven predictor variables, derived via PCA (97.5% variance retained) from 199 base predictors, were incorporated to capture the complex interactions influencing groundwater depth. The model demonstrates strong performance, with ROC-AUC values ranging from 0.823 to 0.962, and accuracy improves with increasing depth. This categorical framework addresses challenges associated with data imbalance and enhances uncertainty quantification compared to traditional regression methods. Probabilistic predictions allow stakeholders to set customizable risk thresholds and manage acceptable error levels based on specific coastal management contexts. By bridging the gap between advanced numerical modeling and practical adaptation planning, the approach provides a robust tool for evidence-based decision making in the face of rising sea levels.

The concentration of fluoride (F) in many geothermal waters worldwide exceeds the World Health Organization's drinking water guideline of 1.5 mg/L, highlighting a widespread water quality issue in geothermal systems. It is essential to predict the mobility of F in geothermal reservoirs compared to that of chloride (Cl). Column experiments were conducted at 45°C under hydrostatic pressures from atmospheric pressure to 12 MPa to simulate geothermal reservoir conditions, with the computer code CXTFIT 2.1 used to fit the data and determine transport and adsorption parameters of Cl and F. The effect of hydrostatic pressure on F transport was evaluated. Results showed that the Convection-Dispersion Equation (CDE) accurately fits the breakthrough curves of Cl and higher pressures enhance its dispersion transport. In contrast, F transport under atmospheric pressure shows rate-limited non-equilibrium, with the Two-Site Model (TSM) fitting F curves better. However, when pressure exceeds 6 MPa, equilibrium adsorption becomes more pronounced, and the equilibrium CDE better characterizes F transport. Higher pressure also increases the retardation factor (R) of F, suggesting more instantaneous adsorption sites under these conditions. The release of hydroxide (OH) during F adsorption increases pH values in the effluent notably. Moreover, greater pH fluctuations occur with increasing pressure due to the enhanced adsorption of F by the packed matrix. The findings would be helpful for understanding the solute transport and adsorption mechanism of fluoride under varying hydrostatic pressures in groundwater and geothermal systems.

Groundwater quality changes in wells and streams lag behind changes to land use due to groundwater travel times. Two contaminant transport methods were compared to assess differences in their simulated travel time distributions (TTDs) to streams and wells in the Wisconsin Central Sands. MODPATH simulates advective groundwater flow with particle tracking, while MT3D simulates age-mass using a finite difference solution without dispersion to allow for direct comparison of the two methods. MODPATH appropriately simulates groundwater TTDs from the water table to surface discharge but is subject to inaccuracies at weak-sink well cells due to the flow-model grid discretization and imprecise location of well discharge within well cells. MT3D better represents weak-sink well cells since it removes mass in proportion to the prescribed pumping rate, although travel time within well cells is neglected. Conversely, MT3D's treatment of surface water boundary cells is not as accurate as MODPATH because mass should be removed from the water table rather than the full cell volume. MT3D simulations of TTDs can also be confounded by the instantaneous vertical distribution of mass introduced throughout recharge cells instead of at the water table, which initiates mass along deeper flow paths. We evaluated 9 MODPATH and 13 MT3D implementations, generating differences in median travel times of up to 18 years. Both methods have strengths and weaknesses, with MT3D better representing weak-sink well cell behavior and MODPATH better representing surficial recharge and discharge. The effect of these characteristics on simulated TTDs, along with ideas for ameliorating method weaknesses, is discussed.

Groundwater models are important and useful tools for answering scientific and technical questions about the quantity and quality of groundwater, as well as for making critical management decisions. However, the heterogeneity of subsurface properties, such as hydraulic conductivity, is known to play a central role in groundwater flow and transport; therefore, its accurate quantification and incorporation into the groundwater workflow are critical. This paper presents a novel tool, ArchPy2Modflow, that efficiently combines a stochastic geological generator, ArchPy, with a groundwater flow software, MODFLOW. ArchPy2Modflow provides a rapid and practical way to convert and link any ArchPy model to a new (or existing) MODFLOW model, where any MODFLOW spatial parameter (such as porosity, hydraulic conductivity, or storativity) can be obtained from an ArchPy property, which is then upscaled according to the MODFLOW grid. ArchPy2Modflow offers several different options for selecting the appropriate MODFLOW grid: using the same grid as in the ArchPy model, defining each ArchPy geological unit as a MODFLOW layer, coarsening the grid by a certain factor, or directly using an existing MODFLOW grid. This flexibility enables users to adapt their models to suit their needs and constraints. The usefulness and practicality of the new tool are demonstrated by a synthetic example considering flow and transport in a heterogeneous aquifer, while the impact of a particular grid selection on the simulations is demonstrated.

Increasing water demands and declining groundwater levels have led to rising interest in managed aquifer recharge. That interest is growing in the United States-the focus of this article-and elsewhere. Increasing interest makes sense; managed aquifer recharge can reduce water-supply challenges and provide environmental benefits, sometimes with lower costs than alternative water-management approaches. But managed aquifer recharge also faces growing pains, which will make it difficult for projects to scale up and may limit the benefits provided by those projects that do go forward. Some of the problems arise from the challenges of finding physically suitable locations for managed aquifer recharge; many derive from economics, public policy, and law; and some derive from ways in which managed aquifer recharge could exacerbate traditional equity challenges of water management. But as we explain, there also are potential solutions to these challenges, and the future success of managed aquifer recharge will likely depend on the extent to which these solutions are adopted.

The tracer flowmeter and depth-dependent sampler (TFDDS) has been used to characterize flow and chemistry in hundreds of public supply wells. The TFDDS surveys address production and water quality compliance issues for utilities facing mounting treatment costs. Surveyed wells often show a stratified distribution for constituents of concern that were underestimated and poorly discretized using conventional zone testing in the pilot hole. Elevated concentrations of metals, semimetals, and radionuclides typically occurring at or near lithologic boundaries between fine- and coarse-grained sediments are missed since the primary focus of the zone test is to estimate the yield of the permeable sediments that are centered within a coarse unit. To minimize the risk of constructing new wells that fail to produce compliant drinking water relative to regulatory standards, we propose a new approach that vertically stacks flow and mass balance chemistry profiling of long screened test wells (LSTWs) under pumping conditions prior to installing the more costly public supply well. This approach integrates the TFDDS with a spectrum of drilling methods to rapidly provide hydraulic and chemistry data. A detailed distribution of anthropogenic and geogenic constituents within the saturated zones is produced, characterizing both permeable sediments and contaminant-laden boundaries between finer- and coarser-grained deposits. Applying this high-resolution data to well design leads to informed predictions on whether the new well will be compliant. Two case studies are presented, one for public supply well design and the other for managed aquifer recharge.

Climate change induced aridity and Euro-American settlement have altered the historical disturbance and flow regimes of large portions of the ponderosa pine forests of northern Arizona. The increased occurrence of high-severity wildfires due to these changes has led to the establishment of various forest restoration programs to protect the region's forests and their watersheds. In 2014, a paired-watershed monitoring project was implemented to compare the impacts of differing levels of forest thinning to watershed hydrology in seven experimental watersheds nested within the Upper Lake Mary (ULM) watershed in Arizona. This study expands the calibration phase of the ULM paired-watershed by synthesizing historic precipitation, surface runoff, groundwater recharge, soil moisture data, and evapotranspiration (ET) data to perform regression analyses and create a holistic water balance for each watershed. The magnitude and timing of seasonal groundwater recharge events were quantified for the first time in this region using a water table fluctuation method. The results showed that recharge did not occur every year and was heavily dependent (P < 0.05) on total winter season precipitation and snowpack duration. On average, recharge composed 9% of the total water budget when present. The results of this study lay the foundation for a greater understanding of how forest restoration alters northern Arizona's forest hydrology and will provide crucial information that should be used in water policy and water resource decision-making as the region plans for future water availability.

Adjoint sensitivity analysis provides an efficient alternative to direct methods when evaluating the influence of many uncertain parameters on a limited number of performance measures in hydrologic and hydrogeologic models. However, most adjoint implementations are "intrusive", requiring extensive modifications of the forward simulation code. This creates significant development and maintenance burdens that limit broad adoption. To address these needs, we present MF6-ADJ, a "non-intrusive" adjoint sensitivity capability for the MODFLOW 6 groundwater flow process that leverages the MODFLOW Application Programming Interface (API) to interact with the forward groundwater flow solution without altering its core code. MF6-ADJ supports both confined and unconfined flow conditions, structured and unstructured grids, and is compatible with both the Standard and Newton-Raphson solution schemes. It computes sensitivities of a wide range of general performance measures, including hydraulic heads, boundary fluxes, and weighted residuals, with respect to key model parameters such as hydraulic conductivity, storage coefficient, injection/extraction rate, recharge rate, boundary head, and boundary conductance. Sensitivities are computed at each node, enabling fine-grained diagnostic and calibration analysis. Validation against analytical solutions and the finite-difference perturbation method confirms excellent agreement, while demonstrating speedups ranging from hundreds to tens of thousands of times depending on grid discretization, since the adjoint state method computes sensitivities efficiently at the grid-block level. This non-intrusive design makes MF6-ADJ highly accessible and maintainable, offering efficient and scalable sensitivity analysis in complex groundwater modeling workflows.