Urbanization effects on the groundwater potential recharge of the aquifers in the southern part of the Basin of Mexico

The growing water pressure caused by the disproportionate increase in population and urbanization in large cities has triggered several adverse effects on groundwater resources. A common effect is the reduction of the natural recharge (produced by rainfall) of the underlying aquifer due to changes in land use/cover, overexploitation, and changes in the precipitation regime. Impact quantification of land use temporal changes is essential to develop groundwater-flow models and generate management programs. For example, Adhikari et al. (2020) analyzed the effects of growing urban Ho Chi Minh City on groundwater recharge, evaluated three future scenarios up to the year 2100 of the low, medium, and high urbanization using an empirical land use projection model and a hydrological SWAT model. They found that an increase in the urban area significantly impacts recharge, reducing it by 30 and 52% in medium and high urbanization, respectively. Another issue is the combined impact of climate change and land use change on groundwater recharge (see, for example, Ghimire et al. 2021; BuhayBucton et al. 2022). Recharge could be any infiltrated water that reaches the water table, but it may return to the atmosphere by evapotranspiration or remain in the unsaturated zone. According to Healey & Scanlon (2010), infiltrating water can be viewed as potential recharge, so for the rest of the document, both terms will be used interchangeably.

Actualized land use cover change (LUCC) maps are needed to estimate potential groundwater recharge (PGR) but are challenging to develop, especially in developing countries. However, technological advancement, satellite images, and software make it possible to integrate remote sensing (RS) techniques and geographic information systems (GIS) to support recharge estimation. Suganthi et al. (2013) delineated potential groundwater zones using combined RS-GIS systems, overlaying lithology, geomorphology, slope, geologic lineaments, drainage, and precipitation maps, to better plan groundwater use in the city of Chennai, India. Oikonomidis et al. (2015) evaluated groundwater recharge potential with RS-GIS in Thessaly, Greece; the recharge potential varied from very high to very low, using and assigning weights to similar thematic layers, including potential recharge and piezometric levels; this allowed them to generate water quality suitability maps. Similarly, Rwanga & Ndambuki (2017) and Wang et al. (2008) used RS-GIS to perform potential recharge calculations. However, they did not address the temporal recharge component caused by the change in land use associated with urban growth.

On the other hand, Guerrero-Morales et al. (2020) evaluated the proportional variation in PGR due to urbanization and climate change in Acapulco, Mexico, where LUCC forecasts, supervised classification of land use images, and a Markov model was used to predict LUCC to 2050. The results indicated recharge losses of 73 and 273 Mm³ in 2017 and 2050, respectively. LUCC is intrinsically related to runoff rates, and several authors have recently applied the runoff coefficient to complete the components in water budgets at unmonitored basins. Peng (2019) introduced a new approach to reconstructing runoff records at large lakes’ ungauged basins, using various components of the water budget. Son & Kwon (2022) analyzed increasing the infiltration rate of rainwater by analyzing the water budget and estimating the runoff coefficient in urban parks.

Mexico City Metropolitan Zone is located in the southern part of the Basin of Mexico (SBM). The urban development has triggered environmental problems, including the upper aquifer overexploitation and reduced recharge zones in the south and southwest of the basin. The Mexican National Commission of Water (CONAGUA, Spanish acronym) has divided the SBM into three administrative aquifers according to hydrogeological and administrative criteria: ‘Mexico City Metropolitan Zone (ZMCM, Spanish acronym)’, ‘Texcoco’, and ‘Chalco-Amecameca’. Development of conceptual models is essential at the beginning of any recharge study. These models can be adjusted as additional data provide new insights into the hydrologic system and may indicate where, when, and why recharge occurs. The main components of conceptual models are influencing factors and water budgets that provide a link between recharge and other processes. Some studies partially analyze the recharge in the SBM. Carrera-Hernández & Gaskin (2008) estimated groundwater recharge through daily water budgets, considering different vegetation and soil types, topography, and climate. Results indicated that recharge occurs mainly in mountainous regions, suggesting that although recharge has diminished in the alluvial plain, urban growth in the SBM has a minimal impact. However, this study covers a short period (1975–1986) and does not consider important local geological aspects such as lithology. Later, Gómez-Reyes (2013) presented a methodology to obtain the water components in hydrographic basins to assess the input and output water volumes using geographical and statistical data. Ruiz (2015) applied water table fluctuation methods and groundwater-flow modeling to evaluate spatiotemporal recharge in the area. The results revealed a variation in mean annual recharge greater than the precipitation that generated it, from 20% in the mountains to 3% in the plains. Meanwhile Peña-Díaz (2019) applied a methodology that accounts for surface water and groundwater to assess recharge availability, using successive approximations to get potential recharge coefficients for the upper aquifer but focusing mainly on high areas.

Other studies that include an integrated conceptual model are those of ‘Sistema de Aguas de la Ciudad de México’ (SACMEX 2005), which determined recharge between 1970 and 2003 considering land use by associating respective runoff coefficients and using a water budget. The results suggested a significant decreasing recharge trend attributed to increased temperatures and reduced precipitation. However, this study was carried out by administrative delegations and municipalities without considering governing parameters that affect potential recharge in the study area. Boyás Martínez et al. (2021) mapped zones of potential artificial recharge using categories from very poor to very good. They used a multiple influence factor (MIF) approach to assign weights to priority variables, such as lithology and land use. The results indicate that mountain ranges such as ‘El Ajusco’ and ‘Chichinautzin’ have the most significant potential recharge; however, essential parameters obtained with remote sensors, such as geological alignments and land use maps, are not included. Furthermore, the technique for assigning weights downplays precipitation, a priority variable. Mautner et al. (2020) carried out a multi-objective analysis to evaluate aquifer management alternatives focusing on flood risk analysis on urbanized ZMCM. They used a groundwater-flow model and incorporated trends in land use, repairs to the water supply distribution network, an increase in wastewater treatment, potential recharge, and imported water potential recharge. However, this study mainly focuses on urban groundwater systems regardless of other essential surface factors.

This work aims to develop a conceptual model of recharge processes and estimate potential recharge volumes considering temporal changes, whose results could be used as an input for a spatially distributed groundwater-flow model. The present study applies a potential recharge estimation methodology in the SBM, which integrates a combined RS-GIS technique with the analytic hierarchy process (AHP) to assign weights to different hydrologic factors as thematic layers. The proposed methodology is divided into three stages. First, geographic data are collected, processed, and analyzed using GIS and RS tools. Second, long-time-changing factors (assumed static) are used to map the study area's potential recharge (PR) regions. The mapping was done by hierarchizing the spatially distributed data such as geological alignments, drainage network, lithology, surface slope, and historical annual precipitation. In the last stage, the dynamic factors whose data is recorded over the 1970–2019 period of analysis were mapped every 5 years to assess the temporal evolution of potential recharge and its spatial distribution. In this way, the temporal dynamic processes that promote PR are considered by employing a surface water budget equation to estimate PR volumes at each PR region. The novelty of the work is to present a methodology that separates the static parameters to produce a potential recharge zonification from those that vary in time used to calculate the evolution of potential recharge, which facilitates the calibration of the recharge in a groundwater model.

Recharge conceptual models are developed using as much input data as is available, and revising and adjusting these models as additional data and analyses provide new insights into the hydrological system. The factors and data that can influence the conceptual model are climate, geology, topography, hydrology, vegetation, and land use. According to Healey & Scanlon (2010), climate variability is often the most critical factor affecting variability in recharge rates. Since precipitation is the source of natural recharge, it is the dominant component in the water budget in most cases. Second, geological aspects such as permeability, soils, and lithologies that can promote recharge must be considered since water quickly infiltrates and drains through the root zone. After that, the surface topography may affect diffuse and local recharge since steep slopes tend to have low infiltration and high runoff rates. Also, hydrology is important because of the natural interaction between surface water and groundwater-flow systems. Not less important is the vegetation and land use since irrigation (for example) may constitute a significant amount of recharge. Lastly, the urbanization component caused by converting natural land to impervious areas that can inhibit recharge must be considered.

The data were collected through public government agencies distributed through official websites and were visualized, analyzed, and pre-processed using ‘geographic information systems (GIS)’ tools. Table 1 summarizes the collected data and their sources.

The proposed methodology for this work consists mainly of collecting, processing, and analyzing geographic data using GIS tools. This information is handled as layers of thematic information, specifically in shape format (⁠⁠), through GIS tools. This information corresponds to the factors that most influence the potential recharge process in the study area, such as geological alignments, drainage network, lithology, land use, topography, and precipitation. The methodology can be grouped into three stages as follows:

• Stage 1. Identification and assignment of weights by zones in each layer of information, according to their physical characteristics favoring recharge.

• Stage 2. Regionalization of ⁠. First, defining weights for each thematic layer using AHP based on importance value. Second, those values are used to apply the ‘GIS-weighted overlay’ tool to map the regions.

• Stage 3. Based on the surface water budget equation, a computation for the infiltrated volumes is performed for each region. regions are areas where groundwater recharge is most likely.

The methodology considers spatial heterogeneity by including long-time-changing factors, defining average surfaces for calculating the water budget. These static factors are used to define the PR regions, which are surfaces that allow qualitative analysis of the potential recharge by classifying them into five categories: very poor, poor, moderate, good, and very good. Furthermore, the dynamic factors are subsequently used for the quantitative analysis that comes in phase three, where climate data and land use change are employed to calculate the associated potential recharge amounts for each defined PR region.

Different thematic layers are developed by processing the geographical data described in Table 1. First, a quality analysis of the collected information was carried out before the primary processing. This analysis ensures spatial and temporal precipitation records consistency with the analysis periods over the entire SBM. In addition, anomalies, inconsistencies, missing data, and atypical values were detected by analyzing time series. Consecutively, thematic data layers were mapped within the study area using GIS tools.

Geology is one of the principal factors controlling groundwater's occurrence and movement. The soil texture can be defined as the relative proportion of silt, sand, and clay present in the soil. Clay's smaller pores induce more runoff and less potential recharge, while sandy soils with larger pores induce more potential recharge. Different geological formations can generally be identified according to their lithology and soil type. Also, the terrain slope is crucial because it determines the time required for water to accumulate in each location. The lower the slope, the lower the runoff coefficient, and the higher the potential recharge.

Regarding LUCC, the first refers to ‘human activities and various uses that are carried out on the land’, while land cover refers to ‘vegetation, water bodies, rocks/soil, artificial cover and others resulting due to transformations’. LUCC defines the potential recharge rate of rain and water. On the other hand, drainage density depends on the slope, the nature, and attribute of the parent rock, and the regional and local fracture patterns. The drainage network reflects the lithology and structure of a particular area and is studied according to two parameters: (a) drainage pattern, associated with the nature and structure of the substrate and (b) the drainage texture, which is related to the permeability of the rock/soil. The less permeable rock, the lesser potential recharge, and minor weights are assigned to those areas. While geological alignments are linear or curvilinear structures on the Earth's surface and represent the weakest zone of the bedrock. Also, alignments may be correlated with faults, fractures, joints, bedding plans, and lithological contacts.

Essential components of the short-time-changing factors are precipitation and temperature since a proportion of the total precipitation volume infiltrates, but the rest is removed by runoff or evapotranspiration processes. In the area, 132 climatic stations are present, but following the criteria proposed by Salas et al. (1980) and Sánchez-Quispe et al. (2021), only 85 were used to interpolate the precipitation and temperature raster files for the established periods of analysis, generated as the average of 5 years, between 1970 and 2019. Next, Table 2 describes the procedure to generate each thematic layer. Thematic maps such as geology, terrain slope, land use and land change, drainage network, and alignment density were generated using satellite images and conventional data to assess PR regions.

Like any other balance, the water budget proposed here is only a basic tool that can be used to assess the occurrence and movement of water in the natural environment. Average state conditions are generally analyzed, and the representing equation can be very simple or complex, depending on the studied environment and its objectives. The water budget analysis would provide an initial understanding of the system that can be later refined as more data becomes available. As a complement, a numerical model would provide a more refined understanding of the subsurface and surface flow system, at different scales and levels. Nevertheless, it is also important to recognize uncertainty in all water budget components associated with the natural variability of the hydrological cycle and measurement errors (Healy et al. 2007).

The AHP is a simple tool to help decision-makers; it is used in decision-making problems that involve very complex processes. It was developed in the late 1960s by Saaty (1990), who based his research on the military field and his teaching experience. In this way, the AHP allows discerning between different processes to determine a solution that best satisfies the combination of possible alternatives (Osorio & Orejuela 2008). The first uses of the AHP were given to solve decision problems in multicriteria environments.

Goepel (2013; 2018) implements the method based on the solution of an eigenvalue problem. This way the results of the pairwise comparisons are arranged in a matrix. The matrix's first (dominant) normalized right eigenvector gives the ratio scale (weighting), and the eigenvalue determines the consistency ratio (CR). The process follows the next steps: (1) get pairwise comparisons of the factors; (2) apply selected AHP judgment scale; (3) fill the decision matrix from pairwise comparisons; (4) find the eigenvector using the power method; (5) calculate the dominant eigenvalue from the eigenvector; (6) calculate the consistency ratio (CR); (7) calculate the inconsistency matrix; (8) identify and highlight the top three inconsistencies; and (9) go back to step 1 until the user submits his judgments.

According to Goepel (2013; 2018), the AHP process is performed by filling out questionnaires which may require further clarifications about the completion procedure or pairwise comparison. In this way, AHP hierarchies are derived from the most consistent matrices. A measure of consistency is the consistency index (CI). From this, a consistency ratio is derived, using a randomized index (RI), the average CI for randomly filled matrices. Also, a parameter α is needed as a level of consistency to adapt the cut-off value (CR <α), for declaring the matrix inconsistent. Saaty's rule of thumb is to accept only judgment matrices with CR < 0.1.

A consensus measurement that includes an aggregated group of results is calculated for a group decision-making factor is calculated based on Goepel (2018). This process includes Shannon entropy in two independent components (alpha and beta diversity) to derive the AHP consensus indicator. The similarity measure S (or consensus indicator) depends on the number of criteria and uses a linear transformation to map it to a range from 0% (no consensus) to 100% (full consensus). This range is then categorized into five categories: very low (<50%), low (50–65%), moderate (65–75%), high (75–85%), and very high (>85%).

The drainage pattern is one of the most important indicators of the potential recharge process because the underlying lithology mainly controls drainage patterns and density. Furthermore, the stream pattern reflects the potential recharge rate of precipitation compared with surface runoff. The potential recharge/runoff rate is mainly controlled by permeability, which is a function of rock type and the underlying bedrock fractures. Comparing mountainous terrains and valleys, the one with the higher drainage density is usually less permeable. Therefore, the denser the drainage network, the lower the recharge, and vice versa. Figure 7(a) and 7(b) shows the resulting density mapped with the spatial analyst GIS software, ranked from 1 to 5.

To map the highest and lowest PR regions, the ‘GIS-weighted overlay tool’ was used. This tool requires analysis layers to be supplied with a certain percentage of influence. This influence was computed using the AHP technique to structure, measure, and synthesize information. In this case, to assign a specific weight, layer importance and hierarchy of the proposed thematic GIS layers were considered. The methodology was applied as in Goepel (2018) using an online AHP calculator. First, ten pairwise comparisons were performed for the five factors using a questionnaire. The answers should include the importance on a scale of 1–9. AHP Scale: 1 – Equal importance, 3 – Moderate importance, 5 – Strong importance, 7 – Very strong importance, and 9 – Extreme importance. The questionnaire answers are highlighted in Table 5; this was done using experience and knowledge of the study area. The consistency ratio was CR = 6.3%, considered satisfactory (<10%).

Table 5

Factors pairwise comparisons questionnaire, with highlighted answers

View Large

The hierarchies of each factor are the resulting weights from the pairwise comparisons. Those weights are assigned to each layer based on the principal eigenvector of the decision matrix. In this case, the principal eigenvalue was 5.28, with 6 iterations and δ = 3.1⁻⁹. Table 6 summarizes the resulting factor hierarchy ranks and priority percentages. Those were used as input to the GIS-weighted overlay tool.

In the present study, PR regions are delineated using GIS techniques. First, thematic layers of long-time-changing factors were included: geology, drainage density, annual precipitation, and geologic alignments. Each layer is a thematic map of regionalized weights (from 1 through 5, where the higher the number, the better PR); after that, total layer weights were assigned to each one using the AHP technique. Then, the ‘GIS-weighted overlay tool PR regions’ were mapped. The PR regions for this work were defined into five categories: very poor, poor, moderate, good, and very good.

The mapped PR regions reflect long-time-changing parameters that are fixed for the rest of the analyses. Notably, the largest PR regions correspond to poor and moderate (orange and yellow areas), with 33.9 and 37.5% of the SBM total area (4,080 km²), respectively. These two areas, together with the very poor PR with just 0.1%, represent the lowest percentage of recharge for all regions; in turn, they present the highest percentage of the urbanized surface along the analyzed 5-year periods. On the other hand, the good and very good PR regions, with 24 and 4% of the SBM surface, respectively, presented the smallest urbanized areas along the periods. These percentages concord with what was found by Boyás Martínez et al. (2021), whose results showed that the medium category predominates over the entire area. However, they also include two additional aquifers in the northern part of the Basin of Mexico. On the other hand, in their study, the good and very good areas represented 19.1 and 5.2% of the SBM total area, respectively, values slightly smaller than the ones reported in this work.

Urban and non-urban areas were remotely detected for the study area developing LUCC layers through Landsat images downloaded from the USGS agency's EarthView website. The images were classified using the free Multispec^® software (Biehl & Landgrebe 2002). Table 7 shows the classified image representative of each 5-year between the years 1970 and 2019.

The 5-year LUCC layers and the PR regions map (Figures 8 and 9) were used to obtain the proportions of urban and non-urban areas. Urban areas include housing developments, roads, sidewalks, and parks, while non-urban areas include bare soil, agricultural land pasture, woodlands, and wetlands. The respective runoff coefficients were selected for each land cover (Table 3), including soil properties and surface slope. A global runoff coefficient for each PR region was computed as an equivalent runoff coefficient (⁠⁠). Table 8 shows the respective computed values of for each of the analyzed 5-year periods by PR regions in the study area.

Potential recharge volumes were calculated using the surface water budget (Equation (4)) for the 1970–2019 period and for each 5-year time interval using data in raster format with the GIS-Map Algebra tool. Then, the potential recharge volumes were computed using Equation (4) for each PR region with the different urban and non-urban land use areas, through their respective associated runoff coefficients ⁠. Afterward, was calculated for the representative or average records of each PR region (Figure 12(b)). A map of the available precipitation layer is shown in Figure 12(a), used to generate the surface runoff and potential recharge, by subtracting the calculated value of the from the total precipitation for each period of analysis. As an illustration, Figures 12(a) and 12(b) show the resulting maps of real evapotranspiration and available precipitation for selected periods, while the rest are found in the Supplementary material.

According to Figure 14, urbanization plays an important role in the percentage of PR. While in the period 2005–2009, where there was a considerable increase in precipitation, its effect on the trajectory of the percentage of PR seems only to level it, in 2015–2019, there was also a slight increase in precipitation, but the percentage of PR is more notably affected by the increase in urbanization. Of course, this is for the percentage of PR in the whole SBM, but precipitation influence might be more significant by analyzing each PR region separately. Also, Table 9 shows the results of estimated infiltrated volumes in Mm³/year. Note that the first period with the highest precipitation infiltrated percentage may suggest a correlation with the urbanization surface percentage for the same period. The precipitation volume varies in a range of 3017.89–3,540.81 Mm³/year corresponding to 1990–1994 and 2005–2009, respectively. While the potential recharge varies from 24% (823.71 Mm³/year) to 19% (630.83 Mm³/year) of precipitation, corresponding to 1975–1979 and 2015–2019, respectively, representing a maximum decrease of 24% in potential recharge. This result is close to the Peña-Díaz (2019) result of 647 Mm³/year, but his study was for the entire Valley of Mexico from 1980 to 2011.

In a previous study (SACMEX 2005), the resulting mean recharge for the area was 12.3 m³/s; this value is close to that obtained in this study for the poor PR region, but the study only accounted for the ZMCM aquifer. The total area used in that study was 1,906.05 km², representing only 57% of the area considered in this work. They also reported an annual volume estimated on average for aquifer recharge of 274.7 Mm³/year. Meanwhile CONAGUA (2015) officially reported 512.8 Mm³/year (which is equivalent to 16 m³/s), for the principal aquifer of the metropolitan zone (ZMCM). The mean value of PR (including all 5-year periods shown in Table 9) has a gross average of 682 Mm³ or 22 m³/s, which is not too far from that officially reported, considering that it was calculated for the three administrative aquifers in the SBM.

The variation observed in the results for the potential recharge shows a dependency on the input data and the hydrological information since they present a high variation instead of unique values, especially if they are analyzed on small scales. In this sense, the balance equation estimates PR as the complement of the other flows involved in the hydrological cycle that are more easily measured or estimated. However, given this uncertainty, the PR is more feasible as a realistic range instead of an exact value that satisfies the objectives of this study. This confirms the importance of accurately estimating a model's parameters, as indicated by other authors, since using more complex models introduces more significant errors and variations in estimating the water balance components (De Silva 2005).

This study developed a conceptual model of recharge processes and a methodology to estimate potential recharge volumes considering its temporal changes. The methodology separates the long-time varying parameters to produce a potential recharge zoning from those that vary in a short time. This separation facilitates using the potential recharge results for the groundwater recharge calibration in a spatially distributed groundwater model. The methodology is applied in the estimation of PR in the SBM aquifer. The methodology mainly consists of collecting, processing, and analyzing georeferenced thematic data using GIS software. Geographical data, such as geological alignments, drainage networks, lithology, land use, topography, and precipitation, define PR regions. Short-time variables and their maps were generated in intervals of 5 years along the period 1970–2019 to evaluate the evolution of potential recharge relative to the dynamic processes over time that generates it. The potential recharge is then estimated using a surface water budget equation at each PR region.

The mapped areas with the most significant recharge potential resulted in the highlands, mainly the Sierra de las Cruces and Chichinautzin, where precipitation and lithology favor potential recharge. On the other hand, the areas with the lowest potential recharge correspond to the valley, where the urban development is placed, especially those areas or neighborhoods with high population density (according to Landsat image reclassification), and the underlying ancient lakebed in the plain areas turned into a clay aquitard with low permeability significantly reduces the potential recharge.

This study proves that potential recharge depends on several critical factors. Although a higher sensibility was observed for changes in precipitation, another sensible parameter was the size of PR surfaces (considered with slightly more impact when analyzing the water budget equation by PR regions separately). In third place, the runoff coefficient depending on LUCC was found sensible, especially using values close to 1. Although several works consider a low urban impact in the total recharge, this work proves that it can be significant in highly urbanized areas. Analyzing the whole SBM surface, urban change significantly influenced PR.

Results also showed slight increases in potential recharge at some periods; this increase may confirm the results reported in SACMEX (2005) regarding the influence of temperatures. With a different behavior in the more recent periods, where a slight decrease in recharge was observed compared to earlier periods, temperatures may also play a role in conjunction with changes in land cover, urbanization, and even global warming.

In this analysis, other non-natural potential recharge sources, such as irrigation returns, leaks from the hydraulic and drainage network, reservoir recharge, and conservation structures, are not included. Future works will analyze these factors. The results can be used as input data to feed a numerical groundwater model for the SBM area. Finally, this study also deals with the sensitivity of the parameters, as shown by the potential recharge results, with the differences between urban and non-urban analyzed regions.

The limitations of the estimations presented in this work are commonly observed in any water budget, such as the precision of the input data and the understanding of the hydrological processes of a study area is defined using information from other studies and introduces uncertainty in the calculations. Also, the time scale with which the variables that change in a short time were analyzed (5-year periods) depended on the available satellite information. The periodicity of satellite information may limit the applicability of the proposed methodology in studies that need to reflect more dynamic changes. These drawbacks constitute a focus point for future work. In addition, the results of this work may be used as a complement for a more complete analysis through the use of numerical models that solve the balance equations that describe the hydrological processes in the saturated zone and in the interface between surface water bodies and the saturated zone.

To project PAPIIT-UNAM: IN110620, ‘Caracterización mediante modelación matemática del origen de nitratos del agua subterránea en la parte sur de la Ciudad de México, cuantificación de incertidumbre y propuestas de medidas de control’, 2020–2022. To Cátedras CONACYT Program.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.