Characterization of Surface Evidence of Groundwater Flow Systems in Continental Mexico

: The dynamics of the underground part of the water cycle greatly inﬂuence the features and characteristics of the Earth’s surface. Using T ó th’s theory of groundwater ﬂow systems, surface indicators in Mexico were analyzed to understand the systemic connection between groundwater and the geological framework, relief, soil, water bodies, vegetation, and climate. Recharge and discharge zones of regional groundwater ﬂow systems were identiﬁed from evidence on the ground surface. A systematic hydrogeological analysis was made of regional surface indicators, published in o ﬃ cial, freely accessible cartographic information at scales of 1:250,000 and 1:1,000,000. From this analysis, six maps of Mexico were generated, titled “Permanent water on the surface”, “Groundwater depth”, “Hydrogeological association of soils”, “Hydrogeological association of vegetation and land use”, “Hydrogeological association of topoforms”, and “Superﬁcial evidence of the presence of groundwater ﬂow systems”. Mexico’s hydrogeological features were produced. The results show that 30% of Mexico is considered to be discharge zones of groundwater ﬂow systems (regional, intermediate, and recharge). Natural recharge processes occur naturally in 57% of the country. This work is the ﬁrst holistic analysis of groundwater in Mexico carried out at a national–regional scale using only the o ﬃ cial information available to the public. These results can be used as the basis for more detailed studies on groundwater and its interaction with the environment, as well as for the development of integrative planning tools to ensure the sustainability of ecosystems and satisfy human needs. Histosol, Lixisol, Planosol, Plinthosol, Solonchak, Solonetz, and Vertisol; (ii) soils associated with recharge zones: Leptosol, Luvisol, Nitisol, and Regosol; and (iii) soils found in indistinct zones, considering their taxonomy: Acrisol, Alisol, Andosol, Arenosol, Calcisol, Cambisol, Chernozem, Durisol, Gypsisol, Kastanozem, Phaeozem, and Umbrisol. rainforest, tamaulipean spiny thicket, coastal rosettophilous thicket, mesquite, high evergreen rainforest, sandy desert vegetation; (iv) throughﬂow–recharge zones: chaparral, low open forest, cultivated forest, medium deciduous forest, mesophilic mountain forest, low deciduous forest, crasicaule thicket, microﬁlo desert thicket, rosetophilic desert thicket, subtropical thicket; (v) recharge zones: pine forest, oak forest, pine–oak forest, oak–pine forest, cedar forest, t á scate forest, oyamel forest, ayarin forest, coniferous thicket, high mountain meadow, high evergreen jungle, low evergreen forest, medium evergreen forest; and (vi) indeﬁnite: induced pastureland, natural grassland and cultivated In this study GDE


Introduction
The prevailing belief is that the water in the hydrological cycle exists in independent phases. Without a comprehensive understanding of its nature, water is simply described as changing in time, space, physical state, and physicochemical composition, without the understanding of its systemic and dynamic functioning. Of all freshwater, 1.2% is surface water (ground ice, permafrost, lakes, atmosphere, soil water, living organisms, swamps, marshes, rivers), 30.1% is under the Earth´s surface, in its underground part, and more than two-thirds is frozen in glaciers and ice caps. Of the nonfrozen part, 97% is in the underground part of the water cycle [1] and "present in the entire porous section of the Earth's upper crust, possibly penetrating by cross formational communication to depths of 15-20 km" [2]. In Tóth's theory, natural effects induced by groundwater and modified by environmental conditions are identified [11]. Some of these effects and/or results can be seen on the surface, generally through the analysis of permanent surface water, soil type, vegetation, and geomorphology. These environmental elements are known as surface indicators of groundwater flow systems [5]. Most of the studies on groundwater focus on estimating its volume rather than understanding its environmental interactions, which in turn define water quantity and quality. This is because studies on groundwater are generally carried out using surface hydrological information or, when studying groundwater, its variations are Groundwater can permanently and simultaneously produce in situ interactions with the environment and can determine the organized spatial distribution of the products of that interaction (heat, solutes, particles of dissolved material, nutrients, etc.). Groundwater flow systems are now recognized as the most important general geological agent responsible for controlling a great number and variety of natural processes and phenomena [6].
In Tóth's theory, natural effects induced by groundwater and modified by environmental conditions are identified [11]. Some of these effects and/or results can be seen on the surface, generally through the analysis of permanent surface water, soil type, vegetation, and geomorphology. These environmental elements are known as surface indicators of groundwater flow systems [5]. Most of the studies on groundwater focus on estimating its volume rather than understanding its environmental interactions, which in turn define water quantity and quality. This is because studies on groundwater are generally carried out using surface hydrological information or, when studying groundwater, its variations are generally associated with local changes, such as rainfall. In general, in México, studies of the dynamic functioning between environmental elements are not included.
Preliminary analyses of regional surface indicators based on Tóth's theory have been carried out for central Mexico by Peñuela-Arévalo and Carrillo-Rivera [29]. These works propose the determination of groundwater flow systems as a basic tool to support recommendations for adequate water and environmental management.
There is enough official and public geographical information concerning the presence of surface water, soil, vegetation, ecosystems, geomorphology, and climate distribution to generate scientific knowledge regarding the possible regional functioning of groundwater in Mexico. The objective of this study is to interpret the natural relations between geomorphology, soil, water, vegetation, and climate that exist in the environment, from a hydrogeological perspective, by applying Tóth's theory. Using these relationships as indicators of groundwater discharge conditions, it should then be possible to describe elements related to recharge and throughflow zones. This study is expected to provide a basis for better understanding of groundwater flow systems in México and may lead to more detailed studies.

Site Description
México is located in the northern hemisphere, stretching from southern USA to Guatemala and Belize. It is home to almost 70% of the biodiversity of flora and fauna on Earth. The biodiversity that occurs in Mexico is due to its intricate geological framework, which results in complex physiography and varied climate distribution [30]. The biodiversity is a product of systemic natural balance, conditioned by the climate, geological context, water, soil, and vegetation of the environment [21]. About two-thirds of Mexico is characterized as arid or semiarid, with annual precipitation of less than 500 mm, while one-third of the country has annual precipitation of over 2000 mm. The rivers constitute 633,000 km and there are 100,000 km 2 of wetlands [31].

Surface Indicators of Groundwater
As there are slight variations in the definitions of some of the indicators related to groundwater characterization, these terms are described here to assist in their understanding as related to the present investigation.

Groundwater and Perennial Surface Water
A surface water body is described as perennial when it is present throughout the year, regardless of seasonal rain or runoff water. The water that is present in low water periods is known as the baseflow. Baseflow is the proportion of stream discharge sustained by groundwater discharge [32]. A water flow that emerges directly from the ground is known as a spring and occurs where there is enough hydraulic conductivity in the rocks, the groundwater table, topoform, and soil intersect to allow this. No matter how small the spring is, perennial water bodies occur when the water flow is protected from evaporation and absorption by nearby plants. In discharge zones, processes of recharge and/or throughflow may also take place as superimposed flow systems, with their discharge occurring in other zones [4]. The principal feature indicating a discharge zone is the permanent presence of water, waterlogging, mud, and/or flooding, and a shallow static water level.

Groundwater Depth
The depth of the groundwater table is essential information for the identification of the groundwater flow system hydraulic regime. Although groundwater depth is a strong indicator, it should be analyzed along with the main hydraulic and hydrological characteristics, as well as other surface indicators [2]. In recharge zones, shallow groundwater depths are rare, while in discharge zones, groundwater is usually on or close to the surface [2].

Groundwater and Soil
Soil is the surface layer of the Earth and defined as a natural body with distinctive, repeatable, and foreseeable properties [33] as a result of the combination and evolution of the soil-forming factors. These factors were defined by Jenny [34] as original material, relief, climate, time, and biota. International systems of edaphic classification, which generally refer to a maximum depth of 2 m below the surface, establish the morphological traits and analytical properties of the edaphic profile in a specific time and space [35]. Soil moisture and temperature regimes are considered and are generally associated with general climate conditions (temperature and precipitation). Likewise, to consider the water regime of the soil, the occurrence of hydromorphic features and the water table depth are analyzed. However, the physicochemistry of the groundwater and its functioning in the soil are not included. If these last aspects were considered, it would be possible to establish the behavior of the local soil in the regional context, and to define sustainable management [36,37]. In the present study, all soil formation agents are considered, as well as the way in which groundwater flows could affect its formation and development, based on Tóth's theory [4,11], which considers soil as a component of local ecosystems, linked to groundwater flow systems [37][38][39][40]. The soil types described in the maps available were analyzed with reference soil groups and principal and supplementary qualifiers [41]. The aspects in the soil that define the relationship with groundwater are taxonomic classification, analytical (pH, salinity, alkalinity, organic matter, base saturation, presence of limestone, etc.), and morphological properties (texture, structure, color, consistency), with special emphasis on the presence or absence of hydromorphic features (redoximorphic, such as colors, spots, Fe and Mn concretions, and water table depth). These features indicate that Fe-Mn oxidation-reduction processes have existed and/or exist (when they are present), depending on the hydric state of the soil (aeration and anoxia).

Groundwater and Vegetation
Worldwide, plants have been classified into xerophytes and phreatophytes. In between these two groups are plants that usually do not directly depend on groundwater, except under certain spatial and temporal conditions. Therefore, they do not easily fit into these classifications, since they only obtain groundwater under very specific conditions, or as a mixture of groundwater and rainwater [42][43][44]. Considering that local, intermediate, and regional flow systems can coexist at the same site, the vegetation can alternate their dependence on different types of groundwater flow (local, intermediate, or regional). Therefore, there is a temporal variation in the proportional dependence on groundwater in terms of volumetric proportion and quality. There are some types of vegetation, or dominant species, which belong to ecosystems that are directly dependent on groundwater, globally recognized as groundwater dependent ecosystems (GDEs).

Groundwater and Topoforms
Some topographical features provide the energy necessary to induce the movement of water, determining the distribution of flow energy and shaping the limits of the dominant flows or "effects of major regional landforms" [2]. Tóth's theory explains that regional relief controls the pattern and direction of flow systems. Flow is generally downward in positive morphological features, while it is upward in topographic depressions. However, this is a relative phenomenon and variable upward and downward movement of groundwater flow in systems is possible [2,11]. The topographical aspects are relative to the spatial scale in which the analysis is performed and the geomorphological aspects are considered.

Groundwater and Geological Framework
The geological context controls the presence, speed, and direction of groundwater flow and its physicochemical characteristics; interaction with tectonic activity, surface water and climate are also relevant. In the geological context, a wide range of effects in the pattern and direction of groundwater movement (vertical ascending, vertical descending, horizontal, inclined) are due mainly to the permeability of the media, in combination with the position of the water table and the depth of the basement rock [2]. The same rock, or groups of rocks, may generate different types of soil and therefore different types of vegetation. Conversely, different types of rocks can produce the same type of soil and associated vegetation. This is because of the combination of geological framework, groundwater flow systems, and climate. Groundwater is a shaping agent, both underground and on the surface. However, it is incorrect to use rock cover in a regional analysis of groundwater flow systems. In this work, it was not considered acceptable to make a simplistic, linear association of surface rock type with respect to permeability and chemical features.

Groundwater and Climate
The climate directly establishes the amount of water in a region, and its distribution over time. Air temperature, precipitation, and rate of evaporation are climatic elements to be evaluated, along with other hydrogeological components. Due to the large area of this study and the cartographic scale of the climatic information, climate is used as a regional reference for the superficial indicators, e.g., to explain the continued perennial presence of water bodies in an arid context.

Discharge, Throughflow, and Recharge Zones
The journey of groundwater will end in discharge zones, where the water emerges onto or close to the surface. Discharge zones are characterized by primary attributes [4]: (i) positive potential slope, (ii) relatively low-lying position, (iii) allochthonous physicochemical composition of water, and (iv) allochthonous water temperature. One or more of these attributes may be present. The appearance and intensity of the primary attributes depends on the combination of hydrogeological conditions: air temperature, slope of the relief, physicochemical composition of the rocks, permeability of the rocks, vegetation, land uses, etc.
Natural features resulting in the discharge of groundwater, either directly onto the surface or in a diffused manner, include perennial surface water bodies, baseflow, springs, geysers, marshland, lands with quicksand, salt accumulations, landslips or landslides, shallow phreatic level, river valleys, swamp, eutrophicated water bodies, topographic depressions, excess humidity in the soil, hydrophilic vegetation, and/or halophyte vegetation. The permanent presence of water, waterlogging, mud, and/or flooding, and shallow static water level depth are the first and main indicators that a site is in a discharge zone of a groundwater flow system.
In recharge zones, the phreatic level varies throughout the year and is usually at considerable depth. Depending on the size of the study area, there may be no sign of permanent water on the surface in relation to the local or regional conditions. The size of the study area is usually meant to be related to the geological framework, i.e., depth of the basement rock, the water bearing and transmissive rocks, rainfall regime, and soil-vegetation associations [11]; if possible, the area should include discharge zones. Throughflow zones are those between the discharge and recharge zones of a groundwater flow system. In surface indicator analysis, throughflow zones are those in which there are no evident characteristics of a discharge or recharge zone.

Data Processing and Interpretation
The methodology used included data collection, definition and analysis of spatial attributes, analysis of surface indicators, and uploading of a GIS (geographical information system) database.
Data collection involved the direct request for data from the National Water Commission (CONAGUA) and National Institute of Statistics and Geography (INEGI), reviewing, and systematizing databases via Excel Office 365 and ArcMap 10.4.1 software.
To analyze the attributes of the spatial database evaluation and selection of the INEGI spatial attributes, the respective guides for cartographic interpretation and data dictionaries were used [61][62][63][64][65][66][67]. For the analysis of the edaphic properties of soils, the criteria from the Global Soil Resource Reference Base [41] and the Keys to Soil Taxonomy [35] were used.

Analysis of Surface Indicators
With an interdisciplinary, wide-view system interpretation, a systematic, nonlinear analysis of surface indicators of the effects and/or manifestations of groundwater flow systems took the following order of importance and scientific validity: (i) perennial surface water, including springs; (ii) depth of static water level; (iii) soil cover; (iv) vegetation and land use; and (v) topoforms.
For each dataset theme, based on the previously defined criteria and the characteristics of the original thematic charts, Excel matrices were made. For the systematic classification of these criteria and a groundwater flow zone classification was associated to every feature into a new attribute in each dataset. Conditional statements were programmed in Python and implemented through ArcGIS (v. 10.4.1).
After associating each of the thematic datasets, they were combined into a single feature class using the ArcGIS Identity tool (v10.4.1) which computes a geometric intersection of the input features inheriting all the attributes of those intersecting features. This process has limitations, as it combines different thematic datasets with different spatial resolution and accuracy, and should be interpreted with caution and only in a general, regional approach. From the attribute table of the feature class obtained, a comprehensive analysis of all the conditions defined for each of the surface indicators was carried out. Then, a final groundwater flow association was assigned to each feature using the same method as with the thematic datasets. Zones associated with discharge on the natural surface water and soil coverage datasets were associated with discharge; the rest of the combinations were analyzed separately in a systematic manner. In all, 2335 unique conditions were analyzed and classified, defining the groundwater flow systems zones, from which the final map "Superficial evidence of the presence of groundwater flow systems" was obtained (Figure 7). The aspects considered relevant to establish the associations between natural elements and groundwater are described as follows.

Natural Surface Water
The first surface indicator examined was the presence of natural, perennial water on the surface. From the 15 maps that contain spatial information regarding surface water, the spatial attributes that are directly associated with discharge zones are perennial freshwater bodies, perennial water currents, water body central line, springs, muddy terrains, wetlands, flood plains, natural water bodies (such as perennial lakes, lagoons and pools, marsh, swamp, humidity-type H 2 O, and perennial rivers of Horton-Strahler 5-7 order, at a scale of 1:250,000). In different official thematic charts related to surface water, the spatial attributes called saline soils and moisture agriculture are included, which, due to their definition in the respective dictionaries, are associated to discharge zones in this study [45,[48][49][50][51][52][53][54][55][56]58,59,68]. The spatial scale at which this study is developed suggests that the process of discharge is the most important of the above-mentioned spatial attributes.

Groundwater Depth
The depth at which groundwater is found was analyzed from CONAGUA and INEGI [60,69]. Boreholes and wells with static water level depth data were classified according to water depth (less or equal to 3 m, between 3 and 11 m, greater than 11m) and temperature. The depth of 3 meters was assigned because it is a value used as the lower limit in evapotranspiration processes. The range of between 2 and 11 meters is an interval in which natural temporal fluctuation is observed. The boreholes were grouped into seven categories and analyzed in relation to the topoforms associated with discharge and swamping-discharge: (1) thermal water with static water level ≤3 m in discharge topoform, (2) thermal water with static water level between >3 m and ≤11 m in discharge topoform, (3) thermal water with static water level between >3 m and ≤11 m not in discharge topoform, (4) no thermal water with static water level ≤3 m in discharge topoform, (5) no thermal water with static water level between >3 m and ≤11 m in discharge topoform, (6) no thermal water with static water level between >3 m and ≤11 m not in discharge topoform, and (7) static water level >11 m. The CONAGUA database has no temperature record, so their data are displayed as nonthermal. When the depth of the static level is greater than 11 meters, a mixture of flows of different order may be occurring, and therefore it is incorrect to assume a single thermal value as an indicator of a flow order (local or regional). The INEGI database contains only one value per well, but the average was used for the CONAGUA database records with more than one value.
Each spatial polygon that forms the edaphic map has three main attributes [64]: (i) dominant soil, covers ≥ 60% of the edaphic polygon; (ii) main qualifier of the dominant soil, indicating its quality; and (iii) secondary soil, which is estimated to cover 20-40% of the edaphic polygon. The main qualifiers associated with discharge processes were: albic, alkali, endogleyic, endopetrocalcic, endopetroduric, endoplinthic, endosalic, endosodic, epigleyc, epiplinthic, episalic, episodic, ferric, fibric, fluvic, folic, gelic, gypsic, gypsyric, gleyc, grumic, hypersalic, hypersodic, hypogleyc, hyposalic, hyposodic, mazic, mesotrophic, natric, petrosalic, plinthic, rheic, rhodic, rubic, salic, stagnic, takyric, and vertic. The rest of the qualifiers are edaphic characteristics that may occur due to both recharge and discharge processes, so they were assigned as indistinct. The analysis of the soil attributes gave nine soil classifications that show the processes of groundwater flow in each polygon of the soil cover ( Table 1). The last two classifications result from polygons with two hydrogeologically contrasting soil types. In the original soil map this polygon is recorded with ≥60% of the soil of one type and 20-40% of another soil type [64]. The category proposed as major discharge with minor recharge surface, indicates that 60% of the polygon is discharge soil, and 20-40% is a soil with conditions associated with groundwater recharge. The major recharge with minor discharge surface category is an edaphic polygon with ≥60% covered by a soil type associated with recharge zones and that 20-40% of the polygon has a typical groundwater discharge soil.

Relief Analysis
The relief of Mexico has ten regional geoforms [51], which were grouped into (i) regional discharge zones, valleys, plains, depressions, sinkholes, beach/bar, plateau, and dune fields; (ii) regional recharge zones, mountain ranges, mountains, and hills. The spatial scale of topoform mapping is 1:1,000,000, so the hydrogeological association is of a regional nature where natural discharge processes of a different order may be taking place (regional, local, intermediate).
The association of these geoforms was refined by classifying deposition material, location, and phase [51], associated with the discharge processes: flood prone, saline, barrier, coastal, deltaic, alluvial, and lacustrine deposits. In a regional context, the topoforms associated with discharge, swamping-discharge and discharge-throughflow zones, are found in the regional relief (scale of 1:1,000,000) with little or no slope. In order to ensure the degree of influence on the subsurface of the presence of shallow perennial natural water (springs, perennial water bodies, lakes, lagoons, wetlands, perennial rivers with Horton-Strahler 5-7 order, swamping areas, pond, mud, marshland, and locations with a static water level depth of ≤3 m), a buffer of 1 km was applied to those attributes when in topoforms associated with discharge. Figure 2 shows the hydrogeological surface landforms that are associated surface indicators of the groundwater flow systems analyzed in the present study. Based on the geomorphological features, 339 types of topoforms, in six classifications, were found: discharge, 11.3% (48 types); swamping discharge 3.26% (24 types); discharge-throughflow, 10.8% (46 types); throughflow, 11.4% (63 types); throughflow-recharge, 24.5% (49 types); and recharge, 38.7% (108 types). According to the definition in the data dictionary, the term swamping discharge (floodable) refers to intertidal zones or seasonally continental water inundation [51].   Hydrogeological association of regional topoforms. Figure 3 shows permanent surface water. The discharge zones of regional groundwater flow systems are shown by blue shading. Approximately 15.5% of the surface is covered by natural perennial water bodies. Nonthermal springs located in topoforms that are not associated with discharge zones indicate discharge points. They indicate discharge zones when found in discharge topoforms. The presence of thermal springs in, or near, confirms that the discharge zones have regional flow system characteristics, because of the range of water temperatures; however, isotopic stable analysis should be carried out. Where regional flow systems occur, local and intermediate flows coexist. All rivers travel through recharge, throughflow and discharge sections. In turn, rivers in discharge topoforms, show that the systemic process of discharge is dominant. Figure 4 shows the groundwater depth. The static water level depths were classified in terms of depth ranges, temperature, and their regional topographic context. There are 23,767 boreholes and wells that have static water level depth records. Of them, 53% have an average groundwater depth of >11 m; 33% of between 3 and ≤11 m; and 14% of ≤3 m. Sink holes (known as cenotes) have groundwater depths of 2-6 to 20.5 m [60,69].    Figure 4 shows the groundwater depth. The static water level depths were classified in terms of depth ranges, temperature, and their regional topographic context. There are 23,767 boreholes and wells that have static water level depth records. Of them, 53% have an average groundwater depth of >11 m; 33% of between 3 and ≤11 m; and 14% of ≤3 m. Sink holes (known as cenotes) have groundwater depths of 2-6 to 20.5 m [60,69].

Results and Discussion
Where static water depths of ≤11 m were found in topoforms classified as discharge, flooding discharge and discharge-throughflow attributes, this was considered as an indicator suggesting the location is in a discharge zone. The static depth level of >11 m in topoforms related to recharge may indicate throughflow or recharge zones; when found in discharge topoforms, they may be classified as throughflow zones. Shallow water depths, close to perennial water bodies in lowlands and plains, were also found in many instances. Data without temperature records are shown as nonthermal.  Figure 5 shows the hydrogeological association of soils. Here, the distribution of edaphic covers suggests the edaphic coverage as discharge 11.3% (discharge 10.3%, discharge-recharge 1%); dischargethroughflow 14.3%; recharge 24%; throughflow-recharge 37% (throughflow-recharge, rechargethroughflow, throughflow); recharge-discharge 2.6%; and indistinct 11%. It is interesting to see the similitude that exists, in terms of magnitude of spatial coverage, between the 11.3% of the soil surface associated with discharge processes and the 15.5% of the surface covered with perennial natural water (Figure 3), as well as the 37% of soil coverage associated with recharge processes, which is almost equal to the 38% for relief associated with recharge processes (Figure 5). Where static water depths of ≤11 m were found in topoforms classified as discharge, flooding discharge and discharge-throughflow attributes, this was considered as an indicator suggesting the location is in a discharge zone. The static depth level of >11 m in topoforms related to recharge may indicate throughflow or recharge zones; when found in discharge topoforms, they may be classified as throughflow zones. Shallow water depths, close to perennial water bodies in lowlands and plains, were also found in many instances. Data without temperature records are shown as nonthermal. Figure 5 shows the hydrogeological association of soils. Here, the distribution of edaphic covers suggests the edaphic coverage as discharge 11.3% (discharge 10.3%, discharge-recharge 1%); discharge-throughflow 14.3%; recharge 24%; throughflow-recharge 37% (throughflow-recharge, recharge-throughflow, throughflow); recharge-discharge 2.6%; and indistinct 11%. It is interesting to see the similitude that exists, in terms of magnitude of spatial coverage, between the 11.3% of the soil surface associated with discharge processes and the 15.5% of the surface covered with perennial natural water (Figure 3), as well as the 37% of soil coverage associated with recharge processes, which is almost equal to the 38% for relief associated with recharge processes ( Figure 5). The national coverage of land and vegetation uses in [45] is the reference with the oldest original conditions. At that time, urban zones covered 0.1% of Mexico, 13% was agricultural, and 7.4% induced pastureland. Agriculture and induced pastureland were classified as indistinct since they can occur in many environments. Figure 6 shows the hydrogeological association of vegetation and land use. From a groundwater flow system perspective, the vegetation and land use attributes were associated as zones of (i) discharge 3.8%; (ii) discharge-throughflow, 6.4%; (iii) throughflow, 14.5%; (iv) throughflow-recharge 31.6%; and (v) recharge 22.4%.
It was found that the relationship between the vegetation associated with discharge and throughflow zones was 10.2% in the 1980s. In discharge zones with perennial water on the surface, the relation was 15.5%, and in soils related to discharge, 11.3% (Figures 3 and 5). A correlation was carried out for recharge zones based on vegetation and soil analysis, yielding 22.4% and 24.0%, respectively.
The conceptual criteria of Tóth's theory for soils and vegetation must be considered as a first reference and a more detailed multifactorial procedure should be carried out, avoiding a linear and direct process. Furthermore, it should be remembered that groundwater flows in multidimensional ways (Figure 1). From this perspective, the vegetation types with the highest degree of certainty are those related to discharge and recharge environments, as explained previously. The national coverage of land and vegetation uses in [45] is the reference with the oldest original conditions. At that time, urban zones covered 0.1% of Mexico, 13% was agricultural, and 7.4% induced pastureland. Agriculture and induced pastureland were classified as indistinct since they can occur in many environments. Figure 6 shows the hydrogeological association of vegetation and land use. From a groundwater flow system perspective, the vegetation and land use attributes were associated as zones of (i) discharge 3.8%; (ii) discharge-throughflow, 6.4%; (iii) throughflow, 14.5%; (iv) throughflow-recharge 31.6%; and (v) recharge 22.4%.
It was found that the relationship between the vegetation associated with discharge and throughflow zones was 10.2% in the 1980s. In discharge zones with perennial water on the surface, the relation was 15.5%, and in soils related to discharge, 11.3% (Figures 3 and 5). A correlation was carried out for recharge zones based on vegetation and soil analysis, yielding 22.4% and 24.0%, respectively.
The conceptual criteria of Tóth's theory for soils and vegetation must be considered as a first reference and a more detailed multifactorial procedure should be carried out, avoiding a linear and direct process. Furthermore, it should be remembered that groundwater flows in multidimensional ways (Figure 1). From this perspective, the vegetation types with the highest degree of certainty are those related to discharge and recharge environments, as explained previously. Water 2020, 12, x FOR PEER REVIEW 13 of 20   Figures 2-6, based on the methodology described in the methods. This figure shows that in at least 21.4% of Mexico there is strong evidence of discharge from groundwater flow systems. The classification of discharge-throughflow (~9.3%) indicates characteristics, phenomena, and/or manifestations of groundwater discharge into the soil, vegetation, and/or topoforms. The recharge processes may be occurring naturally in ~23% of the territory (recharge zones); this could be higher if it is considered that ~34.7% of the territory was classified as throughflowrecharge zones. The throughflow zones are those where, at the time the spatial information was recorded, there were no clear discharge or recharge zone features. They are found in ~11.5% of México.
In the discharge zones shown in Figure 7, water is present permanently, either at the surface or at a shallow depth. However, there is no confirmation of groundwater at a point here, for two reasons, (i) the origin and temporality of the data for each of the surface indicators, as pointed out in the methodology, and (ii) due to the constant dynamics to which the groundwater flow systems are subject, mainly from human activity, directly or indirectly. Groundwater discharge and recharge are affected mainly by water extraction through boreholes and dug wells, evapotranspiration, in soil and agricultural activity, loss/alteration of natural vegetation, alteration of the natural channels of perennial runoff, loss of soil permeability, and compaction-collapse of soil.   Figures 2-6, based on the methodology described in the methods. This figure shows that in at least 21.4% of Mexico there is strong evidence of discharge from groundwater flow systems. The classification of discharge-throughflow (~9.3%) indicates characteristics, phenomena, and/or manifestations of groundwater discharge into the soil, vegetation, and/or topoforms. The recharge processes may be occurring naturally in~23% of the territory (recharge zones); this could be higher if it is considered that~34.7% of the territory was classified as throughflow-recharge zones. The throughflow zones are those where, at the time the spatial information was recorded, there were no clear discharge or recharge zone features. They are found in~11.5% of México.
In the discharge zones shown in Figure 7, water is present permanently, either at the surface or at a shallow depth. However, there is no confirmation of groundwater at a point here, for two reasons, (i) the origin and temporality of the data for each of the surface indicators, as pointed out in the methodology, and (ii) due to the constant dynamics to which the groundwater flow systems are subject, mainly from human activity, directly or indirectly. Groundwater discharge and recharge are affected mainly by water extraction through boreholes and dug wells, evapotranspiration, in soil and agricultural activity, loss/alteration of natural vegetation, alteration of the natural channels of perennial runoff, loss of soil permeability, and compaction-collapse of soil.
Understanding the links between the discharge, throughflow, and recharge zones of a groundwater flow system requires more detailed analysis. Indeed, the definition of the order of groundwater flow (i.e., local, intermediate, regional) should be investigated more fully. Using systematic analysis of the physical parameters, piezometric information (water static level) and the nature of the chemical and isotope water fingerprint, to understand recharge altitude, flow path through different rock formations, water-rock interaction, depth of circulation different groundwater flows can be separated and discharge conditions established. These characteristics might assist in relating groundwater flow with the environment and accurately define its links to ecosystems.
If the reader wishes to reproduce Understanding the links between the discharge, throughflow, and recharge zones of a groundwater flow system requires more detailed analysis. Indeed, the definition of the order of groundwater flow (i.e., local, intermediate, regional) should be investigated more fully. Using systematic analysis of the physical parameters, piezometric information (water static level) and the nature of the chemical and isotope water fingerprint, to understand recharge altitude, flow path through different rock formations, water-rock interaction, depth of circulation different groundwater flows can be separated and discharge conditions established. These characteristics might assist in relating groundwater flow with the environment and accurately define its links to ecosystems.
If the reader wishes to reproduce Figures 2-7, the necessary files and vectors have been provided at https://doi.org/10.5281/zenodo.3928750 as supplementary material.
The dynamics of groundwater do not stop at the coast; there are also maritime discharge zones [70]. The discharge zones on the Mexican coast show the importance of the fresh groundwater component in the submarine groundwater discharge (SGD). In agreement with [10], "The total flux of submarine groundwater discharge, SGD, to the Atlantic Ocean is similar in volume to the riverine flux", SGD is a very important process of the hydrological cycle in México. Certainly SGD is related to the nature of the geological framework; for instance, the presence of submarine springs has been documented offshore from the limestone rocks of the Yucatán Peninsula, as well as those of rhyolitic nature in Baja California Sur, contrasting with the dispersed submarine discharges of groundwater in the Gulf of Mexico [71] or elsewhere in the Mexican Pacific and the Gulf of Mexico.
More than half of Mexico has dry or very dry climates, with an average annual precipitation range of 100-600 mm, and an average annual temperature in the range of 20-26 °C [54]. Dry climates are commonly associated with environments where there is little perennial water on the surface. Of the total area classified as discharge zones, 39.7% have a dry climate. This is consistent with about two-thirds of Mexico having dry conditions, suggesting that the main source of groundwater recharge occurs, as usual, in the highlands, where most of the rainfall occurs. The dynamics of groundwater do not stop at the coast; there are also maritime discharge zones [70]. The discharge zones on the Mexican coast show the importance of the fresh groundwater component in the submarine groundwater discharge (SGD). In agreement with [10], "The total flux of submarine groundwater discharge, SGD, to the Atlantic Ocean is similar in volume to the riverine flux", SGD is a very important process of the hydrological cycle in México. Certainly SGD is related to the nature of the geological framework; for instance, the presence of submarine springs has been documented offshore from the limestone rocks of the Yucatán Peninsula, as well as those of rhyolitic nature in Baja California Sur, contrasting with the dispersed submarine discharges of groundwater in the Gulf of Mexico [71] or elsewhere in the Mexican Pacific and the Gulf of Mexico.
More than half of Mexico has dry or very dry climates, with an average annual precipitation range of 100-600 mm, and an average annual temperature in the range of 20-26 • C [54]. Dry climates are commonly associated with environments where there is little perennial water on the surface. Of the total area classified as discharge zones, 39.7% have a dry climate. This is consistent with about two-thirds of Mexico having dry conditions, suggesting that the main source of groundwater recharge occurs, as usual, in the highlands, where most of the rainfall occurs.
One of the limitations of this study is the low resolution of the spatial information used. This is due to the large size of the country. From the correlation of results of the systematic analysis of the groundwater flow system indicators at the surface, the discharge zones can be identified with a high degree of certainty. However, this certainty decreases with the classification type: discharge-throughflow, throughflow, throughflow-recharge, recharge, and indistinct. As was mentioned earlier, in smaller recharge/throughflow zones, discharge processes also occur, and in discharge zones processes of recharge/throughflow may also occur locally.
Hydrogeological studies are often approached from the local characteristics of an area without full consideration of the regional context and interactions to and from the local and regional hydrogeological environments. This study provides the regional context for local understanding, in order to preserve, recover, or evaluate a site or landscape and promote sustainable interventions.
The results of studies of greater detail are compared with those obtained in the present study, at regional level, and are worth noting. A study of the importance of coastal wetlands was carried out on a coastal floodplain, in central Veracruz state, by Neri-Flores et al. [72]. The study area has coastal dunes and diverse wetland vegetation, riverine vegetation, salt marshes, and an elevation range of <10 m a.s.l. With 31 piezometers installed near the Jamapa and Cotaxtla Rivers and in surrounding urban areas in 2007, 2012, and 2017, the depth of groundwater was registered as between 0.9 and 7.4 m, with an average of 3.1 m. The groundwater table naturally rose by an average 1.3 m in the rainy season, causing groundwater flooding in wetlands and dune lakes. Figure 7 shows this study area to be in a discharge zone (Figure 8), where Neri-Flores et al. [72] conclude that "groundwater and surface water in coastal floodplains are linked and are important to wetland ecosystem functioning" and that "groundwater floods occur every year in the lower-lying areas of the floodplain, which are discharge zones covered with wetlands". In practical terms, it was found that the results of this independent study, with direct field measurements, agrees with the discharge conditions found with the indirect data obtained for the present work from inexpensive published official data.
throughflow, throughflow, throughflow-recharge, recharge, and indistinct. As was mentioned earlier, in smaller recharge/throughflow zones, discharge processes also occur, and in discharge zones processes of recharge/throughflow may also occur locally.
Hydrogeological studies are often approached from the local characteristics of an area without full consideration of the regional context and interactions to and from the local and regional hydrogeological environments. This study provides the regional context for local understanding, in order to preserve, recover, or evaluate a site or landscape and promote sustainable interventions.
The results of studies of greater detail are compared with those obtained in the present study, at regional level, and are worth noting. A study of the importance of coastal wetlands was carried out on a coastal floodplain, in central Veracruz state, by Neri-Flores et al. [72]. The study area has coastal dunes and diverse wetland vegetation, riverine vegetation, salt marshes, and an elevation range of <10 m a.s.l. With 31 piezometers installed near the Jamapa and Cotaxtla Rivers and in surrounding urban areas in 2007, 2012, and 2017, the depth of groundwater was registered as between 0.9 and 7.4 m, with an average of 3.1 m. The groundwater table naturally rose by an average 1.3 m in the rainy season, causing groundwater flooding in wetlands and dune lakes. Figure 7 shows this study area to be in a discharge zone (Figure 8), where Neri-Flores et al. [72] conclude that "groundwater and surface water in coastal floodplains are linked and are important to wetland ecosystem functioning" and that "groundwater floods occur every year in the lower-lying areas of the floodplain, which are discharge zones covered with wetlands". In practical terms, it was found that the results of this independent study, with direct field measurements, agrees with the discharge conditions found with the indirect data obtained for the present work from inexpensive published official data.

Conclusions
From the available environmental information new maps were produced that link components of the environment with groundwater function by following Toth's theory.
There is good correlation between the maps that allows the identification of the natural hydrologic functioning in most Mexican environments at a regional level. When all the elements of the groundwater regional flow system in Mexico are connected, 30% are regarded as discharge zones, and 57% as recharge zones. From the results, only <0.05% of the territory was classified as indistinct in terms of hydrogeological context. These findings are important tools in understanding the interaction of ecosystems and the environmental dynamics related to the groundwater.
Another direct application of the cartography developed is the identification of regional discharge and recharge zones in relation to the stability of protected natural areas and their ecosystem services; the regional discharge zones are zones vulnerable to flooding and inundation, collapse, landslides, and subsidence. On the other hand, discharge zones located on the coast could be potential sites for the development of salinity-gradient energy projects. At the regional level, while areas of recharge, throughflow, or discharge can be identified, more detailed studies will make more evident the coexistence of areas and flow systems. Ideally, this work should be continued, to identify the recharge and throughflow zones with their specific links to discharge zones. This would make it feasible to determine the length and the depth of groundwater flows and define flows of different hierarchies (regional, intermediate or local), as well as identify their possible interactions with valuable environmental components [73]. The methodology used can easily be implemented in other regions like Mexico, where monitoring is scarce and sources of information are not always universally accessible.
The cartography elaborated offers base material in defining environments to be studied in more detail in order to encourage sustainable management. It is hoped that these maps will become the reference for the national-regional context of any study on groundwater and its land-based environmental interaction.