Groundwater Model for Karst and Pelitic Aquifer Systems from a Semi-Arid Region Under Climate Change Scenarios: A Case Study in the Vieira River Watershed, Brazil

Abstract

Water scarcity is a global issue, especially in semi-arid and arid regions where precipitation is irregularly distributed over time and space. Predicting groundwater flow in heterogeneous karst terrains, which are essential water sources, presents a significant challenge. This article integrates geology, hydrology, and water monitoring to develop a pioneering conceptual and numerical model of groundwater flow in the Montes Claros Region (Vieira River Watershed, Brazil). This model was evaluated under various climate change scenarios, considering changes in rainfall, groundwater consumption, and population growth over the current century. The results indicate that a decline in water table levels is inevitable, primarily driven by population growth and high pumping rates rather than rainfall fluctuations. This underscores the urgent need for improved monitoring, model upgrading, and more importantly, targeted water resource management for Montes Claros.

1. Introduction

Global population growth is intensifying the environmental impacts from agricultural and industrial activities necessary for human life. Water is the primary natural resource essential for both rural and urban activities to continue expanding [1]. With increasing issues related to water availability, the effective management of water resources has become crucial [2,3]. Water scarcity arises if water demand exceeds water availability [4], initially putting pressure on surface water sources and eventually, on aquifer systems, which are increasingly used as fresh water sources [5].

Regions with semi-arid and arid climates and/or irregular precipitation distribution throughout the year tend to face significant water scarcity issues. In the context of climate change, this challenge is expected to worsen [6]. The distribution of water resources across continents is uneven, with some countries like Canada and Brazil possessing extensive and water abundant river basins [7]. However, even within these nations, water distribution can be highly irregular, with some regions presenting abundant water resources, while others face severe water shortages.

Karst terrains cover approximately 20% of the Earth’s surface [8], and its aquifer systems are increasingly vital for supplying populations. It is estimated that around 25% of the world’s population relies on the exploitation of these aquifer systems, which feature typical structures such as caves, sinkholes, poljes, conduits, dolines, and springs [9]. Karst aquifers are marked by heterogeneous underground water flow, and many analysis and modeling methods have been discussed in the literature (e.g., [10,11,12,13]).

The northern region of the state of Minas Gerais, bordering the Brazilian semi-arid zone, has faced water availability challenges for decades [14,15]. The region’s scarce rainfall and high evapotranspiration rates make it highly susceptible to drought [16,17]. Considering natural geological complexities and climate control in assessing water availability and its underground flow, this study compiles information from the Vieira River Watershed, Montes Claros Region (Minas Gerais), aiming to achieve the following goals: (i) to detail the methodology used to propose, for the first time, a groundwater flow conceptual and numerical model for the aquifer systems and (ii) to analyze potential water availability/scarcity related to exploitation under different climate change scenarios.

2. The Vieira Watershed

2.1. Local Settings

The Vieira River Watershed is situated in Montes Claros, southeastern Brazil (Figure 1). Climatically, the watershed is on the edge of the Brazilian semi-arid region, which has historically faced water scarcity issues. Montes Claros is the largest city in the northern part of Minas Gerais state, with over 400,000 residents in both urban and rural areas, serving as a vital economic hub. The region’s economic development surged in the 1960s due to specific public policies targeting cities in the Brazilian semi-arid zone [18,19].

Montes Claros experiences concentrated rainfall, primarily during the summer, from November and March [20]. Due to geographic and climatologic conditions, coupled with population and economic growth, Montes Claros has faced significant challenges regarding water supply, with notable conflicts recorded since the 1980s [14,15].

2.2. Physiographic Features

The Vieira River Watershed is part of the Verde Grande River Watershed, a tributary of the São Francisco River. The entire Montes Claros urban area lies within this Watershed. This region features rugged terrain, with elevations below 600 m in the Vieira River Valley, and altitudes exceeding 1000 m a.s.l. (above sea level) to the west of the urban area (Figure 2). Over the past decade, the average annual rainfall has been around 920 mm [20].

In terms of land area, the urban region is not significant compared to the total area of the watershed, which also lacks significant agricultural activities [14].

2.3. Geological Context

The Vieira River basin is geologically situated within the São Francisco Craton [21], predominantly composed of Archean and Paleoproterozoic rocks [22]. The study area is comprised of Proterozoic metasedimentary rocks from to the Bambuí Group outcrop, spanning along the São Franciscana Basin with substantial thicknesses, reaching up to 3000 m [23,24]. The Bambuí group is divided into the Sete Lagoas, Serra da Santa Helena, Lagoa do Jacaré, Serra da Saudade, and Três Marias areas [25,26,27].

In the western part of the watershed, the carbonatic Lagoa do Jacaré Formation and the pelitic Serra da Saudade Formation are predominant, while in the east, the pelitic Lagoa do Jacaré and Serra da Santa Helena Formations outcrop prevail.

The urban area is located over the pelitic domains of the Serra da Santa Helena and Lagoa do Jacaré Formations, but with reduced thickness. A few meters below the surface, the Sete Lagoas Formation can be observed through geophysical deep wells’ data and is characterized by the intercalation of pelitic and carbonatic domains (Figure 3). Regionally, these stratigraphic layers can exceed 400 m in thickness, except for the Serra da Saudade Formation, which has a maximum thickness of 120 m. Structurally, the rocks exhibit subhorizontal to slightly inclined banding, kilometer-long open folds, and some brittle structures, observed on satellite images and confirmed in the field. The perpendicular NNE–SSW and WNW–ESE trends control the water drainage network and various karst features [28,29,30,31].

2.4. Hydrogeology

Bhering et al. [20] have conducted a hydrogeological study of the Vieira River Watershed, subdividing it into three aquifer systems: karstic fissured aquifers, fractured formations, and non-expressively, local granular aquifers (Figure 4a). Of the two most important systems, the karstic fissured aquifers are predominant and are found in limestone from the Lagoa do Jacaré Formation. This domain features high variability in regards to its hydrodynamic parameters, including primary and secondary porosity, in addition to fractures. The main recharge zone coincides with the significant karstic system known as Lapa Grande State Park. On the other hand, poorly fractured aquifers developed in the metasiltstones and metamudstones of the Serra da Santa Helena and Serra da Saudade Formations, exhibiting low transmissivity and often behaving like aquitards, with higher productivity in areas with carbonate intercalations.

The urban area of Montes Claros lies primarily within the poorly fractured aquifers, resulting in generally low groundwater productivity. Specific capacity data indicate that higher productivity zones are associated with regional structural features, such as intersections of NNE–SSW alignments and fractures. The regional groundwater flow is mainly influenced by differences in hydraulic gradients and structural features, with both local NW/SE flows and a regional SW/NE flow pattern observed (Figure 4b). Wells are, on average, 86 m (depth enough to reach the lower limestones from the Sete Lagoas Formation; see Figure 3b), reaching up to 230 m deep, and the average groundwater level is ca. 14.5 m deep, with a maximum around 60 m.

3. Methodology

Groundwater flow and numerical models utilize regional data and specialized software, considering the particularities of the aquifer systems. For the aquifer system in the Vieira River Watershed, the methodology involves regional and hydrogeological characterization to parameterize the conceptual and numerical model applied to various climate change scenarios.

3.1. Topography

Digital terrain models from Brazil’s geomorphometric database, provided by the INPE (National Institute for Space Research) and based on SRTM (Shuttle Radar Topography Mission) data [32], were used. These publicly available models were processed in a the GIS environment (QGIS 3.10.2) for subsequent stages.

3.2. Geological Model

Satellite images were used to define the morphology, textural features, major lineaments, pattern of discontinuities, and structural domains. After analyzing the drilling database (wells) of the Geological Survey of Brazil [27], which contains data for more than 700 wells in the area, 125 logs with lithological descriptions were selected. Geological mapping campaigns were conducted and integrated to build the 3D geological model using the Leapfrog Geo software (2022.1.0).

3.3. Hydrogeological Information

After analyzing existing data and previous work, new information was added to define the main hydrogeological compartments of the watershed. Using the SIAGAS database [27], new interpretations were generated regarding the regional potentiometric surface (Figure 4b), specific capacity, preferred flow directions, and discontinuities of the massif (Figure 4c). The caves database from the region [33] was also used.

3.4. Groundwater Monitoring

More than 700 wells were used to build the models, and new measurements were captured in the field (Figure 4d). Using seven water level loggers, 11,000 measurements were automatically generated between December 2020 and March 2022 in three selected wells, at the source and downstream of the Vieira River, close to the Verde Grande River. Additionally, eleven wells were manually monitored, recording water levels during dry and rainy periods (Figure 4d).

In addition to the use of public monitoring data provided by the Geological Survey of Brazil, other new monitoring data, including automated instruments, were presented in previous works, as shown in the Figure 5, and used in the interpretations, simulations, and calibrations of the model.

3.5. Conceptual and Numerical Model and Parametrization

The FEFLOW platform uses the FEM (finite element method). In summary, this method discretizes the problem into a finite number of elements to approximate the solution of partial differential equations as a set of algebraic equations that can be solved computationally. The model was configured as a variably saturated medium using the Richards’ equation. In this study, groundwater flow modeling was chosen under the assumption that it occurs exclusively in a saturated zone, where moisture is considered constant. This approach simplifies the model and has proven particularly useful, especially given the large area of the watershed under analysis.

The proposed model is based on Darcy’s law, which defines equations governing groundwater flow in fully saturated porous media. Mathematical models based on Richards’ equations (Equation (1)), which define flow in partially saturated zones, are simplified into equations based on Darcy’s law (Equation (2)), which describe groundwater flow under fully saturated conditions [34]. This simplification can increase model efficiency and reduce processing time [35], which is important for simulations in large-scale watersheds.

$$s{S}_o{\displaystyle \frac{𝜕\psi }{𝜕t}}+\epsilon {\displaystyle \frac{𝜕s}{𝜕t}}+\nabla \cdot q=Q+Q_{EOB}$$

(1)

$$q=-k_{r}K{f}_{\mu }\left(\nabla \psi +\left(1+\chi \right)e\right)$$

The governing equation for the fully saturated model simplifies to the following:

$$S_o{\displaystyle \frac{𝜕h}{𝜕t}}+\nabla \cdot q=Q$$

(2)

$$q=-K\cdot \nabla h$$

where

• s: medium saturation;
• S_o: specific storage;
• ψ: hydraulic head;
• ε: saturation adjustment coefficient over time;
• ∂ψ/∂t: variation of hydraulic head;
• ∂s/∂t: variation of saturation;
• ∇⋅q: gradient of the flux vector q;
• Q: generic source or sink function;
• Q_EOB: source/sink correction term of the extended Oberbeck–Boussinesq approximation.

After constructing the 3D geological model using Leapfrog software, Leapfrog Hydro (2023.1.0) was used to generate the final mesh. The solids were exported to FEFLOW software (8.0) to simulate groundwater flow by discretizing the problem into finite elements.

The first layer of the model was established as a phreatic, with subsequent units dependent on the upper layers. A simplification was made in which the entire model volume was considered saturated, with the water level equivalent to the surface, where the pressure is zero and negative charges represent the aerated zone. Simple triangular-type prismatic elements (six nodes per element) were used.

3.5.1. Mesh

After incorporating detailed topography, database, and geology information, a finite element model was generated with the following characteristics:

• The upper limit corresponds to the topographic surface, covering an area of approximately 579 km²;
• The model assumes a thickness of approximately 840 m (i.e., groundwater depth);
• The calculation grid represents the geometry of the hydrogeological system, as well as fault structures and drilled wells.

The discretization and extension of the model’s finite element mesh were defined, as shown in Figure 6. The model’s mesh was refined around watercourses, geological structures, and pumping wells to ensure computational feasibility and detailed representation.

After defining the discretization zones, Leapfrog Hydro software (2023/1) generated regular triangular prisms measuring approximately 130 m. A total of 658,987 elements were discretized for each of the 7 layers, totaling 379,992 nodes. The mesh quality was analyzed using the Delauney criterion violation method, ensuring the equilaterality of the internal angles (Figure 7).

3.5.2. Boundary Conditions, Parametrization, Groundwater Recharge, and Calibration

The hydrographic watershed limit was set with a zero-flow boundary condition, equivalent to a hydrogeological watershed, to eliminate potential flow between neighboring basins.

The model’s boundary conditions were set as first-order (Dirichlet) for the drainages, rivers, and pumping areas. The wells were inserted as multilayers, according to the original database. The model boundaries were defined with a no-flow condition, meaning there is no inflow or outflow of water to adjacent watersheds. Thus, as a premise, the hydrogeological watershed was assumed to be equivalent to the hydrological watershed.

The type 1 condition (Dirichlet), when applied to a differential equation, specifies the values that a solution must take on the boundary of the domain. It determines the value that the solution must assume at a given point as a specified variable. In other words, the head in the cell where the boundary condition is assigned is equal to the specified head value, while the flow entering or leaving this cell is uncontrolled, potentially tending to infinity.

Regarding the hydrodynamic parameters,

Table 1. Hydrodynamic parameters used in the numerical model for the geological units—K (hydraulic conductivity), S_o (specific storage coefficient), and ε (porosity).

Geological Unit	KX	KY	KZ	So	ε
Alluvial deposits	9.00 × 106	1.08 × 105	9.00 × 107	1.00 × 104	2.50 × 101
Soil	1.04 × 106	1.10 × 106	5.79 × 108	5.00 × 105	1.00 × 102
Lower limestone	1.16 × 106	2.31 × 106	1.16 × 107	1.00 × 103	2.50 × 101
Siltstone	9.26 × 107	1.10× 106	1.04 × 107	1.00 × 106	6.00 × 103
Limestone	2.80 × 106	3.36 × 106	2.80 × 107	8.00 × 105	1.60 × 101
Contact zones	3.47 × 105	3.47 × 105	3.47 × 105	5.00 × 104	2.00 × 101

summarizes the main values applied to the geological units described in item 2.3. Values of the same order of magnitude were found in similar geological units of the Bambuí Group, specifically in the Serra da Santa Helena, Lagoa do Jacaré, and Serra da Saudade Formations [33].

Karst aquifers are a typical example of a pore-fractured medium. However, the dual porosity flow, classified by Ref. [36] as triple permeability, can be approximated using equivalent porous aquifer models. Multiple authors [37,38] have applied numerical modeling in karst aquifers as equivalent porous media in areas larger than 150 km². Thus, an equivalent porous aquifer model was used in this study due to the size of the study area. Despite the scale effect in a karst aquifer, where the variability of hydrodynamic parameters can be significant in small areas, the flow system in relatively large areas of a karst aquifer follows a pattern in which the mean flow approximates the median flow. In such cases, the overall system behaves similarly to a porous medium. Therefore, the continuity equation of flow can be applied, meaning that the variation of mass over time in a scalar volume can be used in Darcy’s equation. This approach approximates the continuity equation of flow in three-dimensional Laplace systems, assuming that hydraulic conductivity data do not vary within the same lithotype, and that the flow is homogeneous. This assumption is supported by the average specific yield of the area, which can be linearized and correlated with a constant transmissivity for the karst region, as shown by Refs. [39,40].

Regional structures impacting the landscape were identified through aerial imagery analysis and mapping. These structures, which are important conductors of groundwater, especially when they show more advanced karstification, were incorporated into the model as 2D discrete features. Figure 8 illustrates the boundary conditions and the 2D discrete features implemented in the model.

Estimating recharge in a karst environment is challenging, and various methodologies have been applied. The APLIS method [41] and its modified version [42] were developed for measuring average recharge in Mediterranean karst aquifers [43]. Both methods have also been applied in other karst system environments, such as Cuba [44] and the Peruvian Amazon region [45]. This study also exemplifies how the APLIS method can be adapted and applied in different geological and hydrological contexts to estimate groundwater recharge. For the present work, scoring was initially assigned based on the original method and adjusted during numerical model calibration. Below are descriptions of the parameters:

• Altitude (A): Refers to the elevation of the terrain. The lowest point of the basin was used as a reference, as the Vieiras River is the base level of this basin and is located in a low-energy region in terms of surface flow. The minimum elevation of this point is approximately 585 m. Thus, two altitude zones were considered relative to the 585 m reference, following the original methodology. The first zone ranges from 0 to 300 m (585 to 885 m), and the second from 300 to 600 m (885 to 1185 m).
• Slope (Pendiete, P): Refers to the incline of the terrain. In the area, the maximum slopes reach 70%. In this study, the slope was the first parameter considered differently than in the original method. In the Vieiras River Basin, a significant geomorphological difference is observed between the river plain and the flow regions of its tributaries. The Vieiras River flows through a valley that has been leveled by erosive processes, with low slopes and differing geology. However, this region consists of metapelitic rocks of the Lagoa do Jacaré Formation and mainly from the Serra de Santa Helena Formation, being very close to the local base level. This could indicate a low infiltration potential, mainly due to the permeability characteristics of this lithology. Therefore, adaptations were made in the region, associating a correction factor. To delineate the topographic zones of higher from those with lower altitude, the Terrain Segmentation tool of SAGA GIS was used (Figure 9). This tool divides an area into segments or regions with similar characteristics and is frequently applied in GIS analyses.

• Lithology (L): Refers to the geological composition of the terrain. The classification is based on local geological, geomorphological, and hydrogeological maps, as well as field inspections and observations. In this study, values corresponding with or similar to those proposed in the original method were used.
• Infiltration (I): Determined by observing preferential infiltration pathways in the field, such as major faults, fault-controlled valleys, and karst features, among others. In addition to mapping the main lineaments and features and interpreting them in a GIS environment, the Peak and Depressions algorithm from Surfer^® software version 26 was used. This tool analyzes and identifies the highest (peaks) and lowest (depressions) points on a surface map, such as a DEM or another type of grid map, and is widely used in geology, geomorphology, environmental studies, and territorial planning. Polygons identified as depressions, generated from this tool, were used as preferential recharge points. Another important factor identified in the preferential infiltration zones was the presence of lineaments, which were inserted by constructing a density map (Figure 10). Kernel analysis was performed on the DEM (QGIS), with a search radius of 300 m, considering three average thickness values of 10, 80, and 200 m. These values were used to assign weights to the structures when constructing the lineament density map, with a weight of 1 for 10 m, 8 for 80 m, and 20 for 200 m drainages. Karst features extracted from Ref. [46] were also assigned a score of 20 within a 200 m radius for cavity protection areas, while non-carbonate rocks (mainly metapelites from the Serra da Saudade, Lagoa do Jacaré, and Serra da Santa Helena Formations) were given reduced scores that were gradually increased for the rocks of the Lagoa do Jacaré Formation, especially for the limestones. The estimated density analysis varied between 0 and 0.07 m/m². This lineament density data was aggregated with the depressions map and normalized on a scale from 1 to 10.
• Soil Type (S): Parameters were assumed based on the original values of the method. The soil map presented in the area’s characterization was used as the base, obtained from studies by Ref. [47]. Some soil classes, such as haplic nitisols, are not found in the original method [41]. In this case, data from Ref. [15] were used and correlated based on the average K value.

Figure 10. Infiltration (I) for the APLIS method after using GIS system interpretation and applying tools to determine topographic highs and lows (peaks and depressions) and fracture and to conduct karst feature density analysis (kernel).

Figure 10. Infiltration (I) for the APLIS method after using GIS system interpretation and applying tools to determine topographic highs and lows (peaks and depressions) and fracture and to conduct karst feature density analysis (kernel).

Therefore, some adjustments were made, such as considering “non-aquiferous” regions due to geological formations with hydrogeological characteristics different from those of carbonate rocks. Additionally, a correction was applied by considering two topological regions: one with higher altitudes and slopes, and another with lower altitudes and slopes. The regions were differentiated using a clustering methodology using SAGA GIS 7.8.2 software, which divided the topography into two sets. For the lower topographies, a slope value of approximately 10%, estimated by the applied methodology and weighted by back-analysis, was used (Figure 11).

The average recharge rates were approximately 11%, with the Lagoa do Jacaré Formation showing a recharge of 14% and the Serra da Santa Helena Formation displaying the lowest value at 5%. Similar average values (11.4%) were found in studies of the same river basin [14] using a different method of recharge calculation, as well as in other basins with the same geological formations of the Bambuí Group [48].

The steady-state model calibration (Figure 12) was based on data from the SIAGAS database of the Geological Survey of Brazil [27]. Available static water level data from wells were used, assuming that these tests were carried out during well installation. Therefore, test dates were used when available; otherwise, the well installation dates were considered.

The transient calibration (Figure 13) covered the period from 10 January 1949 to 4 February 2022. The initial hydraulic head was derived from the steady-state model, which represents the initial stage of the aquifer before the start of water withdrawal via pumping. Automated instruments provided calibration and recorded data between late 2020 and early 2022, as detailed in Section 3.4 (groundwater monitoring).

3.6. Climate Change Scenarios

The increasing water scarcity conflicts, particularly in the semi-arid region of northeastern Brazil, exacerbated by climate change and uncertainties, have prompted a series of studies. In the states of Ceará and Piauí, an integrated modeling study projecting climate change scenarios revealed that some areas may experience drastic reductions in precipitation, while neighboring regions could see increases in annual rainfall over the next 50 years [50].

In Pernambuco state, northeastern Brazil, similar scenarios have been observed in recent studies. Ref. [51] analyzed the daily rainfall, flow, and meteorological data collected since 1970 in the Tapacurá River Basin. Simulated impacts of climate change under the low-emission scenarios (B1) predicted an increase in temperature and a decrease in rainfall. Projected river flow is expected to decrease by 4.89%, 14.28%, and 20.58% for the periods of 2010–2039, 2040–2069, and 2070–2099, respectively, with a corresponding reduction in groundwater recharge of 13.90%, 22.63%, and 32.91%. Controversially, high-emission scenarios (A2) predict increases in both temperature and precipitation, with water flow expected to rise by 25.25%, 39.48%, and 21.95% for the same periods, respectively, and groundwater recharge increasing by 14.93%, 26.68%, and 11.49%.

In Araripina, Pernambuco, a 50-year historical series of precipitation and temperature data was analyzed. Projections indicate an increase in both maximum and minimum temperatures, coupled with a reduction in rainfall of approximately 2.54 mm per year [52].

Global climate models (GCMs) are essential tools for studying climate change and are widely used to simulate and project changes on global and regional scales. In response to the need for more comprehensive modeling, a federated structure was adopted for the sixth phase of the Coupled Model Intercomparison Project (CMIP6), significantly expanding the number and scope of experiments [53].

Dantas et al. [54] analyzed temperature and rainfall projections using 15 GCMs included in CMIP6 across the following four socioeconomic scenarios (SSPs):

• SSP1-2.6 (optimistic scenario, with low greenhouse gas emissions);
• SSP2-4.5 (stabilization scenario);
• SSP3-7.0 (intermediate scenario);
• SSP-8.5 (pessimistic scenario, with high greenhouse gas emissions).

Projections for the northeast indicate an irreversible increase in average air temperature of at least 1 °C over the course of the 21st century, with up to a 30% reduction in annual precipitation under the regional rivalry and high emissions scenarios. To distinguish the region’s highly variable climate, the Northeast was divided into three subregions: North, East Coast, and Interior, where the study area in northern Minas Gerais is located. Across all scenarios, the temperature is projected to increase. For precipitation, the Interior shows a variation of 10% for all scenarios.

For this research, four adapted scenarios were defined, primarily based on Scenarios 3 and 4 [51] and Scenarios 1 and 2 [54]. For Scenarios 1 and 2, the adaptation involved applying a negative bias, as projections suggest variations in the range of ±10% by 2100. Given the context of the study area, characterized by significant water availability conflicts and its location in the Brazilian semi-arid region, the focus was on assessing the impacts of these reduction rates, as a 10% increase would not significantly alter the current situation. For Scenarios 3 and 4, the adaptation involved averaging the projected increases and decreases in annual rainfall. This approach was chosen to simplify the computational process using annual outputs.

The four scenarios are as follows:

• Scenario 1: 10% reduction in rainfall by 2100, with no increase in groundwater consumption. Adapted from SSP5-8.5 [54];
• Scenario 2: 10% reduction in rainfall by 2100, with an increase in groundwater consumption proportional to that in recent decades. Adapted from SSP5-8.5 [54];
• Scenario 3: 9.5%, 16.5%, and 25% reduction in rainfall for the periods up to 2039, between 2040 and 2069, and from 2070 to 2100, respectively, considering population growth. Adapted from Scenario B1 [51];
• Scenario 4: 16%, 28%, and 12.5% increase in rainfall for the periods up to 2039, between 2040 and 2069, and from 2070 to 2100, respectively, considering population growth. Adapted from Scenario A2 [51].

4. Results and Discussion

The primary distinction between the proposed climate change scenarios and those outlined in the previous section lies in the variation in rainfall and its subsequent impact on aquifer recharge. To simulate these future scenarios, recharge values were adjusted to observe the aquifer’s water level behavior and any potential drawdowns.

One key observation was that the flow pumped from the watershed, as represented in the numerical model, is expected to stabilize in the coming years. This stabilization is partly due to the lowering of the water level, which may drop below the depth of many wells, rendering them unable to pump water, as they no longer reach the water level. Additionally, it is also important to note that the most karstified zones are located within the first few tens of meters, and that the permeability coefficients tend to decrease with depth. This can further contribute to the stabilization of flows when additional wells with similar average depths are introduced to an already lowered water table. Even with an increase in the number of wells (with average depths similar to the current standard), the overall efficiency of the pumping system diminishes and stabilizes.

To address this, a fifth scenario was proposed, involving an increase in the depth of some wells. In this scenario, rates like those in Scenario 2 were applied, with a 10% reduction in rainfall by 2100 and population growth factored in. Figure 14 compares the results of the five scenarios for the years 2039, 2069, and 2100.

Even with flow stabilization at 3000 m³/h, the local exploitation rates exceed effective recharges, creating a groundwater deficit, particularly in urban areas, as confirmed by the drawdown cone.

In all the scenarios, the potentiometric surfaces related to the water level declines are concentrated in urban areas, with no expansion into the location where the Lagoa do Jacaré Formation limestones outcrop. This is significant because the main recharge zone in the watershed, as indicated by the APLIS method calculations, is not experiencing significant dewatering. The aquifer shows resilience in the short term across all scenarios, except for Scenario 5, with minimal expansion of the drawdown cone until 2039. However, by 2069, the drawdown cone reaches 30 m in all scenarios, with deeper declines in Scenarios 2 and 3 due to reduced rainfall and recharges, alongside population growth. In Scenario 5, the groundwater level drops more than 100 m by this period.

Scenario 1, with a moderate 10% reduction in rainfall by 2100, does not indicate severe aquifer depletion. However, by the 2070s, the groundwater level is estimated to drop by 50 m. It is important to note that the city of Montes Claros, a regional hub, is likely to experience both population and industrial growth, which is not accounted for in this scenario, making its results less probable.

The comparison between Scenarios 1 and 4 highlights how sensitive aquifer sustainability is to population growth. Regardless of rainfall increases or decreases, population growth and increased groundwater demand have more significant impact than rainfall variations at the simulated levels. Scenario 1, without population growth or increased groundwater consumption, shows less depletion by the century’s end than Scenario 4, despite its higher annual rainfall and increased water use.

Scenario 2, with only a 10% reduction in rainfall by 2100, and Scenario 3 both exhibit deep drawdown cones, which significantly affect water levels in urban areas due to reduced recharge and high pumping rates. Even Scenario 4, with increased rainfall at different intervals, shows a water table decline similar to those in Scenarios 2 and 3, confirming that population growth and sustained high pumping rates have a greater impact on groundwater sustainability than rainfall fluctuations.

By 2100, the difference in water level decline between Scenarios 1 and 3 is around 15 m, as shown in Figure 15. Additionally, there is a ca. 50% reduction in the average flow of the Vieira River, leading to zero flow during prolonged droughts. Scenario 5 exhibits the most significant depletion, with groundwater levels dropping as low as 140 m in urban areas.

From the perspective of climate change, variations in rainfall are highly complex to project, as they fluctuate greatly across different scenarios and regions. However, all IPCC scenarios agree on the prediction of substantial temperature increases. In this study, recharge was calculated using APLIS, which, as demonstrated in regional studies (e.g., [14]) and similar geological contexts (e.g., [49]), has shown strong consistency. However, APLIS does not account for evapotranspiration in its recharge calculations, potentially underestimating future recharge reductions, which could accelerate with rising temperatures and consequently increasing evapotranspiration.

In Scenario 5, the deepening of the water table and the expansion of the drawdown cone are evident, driven by increased water extraction, with an average of 5000 m³/h of pumping. Although the scenarios are approximate, the deepening of wells in the region is inevitable unless external water sources are provided. Thus, Scenario 5 is a highly likely outcome, with flows potentially reaching 5000 m³/h m³/h, compared to 2000–3500 m³/h in Scenarios 1 to 4.

The effective recharge in the Vieira River Basin is 70 hm³/year, while groundwater exploitation in the proposed scenarios ranges from 26 hm³/year to 44 hm³/year in Scenario 5. Although the watershed does not currently have a negative overall water balance, there is a clear deficit between recharge and groundwater extraction in the urban area of Montes Claros, as evidenced by the growing drawdown cones in all scenarios.

Given these impacts, it is essential to strengthen and enforce stricter regulations, limiting any significant increases in pumping and curbing the operation of illegal wells, which are common in the region. Additionally, the watershed needs a robust and specific water resource management system to ensure that new monitoring data can feed back into the model, defining sustainable extraction limits.

Land use and occupation in karst terrain is a topic of significant societal and academic debate due to the inherent risks associated with carbonate rocks, such as uncontrolled surface runoff and sinkholes, which pose threats to infrastructure and public safety. Lowering the water table in karst areas can exacerbate these risks, potentially leading to subsidence, collapse, surface damage, and even minor seismic events. Similar situations have been observed in regions such as China [55], Italy [56], and Brazil in similar geological settings [57]. This underscores the urgent need for new monitoring strategies and most importantly, targeted water resource management for the city of Montes Claros.

5. Conclusions

This study presents an integrated approach to developing a conceptual and numerical groundwater flow model for the sensitive Vieira River Watershed, a karstic region in a semi-arid climate. It assesses the watershed’s response to various climate change scenarios, including changes in rainfall, groundwater consumption, and population growth.

The methodology emphasizes two key aspects: (i) the importance of considering multiple factors that influence the basin, such as topography, geology, and hydrogeology, and (ii) the value of integrating all available data, including both publicly accessible and newly acquired field measurements. While the model calibration shows good correlation with the observed points, there are still opportunities for further refinement, particularly as model calibration is one of the most challenging aspects in evaluating karst systems. Given the large area, the lack of data, and the significant complexity of aquifer permeability in this basin, this research represents an important initial contribution to ongoing studies, especially as it is the first model ever developed for this region.

In terms of the watershed’s response to climate change, the findings indicate the following: (i) Population growth and sustained high pumping rates have a more significant impact on groundwater sustainability than the forecasted rainfall variations; (ii) there is a deficit between recharge and groundwater extraction in the urban area of Montes Claros, demonstrated by growing drawdown cones; and (iii) deeper wells will be required unless water demand changes or external sources are provided. These findings underscore the urgent need for stricter regulations, limitations on pumping, curbing illegal wells, and implementing a robust water resource management system to monitor and enforce sustainable extraction limits. Therefore, effective water resource management is crucial to ensure water availability and public safety in changing dynamic karst environments.

In addition to enhancing the local and regional monitoring system to advance the presented numerical models and future projections, it is important to evaluate the underground connectivity of watersheds. This will help define and monitor the hydrogeological basin more accurately. Furthermore, applying different recharge methodologies is crucial to assess the impact of rising temperatures, as projected by various climate change models, which also influence evapotranspiration and consequently, the water balance. At last, it is recommended for future simulations to apply linear models of temperature evolution and consequently, the relative increase in actual evapotranspiration, as a decreasing rate of recharge in the mass balance for better representation of climate effects up to the year 2100.