Water Supply Source Evaluation in Unmanaged Aquifer Recharge Zones : The Mezquital Valley ( Mexico ) Case Study

The Mezquital Valley (MV) hosts the largest unmanaged aquifer recharge scheme in the world. The metropolitan area of Mexico City discharges ~60 m3/s of raw wastewater into the valley, a substantial share of which infiltrates into the regional aquifer. In this work, we aim to develop a comprehensive approach, adapted from oil and gas reservoir modeling frameworks, to assess water supply sources located downgradient from unmanaged aquifer recharge zones. The methodology is demonstrated through its application to the Mezquital Valley region. Geological, geoelectrical, petrophysical and hydraulic information is combined into a 3D subsurface model and used to evaluate downgradient supply sources. Although hydrogeochemical variables are yet to be assessed, outcomes suggest that the newly-found groundwater sources may provide a long-term solution for water supply. Piezometric analyses based on 25-year records suggest that the MV is close to steady-state conditions. Thus, unmanaged recharge seems to have been regulating the groundwater balance for the last decades. The transition from unmanaged to managed recharge is expected to provide benefits to the MV inhabitants. It will also be likely to generate new uncertainties in relation to aquifer dynamics and downgradient systems.


Introduction
Aquifer exploration and water supply evaluation is a key challenge in low-and mid-income countries [1][2][3], where the characterization of potential groundwater sources is still needed and where a spectacular increase in groundwater use has occurred in the past decades [4].
Groundwater exploration is usually most needed in arid regions.This is because scanty rainfall and high evapotranspiration rates tend to make surface sources unreliable.In these areas, aquifers play an essential role in ensuring sustainability, particularly whenever unregulated recharge schemes are Evaluating groundwater resources in the vicinity of the world's largest unmanaged recharge region poses a major technical challenge.To address it, we develop an inclusive methodology based on geological reservoir modeling techniques.This encompasses all aspects related to stratigraphic, lithological and petrophysical properties of subsurface rocks, and serves the purpose of estimating the spatial distribution of target units and the volume of economic fluids [5].Available information for integrative modeling combines static and dynamic data at different scales.Static records include geophysical surveys, well logs and geological data [6] whilst dynamic information provides evidence of the interactions between the static model and the fluid over time.This approach delivers an appropriate means to deal with incomplete datasets [7,8] and has been the standard in oil and gas reservoir modeling for at least three decades [9,10].
The aim of reservoir characterization is to understand and identify regional flow units, i.e., hydrofacies, as well as to assess the distribution of relevant properties in the vicinity of a borehole [11].Static reservoir modeling is frequently coupled with complementary geophysical techniques.This is because drilling data is characterized by high vertical resolution-low horizontal resolution, whereas geophysical soundings present the opposite pattern [8].In the case of the former, well logs provide a quantitative method to assess near-borehole petrophysical signatures.Within the latter, Evaluating groundwater resources in the vicinity of the world's largest unmanaged recharge region poses a major technical challenge.To address it, we develop an inclusive methodology based on geological reservoir modeling techniques.This encompasses all aspects related to stratigraphic, lithological and petrophysical properties of subsurface rocks, and serves the purpose of estimating the spatial distribution of target units and the volume of economic fluids [5].Available information for integrative modeling combines static and dynamic data at different scales.Static records include geophysical surveys, well logs and geological data [6] whilst dynamic information provides evidence of the interactions between the static model and the fluid over time.This approach delivers an appropriate means to deal with incomplete datasets [7,8] and has been the standard in oil and gas reservoir modeling for at least three decades [9,10].
The aim of reservoir characterization is to understand and identify regional flow units, i.e., hydrofacies, as well as to assess the distribution of relevant properties in the vicinity of a borehole [11].Static reservoir modeling is frequently coupled with complementary geophysical techniques.This is because drilling data is characterized by high vertical resolution-low horizontal resolution, whereas geophysical soundings present the opposite pattern [8].In the case of the former, Water 2017, 9, 4 3 of 25 well logs provide a quantitative method to assess near-borehole petrophysical signatures.Within the latter, deep horizons (usually derived from seismic data) deliver the geometric, structural and stratigraphic properties of subsurface units [12][13][14][15][16][17][18].Moreover, advances in three-dimensional (3D) mapping provide alternative means to developing visual representation of subsurface reservoirs and quantifying some of their key properties [7,[19][20][21].
Static modeling frameworks have been successfully used to deal with a variety of practical problems.These include geological CO 2 sequestration and storage [22,23], uncertainty analysis in enhanced oil recovery processes [24], land subsidence related to subsurface gas development [21], correlation between oil field analogues [25] or the understanding of geothermal systems [26], among others.
Mainstream static reservoir modelling techniques, including 3D mapping or well logging, have gained recognition in the field of groundwater exploration in recent years [27][28][29][30][31].However, integrated approaches (e.g., [32]) such as the one we present are comparatively scarce in the groundwater literature.More specifically, these have seldom been used for the purpose of characterizing groundwater bodies, and their application in the context of unmanaged aquifer recharge is practically non-existent.
A combination of geological, geoelectrical, petrophysical and hydraulic methods is used to gain an understanding of the interaction between lithological factors, hydrogeological properties and the spatial distribution of groundwater.Time domain electromagnetic (TDEM) techniques were used for subsurface modeling.TDEM has been widely used to characterize potential groundwater zones [33][34][35].The TDEM-resistivity model was correlated with drilling data, and then interpreted jointly with well logs and pumping tests so as to deduce hydrogeological properties and borehole performance.
Within this context, the main contribution of this paper to the literature consists in demonstrating how these techniques may lead to locating and assessing reliable long-term sources of drinking water in unmanaged regions.This research presents additional value by dealing with a rather unique hydrogeological setting in the world.

Geological and Hydrogeological Setting
The study area corresponds to the Alfajayucan aquifer system (Figure 1b,c), which in turn is located downstream from the Mezquital Valley, central Mexico (Figure 1c).Climate is predominantly semiarid.Rainfall amounts to 450 mm/year, while potential evapotranspiration is 2100 mm/year [36].
Mexico City discharges ~60 m 3 /s of raw wastewater without formal treatment into the Mezquital Valley (Figure 1b,c).Wastewater is transported via a complex network that includes the Central Emitter, the Grand Canal and the Western wastewater outlets (Figure 1c).These are all part of a ~90,000 ha irrigation district that straddles across the Tula and Alfajayucan districts [36] (Figure 1c).
Artificial recharge within the Mezquital Valley is controlled by unmanaged wastewater discharge and estimated at 25 m 3 /s, that is, over 13 times the natural recharge rate [36].This has caused the water table to rise since the 1950s.New artesian wells and springs have appeared as a result.Some showcase substantial flowrates.Take for instance the Cerro Colorado spring (~600 L/s) [37].
Overall, the Mezquital Valley is the largest raw wastewater-based agricultural irrigation region in the world [37].For nearly 110 years, raw wastewater has been used as a source of irrigation water and nutrients for crops.This has allowed local farmers to minimize production costs [38].Numerous authors deal with the environmental impacts of wastewater irrigation [38], organic micropollutants [39], leaching potential of pharmaceutical components [40], health risk assessments [41] and groundwater pollution [42].These and other environmental concerns led the Federal Government to build the Atotonilco water treatment plant in 2013 [43].The plant was designed to treat ~35 m 3 /s of wastewater by means of biological and physicochemical processes.Operations are expected to begin in 2016-2017.
The Mezquital Valley provides a unique example of unmanaged aquifer recharge (UAR).UAR is described as the disposal of water that results in unintentional aquifer replenishment [44], in contrast with managed recharge schemes [45,46].There have been examples of longstanding Water 2017, 9, 4 4 of 25 unmanaged recharge practices in countries such as Australia [47].More recent ones pertain to Italy [48] or Oman [49].In other regions of the world UAR has received less attention.
Previous work suggests that the Mezquital Valley might present hydraulic continuity with downgradient permeable formations, such as the volcanic deposits of the Alfajayucan aquifer [50].Population nuclei located in these areas (San Agustín Village, Ixmiquilpan, Alfajayucan and others, see Figure 1c) are thus likely to be influenced by unregulated activities in terms of water quality and quantity.
Regional geology includes volcanic rocks and calcareous sediments (Figure 2).Lower Cretaceous rocks (El Doctor, Mexcala and Soyatal formations) are made up of limestone with interbedded clay horizons.A dolomite formation caps the sequence.Tectonic activity presents a NW-SE sequence of anticlines and synclines.
Middle and Upper Tertiary deposits such as the Tarango Formation (the main supply aquifer for Mexico City) are made up of lava flows and volcanic ash of basaltic-andesitic composition, with interbedded lacustrine deposits [37].
Tertiary volcanic rocks of the Donguinyó and Huichapan formations, related to the eruption of the Donguinyó-Huichapan caldera complex, predominate across the region.The eruptions released discrete pulses of intermediate flows and minor rhyolitic, ignimbrites, pumice fall and surge deposits [51].
From the regional standpoint (Mezquital Valley and surroundings areas), a multi-layer aquifer system has been proposed as a valid conceptual groundwater flow model.However, there are different conceptualizations.The British Geological Survey (BGS) and CONAGUA [37] delineated three overlapping aquifers: (a) the lava flows and coarse sediments of the Tarango Formation; (b) the Tertiary lava-pyroclastic deposits; and (c) the Cretaceous limestones.Other authors reported the existence of two aquifers based on water level measurements [52], while the National Water Commission [36] describes the Alfajayucan aquifer as a two-layer system.
Previous work suggests that the Mezquital Valley might present hydraulic continuity with downgradient permeable formations, such as the volcanic deposits of the Alfajayucan aquifer [50].Population nuclei located in these areas (San Agustín Village, Ixmiquilpan, Alfajayucan and others, see Figure 1c) are thus likely to be influenced by unregulated activities in terms of water quality and quantity.
Regional geology includes volcanic rocks and calcareous sediments (Figure 2).Lower Cretaceous rocks (El Doctor, Mexcala and Soyatal formations) are made up of limestone with interbedded clay horizons.A dolomite formation caps the sequence.Tectonic activity presents a NW-SE sequence of anticlines and synclines.
Middle and Upper Tertiary deposits such as the Tarango Formation (the main supply aquifer for Mexico City) are made up of lava flows and volcanic ash of basaltic-andesitic composition, with interbedded lacustrine deposits [37].
Tertiary volcanic rocks of the Donguinyó and Huichapan formations, related to the eruption of the Donguinyó-Huichapan caldera complex, predominate across the region.The eruptions released discrete pulses of intermediate flows and minor rhyolitic, ignimbrites, pumice fall and surge deposits [51].
From the regional standpoint (Mezquital Valley and surroundings areas), a multi-layer aquifer system has been proposed as a valid conceptual groundwater flow model.However, there are different conceptualizations.The British Geological Survey (BGS) and CONAGUA [37] delineated three overlapping aquifers: (a) the lava flows and coarse sediments of the Tarango Formation; (b) the Tertiary lava-pyroclastic deposits; and (c) the Cretaceous limestones.Other authors reported the existence of two aquifers based on water level measurements [52], while the National Water Commission [36] describes the Alfajayucan aquifer as a two-layer system.
The upper layer is made up of granular and fractured volcanic rocks up to ~400 m thick, while the lower (confined) layer is made up of limestone sediments and calcareous clays.Target units within the Alfajayucan aquifer correspond to the Donguinyó and Huichapan volcanic deposits.The upper layer is made up of granular and fractured volcanic rocks up to ~400 m thick, while the lower (confined) layer is made up of limestone sediments and calcareous clays.Target units within the Alfajayucan aquifer correspond to the Donguinyó and Huichapan volcanic deposits.
The Mezquital Valley relies on groundwater for domestic supply.There are 456 wells, which extract ~98 Mm 3 /year [52].Approximately 17% of this volume corresponds to human consumption, 38% to agriculture, 33% to the industry and 12% to other uses [54].
On the other hand, 76 water sources were documented in 2012-2013 for the Alfajayucan regional aquifer.These include 45 boreholes, nine dug wells, 21 springs and one gallery.Groundwater abstraction is ~8 Mm 3 /year, out of which 64% is for human consumption, 23% for agriculture and 13% for other uses [36].Within the study area (Figure 2), 11 pumping wells were identified.Out of these, nine are currently inactive.There are significant uncertainties in relation to their depth, screening and hydraulic behavior.

Hydrogeophysical Investigation
Geophysical data was performed in two stages (Figure 2).The first one consisted of 26 soundings.These were acquired with a single loop configuration, 300 m by side.Five additional TDEM were acquired during the second stage with the same loop configuration, but with a size of 150 m side.All measurements were taken using a TerraTem equipment by Alpha Geoscience Co. (Clifton Park, NY, USA) maintaining an intensity of 7-8 amperes.After editing and averaging the decay curves, apparent resistivity curves were obtained for each TDEM station.One-dimensional (1D) models were computed using an Occam inversion scheme [55].

Borehole Drilling and Well Log Analysis
A 200-m, 12 -diamater borehole was drilled using a Gardner-Denver 14 W rig. The location was established based on the integration of geological and geophysical data (Figure 2).A static water table was found at a depth of 81 m.Geophysical logs were gathered so as to obtain the petrophysical properties of subsurface units.An 8144 Century Geophysical Co device was used for this purpose, with a voltage of 36 VDC and a logging speed of 9 m/min.This allowed to log variables such as temperature, spontaneous potential (SP), short and long normal resistivity, resistivity of the drilling fluid and gamma ray (GR).
Rock physics is based on sedimentary sequences [56], where oil and gas reservoirs usually take place.In volcanic environments, however, the radioactive nature of rocks makes the interpretation of these curves more challenging [57].While in sedimentary materials the GR curve is controlled by clay content, the magma composition plays a major role in volcanic rocks.Based on the analysis of 254 samples, Stefansson et al. [58] concluded that there is a linear relationship between GR and silica content in volcanic materials.Though unimportant in the case at hand, it must be noted that this relationship needs to be handled with care in instances where silica content may be altered by the presence of clays.Hence, the borehole log could be completed and analyzed from a quantitative perspective.Completion was necessary because it was only possible to recover material from the upper 66 m of the borehole.
Lithology was inferred from silica content.Criteria established by [57] was used for this purpose.Basalt corresponds to a silica content of <52%, andesitic-basalt to 52%-57%, andesite to 57%-63% and dacite-rhyolite to >63%.Silica content was then checked against standard igneous petrology classifications, to ensure that the inferred lithologies were consistent with chemical composition.
Volcanic facies (lithofacies) were deduced based on silica content and GR data.These were in turn associated to their hydraulic properties (hydrofacies).On the other hand, groundwater salinity was estimated from the SP curve.Its response is controlled by electrochemical and electrokinetic variables, although the latter are generally negligible [56].In-situ measurements were obtained using a Solinst 428 PVC bailer.Electric conductivity and total dissolved solids were measured by means of a Hannah HI 9828 multiparametric device.Field readings were then compared to theoretical salinity contents.Total porosity is estimated by means of Archie's Law [59].Archie's equation remains a method of choice for the estimation of porosity from petrophysical data [60].Archie's equation is valid for saturated non-clayey media, and is expressed as follows [61]: where ρ corresponds to the resistivity of saturated rock [Ω•m], ρ w is the resistivity of pore water [Ω•m], m the cementation coefficient [-] and a the tortuosity factor [-]. Equation ( 1) was used to develop a porosity curve (φ) for the saturated zone.Input data included the long normal resistivity curve and the estimated groundwater resistivity (as inferred from the SP log).The cementation and tortuosity coefficients were modified as per [62].Four configurations were analyzed: (1) sandy media (m = 2.5, a = 0.62); (2) carbonate media (m = 2, a = 1); (3) vitreous media (m = 2, a = 1.85); and (4) m = a = 1.The lowest porosity values were obtained with configuration (4).This is the most consistent estimate with reported data for volcanic environments [63].Finally, all information was integrated within a Schlumberger Petrel ® platform [64] to perform petrophysical data analyses.

Well Yield and Aquifer Parameters
The borehole was equipped as a groundwater well and tested for aquifer parameters and expected yields.A 10 -diameter stainless steel screen casing (grid size 0.125 ) and a gravel pack was fitted to the 100-200 m interval.The remainder of the outer annulus was sealed with a 4:1 cement-bentonite proportion.
The well was developed to remove all remaining drilling fluid and settle the gravel pack (procedures in [65]).Pressurized air was injected for 17 h, until a stable flow with negligible fine-particle content was obtained.
A 350-Horse Power submersible pump was installed at a depth of 100-150 m and used to carry out pumping tests.The pump was connected to a 150-m long 4 -diameter pipe.A 42-h 7-step pumping test was carried out to determine well yield, efficiency and aquifer parameters.Steps included: (1) flow rate = 534 m 3 /day from 0 to 360 min; (2) 817 m 3 /day from 360 to 720 min; (3) 1130 m 3 /day from 720 to 1080 min; (4) 1332 m 3 /day from 1080 to 1440 min; (5) 1532 m 3 /day from 1440 to 1800 min; (6) 1714 m 3 /day from 1800 to 2160 min; and (7) 1904 m 3 /day from 2160 to 2520 min.Field measurements (elapsed time vs. drawdown vs. flow rate) are shown in Table 1.Results allowed for the evaluation of head losses associated to the aquifer properties and to the borehole itself.Classic interpretation models were used for this purpose [66,67].These assume that theoretical drawdown and head losses can be estimated by means of Equation (2): where s w corresponds to the total drawdown inside the well [L], Q is well flow rate for each step and CQ n is the nonlinear component of head losses Water 2017, 9, 4

of 25
Hydraulic parameters were estimated by curve-matching procedures, using a nonlinear weighted least-squares parameter estimation and the Gauss-Newton linearization method.Field data was fitted to the Theis model [68] for a pumping well in a confined aquifer, modified to include linear and nonlinear well losses in a step-drawdown test [69] (Equation (3)): where w is the well function, i.e., the exponential integral function, r represents the well radius [L], S and T correspond to the storage coefficient [-] and aquifer transmissivity [L 2 T −1 ], respectively; t is the elapsed time [T] and σ the skin factor [-], estimated by means of Equation ( 3).The skin factor is an indicator of well damage and has been widely evaluated in the oil reservoir practice [70,71].
The skin factor causes an additional pressure drop, ∆p σ , across the altered zone in hydrocarbon wells [72]: Equation ( 4) was adapted to cater for the conditions of groundwater systems.Considering that the oil-formation volumetric factor, B o , and viscosity, µ, equals 1 in a single-flowing phase, and that matrix permeability (k)-reservoir thickness (h) product is equivalent to T, Equation ( 4) can be rearranged as follows: where ∆h σ represents the additional drawdown caused by the skin factor [L].

Hydrogeoelectrical Model
TDEM resistivity models were interpolated in 10 m steps for the upper 100 m, 20 m between 100 and 200 m and 40 m for depths up to 400 m.A digital elevation model (DEM) was integrated in Petrel [64] as a reference for subsurface layers.Resistivity curves were upscaled and interpolated by means of a moving averages model with a fourth-power inverse distance algorithm.The hydrogeophysical resistivity model was used to correlate hydrofacies based on horizons and surfaces, as commonly done in the oil industry [73].A square 800 × 800 m mesh was used to devise the hydrogeological units.Additionally, a vertical resistivity gradient was developed and used to identify subsurface electrofacies, associated to hydrostratigraphic changes.

Total Porosity Model
Borehole data was integrated with TDEM results using Archie's equation [59].In the absence of predominantly clayey minerals, water resistivity is assumed homogeneous across the study area.Conductivity is assumed to be ionic at constant temperature and the effect of chemical thermalism is considered negligible.

3D Model
Petrel [64], a 3D oil and gas reservoir modeling software, was used to integrate petrophysical data with TDEM results into a 3D visualization model.The model relies on 31 TDEM inverted resistivity curves, geophysical well logs, porosity information, groundwater depth in nine wells, surface geology, a digital elevation model (DEM), the inferred geometry of electrofacies and the geological information from the borehole drilling.Surface geology [51,53] was only useful to define the uppermost part of the top layer.All other layers were inferred from TDEM information.In other words, the model can only be used to infer subsurface geoelectrical properties.
The DEM was based on four DEM scenes from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER GDEM v.2, NASA, California Institute of Technology, Pasadena, CA, USA) and then imported to Petrel as a shape file.DEM scenes present a native pixel resolution of 3 × 3 arcsec.Individual scenes were re-projected to a UTM coordinate system (zone 14 N) with a 90 m × 90 m pixel size.Bilinear resampling was used to account for realistic spatial deformation.This process was carried out using QGIS 2.12 Lyon (QGIS Development Team, Open Source Geospatial Foundation) [74].
The development of the 3D model consisted in the following steps: (1) Grid discretization and creation of geological horizons, based on the TDEM interpretation model.For that purpose, we built a three-layer model (layer-1 from 0 to 100 m; layer-2 from 100 to 200 m; layer-3 from 200 to 400 m), with a grid size of 2 km × 2 km, covering a surface area of ~300 km 2 .Each layer was divided into 10 to 40 m-thick sublayers (e.g., geocellular grids) to gain vertical resolution; (2) Geoelectrical input data using resistivity population from each of the 31 TDEM inverted curves (see Section 2.2); (3) Creation of a 3D electrofacies model of the subsurface that enables the population all geocellular grids with discrete resistivity values of the electrofacies EF-1 (vadose zone), EF-2 (aquifer) and EF-3 (basement).For that purpose, a deterministic modeling approach was used by means of a moving averages algorithm, taking advantage of the Petrel Facies Modeling Module [64].
Since the model is built upon geoelectrical information, uncertainties can be expected to be higher in those areas which are located further away from the 32 control points (e.g., 31 TDEM soundings plus the borehole).

Groundwater Flow Evolution in the Mezquital Valley
Groundwater flow behavior in unmanaged recharge regions is a key factor when unregulated water disposal is a major recharge component for downgradient supply aquifers.In the case at hand, an investigation conducted by Geocalli [50] proposed that the Mezquital Valley might present hydraulic continuity with downgradient aquifers, including the study area (Figures 1c and 2).
In order to assess potential hydraulic connectivity and to evaluate if a decrease in the unmanaged wastewater disposal may compromise downgradient aquifer replenishment, piezometric records were analyzed.These are considered a suitable proxy for temporal changes in the wastewater input from the Mezquital Valley.
Datasets from 1982 [50], 1998 [37] and 2007 [52] were taken into consideration.Water table depth contours were digitized from the original sources and the temporal changes of groundwater levels between 1982 and 1998, 1998 and 2007 and 1982 and 2007 were assessed by contouring the hydraulic evolution.Three new maps were created as a result, based on level measurements in 152, 59 and 30 production wells, respectively.
The groundwater evolution of each period of time was analyzed using conventional statistical tools.The spatial and statistical analysis was carried out in QGIS 2.12 Lyon [74].

Hydrogeoelectrical Model
Figure 3 presents the hydrogeoelectrical 2D model (cross sections are shown in Figure 2).Electrofacies 1 (EF-1) corresponds to the uppermost unit.EF-1 presents a relatively homogeneous geometry.The average thickness is ~70 m, and resistivity ranges between 40 and 300 Ω•m, with a bimodal behavior (modes in 50 and 110 Ω•m).Gradients are negative and observed to change laterally.This unit is interpreted as the vadose zone of the system and correlated with a heterogeneous mix of lacustrine, alluvial and volcanic deposits.Electrofacies 3 (EF-3) provides a hydrogeophysical contrast with EF-1 and EF-2.Resistivity is comparatively low (4-20 Ω•m) at depths of 250-400 m.EF-3 consists of low-permeability clay sediments, marls and clayey limestones, probably related to the Soyatal Formation.Thus, unit EF-3 might represent the system basement.Our geoelectrical interpretation is consistent with those provided by other authors, which suggested that resistivity tends to diminish with increased clay content [47].
TDEM surveys allowed to identify EF-2 as the unit with a better aquifer potential in the region.Resistivity values (5-170; 70-90; 20-80 Ω•m) are consistent with the existing knowledge about the area's subsurface geology, as well as with standard geoelectrical values for similar geological formations [75,76].Correlated materials (Huichapan and Donguinyó tuffs) include lithological contrasts associated to the eruption of lava and pyroclastic flows [51].Unit EF-2 presents a marked resistivity contrast with the basement unit, EF-3.

Reconstructed Borehole Lithology
Figure 4 shows the results of geophysical logs carried out in the borehole.All readings except GR begin at the static groundwater level.The incomplete geological record (from 66 to 200 m) was reconstructed primarily from the GR curve.GR yielded values of ~60°-90° API (API = American Petroleum Institute) for the upper 40 m, growing to 90°-110° API between 40 and 60 m.The former is controlled by the presence of radioactive minerals in andesites and dry basalts.Increase with depth can be explained by the presence of a 20-m thick dacite horizon where an intensification in radioactivity was observed.From 60 to 200 m, values remain uniform at ~30°-60° API.Some radioactive peaks were also identified (e.g., 164-174 m).These were attributed to the presence of silica-rich rocks.
High GR readings usually hint at the presence of clay-rich materials.However, this does not appear to occur.A number studies in igneous terrains have demonstrated that a relationship exists Electrofacies 3 (EF-3) provides a hydrogeophysical contrast with EF-1 and EF-2.Resistivity is comparatively low (4-20 Ω•m) at depths of 250-400 m.EF-3 consists of low-permeability clay sediments, marls and clayey limestones, probably related to the Soyatal Formation.Thus, unit EF-3 might represent the system basement.Our geoelectrical interpretation is consistent with those provided by other authors, which suggested that resistivity tends to diminish with increased clay content [47].
TDEM surveys allowed to identify EF-2 as the unit with a better aquifer potential in the region.Resistivity values (5-170; 70-90; 20-80 Ω•m) are consistent with the existing knowledge about the area's subsurface geology, as well as with standard geoelectrical values for similar geological formations [75,76].Correlated materials (Huichapan and Donguinyó tuffs) include lithological contrasts associated to the eruption of lava and pyroclastic flows [51].Unit EF-2 presents a marked resistivity contrast with the basement unit, EF-3.API.Some radioactive peaks were also identified (e.g., 164-174 m).These were attributed to the presence of silica-rich rocks.

Near-Borehole Petrophysical Response
High GR readings usually hint at the presence of clay-rich materials.However, this does not appear to occur.A number studies in igneous terrains have demonstrated that a relationship exists between the intensity of GR response and the geochemical composition of volcanic rocks [57,77], instead of clayey minerals in sedimentary environments.
Criteria by [58] was used to obtain the SiO 2 content of volcanic rocks as a response of GR values.Lithology was then inferred from silica content.The reconstructed borehole lithology is shown in Figure 4.
between the intensity of GR response and the geochemical composition of volcanic rocks [57,77], instead of clayey minerals in sedimentary environments.
Criteria by [58] was used to obtain the SiO2 content of volcanic rocks as a response of GR values.Lithology was then inferred from silica content.The reconstructed borehole lithology is shown in Figure 4.

Hydrofacies
Figure 5 shows the correlation between different petrophysical variables.These were used to understand the relationship between petrophysical response, groundwater occurrence and aquifer properties.Resistivity depends on compaction and fracturing, therefore, it is assumed that both can be obtained, at least qualitatively, from the resistive range of each lithological unit.

Hydrofacies
Figure 5 shows the correlation between different petrophysical variables.These were used to understand the relationship between petrophysical response, groundwater occurrence and aquifer properties.Resistivity depends on compaction and fracturing, therefore, it is assumed that both can be obtained, at least qualitatively, from the resistive range of each lithological unit.(a) Vadose zone (thickness ~81 m), made up of massive, scarcely fractured basalt.This is not only consistent with collected data, but also with direct observation of rock bits obtained from a depth of 64-66 m.(b) Fractured lava aquifer (thickness ~40 m).Fracturing is expected to be less intense towards the upper boundary.Total porosity is estimated at 7%-24%, resistivity is 40-150 Ω•m, radioactivity 10°-25°API and silica content <52%.(c) Granular aquifer in andesitic breccia and tuffs (thickness ~18 m).The granular aquifer presents some degree of fractured behavior.Total porosity is 12%-18%, resistivity is 25-55 Ω•m, radioactivity 45°-63°API and silica content 48%-62%.(c) Granular aquifer in andesitic breccia and tuffs (thickness ~18 m).The granular aquifer presents some degree of fractured behavior.Total porosity is 12%-18%, resistivity is 25-55 Ω•m, radioactivity 45 • -63 • API and silica content 48%-62%.(d) Aquitard, made up of andesite and acidic rocks (thickness ~11 and 5-7 m, respectively).
The vertical distribution of hydrofacies and lythofacies is shown in Figure 4. Petrophysical features are presented in Figure 5a,b.

Porosity and Qualitative Permeability Estimates
The permeability of hydrofacies was estimated on qualitative terms.Well logs establish that the resistivity of the drilling fluid (4-6 Ω•m) is much lower than that of the hydrofacies (25-200 Ω•m).Under these conditions, positive deflections in the SP curve could be interpreted as permeable sectors.In the absence of clay, the larger the deflection, the more permeable the rock would be.SP and resistivity curves should also present a positive correlation [78].
The SP curve presents a first deflection at ~84-96 m (Figure 4), within the sector identified as massive basalt.There is no clear correlation to establish whether this interval qualifies as permeable.Positive deflections were also found around 174 m.Here, there is a better correlation with resistivity curves.Hence, the SP hints at the presence of permeable rocks between 174 and 200 m.These correspond to the fractured breccia aquifer.
A porosity vs. SP plot was built to evaluate these results in a more systematic manner (Figure 5b).In general terms, porosity and SP correlate well.In volcanic rocks, this could be assimilated to a correlation between porosity and permeability [62], which in turn implies that SP is an indicator of permeability.
Permeability was also analyzed by means of resistivity curves.Previous studies establish that the gap between short and long resistivity curves can be used as an indicator of permeability [78].However, in high permeability rocks, short and long resistivity curves will present very similar values as a result of deep-mud filtrate invasion.This may lead to spurious conclusions.Uncertainty was evaluated by studying the relationship between SP and the difference between long and short normal resistivity curves.The analysis aims at identifying a negative correlation between the SP curve and the other two.This would suggest that SP increases as the separation between the resistivity curves decreases (as a result of mud filtration).In turn, this could lead to identifying high permeability areas.
Figure 5c presents the results.A negative correlation comes through clearly.Maximum SP values (~15-20 mV) occur when the separation between the resistivity curves is minimal (0.1-2 Ω•m).This suggests that deep-mud filtrate invasion is actually taking place as a consequence of high permeability formations.Invasion areas can be identified visually in Figure 4.These are located at 150-163 m and 174-200 m.Deep-mud filtrate invasion has been previously explored in the oil industry [79], but little research has been carried out on the groundwater front.
Similarly, separation between the resistivity curves decreases with depth, while porosity increases.Therefore, a high permeability could be attributed to both sectors, corresponding to the fractured breccia aquifer (hydrofacies (e)).The porosity of this unit (17%-26%) is consistent with reports for similar geological formations, i.e., 10%-32% in [62].
Porosity is obtained from Archie's law, which in turn relies on resistivity logs.As a result, porosity may be overestimated.Although values from geophysical logs are consistent with the lithology, at the regional level there appears to be an overestimation of the subsurface porosity model.The model established a ~6%-42% interval for pore size distribution in the subsurface (Figure 6).
Previous works point out some of the limitations of estimating porosity from Archie's law without sonic logs or lab measurements.Alao et al. [80], for instance, explain that Archie's approach was originally developed for clay-free sedimentary rocks, and it fails to consider additional petrophysical features such as the presence of shale, fresh formation waters or thin-bed layers.Not considering these can lead to overestimate porosity and saturation.In turn, Markov et al. [81] reported that in carbonate rocks with a percolation threshold (connectivity and conductivity loss as the rock reaches critical void space) and chaotic microstructure (various pore systems of different types and scales), porosity cannot be estimated within the framework of Archie's Law.
Finally, Wright et al. [62] argue that permeability, porosity and electrical conductivity may vary by several orders of magnitude within extrusive rocks.This is due to variable tortuosity, vesicle sizes, pore elongation and connected pathways during bubble expansion, eruption conditions and crystallization history.These authors conclude that Archie's two-constant scheme (m and a) could be inappropriate under these circumstances.All this means that the porosity model in this paper should be used with caution.Nevertheless, it is advocated as a useful tool to obtain first-approach estimates of porosity and further groundwater reserves.paper should be used with caution.Nevertheless, it is advocated as a useful tool to obtain firstapproach estimates of porosity and further groundwater reserves.

Groundwater Salinity
Groundwater salinity for the upper permeable sector (174-200 m depth, within hydrofacies (e)) was estimated by means of petrophysical techniques [56].Well logs present a positive SP of ~20.5 mV, while groundwater temperature is ~30 °C, the resistivity of the drilling fluid is ~5 Ω•m and the resistivity of the aquifer (Rw) is 9.6 Ω•m.Based on Shlumberger's log interpretation charts [82], a salinity value (NaCl) of ~400 ppm was calculated.
Direct measures of TDS during borehole construction yield 602-614 ppm.These values are representative of hydrofacies (c) and (d).The BGS and CONAGUA [37] reported TDS concentrations in springs and dug wells in Irrigation District 100, near Ixmiquilpan (9 km to the northeast of the borehole, Figure 1c).These range between 656-898 and 973-1280 ppm.More recently, Chávez et al. [83] reported TDS values within the Tula and Mezquital valleys of 865-1586 (wet season) to 861-1488 ppm (dry season) in abstraction wells, 1061-1488 (wet season) to 1111-1513 ppm (dry season) in springs, and 1434-1702 (wet season) to 1402-1597 (dry season) ppm in dug wells, respectively.Thus, the salinity estimates obtained during fieldwork suggest that wastewater infiltration has had no major impacts on aquifer salinity.

Expected Flow Rates, Well Efficiency and Aquifer Parameters
Step-drawdown test results are shown in Figure 7 and Table 2.These correspond to hydrofacies (b) and (c), that is, to the fractured lava flows and the granular breccia aquifer (100-150 m).The borehole yielded 6-22 L/s with a small drawdown (sw = 1-4.29 m).Head losses attributable to the aquifer (BQ) account for most of the drawdown at each step (0.69 m ≤ BQ ≤ 2.22 m; avg.= 1.66 m).This means that 58%-70% of sw is due to the aquifer component (avg.= 60.9%), while ~12%-19% is due to the well loss component (avg.= 15.9%).Additional nonlinear head losses at each step amounted to 0.12, 0.23, 0.37, 0.47, 0.58, 0.69 and 0.81 m, respectively.

Groundwater Salinity
Groundwater salinity for the upper permeable sector (174-200 m depth, within hydrofacies (e)) was estimated by means of petrophysical techniques [56].Well logs present a positive SP of ~20.5 mV, while groundwater temperature is ~30 • C, the resistivity of the drilling fluid is ~5 Ω•m and the resistivity of the aquifer (R w ) is 9.6 Ω•m.Based on Shlumberger's log interpretation charts [82], a salinity value (NaCl) of ~400 ppm was calculated.
Direct measures of TDS during borehole construction yield 602-614 ppm.These values are representative of hydrofacies (c) and (d).The BGS and CONAGUA [37] reported TDS concentrations in springs and dug wells in Irrigation District 100, near Ixmiquilpan (9 km to the northeast of the borehole, Figure 1c).These range between 656-898 and 973-1280 ppm.More recently, Chávez et al. [83] reported TDS values within the Tula and Mezquital valleys of 865-1586 (wet season) to 861-1488 ppm (dry season) in abstraction wells, 1061-1488 (wet season) to 1111-1513 ppm (dry season) in springs, and 1434-1702 (wet season) to 1402-1597 (dry season) ppm in dug wells, respectively.Thus, the salinity estimates obtained during fieldwork suggest that wastewater infiltration has had no major impacts on aquifer salinity.

Expected Flow Rates, Well Efficiency and Aquifer Parameters
Step-drawdown test results are shown in Figure 7 and Table 2.These correspond to hydrofacies (b) and (c), that is, to the fractured lava flows and the granular breccia aquifer (100-150 m).The borehole yielded 6-22 L/s with a small drawdown (s w = 1-4.29 m).Head losses attributable to the aquifer (BQ) account for most of the drawdown at each step (0.69 m ≤ BQ ≤ 2.22 m; avg.= 1.66 m).This means that 58%-70% of s w is due to the aquifer component (avg.= 60.9%), while ~12%-19% is due to the well loss component (avg.= 15.9%).Additional nonlinear head losses at each step amounted to 0.12, 0.23, 0.37, 0.47, 0.58, 0.69 and 0.81 m, respectively.
As the borehole will supply downgradient towns in the short term (e.g., San Agustín Village, Figure 1c), test results appear encouraging.Well yield (534-1904 m 3 /day, or 6-22 L/s; mean = ~1280 m 3 /day, or ~15 L/s) exceeds by far San Agustin's water demand (~225 m 3 /day, or 2.6 L/s).In larger urban cores, such as Alfajayucan (~19,000 inhabitants, according to INEGI (2013, unpublished), Figure 1c), the borehole will be able to roughly assist with ~25% of the current municipal demand.Domestic demands are estimated by multiplying the number of inhabitants by a daily requirement of 250 L/person/day.This estimate is considered conservative and compliant with international standards.For instance, the WHO/UNICEF Joint Monitoring Programme, established a minimum dosage of only 20 L/person/day [84], to cover basic survival needs.This is consistent with estimates by several authors [85,86], although others suggest that the threshold should be raised to 50 [87].In turn, World Health Organization (WHO) [88] cites minimal domestic requirements in the order of 50-100 L/person/day.From the volumetric standpoint, all this suggests that the newly-identified groundwater source (e.g., new groundwater well) provides a long term solution for a substantial share of the downgradient population of the Mezquital Valley.However, it is mandatory to carry out a detailed hydrogeochemical study to further assess water quality aspects.
Water 2017, 9, 4 14 of 24 50 [87].In turn, World Health Organization (WHO) [88] cites minimal domestic requirements in the order of 50-100 L/person/day.From the volumetric standpoint, all this suggests that the newlyidentified groundwater source (e.g., new groundwater well) provides a long term solution for a substantial share of the downgradient population of the Mezquital Valley.However, it is mandatory to carry out a detailed hydrogeochemical study to further assess water quality aspects.Notes: Q: Flow rate; q: specific flow rate; BQ: theoretical drawdown; CQ P : nonlinear well loss; ∆h σ : additional drawdown in the well related to the skin factor; w eff : well efficiency; SD: standard deviation.
Regarding the hydraulic behavior of the well, the equation 1.86Q + 0.65Q 1.5 (B = 1.86 [min/m 2 ]; C = 0.65 [min 2 /m 5 ]; P = 1.5) allows for the estimation of the linear aquifer loss component, BC, and the nonlinear well loss, CQ P , for any flow rate, Q.In this case, P is the order of nonlinear well losses.Jacob [89] suggested that well loss is approximately proportional to the square of the discharge.Thus, P can be assumed to be ~2 in most situations.Rorabaugh [67] showed that P may range from 1.5 to 3 in cases such as the one at hand.Walton [90] proposed an interpretation criteria for the nonlinear well coefficient, C. Values of 1900-3800 and >3800 [s 2 /m 5 ] are indicators of mild and severe well deterioration, respectively.
A coefficient C = 2365 [s 2 /m 5 ] was obtained for the borehole, suggesting suboptimal (but acceptable) behavior.
From the hydraulic point of view, the main implication is an additional drawdown within the well.This is related to turbulent flow at the aquifer-gravel pack interface, as well as to the skin effect.The latter, σ, refers to an alteration in near-borehole permeability, resulting from well completion effects, formation erosion, acidification, well development and mud filtrate invasion [91].Renard [92] pointed out that σ is positive if the well is clogged and negative for the opposite case (stimulated well), while [93] reports that the skin factor may adopt any value, but that it rarely exceeds 20 in practice.In this case, a skin coefficient of +1.48 was obtained, leading to an additional drawdown of 0.4-1.6 m (avg.≈ 1 m).This amounts to 16% of the total drawdown.Since mud filtrate invasion was earlier identified through the use of petrophysical signatures (Figures 4 and 5c), these effects are best described as mild clogging due to incomplete well development.
Well efficiency (w eff ) was estimated by means of Equation ( 2) in the form w e f f = BQ/BQ + CQ P .Table 2 shows the relative efficiencies of the seven discharge rates.As shown in Figure 7a, w eff exceeds 85% for a rate of 534 m 3 /day (6 L/s).It decreases as flow rate increases, up to ~75% for a rate of 1904 m 3 /day (22 L/s).Although this efficiency is acceptable, hydraulic behavior and productivity can be improved by completing the well development using common techniques such as overpumping, surging, rawhiding or airlift pumping [65].More novel methods, such as hydrojetting with coiled tubing and simultaneous pumping, have also been demonstrated to be a cost-effective cleaning procedure [94].
The Theis Model [68] was used to infer aquifer properties.A transmissivity (T) of 1700 m 2 /day, a hydraulic conductivity (k) of 8-40 m/day and a storage coefficient (S) of 8.8 × 10 −3 were obtained, with a mean residual sum squares (RSS) of 0.0023 m between measured and modeled drawdowns.Hydraulic conductivity was estimated by considering a saturated thickness range between 40 and 220 m.The former is referred to as the thickness of hydrofacies (b) while the latter is referred to the average thickness of electrofacies EF-2.In addition, the storage coefficient estimate must be handled with care because it was obtained from drawdown data within the pumping well.Specific capacity is 5-6 L/s/m drawdown.Provided that the well is further developed, these values suggest that flow rates may be increased considerably in case of need.
Aquifer properties are representative of hydrofacies (b) and (c).These correspond to fractured lava flows and volcanic breccia.Previous works report T = 864-4320 m 2 /day, k = 0.035-103.7 m/day and S = 2 × 10 −4 -0.2 for basaltic rocks within the Alfajayucan aquifer [36,95], which are all consistent with the above results.These also show that hydrofacies are highly heterogeneous from hydraulic standpoint.

3D Visualization and Integration Tool for Management Purposes
Field information was compiled into a Petrel 3D model [64].Figure 8 shows several instances of the 3D visualization model.EF-2 represents a regional aquifer.As such, it may be developed at a larger scale.This aquifer is made up of three main hydrofacies, namely, a fractured lava aquifer, a granular breccia aquifer and a fractured breccia aquifer.These are partially isolated from one another due to the presence of intercalated aquitards.The spatial variation of hydrofacies cannot be inferred from a single borehole.Thus, these were only characterized around the drilling site.The model can be used to calculate the estimated reserves of the system.As it stands, the model spans an area of 301.62 km 2 .EF-2 unit presents an average thickness of 220 m (Figure 3), with a 3D methods are useful to appraise uncertainties on the physical attributes of the aquifer systems and greatly improve conceptual understanding of their behavior [29].In this case, the 3D model provides a welcome visual aid to scientists and water managers.The model is useful for geostatistical analyses.It provides a reference for the study of the spatial variability of electrofacies or the geometry of the aquifer system, as well as to estimate the expected water table depth at any point within the model domain.It is also useful to inform further fieldwork, including new extraction wells.
The model can be used to calculate the estimated reserves of the system.As it stands, the model spans an area of 301.62 km 2 .EF-2 unit presents an average thickness of 220 m (Figure 3), with a porosity range between 15% and 21%.If the 8.8 × 10 −3 storage coefficient is taken as a lower-bound estimate of drainage porosity, the aquifer could be expected to hold ~12,200 Mm 3 of groundwater, as well as to release a further ~580 Mm 3 , as a result of physical drainage of interconnected pores.This means that the system's specific yield amounts to 4% of the total storage, an estimate which is roughly consistent with figures for the whole Mezquital Valley system (7%) [96].While still rough, these figures may be taken as a baseline for detailed studies based on groundwater balance calculations [97] or numerical models [98].
Much like the oil industry, the water sector can benefit from the dynamic integration of data in 3D databases.This will contribute to a more robust understanding about aquifer systems related to unmanaged regions, thus providing a sound scientific base to underpin management decisions.porosity range between 15% and 21%.If the 8.8 × 10 −3 storage coefficient is taken as a lower-bound estimate of drainage porosity, the aquifer could be expected to hold ~12,200 Mm 3 of groundwater, as well as to release a further ~580 Mm 3 , as a result of physical drainage of interconnected pores.This means that the system's specific yield amounts to 4% of the total storage, an estimate which is roughly consistent with figures for the whole Mezquital Valley system (7%) [96].While still rough, these figures may be taken as a baseline for detailed studies based on groundwater balance calculations [97] or numerical models [98].

Groundwater Flow Evolution in the Unmanaged Region of the Mezquital Valley
Much like the oil industry, the water sector can benefit from the dynamic integration of data in 3D databases.This will contribute to a more robust understanding about aquifer systems related to unmanaged regions, thus providing a sound scientific base to underpin management decisions.There are three potential explanations for these large fluctuations.The first one has to do with the presence of spatially-isolated data points (Figure 9).Secondly, fluctuations may be attributed to the fact that, in some sectors of the Mezquital Valley, groundwater levels are controlled by abrupt topographic changes.This effect is common in shallow aquifers [99].Finally, the Mezquital and Alfajayucan units have been described as multilayer systems with mutual interconnections [36].This suggests that the original datasets [37,50,52] mix shallow and deep levels.

Groundwater Flow Evolution in the Unmanaged Region of the Mezquital Valley
The latter point is difficult to confirm.This is because there are significant uncertainties in terms of basic borehole attributes such as their depth or the location of the screen.However, field data There are three potential explanations for these large fluctuations.The first one has to do with the presence of spatially-isolated data points (Figure 9).Secondly, fluctuations may be attributed to the fact that, in some sectors of the Mezquital Valley, groundwater levels are controlled by abrupt topographic changes.This effect is common in shallow aquifers [99].Finally, the Mezquital and Alfajayucan units have been described as multilayer systems with mutual interconnections [36].This suggests that the original datasets [37,50,52] mix shallow and deep levels.
The latter point is difficult to confirm.This is because there are significant uncertainties in terms of basic borehole attributes such as their depth or the location of the screen.However, field data combined with the experience of the authors in the study area suggest that most boreholes are screened from the water table to the bottom (some appear to be screened from throughout their whole length).This means that water table readings represent the average head of all aquifers.As a result, the piezometric analysis must be handled with some care.
Overall, the 1982-1998 interval histogram presents a bimodal behavior, modes being −5 and −30 m (Figure 9a).The boxplot interquartile range (Q 3 -Q 1 ) which is a measure of statistical dispersion, provides a more representative range.This suggests that the water table fluctuated mostly between −8 m and +7 m over those 16 years.
Similar considerations apply to the other dispersive histograms (Figure 9b,c).However, the interquartile range oscillates within a considerable narrower interval, in the order of −1 to +22 for 1982-2007 and between −1 and +12 m for 1998-2007, respectively (consider that positive values refer to a rise in water table and opposite for negative values).Furthermore, the median of the temporal groundwater evolution is −4, +6 and +5 m, respectively, for the 1982-1998; 1998-2007 and 1982-2007 cases.This implies that the aquifer level has fluctuated by only ±20 cm/year, which is close enough to zero.In other words, the Mezquital Valley aquifer (MVA) is closed to steady-state conditions for the 25-year interval.
This argument is supported by independent evidence.Garfias-Quezada [100] analyzed monthly piezometric data from industrial wells located towards the SW of the MVA between 1984 and 2010.Hydrographs from five wells show that the water table depth has remained roughly constant.In some cases, these present a slightly upward trend.
For instance, one of the wells presented a water table depth of ~22 m in early 1984, and its fluctuation over 26 years was just 3% (21.9 ± 0.68 m, as referred to the median ± the standard deviation).Other wells registered similar behaviors (17.21 ± 2.29 m; 14.67 ± 0.44 m; 43.82 ± 1.12 m; 44.67 ± 1.02 m).Besides, linear trends for each hydrograph present near-zero gradients (10 −4 -10 −5 ).This clearly implies that both the shallow (15-25 m) and deep piezometric heads (40-50 m) have remained stable.This is remarkable because the MVA is a highly dynamic system in terms of the water balance.The National Water Commission [96] reported the following recharge and discharge volumes for the Mezquital Valley (2007-2012): abstractions from wells amount to ~137.7 Mm 3 /year, springs yield ~85.1 Mm 3 /year, baseflow through Tula River is estimated to be ~280 Mm 3 /year, while direct evapotranspiration from shallow groundwater is in the order of 9.8 Mm 3 /year.Conversely, induced recharge from irrigation, wastewater infiltration and network losses is about ~343.2Mm 3 /year, while rainfall infiltration amounts to ~49.5 Mm 3 /year.These figures shows an aquifer storage variation of roughly +2.4 Mm 3 /year.Given the size of the system and the magnitude of the different elements of the water budget, this figure is sufficiently close to zero.Thus, it is consistent with the outcomes of the piezometric analysis.
Whereas dropping water tables are relatively common across the world [101], long-term pseudo steady-state conditions such as these are relatively rare.Unmanaged recharge seems to regulate the system, thus preventing major changes in the aquifer balance.This provides a stark contrast with neighboring aquifer systems where there is no unmanaged recharge.For instance, the Cuautitlán-Pachuca aquifer presents a severe storage decrease over the same interval of −58 Mm 3 /year [102].This corresponds to a regional water table drop in the order of −0.5 to −2 m/year.
Based on 2007 records [52], groundwater flows from the south to the north of the Mezquital valley.Equipotential lines range from 2200 to 1980 m. a.s.l.(Figure 10b).These conditions are similar in the study area for 2012 [53], where flow gradients and dynamics are equivalent, albeit with a different range of equipotential values (1900-1720 m. a.s.l., Figure 10a).This suggests that there is hydraulic continuity between both regions.Thus, the Zhinte range (Figures 9 and 10) does not constitute a no-flow boundary.
In summary, pseudo steady-state conditions in the MVA, together with the spatial distribution of hydraulic gradients, suggests that the Alfajayucan aquifer is recharged with lateral flows from the wastewater irrigation of the Mezquital Valley.Thus, newly-developed water supplies within the study area are likely to be affected from the volumetric standpoint by upgradient treated wastewater disposal.In addition, the MVA is likely to evolve into transient conditions, a state which it has not experienced for at least 30 years.The water treatment plant is expected to greatly improve the water quality for municipal use in the MVA and surrounding areas.[52]).ID: district, AA: Alfajayucan aquifer.

Conclusions
Hydrogeological exploration and characterization is still needed in many regions across the world, particularly in low-income countries and emerging economies.Within a context of unregulated wastewater-based recharge schemes, the search for new groundwater sources is essential to underpin the present and future life of local communities.Current water supply challenges call for integrated approaches.Interdisciplinary research combining geological, hydrogeophysical, petrophysical and hydraulic techniques, adapted from oil and gas reservoir modeling frameworks, need to be used jointly to evaluate the complexities of aquifer systems and optimize the development of new sources.[52]).ID: Irrigation district, AA: Alfajayucan aquifer.

Conclusions
Hydrogeological exploration and characterization is still needed in many regions across the world, particularly in low-income countries and emerging economies.Within a context of unregulated wastewater-based recharge schemes, the search for new groundwater sources is essential to underpin Water 2017, 9, 4 20 of 25 the present and future life of local communities.Current water supply challenges call for integrated approaches.Interdisciplinary research combining geological, hydrogeophysical, petrophysical and hydraulic techniques, adapted from oil and gas reservoir modeling frameworks, need to be used jointly to evaluate the complexities of aquifer systems and optimize the development of new sources.
The proposed methodology contributes not only to inferring the hydrogeological parameters controlling aquifer dynamics, but also to increasing knowledge about the behavior of key physical properties in different lithologies.Thus, it can be exported to similar settings, so as to enhance our understanding of subsurface processes in unmanaged and managed recharge regions.
Regional and local-scale assessments are needed to adequately frame the outcomes and guarantee sustainable groundwater use in the long run.Practices such as unmanaged aquifer recharge constitute a mixed blessing.Cases like the Mezquital Valley show how it provides additional aquifer replenishment, and how this contributes to the equilibrium of groundwater systems.On the other hand, unmanaged aquifer recharge is likely to induce lasting hydrochemical changes in aquifers and ecosystems.The transition from unmanaged to managed recharge is expected to provide benefits to the Mezquital Valley inhabitants, and furthermore to generate new uncertainties in relation to aquifer dynamics and downgradient systems.

Figure 2 .
Figure 2. Geological map of the study area, adapted from [53].AA: Alfajayucan aquifer.The 3D model boundary (in yellow) is shown in Section 3.6.Figure 2. Geological map of the study area, adapted from [53].AA: Alfajayucan aquifer.The 3D model boundary (in yellow) is shown in Section 3.6.

Figure 2 .
Figure 2. Geological map of the study area, adapted from [53].AA: Alfajayucan aquifer.The 3D model boundary (in yellow) is shown in Section 3.6.Figure 2. Geological map of the study area, adapted from [53].AA: Alfajayucan aquifer.The 3D model boundary (in yellow) is shown in Section 3.6.

3. 2 . 1 .
Figure 4 shows the results of geophysical logs carried out in the borehole.All readings except GR begin at the static groundwater level.The incomplete geological record (from 66 to 200 m) was reconstructed primarily from the GR curve.GR yielded values of ~60 • -90 • API (API = American Petroleum Institute) for the upper 40 m, growing to 90 • -110 • API between 40 and 60 m.The former is controlled by the presence of radioactive minerals in andesites and dry basalts.Increase with depth can

Figure 6 .
Figure 6.Subsurface porosity model.EF: electrofacies.AA , BB and CC refer to geophysical cross-sections shown in Figure 2.

Figure 7 .
Figure 7. Step test analysis.(a) Drawdown and discharge rate vs. time; (b) well efficiency vs. discharge rate.

Figure 7 .
Figure 7. Step test analysis.(a) Drawdown and discharge rate vs. time; (b) well efficiency vs. discharge rate.

Table 1 .
Field measurements derived from the step-drawdown test.

Table 2 .
Results from the step-drawdown test analysis.

Table 2 .
Results from the step-drawdown test analysis.