Regional Groundwater Modeling of the Guarani Aquifer System

The Guarani Aquifer System (GAS) is a strategic transboundary aquifer system shared by Brazil, Argentina, Paraguay and Uruguay. This article presents a groundwater flow model to assess the GAS system in terms of regional flow patterns, water balance and overall recharge. Despite the continental dimension of GAS, groundwater recharge is restricted to narrow outcrop zones. An important part is discharged into local watersheds, whereas a minor amount reaches the confined part. A three-dimensional finite element groundwater-flow model of the entire GAS system was constructed to obtain a better understanding of the prevailing flow dynamics and more reliable estimates of groundwater recharge. Our results show that recharge rates effectively contributing to the regional GAS water balance are only approximately 0.6 km3/year (about 4.9 mm/year). These rates are much smaller than previous estimates, including of deep recharge approximations commonly used for water resources management. Higher recharge rates were also not compatible with known 81Kr groundwater age estimates, as well as with calculated residence times using a particle tracking algorithm.


Introduction
The Guarani Aquifer System (GAS) is one of the largest groundwater reservoirs in the world, representing a transboundary system that encompasses four countries in Latin America: Argentina, Brazil, Uruguay and Paraguay ( Figure 1). Covering an area of about 1.2 million km 2 , the GAS system is regulated by an important international treaty, the Guarani Aquifer Agreement. While the average thickness of GAS is nearly 250 m, its maximum thickness exceeds 800 m [1]. The GAS comprises a sequence of sandy layers of Triassic-Jurassic age, deposited in continental, aeolian, fluvial and lagoon environments, above a regional erosional surface (dated to 250/Ma) and below an extensive layer of Cretaceous basalts (dated to 145-30/Ma) in the basins of Paraná (Brazil and Paraguay), Chaco-Paraná (Argentina) and North (Uruguay) [2]. The Guarani Aquifer System is confined by underlying pre-GAS deposits and overlying post-GAS deposits (Figure 1), except in outcropping areas where the aquifer ranges mostly from confined to semi-confined.
Despite the continental scale of GAS, groundwater recharge is restricted to narrow outcrop zones, mainly along its western and eastern borders ( Figure 2). Estimates of groundwater recharge in the outcrop zones were calculated by several authors (e.g., [3][4][5][6]). However, most groundwater recharge in the outcrop zones is discharged into local watersheds, with only minor amounts reaching the confined part of GAS to form deep recharge. Other studies [7,8] showed that deep recharge comprises only 10 to 15 mm/year, which represents a scant 1 to 2% of annual rainfall. Simplified geologic map showing pre-GAS and post-GAS sediments, GAS outcropping areas and overlying basalts (modified from [9]).
Because of its strategic importance, several studies have focused on possible actions to improve GAS water management, sustainable use and groundwater protection [1,2,[10][11][12][13]. The extensive area underlain by GAS has a population of about 15 million, whereas 90 million inhabitants are indirectly benefitting from GAS exploitation [10]. Thus, the need for reliable potable water and industrial supplies (of low-cost treatment) is likely to grow significantly, especially when considering also climate change [10]. According to [10], GAS exploitation exceeds 1.0/km 3 /year, being 93.6% in Brazil (of which about 80% occurs in São Paulo State), 2.8% in Uruguay, 2.3% in Paraguay and 1.3% in Argentina. Approximately 80% of total water extraction is used for public water supply, 15% for industrial processes and 5% for geothermal spas [10].
Geological models and structural implications in terms of groundwater flow and hydrochemistry of GAS were studied by [14][15][16][17], among others. Numerical models are crucial tools for water resources management of continental scale aquifers. Despite the importance of GAS as a freshwater resource and the capability of numerical models, only a few studies were previously carried out of the aquifer system, ranging from local scale [3,4,18] to regional scale [7,17] investigations. The modeling studies commonly calibrated the recharge fluxes by considering direct recharge in outcropping areas [3,19]. To enhance our understanding of the precise flow dynamics of the Guarani aquifer system and obtain more reliable estimates of prevailing groundwater recharge rates, we created a three-dimensional numerical flow model of the aquifer. The model used for this purpose was obtained using the finite-element based FEFLOW groundwater flow simulator of [20]. In the following we first describe our conceptual model of the GAS system, the adopted numerical model, and results of the calculations in terms of piezometric levels, recharge estimates, and comparisons with travel time estimates within the aquifer system using the particle tracking option of FEFLOW.

Conceptual Model
Our conceptual model assumes that the GAS comprises layered sandy units and recovering Permian aquitards (the Pre-GAS units), while being confined by basalts of varying thickness, except in outcropping areas. While several studies suggested some hydraulic connection between GAS and confining or underlying layers, this issue remains poorly defined. Since precise data quantifying water exchange by inter-layer vertical flow are lacking, we assume for now that this exchange is negligible. Henceforth, our conceptual model assumes neither recharge from nor discharge to confining basalts.
Despite the large number of studies hypothesizing the structure and compartmentalization of GAS in segmented blocks [15,17,21], we regarded GAS as an entirely connected block based on observed groundwater flow continuity at the regional scale. The depositional context of GAS encompasses aeolian, fluvial and lagoon environments. Despite local differences, the geological history of GAS deposition points to an essentially sandy makeup, while being relatively homogeneous from a regional point of view. We conceptualized that the diagenetic history, mainly the pore filling cement, is the main feature controlling aquifer permeability. Figure 3 shows a schematic of our conceptual model. The hydraulic conductivity (K) is a parameter known to vary with the scale of measurements [22,23], which implies that available local-scale literature values of K may well be inappropriate due to the large scale of the GAS. Based on an evaluation of hundreds of K measurements at different scales and geological contexts, [22] found that as the scale increases, the mean K tends to increase until reaching some plateau, while the variability tends to decrease with scale as the scale increases, more or less following a potency-model. The increase in K occurs until some specific volume (the upper bound of the relationship) is reached, after which K remains constant. An analysis of data presented by [22] showed that K values provided by pumping tests are close to the upper bound. Consequently, we assumed that the pumping tests served as a reliable proxy to upscaled values of K for applications at regional scales.
Tens of pumping tests conducted in several areas of GAS [3,[24][25][26][27] showed a relatively narrow range in K values, varying between 6.4 × 10 −6 m/s and 8.1 × 10 −5 m/s, with the vast majority falling between 1.0 × 10 −5 m/s and 5.0 × 10 −5 m/s. Accordingly, we used these K values as constraints during the model optimization. Furthermore, based on the assumption that GAS potentiometry is controlled by the structural contour of the GAS basement, rather than high variability in the hydraulic conductivity, we conceptualized that K values are relatively homogeneous, with little variation of this parameter at the regional scale.

Numerical Model
To verify the conceptual model assumptions and to calculate GAS regional flow behavior, a three-dimensional steady-state model was constructed using the FEFLOW finite element software [20]. Two-dimensional transient models using hydraulic head fluctuations, such as done by [7], helped to better understand the local impacts of different exchange rates, but this approach did not lead to realistic estimates of long-term regional flow velocities. Our 3D GAS model, on the other hand, focused on recharge rates that are consistent with estimated groundwater ages. Since our timescale is hence millennial (up to 1 Ma), head fluctuations during the last several decades will not control our model. Therefore, a steady-state approach is consistent and suitable for this purpose.
The model boundary extended to the entire aquifer system, encompassing an area of 1,087,560 km 2 . The domain was divided into two layers to reflect the existence of GAS and Post-GAS units (the latter being the younger geological units present in the basin, predominantly confining basalts).
The three-dimensional domain was discretized into 337,392 triangular elements and 273,864 nodes, with irregular spacings in both the horizontal and vertical directions. The element size averaged 6 km 2 , but with mesh refinements along the borders and larger elements in the central region.
The domain bottom was assumed to be a no-flow boundary. In addition, seepage face boundary conditions were introduced along the western and southeastern borders to account for the discharge areas, thus allowing free outflow from the aquifer. Recharge rates were represented using a fluid-flux boundary condition along the top elements over the outcropping areas. The model was calibrated using trial-and-error methods by adjusting hydraulic conductivity values such that simulated levels during steady-state flow agreed with the monitored hydraulic heads from 49 wells [9].
We emphasize that underneath the major rivers (e.g., Paraná River, Figure 3) there are approximately 1000 m of confining basalts with no reliable subsurface information to support the use of the rivers in any calibration. The pumping from wells represents very low rates compared to overall aquifer flow rates, albeit highly concentrated. For this reason, they are expected to locally impact the heads, especially in the unconfined portions, but are negligible for the extremely large confined portion of the aquifer system. In terms of ages (several hundreds of thousands of years), the impact is also negligible. To reduce uncertainties related to non-uniqueness of the model solution, we constrained the model parameterization using available data about groundwater age [28].

Piezometric Levels
The adjusted steady-state flow model was evaluated by comparing observed and simulated piezometric levels ( Figure 4A,C). The final model simulation produced a correlation (r) of 0.9632 and a normalized root mean square error (NRMS) of 6.04% (55.79 m) ( Figure 4B,C). Model results accordingly replicated groundwater flow patterns in terms of the hydraulic gradients and the main flow directions.
Smooth gradients were apparent near the center of the simulated aquifer, consistent with observed potentiometry data, with the potentiometric surface being essentially flat there [29]. Higher hydraulic gradients occurred close to the recharge areas. Regional groundwater flow was directed mainly from the north and northeast radially towards the center and then southwards. Regional flow near the western border was found to be somewhat fragmented by local recharge/discharge systems.

Hydraulic Conductivity
Not enough GAS hydraulic conductivity (K) data were available to construct a regionalized distribution map using interpolation methods such as kriging. Hence, the K zones were distributed according to the described lithology compartments and correlated porosities [29,30]. Hydraulic conductivities varied from about 2.5 × 10 −6 m/s near the outcrops, to 5 × 10 −5 m/s in the central and southern portions ( Figure 5A). Although the K zones were distributed in a relatively simplistic manner (e.g., by using only three zones), they were geologically consistent and sufficient to reproduce the observed regional flow patterns. The vertical K values were one order of magnitude less than the horizontal conductivities. The fillable porosity of GAS was set uniformly at 0.2 and that of the post-GAS at 0.0001.
The estimated hydraulic conductivities were close to reported literature values, which showed little variability between a minimum of 1.2 × 10 −6 m/s and a maximum of 5.3 × 10 −5 m/s [7]. Recent pumping tests in São Paulo State yielded K values of about 2 to 3 × 10 −5 m/s [27], similar to values reported earlier [31]. Calibration of local scale models (100 km 2 to 5000 km 2 ) produced mean K values of 2.8 × 10 −5 m/s (as compiled by [7]). Our calibrated hydraulic conductivity values hence compared well with both pumping tests and local-scale groundwater flow models. This shows that there was no need to implement a much denser discretized K distribution as is often done for larger-scale flow models [32].

Water Budget and Recharge Estimates
The recharge zones in the numerical model were distributed along the GAS outcrop areas ( Figure 5B), whereas the discharge zones were mostly located along the western and southeastern borders ( Figure 5C). The model was calibrated by using four recharge zones: a western zone (W) with 6 mm/year, a northeastern zone (NE) with 2 mm/year, an eastern zone (E) with 6 mm/year, and a southern zone (S) with a recharge of 3 mm/year. Calculated recharge for the entire domain was approximately 0.6 km 3 /year, which represents a mean average recharge rate of 4.9 mm/year. Several studies reported much higher recharge rates, generally close to 2-4% of the mean annual precipitation rate of about 1500 mm/year [7]. For example, [33] estimated a recharge rate of 41 mm/year, whereas [7] obtained 35.2 mm/year using a more complex transient flow model (TRANSIN). Our calculated recharge is similar to that of deep recharge into GAS (i.e., water that effectively reaches the confined portion of the aquifer system). In a discussion of the limitations of estimating this deep recharge, the GAS Project [9] hypothesized a volume close to 0.8 to 1.4 km 3 /year, which translates to about 6.5 to 11.5 mm/year. [34] considered deep recharge to represent less than 1% of annual rainfall (about 1 km 3 /year), but this rate could become as high as 5 km 3 /year or more with increasing pumping. Considering that some GAS portions in our model are actually discharging locally, the deep recharge volume may be even smaller.
Much higher local recharge rates have been reported for the outcrop areas, such as by [8] who also showed large discharge volumes into local drainage systems. Direct recharge at various GAS sites has been estimated mostly using water balance models [3]. Hydrographs separation/analysis techniques have been used frequently to differentiate among water sources (such as by [5] who focused on baseflow), but also the water table fluctuation (WTF) method [35] to estimate recharge from groundwater level fluctuations [36]. However, the WTF relies on having good estimates of the specific yield, which is often problematic, while the recharge rates may be overestimated significantly because of the effects of entrapped air in shallow unconfined aquifers overlaying the areas being investigated [37]. For these reasons it is important to consider only recharge rates that effectively contribute to the aquifer system in terms of the overall groundwater balance. Proper management requires hence a focus on the relatively small water volumes that effectively recharge an aquifer, which may not include water that rapidly discharges again and only influences shallow groundwater level fluctuations and/or river outflows.

Groundwater Age Analysis
We further compared our simulations with results from groundwater dating using 81 Kr as described by [28]. For this we performed forward streamline and travel time calculations using the particle tracking option of the FEFLOW software. Particle seeds were positioned on the recharge areas as shown in Figure 6A. The results revealed ages ranging mostly from 100 ky to 500 ky, with the oldest waters reaching 1 Ma. These estimates are consistent with most of the 81 Kr dating results [28]. Some discrepancies occurred when calculating backward streamlines from the age estimates ( Figure 6B), which may reflect mixing with Pre-GAS units, as suggested by [38], or Post-GAS units (e.g., from flood basalts), as previously discussed by [39].
The particle tracking schemes produced estimates of the advective age (also known as the true groundwater age), which is the time needed for a water particle to move between two points purely via advection. Particle tracking is useful since it reflects the movement of groundwater due to hydraulic forcing as a result of hydraulic gradients, the hydraulic conductivity and porosity, and the recharge and/or discharge rates [40]. The advective age can be used to estimate groundwater flow rates [41], and to more directly infer recharge rates [42,43].
Our analysis of 14 wells showed three krypton ages much older than expected, two in the northern part (728 and 834 ky) and one to the south (777 ky). These results may reflect mixing with Pre-GAS units that contain much older water, and located relatively close to the outcropping zones. The mixing processes in the northern portion is discussed by [38], who also included two other samples (showing ages of 566 and 531 ky) which may have been affected by some mixing. In summary, our regional age pattern agrees very well with most of the krypton ages. It is important to note here that we are comparing general flow patterns within a regional context, and not trying to calibrate precisely the ages. Thus, since krypton dating itself has considerable error margins, while only a few poorly distributed GAS samples where available. For comparison purposes, we found that if we increased the recharge rates ten times (to become similar to previous direct recharge estimates), values of K needed to be increased one order of magnitude to properly calibrate the hydraulic heads. This scenario then resulted in GAS groundwater ages of only about 50 years, which is unrealistic, even when considering mixing and preferential flow. We emphasize again that our model is regional and thus not intended for evaluating local-scale issues in specific areas. Local problems are best addressed using local models that could provide better support for decision-making processes. Inter-formation mixing and groundwater flow through geological structures are certainly present and may play an important role locally [16,38,39]. However, our results show that these exchanges are not requisite to explain overall regional flow patterns of the GAS. Therefore, we assume that the mixing waters do not represent an important fraction of the GAS recharge/discharge volumes.
Based on mixing simulations, [38] demonstrated very complex interactions of the Guarani Aquifer System with both underlying and overlying units in terms of GAS receiving and delivering water. The exact level of GAS interactions with overlying and underlying units remain largely unsolved. Because of the complexity of this issue, we focused in our study on a relatively simplistic scenario by assuming limited interactions of the Guarani aquifer system with other units. The capability of our model to reproduce the hydraulic head distribution and most of the groundwater age estimates provides confidence that the assumed scenario is feasible. Despite the disagreement about the level of mixing with other aquifers, our simulations and those by [38] were surprising similar.

Model Discussion
Because direct measurements of discharge to the Paraná River were lacking, the potential contribution of GAS discharge to the river is difficult to confirm directly. Contrary to the numerical model of [7], our model assumes that the main rivers crossing the basin (e.g., the Paraná, Paraguay and Uruguay rivers) do not represent discharge zones due to the thickness of basalts (>1000 m) beneath them. As described by several authors [44][45][46], the basalt floods consist as persistent alternating thick intervals with extremely low permeabilities (interior flows) and thin permeable intervals (top flows). Given this situation, the significant thick basalt layers serve as very effective barriers preventing seepage flow from GAS to the rivers.
In a practical sense, including the main rivers crossing the basin as boundary conditions may help model calibration by modulating the simulated potentiometry. Our model, however, reveals that the rivers are unnecessary to condition flow in the GAS, and that the model can be calibrated effectively without these features, thereby reinforcing the likely absence of connectivity of GAS with the rivers.
The regional hydraulic head and distribution of groundwater reflects the natural flow conditions stablished by long-term water flow during the last hundred millennia. Our model assumes a steady state flow regime, thereby aiming to reproduce these conditions. Some wells in São Paulo State experienced a noticeable, but localized, drawdown in the last four decades. We did not used this exploitation because the effect of pumping is very restricted, without any evidence that this affected the entire GAS (specially the confined portion). The hydraulic heads calculated with our model, without groundwater extraction, hence will probably not adjust to heads at locations with noticeable drawdown, thus producing a certain error. The potential error ( , in %) of assuming a lack of groundwater abstraction in our model may be estimated using the following expression: where h calc is the calculated hydraulic head [L], h obs is the observed hydraulic head [L], h max is the maximum hydraulic head observed in the studied aquifer [L], and h min is the minimum hydraulic head observed in the aquifer [L]. Assuming an extreme case, with a 50 m-drawdown induced by excessive exploitation, the estimated error will be less than 7% as estimated using Equation (1). This indicates that ignoring groundwater exploitation in our model will not produce a significant error when considering the regional GAS scale. Because the generally non-uniqueness nature of a groundwater flow model, distinct sets of hydraulic conductivities and groundwater recharge rates often produce similar model calibrations, leading to dissimilar water balances and high uncertainties. Our strategy to reduce uncertainties related to the hydraulic conductivity was to carefully compile available K values of the GAS and to force their values in the model to remain close to median of the compiled K values. Since GAS is composed of thick sandy units deposited under aeolian to fluvial environments, the K is expected to be relatively homogeneous at the regional scale. The ability of our model to produce reliable potentiometry using little variability in the hydraulic conductivity field confirmed our assumption that the hydraulic gradient and transmissivity are controlled mostly only by the structural contour of the basement.
The use of groundwater age data as an additional calibration criterion served as an important strategy to reduce the uncertainties related to model parameterization. Because of the direct dependence of advective age on groundwater recharge, the degrees of freedom for producing the calibrated model were strongly restricted. Thus, our model possessed a lower uncertainty as compared with models calibrated exclusively with hydraulic head data.
Several studies of watersheds located in GAS outcrop zones have demonstrated that most of the groundwater recharge discharged along the local rivers, with an exceptionally low fraction contributing to the confined portion as deep recharge [3][4][5]18]. These findings closely match our results that groundwater recharge comprises only about 1% of total precipitation. Consequently, the confirmation of the local outcrop estimates using our 3D model covering the entire GAS has important implications in terms of water management planning. Water extraction by well supplies in a single year comprised 17% of annual groundwater recharge, considering the recharge rates that discharge locally in the watersheds. However, 80% (0.8 km 3 ) of the extracted volume was concentrated in the populous region of São Paulo State, in Brazil. Given this scenario, and despite our model being regional with extraction not affecting the entire GAS, we reinforced local studies that the GAS in the São Paulo area is critically overexploited, leading to untenable drawdown in this particular (and socially important) region. Long-term predictions of the effects of intense groundwater exploitation in this region need to be carried out to stablish the limit of groundwater volumes that can be extracted in a sustainable manner.

Conclusions
Because of the strategic importance of the Guarani Aquifer System and the increased demand for water, more precise operational criteria should be defined for the sustainable use of GAS groundwater within the context of water management. The strategies of management should be supported by a correct diagnose of the flow dynamics and, especially, prevailing recharge rates. Our results demonstrate that the recharge rates that effectively contribute to the regional water balance of GAS (about 0.6 km 3 /year) is much smaller than previous estimates. This is true also when comparing our results with deep recharge approximations commonly used for water resources management. Local discharge rates occurring in the drainage systems along outcrop zones are shown to be very high by looking at the baseflow contribution in river hydrographs, thus only a small portion of the infiltrating water moves deep and reaches the confined GAS.
Comparing Kr-measured groundwater ages with particle tracking and travel times showed that higher recharge rates lead to unreasonable residence times. Moreover, when one considers smaller specific yield (S y ) values, which is realistic and possible, the resulted recharge rates would be even less. This because with less fillable pore volume, less water can be stored, leading to reduced recharge rates. Having krypton ages older as expected is evidence of some inter-formation mixing, especially with older units. Still, this water exchange, as well as aquifer structures influencing preferential flow, were not necessary to reproduce the GAS regional groundwater flow patterns.
For this study we used a hydrogeological model to focus on regional aspects of GAS, with a broad view regarding aquifer homogeneity and continuity associated with age estimates. Our concern was especially to obtain better estimates of recharge from a regional water balance perspective. Local analyses may be needed in critical areas to test some key issues, such as inter-formation mixing, structural controls and pumping impacts. Such studies are needed to confirm our overall findings, as well as to beneficially use our outcomes for local scale modeling and management approaches.