Integration of Isotopic (2H and 18O) and Geophysical Applications to Define a Groundwater Conceptual Model in Semiarid Regions

One-third of the global population depends on groundwater for drinking, which is an even larger proportion for arid regions. The integration of isotopic and geophysical applications has been very useful in understanding the process of groundwater recharge. The aim of this study is to define a conceptual model that describes groundwater functions within an aquifer located in a semi-arid region by identifying recharge patterns based on the isotopic characteristics of: Rainfall, surface water, shallow and deep groundwater, and incorporating regional geophysical data. We demonstrated that rainfall was affected by sub-cloud evaporation and altitude. Shallow and deep modern groundwater samples were clustered and exhibited similar evolution from rainfall. However, different groups recharged from different precipitation sources compared to the local one. In the current study, we analyzed the isotopic evolution of deep groundwater over a 10-year period, which was mainly affected by the incorporation of different flows with different isotopic signatures and the hydrodynamics of the area. We performed two geoelectrical sections in the study area to improve the understanding of the hydrogeological setting and water movement patterns. The new conceptual model should help stakeholders in the context of water management policies for the study area.


Introduction
Arid and semiarid regions represent >30% of the global terrestrial surface area and receive <400 mm yr −1 of rainfall. In these regions, aquifers represent the principal water supply; in many cases, groundwater may be a fossil resource [1][2][3][4]. The study of water resources in arid and semiarid regions has become important due to their increasing demand, and groundwater depletion [5].
Many studies indicate that infiltration through streambeds during flood events is the main form of recharge, whereby water losses by infiltration into alluvial channels during floods initiate the recharge process into groundwater [2,6]. Similarly, recharge likely occurs in only small portions of the basins in arid and semiarid regions, such as depressions and ephemeral stream channels [7][8][9][10]. Identifying the relationship between groundwater-surface water (GW-SW) systems is crucial when developing programs and policies for managing water resources [11]. The main factors controlling groundwater recharge on catchment-scale are: (1) Basin morphology and stream channel position within a landscape; (2) hydraulic conductivity of the porous medium connecting stream channels to electromagnetic (TEM) measurements, was used to: (1) Build a conceptual model of groundwater flow patterns based on the analysis of isotopic water composition, and (2) define a geological cross-section of the study area.

Study Site
The state of Zacatecas is located at the central-northern part of Mexico; it is characterized by low rainfall rates and consequently low surface runoff. Therefore, groundwater resources are vital for the local population; supplying drinking water, an irrigation source, and sustaining the oasis ecosystems in the region. The climate in Zacatecas is mostly semiarid and half of the population lives in rural areas. The study area mostly encompasses three administrative aquifers: Calera, Benito Juarez, and Chupaderos; zone topography presents low relief. Most water samples and geophysical data for this study were taken from the Calera aquifer, which is therefore most relevant. According to the National Commission of Water [36], over 130 thousand hectares of this zone are annually irrigated with 330 hm 3 . This process leads to a decline in annual groundwater levels from 1 to 5 m, largely due to rudimentary and outdated irrigation systems with <45% efficiencies. Figure 1 shows the study area where the water sampling and geophysical survey were performed. Piezometric elevations of deep groundwater with an elevation interval of 20 m were constructed using 2017 static water level data from 84 observation wells, distributed in and around the three aquifers. Important seasonal and long-term changes were documented in the Calera aquifer, water table changes from 1980 to 1994 showed~5 m depletions in areas near the aquifer boundaries. Most irrigation wells are located in central-north of the valley where major depletions of 15 m were recorded in the same period; the average depletion rate was 0.4-1.15 m yr −1 [32].
Water 2019, 11, x FOR PEER REVIEW 3 of 19 flow patterns based on the analysis of isotopic water composition, and (2) define a geological crosssection of the study area.

Study Site
The state of Zacatecas is located at the central-northern part of Mexico; it is characterized by low rainfall rates and consequently low surface runoff. Therefore, groundwater resources are vital for the local population; supplying drinking water, an irrigation source, and sustaining the oasis ecosystems in the region. The climate in Zacatecas is mostly semiarid and half of the population lives in rural areas. The study area mostly encompasses three administrative aquifers: Calera, Benito Juarez, and Chupaderos; zone topography presents low relief. Most water samples and geophysical data for this study were taken from the Calera aquifer, which is therefore most relevant. According to the National Commission of Water [36], over 130 thousand hectares of this zone are annually irrigated with 330 hm³. This process leads to a decline in annual groundwater levels from 1 to 5 m, largely due to rudimentary and outdated irrigation systems with <45% efficiencies. Figure 1 shows the study area where the water sampling and geophysical survey were performed. Piezometric elevations of deep groundwater with an elevation interval of 20 m were constructed using 2017 static water level data from 84 observation wells, distributed in and around the three aquifers. Important seasonal and longterm changes were documented in the Calera aquifer, water table changes from 1980 to 1994 showed ~5 m depletions in areas near the aquifer boundaries. Most irrigation wells are located in centralnorth of the valley where major depletions of 15 m were recorded in the same period; the average depletion rate was 0.4-1.15 m yr −1 [32].

Hydrogeological Setting
The Calera horst-graben structure developed over a long stage during the Late Tertiary and Early Quaternary periods; normal faults had three main trends: (i) North-South; (ii) Northwest-Southeast; and (iii) Northeast-Southwest. These trends uplifted peripheral Tertiary volcanic

Hydrogeological Setting
The Calera horst-graben structure developed over a long stage during the Late Tertiary and Early Quaternary periods; normal faults had three main trends: (i) North-South; (ii) Northwest-Southeast; and (iii) Northeast-Southwest. These trends uplifted peripheral Tertiary volcanic mountains, central area settlements, and the deposition of basin fill sediments (alluvial material interbedded with tuffs), with a maximum thickness of 400 m. From the Late Tertiary period, mountains underwent a rapid uplift, the Tertiary Fractured Volcanic Unit was represented by the following: (i) Rhyolite lava flows with porphyritic texture with quartz, sanidine, and plagioclase in a glassy matrix, biotite and clay as an accessory and secondary mineral; and (ii) tuffs and ignimbrites with a felsic nature [31]. The lithology of the aquifer rocks is polymictic conglomerates from the Quaternary period merged to igneous and metamorphic fractured rocks from the Triassic and Cretaceous periods by tectonic movements [33]. The Calera Aquifer is considered an unconfined aquifer in a granular medium, filling a regional graben, with~2087 km 2 horizontal extent. Saturated thickness estimation of the Calera Aquifer indicated a maximum saturated thickness of 38-570 m from the north to central area, respectively, corresponding to the aquifer's sedimentary fill. Hydraulic conductivity varied from 10 −8 to 10 −5 m s −1 , specific yield from 0.01 to 0.3, and the aquifer's reported average representative specific yield was 0.13 [37]. The abstraction wells within the area (~2035) were identified mainly tapping the Basin Fill Sediments; approximately 67% of the extracted groundwater was used for irrigation and 32% for human consumption. The aquifer basement varied in elevation from~1620 m above sea level (m.a.s.l.) in the South-Central area to~2125 m.a.s.l. in the West-Central area. Three flow patterns were reported in the aquifer: (1) Moving South to North, (2) crossing from West to East, and (3) traveling Northwest to Northeast (possibly regional) [31]. Figure 2 shows the geological map of the study area. mountains, central area settlements, and the deposition of basin fill sediments (alluvial material interbedded with tuffs), with a maximum thickness of 400 m. From the Late Tertiary period, mountains underwent a rapid uplift, the Tertiary Fractured Volcanic Unit was represented by the following: (i) Rhyolite lava flows with porphyritic texture with quartz, sanidine, and plagioclase in a glassy matrix, biotite and clay as an accessory and secondary mineral; and (ii) tuffs and ignimbrites with a felsic nature [31]. The lithology of the aquifer rocks is polymictic conglomerates from the Quaternary period merged to igneous and metamorphic fractured rocks from the Triassic and Cretaceous periods by tectonic movements [33]. The Calera Aquifer is considered an unconfined aquifer in a granular medium, filling a regional graben, with ~2087 km² horizontal extent. Saturated thickness estimation of the Calera Aquifer indicated a maximum saturated thickness of 38-570 m from the north to central area, respectively, corresponding to the aquifer's sedimentary fill. Hydraulic conductivity varied from 10 −8 to 10 −5 m s −1 , specific yield from 0.01 to 0.3, and the aquifer's reported average representative specific yield was 0.13 [37]. The abstraction wells within the area (~2035) were identified mainly tapping the Basin Fill Sediments; approximately 67% of the extracted groundwater was used for irrigation and 32% for human consumption. The aquifer basement varied in elevation from ~1620 m above sea level (m.a.s.l.) in the South-Central area to ~2125 m.a.s.l. in the West-Central area. Three flow patterns were reported in the aquifer: (1) Moving South to North, (2) crossing from West to East, and (3) traveling Northwest to Northeast (possibly regional) [31]. Figure 2 shows the geological map of the study area.

Sampling of Stable Isotopes and Collected Data
In this study, we followed the same criteria used by Gonzalez-Trinidad et al. [7] to define sampling locations, which was in accordance with the IAEA (International Atomic Energy Agency).

Sampling of Stable Isotopes and Collected Data
In this study, we followed the same criteria used by Gonzalez-Trinidad et al. [7] to define sampling locations, which was in accordance with the IAEA (International Atomic Energy Agency). A hundred and forty-four water samples were taken from Gonzalez-Trinidad et al. [7]: 115 groundwater and 29 rainfall samples. We further added 123 newly collected samples: 14 rainfall (the total amount of rainfall in the 4-months rainy season, from June to August, in the same locations of Gonzalez-Trinidad et al. [7]), 37 surface water bodies, 3 stream runoff, and 69 shallow groundwater samples, which were mainly collected from the Chupaderos aquifer for the same season in 2017. The water's isotopic composition was measured using a GV-Isoprime isotope-ratio mass spectrometer, with ±0.33 and ±1.78 precision for δ 18 O and δ 2 H, respectively. 18 O and 2 H were expressed in δ notation relative to the Vienna Standard Mean Ocean Water (VSMOW). In addition, 35 groundwater samples collected by Navarro-Velasco [35] were considered. Groundwater was sampled through wells at 65-125 m below the ground's surface (depth) in all sampling periods, whereas shallow groundwater was sampled through shallow wells at 3-10 m below the surface. This shallow water corresponds to perched aquifers settled in the area's vadose zone above the regional water table shown in Figure 1. All in all, we used a total of 302 water samples and their locations are shown in Figure 3. All 43 rainwater isotope data were used to determine the improved local meteoric water line. Data are presented in the Supplementary File. A hundred and forty-four water samples were taken from Gonzalez-Trinidad et al. [7]: 115 groundwater and 29 rainfall samples. We further added 123 newly collected samples: 14 rainfall (the total amount of rainfall in the 4-months rainy season, from June to August, in the same locations of Gonzalez-Trinidad et al. [7]), 37 surface water bodies, 3 stream runoff, and 69 shallow groundwater samples, which were mainly collected from the Chupaderos aquifer for the same season in 2017. The water's isotopic composition was measured using a GV-Isoprime isotope-ratio mass spectrometer, with ±0.33‰ and ±1.78‰ precision for δ 18 O and δ 2 H, respectively. 18 O and 2 H were expressed in δ notation relative to the Vienna Standard Mean Ocean Water (VSMOW). In addition, 35 groundwater samples collected by Navarro-Velasco [35] were considered. Groundwater was sampled through wells at 65-125 m below the ground's surface (depth) in all sampling periods, whereas shallow groundwater was sampled through shallow wells at 3-10 m below the surface. This shallow water corresponds to perched aquifers settled in the area's vadose zone above the regional water table shown in Figure 1. All in all, we used a total of 302 water samples and their locations are shown in Figure 3. All 43 rainwater isotope data were used to determine the improved local meteoric water line. Data are presented in the Supplementary File.

Bivariate Data Analysis
The BiDASys (bivariate data analysis system) software [38] was used to apply ordinary and uncertainty weighted least-squares linear regression models (OLR and UWLR). Rosales Rivera et al. [38] and Verma [39] reported that the UWLR should be considered as a better alternative because the use of uncertainty has a probability connotation with a strict confidence level of 99%, which was used in the present study. The OLR and UWLR were used to estimate the new LMWL, and several regressions with the different sampled sources of water.

Bivariate Data Analysis
The BiDASys (bivariate data analysis system) software [38] was used to apply ordinary and uncertainty weighted least-squares linear regression models (OLR and UWLR). Rosales Rivera et al. [38] and Verma [39] reported that the UWLR should be considered as a better alternative because the use of uncertainty has a probability connotation with a strict confidence level of 99%, which was used in the present study. The OLR and UWLR were used to estimate the new LMWL, and several regressions with the different sampled sources of water.

Geophysical Data Acquisition
Aquifers are regularly characterized by low electrical resistivity, which correlates with differences between rocks (density, shape, and porosity), water content, and temperature. Geological units are also affected by differences in humidity percentages [25]; therefore, they were separated into different geoelectrical units [21,26]. The transient electromagnetic method (TEM) is sensitive to bulk resistivity (conductivity) of the studied medium. Quantitative interpretation of resistivity is based on the modified Archie's Law [40], which establishes that in saturated geologic materials, with conductive water in the pores, bulk electrical resistivity depends on the porosity and resistivity of the pore fluid.
The TEM data were measured using the TerraTEM conductivity meter of Monex Geoscope PTY LTD. TEMIX XL4 and WinGLINK programs were used for processing and inversion of TEM data. Data were inverted by applying Occam´s inversion technique developed by Constable et al. [41]. The final output was a set of cross-sections describing the electrical parameters of the subsurface medium and their corresponding geological material in the study area. In total, 54 VES (vertical electrical soundings) points were placed in two profiles crossing the Calera aquifer (  into different geoelectrical units [21,26]. The transient electromagnetic method (TEM) is sensitive to bulk resistivity (conductivity) of the studied medium. Quantitative interpretation of resistivity is based on the modified Archie's Law [40], which establishes that in saturated geologic materials, with conductive water in the pores, bulk electrical resistivity depends on the porosity and resistivity of the pore fluid.
The TEM data were measured using the TerraTEM conductivity meter of Monex Geoscope PTY LTD. TEMIX XL4 and WinGLINK programs were used for processing and inversion of TEM data. Data were inverted by applying Occam´s inversion technique developed by Constable et al. [41]. The final output was a set of cross-sections describing the electrical parameters of the subsurface medium and their corresponding geological material in the study area. In total, 54 VES (vertical electrical soundings) points were placed in two profiles crossing the Calera aquifer (

Isotopic Composition of Water
The summary of collected water samples and stable isotopes are presented in Table 1.

Isotopic Composition of Water
The summary of collected water samples and stable isotopes are presented in Table 1.  Tables 2 and 3, respectively. The grouping (for shallow and deep groundwater) and number of discordant values observed, the R, and R 2 values are also presented (Tables 2 and 3). A better fit was obtained using UWLR with lower values of uncertainties compared to those for OLR, which was consistent with reports by Rosales-Rivera et al. [38].   as reported by Ingraham et al. [4] where 5.5-6.5 rainfall slopes in arid regions were not uncommon. Similarly, lower slopes for arid environments have been reported by several authors [43][44][45]. All meteoric sampling sites were higher than 2000 m.a.s.l. Thus, the "altitude effect" [44]-the progressive depletion in heavy isotopes-could be affecting the isotopic composition of meteoric waters; however, this effect is not a direct cause of the LMWL slope.   The LMWL slope, as noted by Gonzalez-Trinidad et al. [7], was lower than the global water meteoric line, δ 2 H = 8 δ18H + 10 [42]; this difference can be attributed to several factors. Since the sampling sites were placed in arid and semiarid areas, the subcloud evaporation was an important factor affecting the stable isotopic composition of liquid precipitation falling through dry air. However, it is possible that at higher elevations, subcloud evaporation is lower than at lower altitude; as reported by Ingraham et al. [4] where 5.5-6.5 rainfall slopes in arid regions were not uncommon. Similarly, lower slopes for arid environments have been reported by several authors [43][44][45]. All meteoric sampling sites were higher than 2000 m.a.s.l. Thus, the "altitude effect" [44]-the progressive depletion in heavy isotopes-could be affecting the isotopic composition of meteoric waters; however, this effect is not a direct cause of the LMWL slope.
The sampled water bodies and streams are ephemeral. Water bodies in arid regions suffer large degrees of evaporation due to long residence times and free-water surface areas in contact with high temperatures, especially during the dry season ( Figure 6). This type of water becomes progressively enriched with high rates of continued evaporation (the estimated potential evaporation for the area is 2227 mm yr −1 ), which renders it isotopically distinct compared to surrounding meteoric water. Isotopic values for water bodies lie on a line (SWL) that probably identifies some evaporation lines (since these do not intersect, the LMWL is not considered as an evaporation line), where the slope = 5.0795 (similar to LMWL), and intercept = −22.9634 (Table 3). In fact, δ 18 O positive values were found in meteoric and surface waters, reflecting larger durations of residence. Similar positive values in water bodies have been reported in the same arid regions [46][47][48]. Stream runoff samples do not fall in the LMWL. In fact, these values lie under the SWL, but not close to positive values. The isotopic-depleted composition of these three samples may be due to their higher elevation with respect to other water bodies sampling sites [4,49].  The three ShGw groups were analyzed using a soil water approach, which provided relevant information on groundwater recharge [53]. Group 2 lay almost on the LMWL, reflecting the fact that shallow wells recharged directly from local precipitation. Figure 7 shows this group to be in the lowest area (see Figure 1), indicating rapid watershed responses to precipitation that allow rapid infiltration. This type of infiltration could also be due to microtopographic features, which allow the Analysis of water sources' isotopic compositions represents a useful tool for identifying GW-SW interactions, which has been applied in alluvial aquifers with relevant results [7,16,[50][51][52]. Shallow groundwater (ShGw) samples, LMWL, and SWL are plotted in Figure 6. ShGw was clustered and plotted following the methodology employed by Gonzalez-Trinidad et al. [7], which consisted of clustering the samples depending on their perpendicular distance from the LMWL. We also performed a hierarchical clustering analysis (HCA); however, it did not explain the recharge process in the study area. In addition, Figure 6 shows UWLR for three groups for comparative analysis against data from previously explained water lines. Figure 7 shows the spatial variation of δ 18 O and deuterium for the ShGw groups in the study area, with soils classified according to FAO (Food and Agriculture Organization).
The three ShGw groups were analyzed using a soil water approach, which provided relevant information on groundwater recharge [53]. Group 2 lay almost on the LMWL, reflecting the fact that shallow wells recharged directly from local precipitation. Figure 7 shows this group to be in the lowest area (see Figure 1), indicating rapid watershed responses to precipitation that allow rapid infiltration. This type of infiltration could also be due to microtopographic features, which allow the rapid precipitation movement downward into deep soil layers, thereby greatly reducing evaporation [53]. Group 2 settled in fluvisols, which allowed water flux depending on its volume, porosity, as well as sand and clay content [54]. Group 3 lay under LMWL and SWL; with a slope = 4.1607 and an intercept = −36.3241. The stable composition of soil water in arid regions is a function of the isotopic composition of rainfall and evaporation losses. Evaporation flux was predominant in Group 3, and the resultant isotopic signature was dependent on the relative magnitude of rainfall or evaporation contributions, which could be explained by the enrichment in Figure 6. Group 1 elements were widespread at both lower and higher elevations, while Groups 2 and 3 settled at lower elevations. Nevertheless, it can be assumed that there is no difference in isotopic signatures as a function of well location. To analyze the effects of surface water on groundwater, two databases previously reported for the study area were considered; the first one by Gonzalez-Trinidad et al. [7] who reported 115 samples collected in 2014-2015, and the second by Navarro-Velasco [35] involving 35 samples collected in 2007. Figure 8 shows the UWLR for the samples obtained in 2007 with respect to the LMWL, all samples were under the LMWL. Depletion of groundwater in arid and semiarid regions was expected in comparison to local precipitation, the difference in the intercept was in the order of 20 δD. This difference suggested that groundwater had undergone some evaporation before or during its underground transit, same evaporated waters were observed in the sampled surface waters and the stream water flow [7,[55][56][57]. Arid and semiarid regions are recharged during humid climate or from weather systems traveling in different atmospheric trajectories, resulting in a more depleted hydrogen isotope composition than present day values, which occurs in the study area [4].  The three ShGw groups were analyzed using a soil water approach, which provided relevant information on groundwater recharge [53]. Group 2 lay almost on the LMWL, reflecting the fact that shallow wells recharged directly from local precipitation. Figure 7 shows this group to be in the lowest area (see Figure 1), indicating rapid watershed responses to precipitation that allow rapid infiltration. This type of infiltration could also be due to microtopographic features, which allow the Groundwater Groups 1 and 2 were less depleted in δD, suggesting that this groundwater was recharged by different sampled precipitation, probably faster traveling from the ocean due to their isotopic signature. The isotopic composition of Group 3 plotted on/or close to LMWL, suggesting a meteoric origin, and implying that modern rainfall is the dominant component of that groundwater. Furthermore, the absence of evaporation during infiltration suggested that the surface recharge towards underground networks was rapid, and that the isotopic effect due to evaporation was unimportant [17,57,58]. This preferential infiltration through the unsaturated zone that reaches the groundwater level has been previously reported for arid regions [53]. Group 4 plotted similar to the 2007 samples, suggesting that groundwater had undergone some evaporation before or during its underground transit (previously discussed). The regression line's slope for Group 4 was lower than that for surface water. Therefore, pools can be considered as a potential source of groundwater. The intercept of LMWL and groundwater lines from Groups 3 and 4 identified the stable isotope composition for the water source prior to recharge. Figure 9 depicts the spatial distribution of the stable isotopic composition of groundwater sampled during 2014-2015.
20‰ δD. This difference suggested that groundwater had undergone some evaporation before or during its underground transit, same evaporated waters were observed in the sampled surface waters and the stream water flow [7,[55][56][57]. Arid and semiarid regions are recharged during humid climate or from weather systems traveling in different atmospheric trajectories, resulting in a more depleted hydrogen isotope composition than present day values, which occurs in the study area [4].
Spatial distribution of the sampling performed in 2014-2015 was practically distributed in the same area where the sampling was performed in 2007, except for those samples settled in Benito Juarez and Chupaderos aquifers (see Figure 2). Sampling performed in 2014-2015 revealed four different clusters (Figure 8), Group 4 and partially group 3 were similar to GW sampled in 2007, while Groups 1 and 2 had similar isotopic signatures ShGw Group 1. Groundwater Groups 1 and 2 were less depleted in δD, suggesting that this groundwater was recharged by different sampled precipitation, probably faster traveling from the ocean due to their isotopic signature. The isotopic composition of Group 3 plotted on/or close to LMWL, suggesting a meteoric origin, and implying that modern rainfall is the dominant component of that groundwater. Furthermore, the absence of evaporation during infiltration suggested that the surface recharge towards underground networks was rapid, and that the isotopic effect due to evaporation was unimportant [17,57,58]. This preferential infiltration through the unsaturated zone that reaches the groundwater level has been previously reported for arid regions [53]. Group 4 plotted similar to the 2007 samples, suggesting that groundwater had undergone some evaporation before or during its underground transit (previously discussed). The regression line's slope for Group 4 was lower than that for surface water. Therefore, pools can be considered as a potential source of groundwater. The intercept of LMWL and groundwater lines from Groups 3 and 4 identified the stable isotope composition for the water source prior to recharge. Figure 9 depicts the spatial distribution of the stable isotopic composition of groundwater sampled during 2014-2015. Similar to ShGw, deep groundwater does not only reflect recharge by local precipitation events, as most samples did not plot close to LMWL, implying that groundwater can be recharged from a mixture of sources [53]. An insignificant contrast was found between shallow (recent) and deep (old) groundwater, reflecting circulation and mixing of ShGw with deeper groundwater [4,59,60].
Different groups of groundwater were presented in 2014-2015 and 2007. In fact, the former exhibited isotopic signatures similar to those of ShGw; which may be due to several reasons, for example hydrodynamic processes in the vadose zone allow meteoric water evolution into ShGw and finally into deep groundwater through unsaturated zone material (water-rock interaction following the recharge). ShGw represents water sampled from wells settled in perched aquifers, this water may be affected by soil properties and hydrodynamic processes in the ground's upper layers. Vertical hydraulic conductivity is heterogeneous in the upper layer, and therefore, the infiltration rates are widely different. The water residence time within the evaporation zone in the vadose zone determines its level of exposure to evaporation impacts arising from high ambient temperatures and low humidity-typical in arid areas and, therefore, different levels of fractionation. The relative heterogeneity of groundwater between 2007 and 2014-2015 suggests that the recharge sources are not largely similar in space and/or time, reflecting mixed conditions (young and old groundwater) [61][62][63]. Deuterium excess "d" also reflects different recharge sources. In accordance with Ingraham Similar to ShGw, deep groundwater does not only reflect recharge by local precipitation events, as most samples did not plot close to LMWL, implying that groundwater can be recharged from a mixture of sources [53]. An insignificant contrast was found between shallow (recent) and deep (old) groundwater, reflecting circulation and mixing of ShGw with deeper groundwater [4,59,60].
Different groups of groundwater were presented in 2014-2015 and 2007. In fact, the former exhibited isotopic signatures similar to those of ShGw; which may be due to several reasons, for example hydrodynamic processes in the vadose zone allow meteoric water evolution into ShGw and finally into deep groundwater through unsaturated zone material (water-rock interaction following the recharge). ShGw represents water sampled from wells settled in perched aquifers, this water may be affected by soil properties and hydrodynamic processes in the ground's upper layers. Vertical hydraulic conductivity is heterogeneous in the upper layer, and therefore, the infiltration rates are widely different. The water residence time within the evaporation zone in the vadose zone determines its level of exposure to evaporation impacts arising from high ambient temperatures and low humidity-typical in arid areas and, therefore, different levels of fractionation. The relative heterogeneity of groundwater between 2007 and 2014-2015 suggests that the recharge sources are not largely similar in space and/or time, reflecting mixed conditions (young and old groundwater) [61][62][63]. Deuterium excess "d" also reflects different recharge sources. In accordance with Ingraham et al. [4], deep groundwater in arid regions is observed to have lower deuterium excess compared to LMWL, average values of groundwater sampled in 2007 and Group 1 sampled in 2007 had lower d than LMWL (see Table 4), which may be due to: (1) Evaporation occurring during rainfall through an atmosphere of lower humidity, and (2) the sampled meteoric water was initially evaporated from the ocean under more humid conditions than the present day.

Geoelectrical Sections
Four geolectrical units were identified in both profiles, West-East and North-South. These units are discussed for each geoelectrical section (Table 5). Results from the West-South profile are plotted in Figure 10 (Section 1). Section 1 resistivity values were <150 Ωm, regions with higher values were located in the section extremes, while 5-25 Ωm values were presented in the central region. This data indicated the presence of an alluvial deposit with low resistivity values limited by two mountain systems named Sierra de Fresnillo (west) and Sierra de Zacatecas (east) (Figure 3), with different geological compositions (Figure 2) compared to alluvial deposits in the middle profile. Unit U1A was defined as predominant clay material mixed with limestone. Unit 1B was the basin fill sediments (alluvial material interbedded with tuffs) with different particle sizes constituted by gravels, sands, clays, and limestones (Q(al)). Unit 2 was a mixture between unit 1B, with altered and fractured volcanic rocks from a conglomerate of tertiary volcanic mountains. Unit 3 corresponded to the tertiary fractured volcanic unit, which was represented by rhyolite lava flows with a porphyritic texture as well as quartz, sanidine, and plagioclase in a glassy matrix, biotite and clay were also present as accessory and secondary minerals.  Figure 10. Two-dimensional time-domain electromagnetic models along West-East Section 1.
Section 2 ( Figure 11) had a North-South direction, 12% of the section had values h > 150 Ωm, reaching almost 6000 Ωm. Section 2 had the highest values in its extremes, indicating boundary conditions for the study area in terms of the geological setting. Similar to Section 1, the low resistivity values were in the center region of the area, reflecting the basin fill sediments. Units 1A and 1B corresponded to the same material explained in Section 1; therefore, units U2 and U3 in Section 2 had different origins. Unit 3 in the left side of the profile corresponded to the rhyolite-acid tuff (Tom(R-Ta)), which has several rocky outcrops in that region. Lower values of resistivity in this region may correspond to the conglomerate (Ts(cg)) which is in the upper section compared to rhyolite. Unit 3 in the right side of the section was analyzed under the same approach, although it lacks outcrops, there is an important area covered by extrusive igneous rock (Tom(R-Ta)) south of Benito Juarez aquifer, a region where the section can be projected (see Figure 2).

Conceptual Model
Hydrogen and oxygen heavy isotope contents of rainwater decrease with increasing altitude [58,64,65]. The relationship between the altitude of meteoric water stations and O ratios is known as the Regional Isotopic Gradient Line (RIGL). The proposed RIGL was elaborated using the annual average isotopic signature of 12 stations at different altitudes [64], between 1660 and 2580 m.s.n.m. (see Supplementary File).
To define the local recharge of groundwater, the O ratios versus altitude were plotted together with the new proposed RIGL (Figure 12). The RIGL equation is δ O = 7.97619 ±3.69416 − Figure 10. Two-dimensional time-domain electromagnetic models along West-East Section 1. Section 2 ( Figure 11) had a North-South direction, 12% of the section had values h > 150 Ωm, reaching almost 6000 Ωm. Section 2 had the highest values in its extremes, indicating boundary conditions for the study area in terms of the geological setting. Similar to Section 1, the low resistivity values were in the center region of the area, reflecting the basin fill sediments. Units 1A and 1B corresponded to the same material explained in Section 1; therefore, units U2 and U3 in Section 2 had different origins. Unit 3 in the left side of the profile corresponded to the rhyolite-acid tuff (Tom(R-Ta)), which has several rocky outcrops in that region. Lower values of resistivity in this region may correspond to the conglomerate (Ts(cg)) which is in the upper section compared to rhyolite. Unit 3 in the right side of the section was analyzed under the same approach, although it lacks outcrops, there is an important area covered by extrusive igneous rock (Tom(R-Ta)) south of Benito Juarez aquifer, a region where the section can be projected (see Figure 2). Section 2 ( Figure 11) had a North-South direction, 12% of the section had values h > 150 Ωm, reaching almost 6000 Ωm. Section 2 had the highest values in its extremes, indicating boundary conditions for the study area in terms of the geological setting. Similar to Section 1, the low resistivity values were in the center region of the area, reflecting the basin fill sediments. Units 1A and 1B corresponded to the same material explained in Section 1; therefore, units U2 and U3 in Section 2 had different origins. Unit 3 in the left side of the profile corresponded to the rhyolite-acid tuff (Tom(R-Ta)), which has several rocky outcrops in that region. Lower values of resistivity in this region may correspond to the conglomerate (Ts(cg)) which is in the upper section compared to rhyolite. Unit 3 in the right side of the section was analyzed under the same approach, although it lacks outcrops, there is an important area covered by extrusive igneous rock (Tom(R-Ta)) south of Benito Juarez aquifer, a region where the section can be projected (see Figure 2).

Conceptual Model
Hydrogen and oxygen heavy isotope contents of rainwater decrease with increasing altitude [58,64,65]. The relationship between the altitude of meteoric water stations and O ratios is known as the Regional Isotopic Gradient Line (RIGL). The proposed RIGL was elaborated using the annual average isotopic signature of 12 stations at different altitudes [64], between 1660 and 2580 m.s.n.m. (see Supplementary File).
To define the local recharge of groundwater, the O ratios versus altitude were plotted together with the new proposed RIGL (Figure 12). The RIGL equation is δ O = 7.97619 ±3.69416 − Figure 11. Two dimensional time-domain electromagnetic models along the North-South Section 2.

Conceptual Model
Hydrogen and oxygen heavy isotope contents of rainwater decrease with increasing altitude [58,64,65]. The relationship between the altitude of meteoric water stations and 18 O ratios is known as the Regional Isotopic Gradient Line (RIGL). The proposed RIGL was elaborated using the annual average isotopic signature of 12 stations at different altitudes [64], between 1660 and 2580 m.s.n.m. (see Supplementary File).
To define the local recharge of groundwater, the 18 O ratios versus altitude were plotted together with the new proposed RIGL (Figure 12). The RIGL equation is δ 18 O = 7.97619 (±3.69416) − 0.00752 (±0.00176) * alt (R 2 = 0.6465, R = −0.80405), the calculated depletion was −0.752 /100 m; even lower values of correlation have been previously reported [58]. even lower values of correlation have been previously reported [58]. Figure 12. Estimation of the recharge elevation using δ 18 O versus elevation, the first RIGL (regional isotopic gradient line) for the region is also presented. Figure 12 shows Group 4 of groundwater sampled in 2014-2015 and groundwater sampled in 2007, which were recharged from local precipitation. According to References [57,66,67], projecting different points on the line determines the mean elevation effect of different aquifers, if we consider similar assumptions (the gradient has not changed during the course of time). The results indicated that water originates from average altitudes of 2000-2600 m.a.s.l.; foothills in the area can be found between both elevations, where notorious changes in resistivity values at the surface of both geoelectrical sections are also found, which may indicate the fractured medium.
A conceptual model for groundwater recharge was proposed based on the isotopic and geoelectrical data ( Figure 13). Geoelectrical Section 1 confirmed water table elevation, starting with ~2180 m.a.s.l. (west) and finishing at ~2000 m.a.s.l. (east), corresponding to the potentiometric lines shown in Figure 1. Section 2 also reflected the water table along the section, with water elevations between 2140 and 2020 m.a.s.l. in the south and north extremes, respectively. Section 2 was used and projected to the Benito Juarez aquifer to build the conceptual model. Groundwater flow followed the slope of the section. Infiltration associated to precipitation, soil water replenishment, and groundwater recharge are complicated processes; most of the recharge from precipitation originates at an average elevation of 2200 m, nevertheless streams and water bodies may be contributing to recharge. Figure 12. Estimation of the recharge elevation using δ 18 O versus elevation, the first RIGL (regional isotopic gradient line) for the region is also presented. Figure 12 shows Group 4 of groundwater sampled in 2014-2015 and groundwater sampled in 2007, which were recharged from local precipitation. According to References [57,66,67], projecting different points on the line determines the mean elevation effect of different aquifers, if we consider similar assumptions (the gradient has not changed during the course of time). The results indicated that water originates from average altitudes of 2000-2600 m.a.s.l.; foothills in the area can be found between both elevations, where notorious changes in resistivity values at the surface of both geoelectrical sections are also found, which may indicate the fractured medium.
A conceptual model for groundwater recharge was proposed based on the isotopic and geoelectrical data ( Figure 13). Geoelectrical Section 1 confirmed water table elevation, starting with 2180 m.a.s.l. (west) and finishing at~2000 m.a.s.l. (east), corresponding to the potentiometric lines shown in Figure 1. Section 2 also reflected the water table along the section, with water elevations between 2140 and 2020 m.a.s.l. in the south and north extremes, respectively. Section 2 was used and projected to the Benito Juarez aquifer to build the conceptual model. Groundwater flow followed the slope of the section. Infiltration associated to precipitation, soil water replenishment, and groundwater recharge are complicated processes; most of the recharge from precipitation originates at an average elevation of 2200 m, nevertheless streams and water bodies may be contributing to recharge.

Conclusions
Our study aimed to use stable isotopes and geophysical applications to define a groundwater conceptual model in a semiarid region of Mexico. To define the correlation between δ O and δ H, two linear regression models (OLR and UWLR) were used. The results showed that UWLR was the better option. Stable isotopes of precipitation were analyzed in the state of Zacatecas, Mexico and a new LMWL is proposed for the region. The new LMWL shows a lower slope compared to GMWL, mainly due to the subcloud evaporation and altitude effect, which are consistent with similar water lines reported for semiarid regions.
The isotopic compositions of shallow and deep groundwater are similar, due to several factors including the hydrodynamic process in the vadose zone allowing the evolution of meteoric water into groundwater. However, this evolution is a complex process. The ShGw isotopic signature may be affected by hydrodynamic processes of the upper ground layers, soil properties, and mainly fractionation by evaporation. Similarly, groundwater heterogeneity reflects the spatial and temporal complexity of the recharge process, leading to potential mixed conditions (young and old groundwater). Fifty-four VES points were placed in two profiles crossing the Calera aquifer to define the geological setting in the study area and build a conceptual model. The integration of both applications allows a better understanding of groundwater recharge.
Several studies do not include a correct understanding of geological setting in aquifers; however, geophysical applications could bridge this gap. In addition, environmental tracers are a novel technique that helps define groundwater recharge sites. Here, we presented a study where both applications are properly integrated for future research. It is expected that the new conceptual model will help stakeholders in defining water management policies for the study area. In order to support the proposed conceptual model, future research should include an analysis of radioactive isotopes such as radiocarbon and tritium as well as hydrogeochemical data.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, File S1: Data.

Conclusions
Our study aimed to use stable isotopes and geophysical applications to define a groundwater conceptual model in a semiarid region of Mexico. To define the correlation between δ 18 O and δ 2 H, two linear regression models (OLR and UWLR) were used. The results showed that UWLR was the better option. Stable isotopes of precipitation were analyzed in the state of Zacatecas, Mexico and a new LMWL is proposed for the region. The new LMWL shows a lower slope compared to GMWL, mainly due to the subcloud evaporation and altitude effect, which are consistent with similar water lines reported for semiarid regions.
The isotopic compositions of shallow and deep groundwater are similar, due to several factors including the hydrodynamic process in the vadose zone allowing the evolution of meteoric water into groundwater. However, this evolution is a complex process. The ShGw isotopic signature may be affected by hydrodynamic processes of the upper ground layers, soil properties, and mainly fractionation by evaporation. Similarly, groundwater heterogeneity reflects the spatial and temporal complexity of the recharge process, leading to potential mixed conditions (young and old groundwater). Fifty-four VES points were placed in two profiles crossing the Calera aquifer to define the geological setting in the study area and build a conceptual model. The integration of both applications allows a better understanding of groundwater recharge.
Several studies do not include a correct understanding of geological setting in aquifers; however, geophysical applications could bridge this gap. In addition, environmental tracers are a novel technique that helps define groundwater recharge sites. Here, we presented a study where both applications are properly integrated for future research. It is expected that the new conceptual model will help stakeholders in defining water management policies for the study area. In order to support the proposed conceptual model, future research should include an analysis of radioactive isotopes such as radiocarbon and tritium as well as hydrogeochemical data.