water

Article

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

Antonio Hernández-Espriú 1,*, Claudia Arango-Galván 2, Alfonso Reyes-Pimentel 3,
Pedro Martínez-Santos 4, Carlos Pita de la Paz 5, Sergio Macías-Medrano 1, Alberto Arias-Paz 1

and José Agustín Breña-Naranjo 6

1 Hydrogeology Group, Earth Sciences Division, Faculty of Engineering, Universidad Nacional Autónoma de
México (UNAM), C.P. 04510 Mexico City, Mexico; sergioenriquemm@gmail.com (S.M.-M.);
ariaspaz@unam.mx (A.A.-P.)

2 Department of Geomagnetism and Exploration, Institute of Geophysics, Universidad Nacional Autónoma
de México (UNAM), C.P. 04510 Mexico City, Mexico; claudiar@geofisica.unam.mx

3 Comisión Nacional de Hidrocarburos, C.P. 03700 Ciudad de México, Mexico; alfonsoreyesp@gmail.com
4 Department of Geodynamics, Faculty of Geological Sciences, Universidad Complutense de Madrid,

28040 Madrid, Spain; pemartin@ucm.es
5 GEOTEM Ingeniería, S.A. de C.V., C.P. 14640 Mexico D.F., Mexico; cpita@geotem.com.mx
6 Institute of Engineering, Universidad Nacional Autónoma de México (UNAM), C.P. 04510 Mexico City,

Mexico; JBrenaN@iingen.unam.mx
* Correspondence: ahespriu@gmail.com or ahespriu@unam.mx; Tel.: +52-55-5622-3297

Academic Editor: David Polya
Received: 4 October 2016; Accepted: 20 December 2016; Published: 24 December 2016

Abstract: 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.

Keywords: aquifer; groundwater exploration; water supply; well logs; geophysical prospection;
TDEM; pumping test; unmanaged recharge; wastewater irrigation; Mezquital Valley

1. Introduction

Aquifer exploration and water supply evaluation is a key challenge in low- and mid-income
countries [1–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

Water 2017, 9, 4; doi:10.3390/w9010004 www.mdpi.com/journal/water

http://www.mdpi.com/journal/water
http://www.mdpi.com
http://www.mdpi.com/journal/water


Water 2017, 9, 4 2 of 25

involved. This paper focuses on developing a comprehensive and transferable approach to assessing
groundwater supply sources. The methodology is illustrated through its application to a semiarid
region subject to unmanaged aquifer recharge practices: the Mezquital Valley in central Mexico
(Figure 1a). More specifically, we assess a representative sector of the so-called Alfajayucan aquifer
(green boundary in Figure 1b,c), which is located downstream to the Mezquital Valley unmanaged
region (Figure 1c).

Water 2017, 9, 4  2 of 24 

 

play an essential role in ensuring sustainability, particularly whenever unregulated recharge schemes 

are  involved.  This  paper  focuses  on  developing  a  comprehensive  and  transferable  approach  to 

assessing groundwater supply sources. The methodology is illustrated through its application to a 

semiarid  region  subject  to unmanaged aquifer  recharge practices:  the Mezquital Valley  in central 

Mexico (Figure 1a). More specifically, we assess a representative sector of the so‐called Alfajayucan 

aquifer  (green  boundary  in  Figure  1b,c), which  is  located  downstream  to  the Mezquital Valley 

unmanaged region (Figure 1c). 

 

Figure 1. Geographical setting. (a) Mexico; (b) Mexico City (MC) and hydrologic watersheds; (c) study 

area. For part (c), MC: Mexico City, ACT: Actopan, TUL: Tula, IMX: Ixmiquilpan, ALF: Alfajayucan, 

ID: Irrigation District, AA: Alfajayucan aquifer (the management boundary is shown in green line, 

parts b and c), (1) Central Emitter; (2) Western wastewater outlet; (3) Grand Canal; (4) Requena Dam; 

(5) Endhó Dam; (6) Zinthe Mountain Range; (7) Rojo Gómez Dam; (8) Tula River. Note that the 3D 

Model boundary (yellow line in parts b and c) is shown in Section 3.6. 

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, 

Figure 1. Geographical setting. (a) Mexico; (b) Mexico City (MC) and hydrologic watersheds;
(c) study area. For part (c), MC: Mexico City, ACT: Actopan, TUL: Tula, IMX: Ixmiquilpan, ALF:
Alfajayucan, ID: Irrigation District, AA: Alfajayucan aquifer (the management boundary is shown
in green line, parts b and c), (1) Central Emitter; (2) Western wastewater outlet; (3) Grand Canal;
(4) Requena Dam; (5) Endhó Dam; (6) Zinthe Mountain Range; (7) Rojo Gómez Dam; (8) Tula River.
Note that the 3D Model boundary (yellow line in parts b and c) is shown in Section 3.6.

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–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–21].

Static modeling frameworks have been successfully used to deal with a variety of practical
problems. These include geological CO2 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–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–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.

2. Materials and Methods

2.1. 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 m3/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 m3/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 m3/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.

Water 2017, 9, 4  4 of 24 

 

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. 

 

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.



Water 2017, 9, 4 5 of 25

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 Mm3/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 Mm3/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.

2.2. 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].

2.3. 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.



Water 2017, 9, 4 6 of 25

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]:

ρ = aρwφ−m (1)

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.

2.4. 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 m3/day from 0 to 360 min; (2) 817 m3/day from 360 to 720 min; (3) 1130 m3/day
from 720 to 1080 min; (4) 1332 m3/day from 1080 to 1440 min; (5) 1532 m3/day from 1440 to 1800 min;
(6) 1714 m3/day from 1800 to 2160 min; and (7) 1904 m3/day from 2160 to 2520 min.
Field measurements (elapsed time vs. drawdown vs. flow rate) are shown in Table 1.

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

534 m3/day 1 817 m3/day 1 1130 m3/day 1 1332 m3/day 1 1532 m3/day 1 1714 m3/day 1 1904 m3/day 1

Time
(min)

s
(m)

Time
(min)

s
(m)

Time
(min)

s
(m)

Time
(min)

s
(m)

Time
(min)

s
(m)

Time
(min)

s
(m)

Time
(min)

s
(m)

60 0.98 420 1.63 780 2.35 1140 2.89 1500 3.43 1860 3.75 2220 4.19
120 1.02 480 1.63 840 2.37 1200 2.87 1560 3.45 1920 3.76 2280 4.21
180 1.01 540 1.65 900 2.37 1260 2.95 1620 3.43 1980 3.76 2340 4.21
240 1 600 1.68 960 2.34 1320 2.96 1680 3.43 2040 3.76 2400 4.21
300 1.04 660 1.68 1020 2.32 1380 2.98 1740 3.44 2100 3.76 2460 4.22
360 0.99 720 1.68 1080 2.34 1440 2.98 1800 3.44 2160 3.76 2520 4.23

Notes: 1 Flow rate; s: Water table drawdown.

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):

sw = BQ + CQP (2)

where sw corresponds to the total drawdown inside the well [L], Q is well flow rate for each step
[LT−3], BQ is the theoretical drawdown component (the B coefficient is expressed in [TL−2] units),
and CQn is the nonlinear component of head losses (coefficient C is measured in [T2L−5] units).



Water 2017, 9, 4 7 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)):

sw =
Q

4πT

[
w
(

r2S
4tT

)
+ 2σ

]
+ CQP (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 [L2T−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]:

∆pσ =
141.2qBoµ

kh
σ (4)

Equation (4) was adapted to cater for the conditions of groundwater systems. Considering that the
oil-formation volumetric factor, Bo, 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:

∆hσ =
Q
T
σ (5)

where ∆hσ represents the additional drawdown caused by the skin factor [L]. Step test results were
analyzed with AQTESOLV Pro, v4.5 (Hydrosolve Inc., Reston, VA, USA) [69].

2.5. Subsurface Models

2.5.1. 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.

2.5.2. 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.

2.5.3. 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.



Water 2017, 9, 4 8 of 25

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 km2. 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).

2.6. 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].

3. Results and Discussion

3.1. 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.



Water 2017, 9, 4 9 of 25

Water 2017, 9, 4  9 of 24 

 

 

Figure 3. Subsurface hydrogeoelectrical model. EF: electrofacies. AA′, BB′ and CC′ refer to geophysical 

cross‐sections shown in Figure 2. 

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. 

3.2. Near‐Borehole Petrophysical Response 

3.2.1. 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 

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

Electrofacies 2 (EF-2) presents negative gradients and variable resistivities of 5–170, 70–90 and
20–80 Ω·m. Thickness ranges between 190 and 250 m (average = 220 m). Geoelectrical data suggests
that geological heterogeneity is controlled by changes in volcanic facies. This unit is correlated to
aquifer formations made up of Quaternary basaltic deposits, pyroclastic flows, volcanic breccia and
sediments from Mio-Pliocene formations Donguinyó and Huichapan.

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.

3.2. Near-Borehole Petrophysical Response

3.2.1. 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



Water 2017, 9, 4 10 of 25

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
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.

Water 2017, 9, 4  10 of 24 

 

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. 

A 48‐m thick basalt unit was identified between 63 and 111 m. Basalt presents an estimated silica 

content of <52%. This unit is underlain by an andesitic breccia (~48%–62% SiO2), andesitic flows (silica 

~50%–65%, depths 111.5–121 m and 130–150 m) and acidic rocks (silica ~62%–95%, depth 163 to 174 

m). This interbedded sequence can be correlated to EF‐2. Thus, the borehole log reflects the presence 

of Huichapan and Donguinyó tuffs. Aguirre‐Díaz et al. [24] report silica contents of 73.7%–74.5% and 

55.6%–64.9% in outcrops of the Huichapan ignimbrite and Donguinyó tuffs, respectively, which are 

consistent with our results for the subsurface analysis. 

 

Figure  4. Geophysical  logs  acquired during drilling,  reconstructed  lithology  and  inferred vertical 

hydrofacies. GR: Gamma Ray; API: American Petroleum Institute.   
Figure 4. Geophysical logs acquired during drilling, reconstructed lithology and inferred vertical
hydrofacies. GR: Gamma Ray; API: American Petroleum Institute.



Water 2017, 9, 4 11 of 25

A 48-m thick basalt unit was identified between 63 and 111 m. Basalt presents an estimated silica
content of <52%. This unit is underlain by an andesitic breccia (~48%–62% SiO2), andesitic flows
(silica ~50%–65%, depths 111.5–121 m and 130–150 m) and acidic rocks (silica ~62%–95%, depth 163 to
174 m). This interbedded sequence can be correlated to EF-2. Thus, the borehole log reflects the presence
of Huichapan and Donguinyó tuffs. Aguirre-Díaz et al. [24] report silica contents of 73.7%–74.5% and
55.6%–64.9% in outcrops of the Huichapan ignimbrite and Donguinyó tuffs, respectively, which are
consistent with our results for the subsurface analysis.

3.2.2. 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.

Water 2017, 9, 4  11 of 24 

 

3.2.2. 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. 

 

Figure 5. Cross‐plots and petrophysical signatures. (a) Silica content vs. Long normal resistivity; (b) 

Porosity vs. Spontaneous potential; (c) Long normal—Short normal resistivity vs. Spontaneous potential 

vs. total porosity. For parts (a,b): (1) Massive basalt; (2) fractured basalt; (3) fractured andesite; (4) 

andesitic breccia; (5) massive andesite; (6) dacite; and (7) rhyolite. 

Seven  lythofacies were  identified attending  to  these  criteria:  (1) massive basalt;  (2)  fractured 

basalt; (3) fractured andesite; (4) andesitic breccia (both granular and fractured); (5) massive andesite; 

(6)  dacite;  and  (7)  rhyolite.  Lythofacies  provide  the  basis  for  the  identification  of  the  following 

hydrofacies: 

(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%.   

(d) Aquitard, made up of andesite and acidic rocks (thickness ~11 and 5–7 m, respectively). Total 

porosity is 7%–15%, resistivity is 40–143 Ω∙m, radioactivity 45°–150°API and silica content 62%–70%. 

(e) Fractured breccia aquifer  (thickness 13–26 m). Total porosity  is 17%–26%,  resistivity  is 30–43 

Ω∙m, radioactivity 30°–50°API and silica content 50%–59%.   

The  vertical distribution  of  hydrofacies  and  lythofacies  is  shown  in  Figure  4.  Petrophysical 

features are presented in Figure 5a,b. 

Figure 5. Cross-plots and petrophysical signatures. (a) Silica content vs. Long normal resistivity;
(b) Porosity vs. Spontaneous potential; (c) Long normal—Short normal resistivity vs. Spontaneous
potential vs. total porosity. For parts (a,b): (1) Massive basalt; (2) fractured basalt; (3) fractured andesite;
(4) andesitic breccia; (5) massive andesite; (6) dacite; and (7) rhyolite.

Seven lythofacies were identified attending to these criteria: (1) massive basalt; (2) fractured
basalt; (3) fractured andesite; (4) andesitic breccia (both granular and fractured); (5) massive andesite;
(6) dacite; and (7) rhyolite. Lythofacies provide the basis for the identification of the following hydrofacies:

(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%.



Water 2017, 9, 4 12 of 25

(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).
Total porosity is 7%–15%, resistivity is 40–143 Ω·m, radioactivity 45◦–150◦API and silica content
62%–70%.

(e) Fractured breccia aquifer (thickness 13–26 m). Total porosity is 17%–26%, resistivity is 30–43 Ω·m,
radioactivity 30◦–50◦API and silica content 50%–59%.

The vertical distribution of hydrofacies and lythofacies is shown in Figure 4. Petrophysical
features are presented in Figure 5a,b.

3.3. 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



Water 2017, 9, 4 13 of 25

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 presented 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.

Water 2017, 9, 4  13 of 24 

 

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. 

 

Figure 6. Subsurface porosity model. EF: electrofacies. AA′, BB′ and CC′ refer to geophysical cross‐

sections shown in Figure 2. 

3.4. 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. 

3.5. 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. 

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 m3/day, or 6–22 L/s; mean = ~1280 

m3/day, or ~15 L/s) exceeds by far San Agustin’s water demand (~225 m3/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 

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

3.4. 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.

3.5. 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



Water 2017, 9, 4 14 of 25

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 m3/day, or 6–22 L/s; mean = ~1280 m3/day,
or ~15 L/s) exceeds by far San Agustin’s water demand (~225 m3/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 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. 

 

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

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

Step  Q (m3/Day)  sw (m)  q (L/s/m)  BQ (m) CQP (m) Δhσ (m) Aquifer Loss (%) Well Loss (%)  weff (%)

1  534  1.00  6  0.69  0.12  0.46  69.23  12.05  85.18 

2  817  1.68  6  1.06  0.23  0.71  62.98  13.55  82.29 

3  1130  2.32  6  1.46  0.37  0.98  63.10  15.97  79.80 

4  1332  3.00  5  1.72  0.47  1.16  57.49  15.79  78.45 

5  1532  3.47  5  1.98  0.58  1.33  57.14  16.83  77.25 

6  1714  3.77  5  2.22  0.69  1.49  58.86  18.34  76.24 

7  1904  4.29  5  2.46  0.81  1.66  57.46  18.87  75.28 

Min  534.00  1.00  5.11  0.69  0.12  0.46  57.14  12.05  75.28 

Max  1904.00  4.29  6.18  2.46  0.81  1.66  69.23  18.87  85.18 

Mean  1280.43  2.79  5.44  1.66  0.47  1.11  60.89  15.91  79.21 

SD  489.44  1.18  0.40  0.63  0.25  0.43  4.48  2.45  3.52 

Notes: Q: Flow  rate;  q:  specific  flow  rate; BQ:  theoretical drawdown; CQP: nonlinear well  loss;  Δhσ: 

additional drawdown in the well related to the skin factor; weff: well efficiency; SD: standard deviation. 

Regarding the hydraulic behavior of the well, the equation 1.86Q + 0.65Q1.5 (B = 1.86 [min/m2]; C 

= 0.65 [min2/m5]; P = 1.5) allows for the estimation of the linear aquifer loss component, BC, and the 

nonlinear well loss, CQP, 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 

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



Water 2017, 9, 4 15 of 25

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

Step Q (m3/Day) sw (m) q (L/s/m) BQ (m) CQP (m) ∆hσ (m) Aquifer Loss (%) Well Loss (%) weff (%)

1 534 1.00 6 0.69 0.12 0.46 69.23 12.05 85.18
2 817 1.68 6 1.06 0.23 0.71 62.98 13.55 82.29
3 1130 2.32 6 1.46 0.37 0.98 63.10 15.97 79.80
4 1332 3.00 5 1.72 0.47 1.16 57.49 15.79 78.45
5 1532 3.47 5 1.98 0.58 1.33 57.14 16.83 77.25
6 1714 3.77 5 2.22 0.69 1.49 58.86 18.34 76.24
7 1904 4.29 5 2.46 0.81 1.66 57.46 18.87 75.28

Min 534.00 1.00 5.11 0.69 0.12 0.46 57.14 12.05 75.28
Max 1904.00 4.29 6.18 2.46 0.81 1.66 69.23 18.87 85.18

Mean 1280.43 2.79 5.44 1.66 0.47 1.11 60.89 15.91 79.21
SD 489.44 1.18 0.40 0.63 0.25 0.43 4.48 2.45 3.52

Notes: Q: Flow rate; q: specific flow rate; BQ: theoretical drawdown; CQP: nonlinear well loss; ∆hσ: additional
drawdown in the well related to the skin factor; weff: well efficiency; SD: standard deviation.

Regarding the hydraulic behavior of the well, the equation 1.86Q + 0.65Q1.5 (B = 1.86 [min/m2];
C = 0.65 [min2/m5]; P = 1.5) allows for the estimation of the linear aquifer loss component, BC,
and the nonlinear well loss, CQP, 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 [s2/m5] are indicators of mild and severe well
deterioration, respectively.

A coefficient C = 2365 [s2/m5] 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 (weff) was estimated by means of Equation (2) in the form we f f = BQ/BQ + CQP.
Table 2 shows the relative efficiencies of the seven discharge rates. As shown in Figure 7a, weff exceeds
85% for a rate of 534 m3/day (6 L/s). It decreases as flow rate increases, up to ~75% for a rate
of 1904 m3/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 m2/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 m2/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.



Water 2017, 9, 4 16 of 25

3.6. 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.

Water 2017, 9, 4  16 of 24 

 

 

Figure 8. 3D visualization and integration model. (a) Vadose zone; (b) fractured lava aquifer; (c) granular 

aquifer in andesitic breccia; (d) aquitard; (e) fractured breccia aquifer. The colors in the uppermost 

layer refer to lithologies shown in Figure 2. 

The model can be used to calculate the estimated reserves of the system. As it stands, the model 

spans an area of 301.62 km2. EF‐2 unit presents an average  thickness of 220 m  (Figure 3), with a 

Figure 8. 3D visualization and integration model. (a) Vadose zone; (b) fractured lava aquifer;
(c) granular aquifer in andesitic breccia; (d) aquitard; (e) fractured breccia aquifer. The colors in
the uppermost layer refer to lithologies shown in Figure 2.



Water 2017, 9, 4 17 of 25

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 km2. 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 Mm3 of
groundwater, as well as to release a further ~580 Mm3, 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.

3.7. Groundwater Flow Evolution in the Unmanaged Region of the Mezquital Valley

Figure 9 provides an overview of the 1982–1998, 1998–2007 and 1982–2007 evolution intervals.
Data shows large dispersions. For instance, changes of −60 to +60 m, −40 to +60 m and −65 to +55 m
are observed for 1982–1998, 1998–2007 and 1982–2007, respectively.

Water 2017, 9, 4  17 of 24 

 

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 Mm3 of groundwater, as 

well as to release a further ~580 Mm3, 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. 

3.7. Groundwater Flow Evolution in the Unmanaged Region of the Mezquital Valley 

Figure 9 provides an overview of the 1982–1998, 1998–2007 and 1982–2007 evolution intervals. 

Data shows large dispersions. For instance, changes of −60 to +60 m, −40 to +60 m and −65 to +55 m 

are observed for 1982–1998, 1998–2007 and 1982–2007, respectively. 

 

Figure 9. Groundwater evolution  in  the Mezquital valley aquifer.  (a) 1982–1998;  (b) 1998–2007;  (c) 

1982–2007. Negative contours (m) refer to a decline in potentiometric levels and opposite for positive 

values. ID: Irrigation district. Red crosses in the boxplots refer to outliers. AA: Alfajayucan aquifer.   

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. 

Figure 9. Groundwater evolution in the Mezquital valley aquifer. (a) 1982–1998; (b) 1998–2007;
(c) 1982–2007. Negative contours (m) refer to a decline in potentiometric levels and opposite for positive
values. ID: Irrigation district. Red crosses in the boxplots refer to outliers. AA: Alfajayucan aquifer.

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



Water 2017, 9, 4 18 of 25

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 (Q3–Q1) 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 Mm3/year,
springs yield ~85.1 Mm3/year, baseflow through Tula River is estimated to be ~280 Mm3/year,
while direct evapotranspiration from shallow groundwater is in the order of 9.8 Mm3/year. Conversely,
induced recharge from irrigation, wastewater infiltration and network losses is about ~343.2 Mm3/year,
while rainfall infiltration amounts to ~49.5 Mm3/year. These figures shows an aquifer storage variation
of roughly +2.4 Mm3/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 Mm3/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



Water 2017, 9, 4 19 of 25

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.Water 2017, 9, 4  19 of 24 

 

 

Figure 10. Groundwater flow net in (a) the Alfajayucan aquifer for 2012 (adapted from [53]); and (b) 

the Mezquital Valley aquifer for 2007 (adapted from [52]). ID: Irrigation district, AA: Alfajayucan aquifer.   

4. 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.   

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 

Figure 10. Groundwater flow net in (a) the Alfajayucan aquifer for 2012 (adapted from [53]); and (b) the
Mezquital Valley aquifer for 2007 (adapted from [52]). ID: Irrigation district, AA: Alfajayucan aquifer.

4. 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.

Acknowledgments: The Earth Sciences Division and the Institute of Geophysics (UNAM) financially supported
this research. The authors are grateful to Schlumberger for the donation of the Petrel® Platform (provisional
academic license number PETREL_151139740_MAAAIACBH/aUA), during March-June 2015. We thank
R. Castrejón for providing us technical guidelines during well logging analysis and to Fátima Espinosa,
Alejandro Gómez, Francisco Miguel, Hermes Rochín, Alberto Garfias and Diego Ruiz, for helping out with
field work and analysis. We thank Santiago Arellano Islas (Comision Estatal del Agua y Alcantarillado de
Hidalgo) and the residents of San Agustín Village for their permission and support, respectively. The authors
thank the reviewers and editor for helpful discussions which considerably improved the quality of the manuscript.
This is the 9th contribution from the Hydrogeology Group (Earth Sciences Division), Faculty of Engineering,
UNAM, Mexico.

Author Contributions: All of the authors contributed extensively to the work. Antonio Hernández-Espriú led
the research. Claudia Arango-Galván and Carlos Pita de la Paz performed the geophysical field work and the
TDEM analysis. Alberto Arias-Paz performed the hydrogeological field work and framed the geology of the
study area. Alfonso Reyes-Pimentel developed the 3D visualization model and the petrophysical evaluation.
Sergio Macías-Medrano developed the geospatial database, GIS analysis and compiled the figures. José Agustín
Breña-Naranjo performed the statistical assessment of groundwater level records. Pedro Martínez-Santos
interpreted the 3D model and the subsurface porosity model, drafted the discussion, edited the format and
developed the English version of the Manuscript. All the authors contributed to the manuscript preparation,
analyzing, discussing and sharing valuable input.

Conflicts of Interest: The authors declare no conflict of interest.

References

1. Oikonomidis, D.; Dimogianni, S.; Kazakis, N.; Voudouris, K. A GIS/Remote Sensing-based methodology for
groundwater potentiality assessment in Tirnavos area, Greece. J. Hydrol. 2015, 525, 197–208. [CrossRef]

2. Panigrahi, B.; Nayak, A.K.; Sharma, S.D. Application of remote sensing technology for groundwater potential
evaluation. Water Resour. Manag. 1995, 9, 161–173. [CrossRef]

3. Dawoud, M.A.; Raouf, A.R.A. Groundwater exploration and assessment in rural communities of Yobe State,
Northern Nigeria. Water Resour. Manag. 2008, 23, 581–601. [CrossRef]

4. Llamas, M.R.; Martínez-Santos, P. The silent revolution of intensive ground water use: Pros and cons.
Ground Water 2005, 43, 161. [CrossRef] [PubMed]

5. Mallet, J.L. Geomodeling; Oxford University Press: New York, NY, USA, 2002.
6. Duca, C.; Abruzzi, D. Fully integrated hydrocarbon reservoir studies: Myth or reality? Am. J. Appl. Sci. 2010,

7, 1477–1486.
7. Warrlich, G.; Bosence, D.; Waltham, D.; Wood, C.; Boylan, A.; Badenas, B. 3D stratigraphic forward

modelling for analysis and prediction of carbonate platform stratigraphies in exploration and production.
Mar. Pet. Geol. 2008, 25, 35–58. [CrossRef]

8. Sacchi, Q.; Salina Borello, E.; Weltje, G.J.; Dalman, R. Increasing the predictive power of geostatistical
reservoir models by integration of geological constraints from stratigraphic forward modeling.
Mar. Pet. Geol. 2016, 69, 112–126. [CrossRef]

http://dx.doi.org/10.1016/j.jhydrol.2015.03.056
http://dx.doi.org/10.1007/BF00872127
http://dx.doi.org/10.1007/s11269-008-9289-x
http://dx.doi.org/10.1111/j.1745-6584.2005.0012.x
http://www.ncbi.nlm.nih.gov/pubmed/15819932
http://dx.doi.org/10.1016/j.marpetgeo.2007.04.005
http://dx.doi.org/10.1016/j.marpetgeo.2015.10.018


Water 2017, 9, 4 21 of 25

9. Dickey, P.A. Petroleum Development Geology; Pen Well Books: Tulsa, OK, USA, 1981.
10. Lake, L.; Carrol, B.H. Reservoir Characterization; Academic Press, Inc.: New York, NY, USA, 1988.
11. Musliu-Kehinde, O. Integrated Reservoir Characterization: A Case Study of an Onshore Reservoir in Niger

Delta Basin. Master’s Thesis, African University of Science and Technology, Abuja, Nigeria, 2011.
12. Guyoton, F.; Grasso, J.R.; Volant, P. Interrelation between induced seismic instabilities and complex geological

structure. Geophys. Res. Lett. 1999, 19, 705–708. [CrossRef]
13. Spikes, K.T. Model-based prediction of porosity and reservoir quality from P- and S-wave data.

Geophys. Res. Lett. 2003, 30, 1–4. [CrossRef]
14. Falahat, R.; Shams, A.; MacBeth, C. Towards quantitative evaluation of gas injection using time-lapse seismic

data. Geophys. Prospect. 2011, 59, 310–322. [CrossRef]
15. Arnulf, A.F.; Singh, S.C.; Pye, J.W. Seismic evidence of a complex multi-lens melt reservoir beneath the 9◦ N

Overlapping Spreading Center at the East Pacific Rise. Geophys. Res. Lett. 2014, 41, 6109–6115. [CrossRef]
16. Arogunmati, A.; Harris, J.M. Quasi-continuous reservoir monitoring with surface seismic data.

Geophys. Prospect. 2014, 62, 117–132. [CrossRef]
17. Grana, D.; Mukerji, T. Bayesian inversion of time-lapse seismic data for the estimation of static reservoir

properties and dynamic property changes. Geophys. Prospect. 2015, 63, 637–655. [CrossRef]
18. Wyllie, M.R.; Gregory, A.R.; Gardner, L.W. Elastic wave velocities in heterogeneous and porous media.

Geophysics 1956, 27, 569–589. [CrossRef]
19. Pringle, J.K.; Westerman, A.R.; Clark, J.D.; Drinkwater, N.J.; Gardiner, A.R. 3D high-resolution digital models

of outcrop analogue study sites to constrain reservoir model uncertainty: An example from Alport Castles,
Derbyshire, UK. Pet. Geosci. 2004, 10, 343–352. [CrossRef]

20. Bueno, J.F.; Drummond, R.D.; Vidal, A.C.; Sancevero, S.S. Constraining uncertainty in volumetric estimation:
A case study from Namorado Field, Brazil. J. Pet. Sci. Eng. 2011, 77, 200–208. [CrossRef]

21. Fokker, P.A.; Visser, K.; Peters, E.; Kunakbayeva, G.; Muntendam-Bos, A.G. Inversion of surface subsidence
data to quantify reservoir compartmentalization: A field study. J. Pet. Sci. Eng. 2012, 96–97, 10–21. [CrossRef]

22. Hosseini, S.A.; Lashgari, H.; Choi, J.W.; Nicot, J.P.; Lu, J.; Hovorka, S.D. Static and dynamic reservoir
modeling for geological CO2 sequestration at Cranfield, Mississippi, U.S.A. Int. J. Greenh. Gas Control 2013,
18, 449–462. [CrossRef]

23. Bunch, M. A live test of automated facies prediction at wells for CO2 storage projects. Energy Procedia 2014,
63, 3432–3446. [CrossRef]

24. Pan, F.; McPherson, B.J.; Dai, Z.; Jia, W.; Lee, S.Y.; Ampomah, W.; Viswanathan, H.; Esser, R. Uncertainty
analysis of carbon sequestration in an active CO2-EOR field. Int. J. Greenh. Gas Control 2016, 51, 18–28.
[CrossRef]

25. Noad, J. The use of field analogues in the correlation and static reservoir methodology used in the Tern Field,
Northern North Sea, UK. Mar. Pet. Geol. 2004, 21, 485–497. [CrossRef]

26. Blöcher, M.G.; Zimmermann, G.; Moeck, I.; Brandt, W.; Hassanzadegan, A.; Magri, F. 3D numerical modeling
of hydrothermal processes during the lifetime of a deep geothermal reservoir. Geofluids 2010, 10, 406–421.
[CrossRef]

27. Ross, M.; Parent, M.; Lefebvre, R. 3D geologic framework models for regional hydrogeology and land-use
management: A case study from a Quaternary basin of southwestern Quebec, Canada. Hydrogeol. J. 2005, 13,
690–707. [CrossRef]

28. Audouin, O.; Bodin, J.; Porel, G.; Bourbiaux, B. Flowpath structure in a limestone aquifer: Multi-borehole
logging investigations at the hydrogeological experimental site of Poitiers, France. Hydrogeol. J. 2008, 16,
939–950. [CrossRef]

29. Gill, B.; Cherry, D.; Adelana, M.; Chiang, X.; Mark, R. Using three-dimensional geological mapping methods
to inform sustainable groundwater development in a volcanic landscape, Victoria, Australia. Hydrogeol. J.
2011, 19, 1349–1365. [CrossRef]

30. Szabó, N.P.; Dobróka, M.; Turai, E.; Szűcs, P. Factor analysis of borehole logs for evaluating formation
shaliness: A hydrogeophysical application for groundwater studies. Hydrogeol. J. 2013, 22, 511–526.
[CrossRef]

31. Turner, R.J.; Mansour, M.M.; Dearden, R. Improved understanding of groundwater flow in complex
superficial deposits using three-dimensional geological-framework and groundwater models: An example
from Glasgow, Scotland (UK). Hydrogeol. J. 2015, 23, 493–506. [CrossRef]

http://dx.doi.org/10.1029/92GL00359
http://dx.doi.org/10.1029/2003GL018450
http://dx.doi.org/10.1111/j.1365-2478.2010.00925.x
http://dx.doi.org/10.1002/2014GL060859
http://dx.doi.org/10.1111/1365-2478.12054
http://dx.doi.org/10.1111/1365-2478.12203
http://dx.doi.org/10.1190/1.1439063
http://dx.doi.org/10.1144/1354-079303-617
http://dx.doi.org/10.1016/j.petrol.2011.03.003
http://dx.doi.org/10.1016/j.petrol.2012.06.032
http://dx.doi.org/10.1016/j.ijggc.2012.11.009
http://dx.doi.org/10.1016/j.egypro.2014.11.372
http://dx.doi.org/10.1016/j.ijggc.2016.04.010
http://dx.doi.org/10.1016/j.marpetgeo.2004.03.006
http://dx.doi.org/10.1111/j.1468-8123.2010.00284.x
http://dx.doi.org/10.1007/s10040-004-0365-x
http://dx.doi.org/10.1007/s10040-008-0275-4
http://dx.doi.org/10.1007/s10040-011-0757-7
http://dx.doi.org/10.1007/s10040-013-1067-z
http://dx.doi.org/10.1007/s10040-014-1207-0


Water 2017, 9, 4 22 of 25

32. Somaratne, N.; Mann, S. Integrated Use of Geological, Geophysical, Radiocarbon and Stable Isotopes
Data for Tracing the Conduit Flow Paths in a Small Karstic Aquifer: Poocher Swamp Freshwater Lens,
South Australia. Environ. Nat. Resour. Res. 2016, 6, 119. [CrossRef]

33. El-Kaliouby, H.; Abdalla, O. Application of time-domain electromagnetic method in mapping saltwater
intrusion of a coastal alluvial aquifer, North Oman. J. Appl. Geophys. 2015, 115, 59–64. [CrossRef]

34. Porsani, J.L.; Almeida, E.R.; Bortolozo, C.A. TDEM survey in an area of seismicity induced by water wells in
Paraná sedimentary basin, Northern São Paulo State, Brazil. J. Appl. Geophys. 2012, 82, 75–83. [CrossRef]

35. Porsani, J.L.; Bortolozo, C.A.; Almeida, E.R.; Santos Sobrinho, N.S. TDEM survey in urban environmental for
hydrogeological study at USP campus in São Paulo city, Brazil. J. Appl. Geophys. 2012, 76, 102–108. [CrossRef]

36. Comisión Nacional del Agua (CONAGUA). Actualización de la Disponibilidad de Agua en el Acuífero
Chapantongo-Alfajayucan (1309), Estado de Hidalgo [Water Budget Update of the Chapantongo-Alfajayucan Aquifer
(1309), Hidalgo State]; CONAGUA: Mexico City, Mexico, 2013.

37. British Geological Survey (BGS); Comisión Nacional del Agua (CONAGUA). Impact of Wastewater Reuse
on Groundwater in the Mezquital Valley, Hidalgo State, Mexico; BGS Final Report WC/98/42; CONAGUA:
Mexico City, Mexico, 1998.

38. Lucho-Constantino, C.A.; Prieto-García, F.; Del Razo, L.M.; Rodríguez-Vázquez, R.; Poggi-Varaldo, H.M.
Chemical fractionation of boron and heavy metals in soils irrigated with wastewater in central Mexico.
Agric. Ecosyst. Environ. 2005, 108, 57–71. [CrossRef]

39. Félix-Cañedo, T.E.; Durán-Álvarez, J.C.; Jiménez-Cisneros, B. The occurrence and distribution of a group of
organic micropollutants in Mexico City’s water sources. Sci. Total Environ. 2013, 454–455, 109–118.

40. Gibson, R.; Durán-Álvarez, J.C.; Estrada, K.L.; Chávez, A.; Jiménez Cisneros, B. Accumulation and leaching
potential of some pharmaceuticals and potential endocrine disruptors in soils irrigated with wastewater in
the Tula Valley, Mexico. Chemosphere 2010, 81, 1437–1445. [CrossRef] [PubMed]

41. Downs, T.J.; Cifuentes-García, E.; Suffet, I.M. Risk screening for exposure to groundwater pollution in
a wastewater irrigation district of the Mexico City region. Environ. Health Perspect. 1999, 107, 553–561.
[CrossRef] [PubMed]

42. Gallegos, E.; Warren, A.; Robles, E.; Campoy, E.; Calderon, A.; Sainz, M.G.; Bonilla, P.; Escolero, O. The effects
of wastewater irrigation on groundwater quality in Mexico. Water Sci. Technol. 1999, 40, 45–52. [CrossRef]

43. Comisión Nacional del Agua (CONAGUA). Planta de Tratamiento de Aguas Residuales Atotonilco [Atotonilco
Wastewater Treatment Plant]; Federal Government, CONAGUA: Mexico City, Mexico, 2013.

44. Ward, J.; Dillon, P. Robust Design of Managed Aquifer Recharge Policy in Australia; The Commonwealth Scientific
and Industrial Research Organisation (CSIRO): Canberra, Australia, 2009.

45. Lacher, L.J.; Turner, D.S.; Gungle, B.; Bushman, B.M.; Richter, H.E. Application of hydrologic tools and
monitoring to support managed aquifer recharge decision making in the upper San Pedro River, Arizona,
USA. Water 2011, 6, 3495–3527. [CrossRef]

46. Megdal, S.B.; Dillon, P. Policy and economics of managed aquifer recharge and water banking. Water 2015, 7,
592–598. [CrossRef]

47. Sinclair Knight Merz (SKM); The Commonwealth Scientific and Industrial Research Organisation (CSIRO).
Progress in Managed Aquifer Recharge in Australia; Waterlines Report No. 73; National Water Commission:
Canberra, Australia, 2012.

48. Rossetto, R.; Ansiati, A.; Barbagli, A.; Borsi, I.; Costabile, G.; Dietrich, P.; Mazzanti, G.; Picciaia, D.; Bonari, E.
Managing Induced Riverbank Filtration (IRF) at the Serchio Rivel Well Field, Tuscany, Italy (Italy); Geophysical
Research Abstracts, 16, EGU2014-7395; EGU General Assembly: Viena, Austria, 2014.

49. Kacimov, A.R.; Zlotnik, V.; Al-Maktoumi, A.; Al-Abri, R. The green-Ampt one-dimensional infiltration from
a ponded surface into a heterogeneous soil. J. Irrig. Drain. Eng. 2015, 136, 68–72. [CrossRef]

50. Geocalli. Actualización del Estudio Geohidrológico del Valle del Mezquital, Hidalgo [Update of the geohydrologic study
of the Mezquital Valley, Hidalgo]; Contract number GZA-81-68-GD; Secretaría de Agricultura and Recursos
Hidráulicos (SARH): San Luis Potosi, Mexico, 1981.

51. Aguirre-Díaz, G.J.; López-Martínez, M. Geologic evolution of the Donguinyó-Huichapan caldera complex,
central Mexican Volcanic Belt, Mexico. J. Volcanol. Geotherm. Res. 2009, 179, 133–148. [CrossRef]

52. Lesser-carrillo, L.E.; Lesser-illades, J.M.; Arellano-Islas, S.; González-Posadas, D. Balance hídrico y calidad
del agua subterránea en el acuífero del Valle del Mezquital, México central. Rev. Mex. Cienc. Geol. 2011, 28,
323–336.

http://dx.doi.org/10.5539/enrr.v6n3p119
http://dx.doi.org/10.1016/j.jappgeo.2015.02.003
http://dx.doi.org/10.1016/j.jappgeo.2012.02.005
http://dx.doi.org/10.1016/j.jappgeo.2011.10.001
http://dx.doi.org/10.1016/j.agee.2004.12.013
http://dx.doi.org/10.1016/j.chemosphere.2010.09.006
http://www.ncbi.nlm.nih.gov/pubmed/20933253
http://dx.doi.org/10.1289/ehp.99107553
http://www.ncbi.nlm.nih.gov/pubmed/10398590
http://dx.doi.org/10.1016/S0273-1223(99)00429-1
http://dx.doi.org/10.3390/w6113495
http://dx.doi.org/10.3390/w7020592
http://dx.doi.org/10.1061/(ASCE)IR.1943-4774.0000121
http://dx.doi.org/10.1016/j.jvolgeores.2008.10.013


Water 2017, 9, 4 23 of 25

53. Miguel-Cortés, F.; Rochín, G.H. Prospección Hidrogeológica y Factibilidad de Extracción de agua Subterránea
en el Poblado de San Agustín Tlalixticapa, Hidalgo [Hydrogeological prospection and groundwater
abstraction feasibility study, within the San Agustín Tlalixticapa Village, Hidalgo]. Bachelor’s Thesis,
Faculty of Engineering, Universidad Nacional Autónoma de México, México City, Mexico, 2010.

54. Jiménez, B.; Chávez, A. Quality assessment of an aquifer recharged with wastewater for its potential use as
drinking source: “El Mezquital Valley” case. Water Sci. Technol. 2004, 50, 269–273. [PubMed]

55. Constable, S.C.; Parker, R.L.; Constable, C.G. Occam’s inversion: A practical algorithm for generating smooth
models from electromagnetic sounding data. Geophysics 1987, 52, 289–300. [CrossRef]

56. Bassiouni, Z. Theory, Measurement, and Interpretation of Well Logs; Society of Petroleum Engineers (SPE):
Woodlands, TX, USA, 1994.

57. Ran, Q.; Wang, Y.; Sun, Y.; Yan, L.; Tong, M. Volcanic Gas Reservoir Characterization; Elsevier: Amsterdam,
The Netherlands, 2014.

58. Stefansson, V.; Gudlaugsson, S.T.; Gudmundsson, A. Silica content and Gamma Ray logs in volcanic rocks. In
Proceedings of the World Geothermal Congress, Kyushu-Tohoku, Japan, 28 May–10 June 2000; pp. 2893–2897.

59. Archie, G.E. The electrical resistivity logs as an aid in determing some reservoir characteristics. Trans. Am.
Inst. Mech. Eng. 1942, 146, 54–62.

60. Zhou, F.; Yao, G. Sensitivity analysis in permeability estimation using logging and injection-falloff test data
for an anthracite coalbed methane reservoir in Southeast Qinshui Basin, China. Int. J. Coal Geol. 2014, 131,
41–51. [CrossRef]

61. Rapti-Caputo, D.; Bratus, A.; Santarato, G. Strategic groundwater resources in the Tagliamento River basin
(northern Italy): Hydrogeological investigation integrated with geophysical exploration. Hydrogeol. J. 2009,
17, 1393–1409. [CrossRef]

62. Wright, H.M.N.; Cashman, K.V.; Gottesfeld, E.H.; Roberts, J.J. Pore structure of volcanic clasts: Measurements
of permeability and electrical conductivity. Earth Planet. Sci. Lett. 2009, 280, 93–104. [CrossRef]

63. Asfahani, J. Basalt Characterization by Means of Nuclear and Electrical Well Logging Techniques.
Case Study from Southern Syria. Appl. Radiat. Isot. 2011, 69, 641–647. [CrossRef] [PubMed]

64. Schlumberger. Petrel E&P Software Platform 2014; Release Notes; Schlumberger: Houston, TX, USA, 2014.
65. Driscoll, F.G. Groundwater and Wells, 2nd ed.; Johnson Screens: St. Paul, MN, USA, 1986.
66. Eden, R.N.; Hazel, C.P. Computer and Graphical Analysis of Variable Discharge Pumping Tests of Wells;

Institution of Engineers Australia: Perth, Australia, 1973; pp. 5–10.
67. Rorabaugh, M.I. Graphical and theoretical analysis of step- drawdown tests of artesian wells. Proc. Am. Soc.

Civ. Eng. 1953, 79, 362.
68. Theis, C.V. The relation between the lowering of the Piezometric surface and the rate and duration of

discharge of a well using ground-water storage. Trans. Am. Geophys. Union 1935, 16, 519. [CrossRef]
69. Duffield, G.M. AQTESOLV™ for Windows, Version 4.5; HydroSOLV, Inc.: Reston, VA, USA, 2007.
70. Agarwal, R.G.; Al-Hussainy, R.; Ramey, H.J. An Investigation of Wellbore Storage and Skin Effect in Unsteady

Liquid Flow. I: Analytical Treatment. SPEJ 1970, 10, 279–290. [CrossRef]
71. Romero, D.J.; Valko, P.P.; Economides, M.J. Optimization of the Productivity Index and the Fracture Geometry

of a Stimulated Well with Fracture Face and Choke Skins. In Proceedings of the International Symposium
and Exhibition on Formation Damage Control, Society of Petroleum Engineers, San Antonio, TX, USA,
29 September–2 October 2002.

72. Bourdet, D. Well Test Analysis: The use of advanced interpretation models. In Handbook of Petroleum
Exploration and Production; Elsevier: Amsterdam, The Netherlands, 2002.

73. Raef, A.E.; Mattern, F.; Philip, C.; Totten, M.W. 3D seismic attributes and well-log facies analysis for prospect
identification and evaluation: Interpreted palaeoshoreline implications, Weirman Field, Kansas, USA. J. Pet.
Sci. Eng. 2015, 133, 40–51. [CrossRef]

74. QGIS Development Team. QGIS Geographic Information System Developers Manual; Open Source Geospatial
Foundation Project, QGIS Version 2.12; QGIS Development Team: Lyon, France, 2015.

75. Han, T.; Best, A.I.; Sothcott, J.; North, L.J.; MacGregor, L.M. Relationships among low frequency (2 Hz)
electrical resistivity, porosity, clay content and permeability in reservoir sandstones. J. Appl. Geophys. 2015,
112, 279–289. [CrossRef]

76. Goldman, M.; Neubauer, F.M. Groundwater exploration using integrated geophysical techniques.
Surv. Geophys. 1994, 15, 331–361. [CrossRef]

http://www.ncbi.nlm.nih.gov/pubmed/15344801
http://dx.doi.org/10.1190/1.1442303
http://dx.doi.org/10.1016/j.coal.2014.05.014
http://dx.doi.org/10.1007/s10040-009-0459-6
http://dx.doi.org/10.1016/j.epsl.2009.01.023
http://dx.doi.org/10.1016/j.apradiso.2010.12.008
http://www.ncbi.nlm.nih.gov/pubmed/21211985
http://dx.doi.org/10.1029/TR016i002p00519
http://dx.doi.org/10.2118/2466-PA
http://dx.doi.org/10.1016/j.petrol.2015.04.028
http://dx.doi.org/10.1016/j.jappgeo.2014.12.006
http://dx.doi.org/10.1007/BF00665814


Water 2017, 9, 4 24 of 25

77. Stefansson, V.; Gudmundsson, A.; Emmerman, R. Gamma Ray Logging in Icelandic Rocks. The Log Analyst,
November-December; SPWLA: Houston, TX, USA, 1982; pp. 11–16.

78. Cannon, S. Petrophysics: A Practical Guide; John Wiley & Sons: Hoboken, NJ, USA, 2005.
79. Chi, S.; Torres-Verdin, C.; Wu, J.J.J.; Alpak, F.O. Assessment of Mud-Filtrate-Invasion Effects on Borehole

Acoustic Logs and Radial Profiling of Formation Elastic Properties. SPE Reserv. Eval. Eng. 2006, 9, 553–564.
[CrossRef]

80. Alao, P.A.; Ata, A.I.; Nwoke, C.E. Subsurface and Petrophysical Studies of Shaly-Sand Reservoir Targets in
Apete Field, Niger Delta. ISRN Geophys. 2013, 2013, 102450. [CrossRef]

81. Markov, M.; Mousatov, A.; Kazatchenko, E.; Markova, I. Determination of electrical conductivity
of double-porosity formations by using generalized differential effective medium approximation.
J. Appl. Geophys. 2014, 108, 104–109. [CrossRef]

82. Schulmberger. Schulmberger Log Interpretation Charts; Schulmberger: Houston, TX, USA, 2009.
83. Chávez, A.; Maya, C.; Gibson, R.; Jiménez, B. The removal of microorganisms and organic micropollutants

from wastewater during infiltration to aquifers after irrigation of farmland in the Tula Valley, Mexico.
Environ. Pollut. 2011, 159, 1354–1362. [CrossRef] [PubMed]

84. World Health Organization (WHO); The United Nations Children’s Fund (UNICEF). Global Water Supply and
Sanitation Assessment 2000 Report; The WHO and UNICEF: Washington, DC, USA, 2000.

85. Carter, R.C.; Tyrrel, S.F.; Howsam, P. The impact and sustainability of water and sanitation programmes in
developing countries. Water Environ. Manag. 1997, 13, 292–296. [CrossRef]

86. Water and Environmental Health at London and Loughborough (WELL). Guidance Manual on Water Supply
and Sanitation Programmes; Water Engineering DC: Loughborough, UK, 1998.

87. Gleick, P.H. Basic water requirements for human activities: Meeting basic needs. Water Int. 1996, 21, 83–92.
[CrossRef]

88. World Health Organization (WHO). The Human Right to Water and Sanitation; Media Brief; WHO: Geneva,
Switzerland, 2011; pp. 1–8.

89. Jacob, C.E. Drawdown test to determine effective radius of artesian well. Trans. Am. Soc. Civ. Eng. 1947, 112,
1047–1064.

90. Walton, W.C. Selected Analytical Methods for Well and Aquifer Evaluation; Bulletin 49; Illinois State Water Survey:
Urbana, IL, USA, 1962.

91. Ramey, H.J. Well-loss function and the skin effect: A review. In Recent Trends in Hydrogeology;
Narasimhan, T.N., Ed.; Special Paper 189; Geological Society of America: Boulder, CO, USA, 1982;
pp. 265–271.

92. Renard, P. Hydraulics of Well and Well Testing. In Encyclopedia of Hydrological Sciences; Anderson, M.G.,
McDonnell, J.J., Eds.; John Wiley & Sons, Ltd.: Chichester, UK, 2006.

93. Horne, R.N. Modern Well Test Analysis: A Computer Aided Approach; Petroway Inc.: Palo Alto, CA, USA, 1995.
94. Rosberg, J.E.; Bjelm, L. Well development by jetting using coiled tubing and simultaneous pumping.

Ground Water 2009, 47, 816–821. [CrossRef] [PubMed]
95. SARH. Servicios de Prospección y Levantamientos Geológicos y Geofísicos en la Zona de Alfajayucan-Chapantongo-Actopan,

[Geological and Geophysical Prospection Studies in the Alfajayucan-Chapantongo-Actopan Area]; Final Report;
SARH: San Luis Potosi, Mexico, 1981.

96. Comisión Nacional del Agua (CONAGUA). Actualización de la Disponibilidad de Agua en el Acuífero Valle del
Mezquital (1310), Estado de Hidalgo [Water Budget Update of the Mezquital Valley Aquifer (1310), Hidalgo State];
CONAGUA: Mexico City, Mexico, 2015.

97. Luijendijk, E.; Bruggeman, A. Groundwater resources in the Jabal Al Hass region, northwest Syria:
An assessment of past use and future potential. Hydrogeol. J. 2008, 16, 511–530. [CrossRef]

98. Jaworska-Szulc, B. Groundwater flow modelling of multi-aquifer systems for regional resources evaluation:
The Gdansk hydrogeological system, Poland. Hydrogeol. J. 2009, 1, 1521–1542. [CrossRef]

99. Haitjema, H.M.; Mitchell-Bruker, S. Are water tables a subdued replica of the topography? Ground Water
2005, 43, 781–786. [CrossRef] [PubMed]

http://dx.doi.org/10.2118/90159-PA
http://dx.doi.org/10.1155/2013/102450
http://dx.doi.org/10.1016/j.jappgeo.2014.07.006
http://dx.doi.org/10.1016/j.envpol.2011.01.008
http://www.ncbi.nlm.nih.gov/pubmed/21316131
http://dx.doi.org/10.1111/j.1747-6593.1999.tb01050.x
http://dx.doi.org/10.1080/02508069608686494
http://dx.doi.org/10.1111/j.1745-6584.2009.00588.x
http://www.ncbi.nlm.nih.gov/pubmed/19473264
http://dx.doi.org/10.1007/s10040-008-0282-5
http://dx.doi.org/10.1007/s10040-009-0473-8
http://dx.doi.org/10.1111/j.1745-6584.2005.00090.x
http://www.ncbi.nlm.nih.gov/pubmed/16323999


Water 2017, 9, 4 25 of 25

100. Garfias-Quezada, A. Modelación Numérica de Acuíferos en Diferencias Finitas Utilizando la Interfaz
Libre PMWIN: Fundamentos Teóricos y caso de Aplicación en el Valle del Mezquital [Finite-Difference
Groundwater Flow Modeling Using the Open Source Interface PWMIN: Theoretical Foundations and
Application Case in the Mezquital Valley Aquifer]. Bachelor’s Thesis, Faculty of Engineering, Universidad
Nacional Autónoma de México, México City, Mexico, 2015.

101. Wada, Y.; van Beek, L.P.H.; van Kempen, C.M.; Reckman, J.W.T.M.; Vasak, S.; Bierkens, M.F.P.
Global depletion of groundwater resources. Geophys. Res. Lett. 2010, 37, L20402. [CrossRef]

102. Comisión Nacional del Agua (CONAGUA). Actualización de la Disponibilidad de Agua en el Acuífero
Cuautitlán-Pachuca (1508), Estado de México [Water Budget Update of the Cuautitlán-Pachuca Aquifer (1310),
Hidalgo State]; CONAGUA: Mexico City, Mexico, 2015.

© 2016 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access
article distributed under the terms and conditions of the Creative Commons Attribution
(CC-BY) license (http://creativecommons.org/licenses/by/4.0/).

http://dx.doi.org/10.1029/2010GL044571
http://creativecommons.org/
http://creativecommons.org/licenses/by/4.0/.

	Introduction 
	Materials and Methods 
	Geological and Hydrogeological Setting 
	Hydrogeophysical Investigation 
	Borehole Drilling and Well Log Analysis 
	Well Yield and Aquifer Parameters 
	Subsurface Models 
	Hydrogeoelectrical Model 
	Total Porosity Model 
	3D Model 

	Groundwater Flow Evolution in the Mezquital Valley 

	Results and Discussion 
	Hydrogeoelectrical Model 
	Near-Borehole Petrophysical Response 
	Reconstructed Borehole Lithology 
	Hydrofacies 

	Porosity and Qualitative Permeability Estimates 
	Groundwater Salinity 
	Expected Flow Rates, Well Efficiency and Aquifer Parameters 
	3D Visualization and Integration Tool for Management Purposes 
	Groundwater Flow Evolution in the Unmanaged Region of the Mezquital Valley 

	Conclusions