The AVA is located in the state of Aguascalientes in the high plains of Central Mexico (Figure 1). The elevation varies from 1,900 m on the valley floor to nearly 3,000 m in the surrounding mountain ranges (INEGI, 2013). The climate is semi-arid (BS1kw Köppen) with an average annual temperature ranging between 18 and 22°C, while the average annual precipitation is limited to 510 mm, due to its location relative to the subtropical high-pressure belt and the general orientation of the mountain ranges that limit and isolate it from the seas, distributed mainly between May and October (García Amaro de Miranda, 2003), and an average annual potential evaporation of 2,010 mm (National Water Commission; CONAGUA, 2020). Administratively, the AVA belongs to the Hydrological Region VIII “Lerma-Santiago-Pacífico” and the Hydrological Region No. 12 “Lerma-Chapala-Santiago.” In this report, we focus on the narrow valley region of the AVA, consistent with prior work on the AVA (CONAGUA, 2020; Guerrero Martínez, 2016, 2020), where data are available and where most of the human activity occurs.

The regional geology associated with the AVA is the result of the events that originated the Sierra Madre Occidental (CONAGUA, 2020), with units of different lithographic characteristics, including metamorphosed strata from the Cretaceous era, volcanic in the Sierra Madre Occidental from the Tertiary, and Quaternary, and basalt/alluvial material from the Quaternary. Based on well boring logs, CONAGUA (2020) describes the AVA as an unconfined, heterogeneous, and anisotropic aquifer, where groundwater flows through three different media: (i) porous media with primary, secondary, intergranular, and fractured permeability; (ii) fractured media with secondary permeability; (iii) double porous media with combined, intergranular and fracture permeability. The AVA was simplified considerably in the groundwater model developed for this work.

The data required for groundwater and transport modeling of the aquifer were obtained from the databases of CONAGUA and the local Citizen Commission for Drinking Water and Sewerage of the Municipality of Aguascalientes (CCAPAMA, 2022), as well as complementary sources from CONAGUA (2020) and previous studies (Hernández-Marín et al., 2018; Guerrero Martínez, 2016; De Figueroa Jesús, 2007). Groundwater piezometric head and water quality data (1,469 and 320 sampling points, respectively) was collected for the period 1985–2014.

To assess AVA groundwater conditions, we obtained data from CONAGUA (2020), which showed 1,830 registered users (1,769 deep-drilled wells and 61 shallow wells), with 1,468 active at the time. The available piezometric data is specified by sampling point (usually a production well) and collected annually (from 1990 through 2014), once the pumping is turned off and static conditions in the well are reached. Scarce and semiannual F- and As concentration data were obtained from 320 wells from 2003 to 2022, using a limited subset (years 2004–2014) to later compare the performance of the model. These concentrations were obtained from the water quality laboratory at the local water utility (CCAPAMA), in accordance with Mexican standards (NMX–AA–077–SCFI–2001, 2002—Test Method, 2016). Piezometric data and F- and As concentrations were spatially distributed through kriging to obtain gridded maps for each available year. Simple statistics (mean, standard deviation, maximum, and minimum values) were determined to describe and visualize the piezometric surface and chemical concentrations through time.

We used MODFLOW (Harbaugh, 2005) within the ModelMuse (Winston, 2022) users’ interface to develop the AVA groundwater model. We employed a single-layer model structure, owing to the preliminary nature of this effort. Given the digital elevation model (DEM) for the region and the well borehole logs, the aquifer model has a variable thickness, from 300 m to the west to 450 m to the East, which was consistent with thicknesses reported in prior work (Pacheco-Martínez et al., 2013). We calculated an effective value of hydraulic conductivity for each well using a weighted average of estimated values from well borehole records (SM-2, Sánchez San Román, 2022). Simple kriging was then used to create a gridded hydraulic conductivity product as the initial MODFLOW layer. We applied a uniform initial value of 0.16 for the aquifer-specific yield parameter (Sy), as proposed by CONAGUA (2020). The resulting model geometry included 335 active 2 × 2 km cells (828 total) arranged in 46 rows × 18 columns in the single, heterogeneous layer.

The model simulations required information on natural and anthropogenic inflow and outflow, such as potential recharge from precipitation (vertical recharge), well extractions based on contracted annual values, and horizontal groundwater flows. In order to estimate vertical recharge, we employed annual infiltration rates obtained from a monthly HBV precipitation-runoff model (Bergström, 1992) calibrated for the area (Hughes Lomelín, 2023). Regional inflows and outflows were distributed within the model boundaries, with values based on an AVA water budget created by CONAGUA (2020).

Starting with the initial piezometric surface, the groundwater model was executed in a transient state for 29 years (1985–2014), driven by annual changes in extraction and vertical recharge. The model was calibrated following a scheme focused on 10 wells located along a longitudinal transect (A-A’ in Figure 2), considering the piezometric levels of the years 1990, 1996, and 2000, and validated using 2007 and 2014. We adjusted the hydraulic conductivity (K) and specific yield (Sy) distributions in the model by applying a scaling factor, optimizing model-fitting metrics for the 11 transect well locations (Moriasi et al., 2007; Lucas Urbina, 2018; Roohollah et al., 2020; Akter and Ahmed, 2021) (Figure 2; Supplementary Tables 1, 2).

In the case of the solute transport model (MT3DMS, Zheng, 2012) utilizes groundwater discharge distributions from the groundwater flow model to distribute chemical solutes by advection and dispersion. Its model geometry was also developed using ModelMuse with the same computational mesh as the MODFLOW model, but reducing the total active cells to 248, to overlap with the region where water quality data were available. Initial fluoride and arsenic concentration conditions were set to the gridded 2003 values, and the source zones were indicated in MT3DMS cells associated with the maximum observed concentrations (Sathe and Mahanta, 2019). The transport simulations were limited to the period 2003–2014, the final year of the MODFLOW simulation.

The transport model calibration and validation were carried out similarly to the flow model, applying scaling factors to the initial longitudinal dispersivity and porosity values. Initial values for longitudinal hydrodynamic dispersion and effective porosity were assigned to two zones with different geology, corresponding to conglomerate sandstone and rhyolitic tuff. These values, both initial and calibrated, are shown in Table 1. In the absence of additional information about these contaminants in the aquifer, both fluoride and arsenic were assumed in this work to behave as nonreactive and nonsorbing solutes (i.e., conservative tracers). While the fluoride anion is not expected to react or sorb to aquifer sediments or confining units, arsenic is known to adsorb to various metal oxides and clays (Stollenwerk, 2003). However, assuming that the current regional dissolved concentrations represent long-term steady-state plumes with respect to the source zones, we can use the conservative tracer assumption to facilitate the simulations; although geogenic sources are expected in the area due to thermalism and similar behavior in other aquifers in the region (Morales-deAvila et al., 2023; Padilla-Reyes et al., 2024). Due to the limited simulation period, fluoride and arsenic observations were applied to model calibration and validation, respectively. The transport model was calibrated by comparing observed and simulated fluoride concentrations in a group of eight additional calibration cells, which differed from those used to calibrate the groundwater model due to the reduced modeled area and the location of the sampling wells. Kriging interpolation method was chosen for both, piezometric and water quality distributions, since correlations between observations in space is better represented by the spatial structure of the variogram, avoiding over-smoothing the spatial distributions.

Groundwater recharge as the main input to the system was obtained in a parallel work developed for the area, then it was later imposed as a lateral flow entering the active domain of the modeled aquifer. Later, the observed AVA piezometric surface was compared with observations, showing a decrease throughout the 1990–2014 period (Figure 2), as expected with increasing population and groundwater extractions during this period. The calibrated groundwater model reflected the slope of the observed surface, ramping downward from north to south in accordance with regional topography (see section A-A’ in part D of Figure 2, years 1990, 1996, 2000; Supplementary Table 1, models fit quality indexes; Supplementary Table 2, classification of the fit quality), but was consistently above the observed surface. For the validation of the groundwater model, we used the results for 2007, were simulated piezometric head is more consistent with the observed water table. Although 2014 simulations yielded acceptable results, they significantly overestimated drawdowns in the urban region of the AVA. Additionally, the hydrological model shows a lower mean precipitation in 2009–2014 (305 mm/year) compared to previous years (418.4 mm/year in 2004–2008), that might have impacted the recharge in the period. The argument of increasing urbanization and impervious surfaces supports the hypothesis of decreased recharge in the AVA area.

In the case of the transport model, the observed arsenic concentrations at the sources were consistently above the maximum permissible levels during the studied period, as were the fluoride concentrations, except for source 3 (Figure 3). For the calibration points, the F- concentrations range from about 1 to more than 4 mg/L; As levels clustered between 0.015 and 0.0175 mg/L for 2004 and 2006 data. Later observations found the arsenic data in the range of 0.009 and 0.017 mg/L, suggesting that some dilution occurred in later years. Recalling that the fluoride observations were used for calibration, the MT3DMS simulations show good agreement with observed concentrations over a range of values, particularly for the earlier 3 years (2004, 2006, 2008). The model slightly overestimates the observed concentrations in these years. However, in the later cases (2011, 2014), there is significantly more variation in the observed fluoride concentrations, and the model underestimates are more evident.

The validation data (arsenic) are difficult to describe, as the earlier outcomes (2004, 2006) exhibit a clustered pattern, with the simulation consistently overestimating the observed concentrations. Nevertheless, the validation improved in the later arsenic data (2008, 2011, 2014), as the clustering disperses.