Conceptual and numerical modeling approach of the Guarani Aquifer System

In large aquifers, relevant for their considerable size, regional groundwater modeling remains challenging given geologic complexity and data scarcity in space and time. Yet, it may be conjectured that regional scale groundwater flow models can help in understanding the flow system functioning and the relative magnitude of water budget components, which are important for aquifer management. The Guarańı Aquifer System is the largest transboundary aquifer in South America. It contains an enormous volume of water; however, it is not well known, being difficult to assess the impact of exploitation currently used to supply over 25 million inhabitants. This is a sensitive issue because the aquifer is shared by four countries. Moreover, an integrated groundwater model, and therefore a global water balance, were not available. In this work, a transient regional scale model for the entire aquifer based upon five simplified, equally plausible conceptual models represented by different hydraulic conductivity parametrizations is used to analyze the flow system and water balance components. Combining an increasing number of hydraulic conductivity zones and an appropriate set of boundary conditions, the hypothesis of a continuous sedimentary unit yielded errors within the calibration target in a regional sense. The magnitude of the water budget terms resulted very similar for all parametrizations. Recharge and stream/aquifer fluxes were the dominant components representing, on average, 84.2 % of total inflows and 61.4 % of total outflows, respectively. However, leakage was small compared to stream discharges of main rivers. For instance, the simulated average leakage for the Uruguay River was 8 m3 s−1 while the observed absolute minimum discharge was 382 m 3 s−1. Streams located in heavily pumped regions switched from a gaining condition in early years to a losing condition over time. Water is discharged through the aquifer boundaries, except at the eastern boundary. On average, pumping represented 16.2 % of inflows while aquifer storage experienced a small overall increment. The model water balance indicates that the current rate of groundwater withdrawals does not exceed the rate of recharge in a regional sense.


Introduction
The Guaraní Aquifer System, hereafter GAS (also known as SAG from its Spanish and Portuguese name), is the largest transboundary aquifer system in South America.It extends for some 1.2 million km 2 over four countries: 70 % in Brazil, 19 % in Argentina, 6% in Paraguay, and 5 % in Uruguay (Fig. 1).The aquifer is formed by sandstones and confined by basalts in about 90 % of its extent.Sandstones outcrop along aquifer edges, deepening toward the center of the basin, where they can reach a maximum thickness of some 600 m and depths of 2200 m.
These data point to several peculiarities of the GAS.First, it contains an enormous volume of water, which makes it appealing for groundwater pumping.But because it is not well known, it is difficult to assess the impact of exploitation.This is a sensitive issue because the aquifer is shared by four Published by Copernicus Publications on behalf of the European Geosciences Union.9) Heine (2008); (10) Rabelo and Wendland (2009).
countries.These peculiarities are not specific to GAS.Other large aquifers relevant for their considerable size are the High Plains Aquifer in the USA (Luckey and Becker, 1999), the Nubean Aquifer shared by Egypt, Chad, Sudan and Lybia (Robinson et al., 2007), the Great Artesian Aquifer in Australia (Habermehl and Lau, 1997), the Yrenda-Toba-Tarijeño Aquifer System shared by Argentina, Bolivia and Paraguay (UNESCO-IHP, 2009), and the Navajo Aquifer System in the USA (Heilweil et al., 2002).
To address these issues and aim towards sustainable management and development of the GAS, the Global Enrironmental Facility (GEF) financially supported the four countries to develop the "Proyecto para la Protección Ambiental y Desarrollo Sostenible del Sistema Acuífero Guaraní -PSAG" (Environmental Protection and Sustainable Development Project for the GAS).The PSAG was a multidisciplinary scientific effort and regional collaboration requiring good cooperation between various disciplines.The project revealed some of the difficulties associated to the study of regional aquifers.
In this article, like in the recent work of Barthel (2011), the term "regional scale" refers to areas of approximately 10 5 -10 6 km 2 in size.Some of the largest aquifers in the world have been studied for diverse purposes and with different modeling approaches.For instance, the US Geological Survey and the Oklahoma Water Resources Board, USA, conjunctively developed a groundwater flow model of the High Plains Aquifer to be used for allocating the amount of water withdrawn from the aquifer (Luckey and Becker, 1999).The Nubean Sandstone Aquifer has been studied for many years.A preliminary modeling effort on this aquifer was reported by Heinl and Brinkmann (1984), who used a finite element model to address various basic questions regarding the dynamics of the aquifer.In a recent work, Gossel et al. ( 2004) presented an integrated GIS-based groundwater flow model for the Nubean Aquifer intended for improving previous modeling efforts.A transient groundwater model was constructed for the Great Artesian Basin for management purposes (Welsh, 2006).One recent contribution to regional-scale groundwater modeling was given by Michael and Voss (2009), who focused their work on the estimation of regional-scale aquifer properties in the Bengal Basin of India and Bangladesh.In order to accomplish their objective, they combined inverse groundwater modeling using measured heads, model calibration using estimated water ages, and statistical analysis of driller logs.
Several local-scale groundwater models were built in the GAS, either to meet local or state requirements or to analyze a particular behavior/characteristic of the aquifer.Here, the term "local scale" refers to models covering areas of approximately 10 3 km 2 .A thorough study of the aquifer within Paraguayan territory culminated with a multi-layer, steady state groundwater flow model (Vassolo, 2007).Within Brazilian territory, Heine (2008) built a MODFLOW model to quantify recharge for management purposes around the city of Ivoti in southern Brazil, while Rabelo and Wendland (2009) assessed groundwater recharge and water fluxes in the state of Sao Paulo, Brazil, through a numerical, finite element model, covering over 5000 km 2 of outcropping sandstones.Within the PSAG project, four pilot areas were selected for detailed study on the basis of their distinct hydrogeologic conditions and potential groundwater exploitation conflicts.In each of them, a local-scale flow model run under MODFLOW intended for management practices was developed (SNC Lavalin, 2008a, b, c, d).Recently, Gómez et al. (2010) implemented MODFLOW (Harbaugh et al., 2000) on a region located along the Brazilian-Uruguayan border in order to validate a newly proposed multi-layer conceptual model, verify previous recharge estimates and test future exploitation scenarios.The location of these and all other modeling sites/boundaries described in this section are indicated in Fig. 1.
On the basis of previous work by Campos (1998), Vives et al. (2001) developed the first regional scale groundwater model of the GAS, known as "pre-model" because it was Hydrol.Earth Syst.Sci., 17, 295-314, 2013 www.hydrol-earth-syst-sci.net/17/295/2013/ built upon scarce geologic and hydrogeologic data.It was a finite element, two-dimensional, steady state model, which extended over 60 % of the currently identified aquifer area.Regional geologic structures that may condition groundwater flow were included in the model, and a handful of short stream reaches in outcropping areas to the east were explicitly simulated.Nonetheless, many questions remained unanswered or partially addressed due to the simplifications set forth during the modeling process, for example: Where is the south-western boundary of the aquifer located?What are the most likely discharge zones in the southern sector?Could reaches along the Paraná and Uruguay rivers be potential discharge zones?What is the dynamics of the system along the aquifer western boundary?How important are stream/aquifer interaction processes in outcropping areas compared to other mass balance components?Are there local recharge/discharge systems identified within a regional context?More importantly, the overall water balance of the aquifer is not known.Field data generated during the PSAG allowed revisiting the conceptual model, proposing a new southern boundary for the aquifer and postulating its interaction between numerous streams (Gastmans et al., 2012).The newly proposed conceptual model was numerically validated by Vives et al. (2008), who also hypothesized and numerically tested that the GAS may discharge through selected reaches along the Uruguay and Paraná rivers.Yet, the magnitude of these discharges is unknown.
All groundwater flow models built so far on the GAS, either at local or regional scale, were run under a steady state regime, which, in turn, limits their use as learning tools as well as management tools.In this work, a transient regional scale model covering the full extent of the GAS is presented.The new model complements its steady state predecessor by Vives et al. (2008), rendering the new modeling approach more informative in the process of enhancing the current hydrogeologic understanding of the aquifer and its potential use as a management device of subsurface resources.
The model was used to test whether the aforementioned questions can be modeled consistently with aquifer head data.It was also instrumental for evaluating water balance components for the entire aquifer, emphasizing the role of processes such as the stream/aquifer interaction as a leading discharge mechanism in outcropping areas, and the feasibility of some aquifer discharge in the southern portion of the aquifer, a hypothesis that was disregarded in previous studies.
The GAS is located between 16°S and 32°S latitude and 47°W and 56°W longitude underlying the Río de La Plata drainage basin in South America (Fig. 1).Ground elevations vary from 1700 m a.s.l.(meters above sea level) in the southeastern border down to approximately 30 m a.s.l.within Argentinean territory.Due to its considerable extent and variations of relief, diverse climates are identified.Mean annual, altitude-dependent precipitation shows a southward gradient, from 2000 mm in the north to 1400 mm in the south while mean temperature is above 20 • C almost everywhere.Mean annual evaporation has been estimated to be around 60 to 70 % of the annual precipitation.
The GAS sedimentary sequence consists of aeolian, and fluvial weakly-cemented sandstones beds of Upper Jurassic-Lower Cretaceous age deposited in parts of the tectonic Paraná Basin and Chaco-Paraná Basin (Araújo et al., 1999).Sandstones range in thickness from a few meters in outcropping areas along western and eastern aquifer boundaries, to more than 600 m at the center of the basin.Upper Cretaceous basalt flows as thick as 1500 m and varying degrees of fracturing/fissuring cover 90 % of sandstone deposits.The stratigraphic sequence completes with Quaternary, non-uniformly distributed sediments.Figure 2 shows a simplified geologic map and transverse and longitudinal geologic profiles (Foster et al., 2009).There has been intense debate as to whether the aquifer can be regarded as a single, continuous unit or is it actually separated in geologic compartments that may or may not introduce regional flow discontinuities (Ferreira, 1982;Campos, 1998;Soares, 2008).At regional scale, the GAS can be conceptualized as a sedimentary formation, spatially continuous, composed of sandstones, confined by underlying pre-GAS deposits, and overlying post-GAS deposits (see Fig. 2), except in outcropping areas.The GAS is assumed to range from unconfined to semi-confined in recharge areas.Towards the center of the tectonic basins, it becomes increasingly confined due to the thickening of overlying basalts, which leads to artesian conditions over large areas.Foster et al. (2009) summarized the hydrogeologic framework as follows: "The aquifer occurs in three main 'hydrogeological domains' delimited by two geological structures that have exerted a control on aquifer thickness and depth, and today influence regional groundwater flow: the Ponta Grossa Arch (in the north of Paraná State-Brazil), which forces groundwater to flow from east to west in São Paulo State, -Brazil, the Asunción-Rio Grande Arch, which divides the portion south of the Ponta Grossa Arch into two semi-independent sedimentary basins -the Central Paraná and the south-western Chaco-Lower Paraná.The GAS is also affected by many tectonic structures and crossed by numerous volcanic dykes, but despite these important discontinuities at local scale it is considered to be  a "continuous groundwater body" across the entire region".
Regional flow is from northeast to southwest.In outcropping areas, predominantly-recharge regions alternate with predominantly-discharge ones.The latter may discharge regional as well as local flows.
Recharge to the GAS occurs by infiltration of excess rainfall in outcropping regions, which cover approximately 10 % of the entire aquifer extent.Estimates of recharge rates range from 10 % (Reboucas, 1976) to 4 % (Chang, 2001) of mean annual precipitation.So far, there is no evidence of recharge from streams.However, it should not be ruled out.Well withdrawals are considered the main source of discharge from the aquifer.Other sinks, though not well known, include seepage to streams and seepage to underlying/overlying formations.
Average salinity in recharge areas is about 50 mg L −1 , and can be as much as 500 mg L −1 on the southwestern region.With increasing depth and confinement, and following the general flow direction, the groundwater temperature also increases from 25 to 65 • C, an increment that has been mostly attributed to the effect of the normal geothermal gradient.A temperature effect on the aquifer's hydraulic conductivity (permeability) may be expected as a result of changing kinematic viscosity.The isotherm map of the GAS developed by Gastmans et al. (2012) is very informative, showing the highest values located within Brazilian territory; nonetheless, about 65 % of the aquifer has temperatures between 25-45 • C.

Conceptual hydrogeologic model
A new hydrogeologic conceptual model and a numerical model for the entire aquifer were end products of the multidisciplinary work within the PSAG.Gastmans et al. (2012) put together a revised, much-improved conceptual model that incorporates many aspects that were overlooked in previous versions.As mentioned previously, from a regional point of view, the GAS was defined as a spatially continuous sedimentary formation, composed of sandstones, confined by underlying pre-GAS deposits, and overlying post-GAS, i.e. basalts and quaternary deposits (see Fig. 2), except in outcropping areas.
The potentiometric map shown in Fig. 2  Regionally, groundwater flows from recharge to discharge areas, presenting a directional trend from the northeast and northwest toward the center of the sedimentary basin and then south.This regional flow pattern is influenced by the tectonics of the geologic basins and its evolution, as pointed out by Araújo et al. (1999).The potentiometric map reflects the presence of mega-structures.Hydraulic gradients are steeper at or near outcropping areas, with values as high as 3 to 5 m km −1 in the northeast, 2 to 3 m km −1 at eastcentral locations and 1.5 to 2 m km −1 at the west-northwest.They decrease toward the center of the sedimentary basin as the aquifer deepens.In spite of the presence of structural discontinuities at different spatial scales, groundwater flow continuity still persists at regional scale (Gastmans et al., 2012;Foster et al., 2009).Hydraulic connection across confining layers has been poorly addressed, though it is likely.At present, available data are limited for assessing and quantifying inter-layer vertical occurrence, its direction and its magnitude.Based on a handful of deep wells that tap underlying Permian deposits, Gastmans et al. (2012) postulated that the GAS lies over an erosive basal surface extended over parts of the Paraná sedimentary basin, putting the GAS in contact with argillaceous, low permeability units in the north; silty-clayed, more permeable formations at the center and silty-sandy formations of moderate permeability in the south.Neuzil (1994) reported argillaceous permeability for both lab test and regional scale studies, relating porosity from different materials to permeability to identify possible trends.For clayed sandstone, this author suggested a hydraulic conductivity range between 8.6×10 −5 and 8.6×10 −11 m d −1 for flow parallel and normal to bedding, respectively, which are very low for a quantifiable GAS/pre-GAS hydraulic connection.Considering that interlayer flow may be controlled by the less permeable material, the magnitude of vertical flows between sandstones and pre-GAS sediments may be assumed negligible.
Water quality data may help to elucidate the interlayer flow connectivity issue.This type of data is scarce and some of dubious quality to extract definite conclusions about the existence of vertical density gradients; however, that possibility is feasible given the great extension of the aquifer.Manzano and Guimaraens (2009) exhaustively analyzed background hydrochemical information and field data generated during the PSAG, concluding that the aerial distribution of solutes such as TDS, Na, Cl and SO 4 , among others, is dominated by a mixture of GAS waters with saline waters from pre-GAS units.These authors also identified two mixture-dominated areas, suggesting the existence of vertical, upward flows from pre-GAS formations controlled by stratigraphic and structural conditions.Nonetheless, they concluded that to date, the magnitude of those fluxes is completely unknown, per-ceived to be rather small compared to other components of the regional water balance.
A number of authors sought evidence of the hydraulic connection between confining basalts and underlying sandstones, showing divergent conclusive results.Reboucas and Fraga (1988) argued that water flow can predominantly occur along horizontal discontinuities surfaces at the top and bottom of lava flows, i.e. at interflow contacts, and at vertical column disjunctions present at the center of flows.Given this type of flow description, restricted to specific areas within the basalt packet, uncertainties regarding hydrogeologic units interconnections are still great (Lastoria et al., 2007).Most gathered evidence regarding GAS/basalts hydraulic connection belongs to the Sao Paulo and Mato Groso do Sul States in Brazil, located to the NE and NW of the aquifer, respectively.Fernandes et al. (2008) intensively studied basalts within Sao Paulo State, proposing a conceptual model for water circulation within basalts around the Ribeirão Preto area.In a recent work, Fernandes et al. (2012) used hydrochemistry and basalt fractures mapping to investigate whether vertical conductive structures might conduct water.So far, they have found vertical fractures present only on dense layers that do not penetrate into vesicular layers preventing a hydraulic connection with sandstones underneath.Even though Lastoria et al. (2007) provided some evidence of ascending/descending flows within Mato Groso do Sul using hydrochemistry, it would be speculative to extrapolate this condition to the entire aquifer.Hence, the current version of the conceptual model assumes neither recharge from nor discharge to confining basalts.
These conceptual simplifications regarding the layer structure of the GAS imply an essentially two-dimensional flow regime at regional scale.Therefore, the model will be treated as two-dimensional.
One of the great uncertainties of the conceptual model is the location of regional discharge zones and the magnitude of discharge fluxes.Even though local recharge/discharge systems have been clearly identified in outcropping areas along the western and southern boundaries (see Fig. 2), regional discharge may occur through selected portions of the boundary or other sinks.In this work it is assumed that short reaches of the Paraná and Uruguay rivers could discharge GAS water.
Net recharge, i.e. effective rainfall minus evapotranspiration, occurs along outcropping areas along portions of the aquifer boundary.If not intercepted by pumping or streams, some of it may become deep recharge.

Methods
The code TRANSIN used in this study allows simulating groundwater flow and solute transport (Medina and Carrera, 1996;Medina et al., 1996).TRANSIN includes an algorithm for automatic calibration of all flow and transport parameters based on the maximum likelihood method for parameter estimation, as explained in Carrera and Neuman (1986).In essence, the values of hydraulic properties variable in space (and sometimes in time), are calculated based on a previous estimate of the parameters and measured values of heads (and concentrations if transport is solved).
It should be pointed out that, besides data errors and conceptual model uncertainties, there are also errors associated with the numerical method.The numerical error is related to the discretization, the numerical procedure and round-off errors.TRANSIN is flexible as to choose a weighting parameter for flow 0 is the Crank-Nicolson scheme of second order in time.More details on the theoretical background of this model can be found in Medina et al. (1996).

Model structure
Following the definition by Carrera et al. (2005), parametrization is one element of model structure.In this work, different parametrizations for hydraulic conductivity are explored to improve model calibration.Five hydraulic conductivity zonations are proposed and evaluated through transient modeling along a 39-yr period based upon an annual time step, the first step being the steady state.

Finite element mesh
The model boundary extends to the entire GAS, as shown in Fig. 1.The two-dimensional domain was discretized into 46 862 triangular elements and 23 890 nodes using the mesh generator 2DUMG (Bugeda, 1990).Mesh refinements were introduced in areas of expected steep hydraulic gradients generally associated with heavy pumping in outcropping areas, and in areas of steep topographic gradients in mountain regions located in the central/southeastern portions of the aquifer boundary.Elements area averaged 25 km 2 ; the largest elements were located in the central region of the domain.The element size used here is similar to the cell size used in models of aquifers of comparable size, for example the Great Artesian Basin in Australia, where a uniform 5 km × 5 km cell size was used (Welsh, 2006).

Spatial zonation of hydraulic conductivity
At the time of the construction of the conceptual and numerical models, there were not enough point hydraulic conductivity (K) data as to construct a regionalized distribution map for K using a tool such as kriging.An alternative approach to produce K maps was used instead.
Zonation is one of the methods to parametrize hydraulic properties needed to solve mathematical equations set forth in inverse modeling (Carrera et al., 2005), each producing an alternative model.Figure 3 shows five hydraulic conductiv-ity zonings defined upon different criteria: (Z1) one zone, i.e. uniform K; (Z2) nine zones, their geometry closely replicating the zoning previously defined by Vives et al. (2001) who defined hydraulic conductivity zones based on the location of main geologic structures; (Z3) seventeen zones, their geometry and boundaries delineated following changes on aquifer thickness; (Z4) nineteen zones, their geometry combining the patterns of the piezometric map, namely transitions of hydraulic gradients (see Fig. 2) and the zoning defined by Vives et al. (2001); and (Z5) thirty one zones, their geometry accompanying hydraulic gradients of the piezometric map more closely than Z4.
Limited and scattered information on hydraulic conductivity from aquifer tests and anecdotal values served as previous estimates.Freeze and Cherry (1979) indicated a maximum threshold value for sandstones of 1-2 m d −1 .Araújo et al. (1999) and Sracek and Hirata (2002) reported values of 8.7 m d −1 and 13 m d −1 , respectively, for mean K in the Brazilian states of Sao Paulo, Paraná and Rio Grande do Sul (see Fig. 1).Within Paraguayan territory, values of 1.6-3.8m d −1 were reported by Vassolo (2007), whilst hydraulic conductivities between 0.12 and 5.76 m d −1 , with an average of 1.5 m d −1 , were published for northern Uruguay (Gómez et al., 2010).
The effect of temperature over the flow field was taken into account by correcting the hydraulic conductivity by means of the Schneebeli formulae (Custodio and Llamas, 1976) with θ the temperature in • C. Each finite element within a given zone was assigned a K-value corrected for temperature.
The temperature was interpolated from the isotherms map provided by Gastmans et al. (2012).Hydraulic conductivity was then converted to transmissivity in TRANSIN by multiplying the corrected initial guess for K for a given element by the mean thickness of the aquifer for that element.Initial K-values were recalculated through automatic calibration.

Boundary conditions, recharge, streams and pumping
According to Gastmans et al. (2012), the GAS limit presents a combination of no flow, outward and inward flow conditions.Figure 4 shows the boundary conditions implemented in the model, resulting from the proposed conceptual model, the piezometric map and the calibration process.Some portions of the boundary were simulated with a mixed or streamlike boundary condition, such as the southern border within Brazilian territory.The magnitude of eastern fluxes was previously estimated multiplying the area of a narrow strip of outcropping GAS along a stretch of the boundary by a recharge rate of 3 % of the mean annual precipitation in the area.In the rest of the outcropping areas (see Fig.  annual precipitation which varied from north to south according to the precipitation gradient previously mentioned; 3 % was imposed along the western area and 1.5 % along eastern and southern areas.Table 1 shows recharge rates and recharge volumes for each simulated recharge area shown in Fig. 4. Rates were not automatically calibrated; however, they were modified annually, multiplying the steady state rate by the value of the temporal function shown in Fig. 5.The value of that function represents anomalies with respect to the mean annual precipitation, i.e. a value of the time function equal to 1.2 for a particular year means that the annual precipitation for that year is 20 % higher than the mean annual precipitation.The precipitation series corresponding to the Rivera-Santana station located on the border between Uruguay and Brazil was used to construct Fig. 5, which was deemed indicative of the precipitation temporal variability. The value of the temporal function corresponding to the first time step is equal to one.The GAS underlies the Río de la Plata Basin, the second largest in South America, characterized by a highly dense drainage network that discharges into the main waterways of the region: Paraná, Paraguay and Uruguay rivers.In the occidental border, numerous streams drain toward the Paraguay River, located outside the model area, while the rest of the streams drain toward the Uruguay and Paraná rivers.Stream-aquifer interactions were simulated along 29 streams reaches (see Fig. 4).For their quick location and analysis, all simulated rivers were identified by a number.As previously described, there are still great uncertainties regarding discharge pathways, without which a sound water balance would be difficult to close.In this work, it was assumed that reaches along the Paraná and Uruguay rivers, located in confined areas but nearby the region where the aquifer is closest to the surface, could interact with GAS waters.
Within TRANSIN, leakage between surface water bodies and the adjacent aquifer is computed as where Q is the stream leakage (L 3 T −1 ); α is the leakance coefficient (L 2 T −1 ); h is the piezometric head (L), and H ext (L) is a reference, external water level.If field data is available H ext is usually the stream stage.In absence of such information, in this work, H ext was interpolated from ground elevations and piezometric levels, in the latter case only for the Paraná and Uruguay river reaches.The leakance coefficient is the ratio between the stream cross sectional area times the hydraulic conductivity of streambed sediments, Kb, and the thickness of these sediments.The coefficient was estimated assuming Kb-values two orders of magnitude smaller than Kvalues for sandstones (around 1 m d −1 according to Freeze and Cherry (1979)), a thickness of streambed sediments equal to 1 m, and a cross-sectional area ranging from 100 to 1000 m 2 for all streams, except for the Uruguay and Paraná rivers, for which higher values were adopted.Leakance coefficients for the Paraná and Uruguay rivers were 50 and 20 m 2 d −1 , respectively, whilst for the rest of the streams were between 1 and 10 m 2 d −1 .Sensitivity runs performed during the early stages of the model development showed no significant changes in model results for the range of α tested.

Pumping
Total pumping was estimated based on pumping rates reported at the time of wells construction; therefore, it could be either overestimated as some wells may not be currently operational, or underestimated as others may not be accounted for.The current groundwater exploitation totals 2 847 013 m 3 day −1 , i.e. 1040 hm 3 yr −1 , distributed as follows: 1.3 % in Argentina, 93.6 % in Brazil, 2.2 % in Paraguay, and 2.9 % in Uruguay.Twenty pumping zones were defined (Fig. 6), and each zone was assigned a different rate based on the geographical distribution of wells.The Sao Pablo State in the northeastern region of the aquifer concentrates the highest amount of wells and water extraction, reaching 63 % of total pumping.

Calibration
The model calibration approach consisted of using a combination of automated and manual methods.The primary objective of the calibration was to minimize the difference between simulated and observed hydraulic heads, while seeking hydrogeologic parameters values consistent with the current knowledge of the aquifer.Five different hydraulic conductivity zonations, i.e. models, were set forth using the criteria explained in Sect.5.2, increasing the number of zones on each zonation pursuing the principle of parsimony whilst relating the zones geometry to different aquifer characteristics.Previous estimates of hydraulic conductivity for each zone were modified through automatic calibration.
Insufficient head data precluded performing a conventional modeling approach, i.e. model calibration for steady state with field data representing pre-development conditions as close as possible, and then use of the steady sate simulated aquifer head as the initial condition for the transient simulation.A different calibration strategy was used instead.
A 39-yr transient simulation with yearly periods, the steady state being the first time year of the time series, was performed.The model was calibrated against 317 observed piezometric levels.Those levels span a 30-40 yr time window starting in the 1970s, therefore some may not be representative of current conditions in areas of intense pumping.Moreover, the reliability of some field data is questionable as wells may not be cased along confining units, resulting in an integrated reading that may be interfered by local-type flows.An additional source of error is the value of well elevation, which is needed for estimating piezometric level, especially in areas of steep slopes.Observations are not evenly distributed across the study area; the majority of measurements are located near or at outcropping areas of the aquifer.Along the central region, where the aquifer reaches its maximum depth, data points are sparse.The absence of data is notorious within Argentinean territory.Observed levels range between a maximum of 1202 m a.s.l. and a minimum of 10 m a.s.l.Head measurements were taken at the time of well drilling; about 70-80 % of the wells were constructed during the last decade when there was a significant increment in groundwater exploitation not evenly distributed across the entire aquifer area.The rest of the measurements are spread in time since the seventies.Therefore, in the absence of transient head data for calibration, all available observations were assigned to the last period.Given the quality of observations and the piezometric levels maximumminimum range, the calibration target was set at ±40 m for mean error statistics.

Results and discussion
Model performance and results were analyzed following various criteria applied to piezometric levels, mass balance com- ponents, hydraulic conductivity zonations and ranges of hydraulic conductivity.

Piezometric levels
The optimized transient models were evaluated with respect to the match between observed and simulated piezometric levels.The root mean square residuals (RMS) were calculated as follows (Zheng and Bennett, 1995): where h sim and h obs are simulated and observed hydraulic head, respectively, and N is the number of observations.The subscript i indicates observation number, while the term in parentheses is called model error.RMS values show an improvement of model calibration as the number of K zones increases (Table 2).It remained outside the calibration target for all but the last zonation reaching a minimum of 36.99 m.The linear correlation coefficient R 2 varied from 0.922 for Z1 zonation to 0.968 for Z5 zonation, indicating a significant linear correlation between calculated and observed levels for all cases (Table 2).Figure 7 shows calibration results for scenario Z5.Model results closely replicate groundwater flow patterns concerning both flow directions and hydraulic gradients.Regional groundwater flow is from northeast and northwest toward the center of the sedimentary basin and then south.Modeling results along the western boundary show that regional flow is disrupted by local recharge/discharge systems, a pattern present in the observed piezometric map of Fig. 2. The eastern boundary and adjacent areas are characterized by steeper simulated hydraulic gradients, in coincidence with observed gradients within highly exploited areas in the northeast and steep terrain in the central-east, while smoother gradients are encountered at the center of the simulated area.
Nonetheless, R 2 is not a good indicator to detect overestimated/underestimated areas.The geographic distribution of errors, with their corresponding sign and magnitude, not only highlights the location and density of calibration data but also helps identifying underestimated/overestimated model areas that would require further modeling efforts, either on the conceptual model or on the calibration process (Fig. 8).The model performs evenly across the modeling domain, with no  identifiable overestimated/underestimated regions.As shown in the histogram at the bottom of Fig. 8, the number of outliers, i.e. points for which the absolute value of model error is greater than 80 m, reduces progressively as the number of K zones increases.By the same token, the number of data points within the calibration target also increases reflecting the improvement of model performance, though some extreme errors persist no matter what K zoning is used.This can be explained by several factors.Firstly, alternative conceptual models based solely in the K parameter may be an acceptable approach at regional scale but it may be questionable at different spatial scales as the presence of geologic structures not explicitly included in the regional model may influence groundwater circulation.Secondly, calibration data are far from ideal, affecting in turn, model fit.Ideally, availability of transient piezometric levels would be desirable.Only a single set of calibration data points was available with no identification of the time of measurement, limiting transient calibration strategies.This situation is particularly critical in areas of intense pumping.Then, observations were assigned to the last simulated year.Since pumping is very low in regional terms, this approach was considered reasonable, although it is recognized that it introduces calibration errors at local-meso spatial scales.

Water budget
An independent water balance for the entire aquifer is not available; nonetheless, qualitative and quantitative analyses are carried out wherever possible to verify the model water budget.Even though only results for scenario Z5 are shown in Fig. 9, water budget components were very similar for all zonations.
Recharge and stream/aquifer fluxes are the dominant input and output flow components, Fig. 9a and d Simulated stream leakage showed sensitivity to the forcing terms set forth by recharge, prescribed eastern flow (Fig. 9b) and pumping (Fig. 9c).On the contrary, flux through constant head boundaries located to the west (Fig. 9e) and western and southern boundary flows (Fig. 9f and g, respectively) were almost invariant over time and of comparable magnitude, ranging between 155.2 hm 3 yr −1 and 622 hm 3 yr −1 (5 to 20 m 3 s −1 , respectively).Storage augments and decreases in response to sink/sources (Fig. 9h), but no clear trend can be identified.However, during the first years of exploitation, the combination of abundant recharge and low pumping produces a rapid increase in cumulative storage (Fig. 9i).As pumping increases, cumulative storage stabilizes to decline rapidly towards the end of the simulation period, in tune with increasing pumping rates and low-recharge years.Declining recharge rates over time are illustrated by the straight line in Fig. 9i, representing the ratio between the recharge rate for the i-th year, REC(i) and the steady state recharge rate REC(SS).
The recent inventory of production boreholes in the GAS (Vives et al., 2008) resulted in a current groundwater exploitation of about 1040 hm 3 yr −1 .In global terms, the model water balance indicates that the current rate of groundwater withdrawals does not exceed the rate of recharge.Notwithstanding, pumping is concentrated in heavily populated and industrialized areas where groundwater withdrawals are expected to continue rising in coming years; consequently, at local scale the situation may be reversed, even at present.The regional model presented in this work did not intent to quantify local-scale issues.Local models already developed in critical areas would serve that purpose.
Table 3 shows budget terms for the years of maximum and minimum recharge as well as averages for the whole transient period.Recharge ranged from 2014 to 6470 hm 3 yr −1 , averaging 3156 hm 3 yr −1 , equivalent to 84.2 % of inflows.On average, pumping totaled 665 hm 3 yr −1 , representing 16.2 % of outflows.Part of the recharge is converted to leakage along streams.For the minimum recharge year, leakage constituted 53 % of outflows, reaching 70 % for the maximum recharge year, with an average of 61.4 %.Water is discharged through the aquifer boundaries, except at the eastern boundary.This result is consistent with the conceptual model.The magnitude of the western flow plus the outward flow through constant heads to the west is comparable to pumping, while southern flows represent less than 10 % of outflows.Hirata et al. (2008) made an attempt to independently quantify boundary fluxes using Darcy's law, assuming a hydraulic conductivity range between 1 and 3 m d −1 and an aquifer thickness between 50 and 300 m, depending on location.They estimated that the southern outward flow would be between 36 and 216 hm 3 yr −1 ; the simulated value for zonation Z5 was 230 hm 3 yr −1 for steady state conditions, remaining almost invariant throughout transient years.Western boundary flow estimated by Hirata et al. (2008)   between 137 and 353 hm 3 yr −1 ; the simulated value for steady state was 335 hm 3 yr −1 , augmenting during the transient period.Notwithstanding the uncertainties and model limitations, this comparison contributes to building confidence in modeling results, helping to progressively close a water balance for the aquifer.
Recent calculations limit recharge to less than 10 % of mean annual precipitation, with values closer to 3-4 %.In a recent study, Rabelo and Wendland (2009) reported 3.5 % of mean annual precipitation of net recharge obtained through a groundwater model calibration in the northeastern region of the aquifer.In this work, the steady state recharge was 3516 hm 3 yr −1 , equivalent to 35.2 mm yr −1 .Considering a mean annual regional precipitation of 1400 mm, modeled recharge amounts to 2.5 % of that value.Average stream leakage resulted 2512 hm 3 yr −1 (81 m 3 s −1 ) (in actuality, it ranged between 81 m 3 s −1 for Z5 zonation to 93.8 m 3 s −1 for Z2 zonation).This result leads to two conclusions: firstly, stream/aquifer fluxes are not very sensitive to the number of K zones; and secondly, total leakage is small in comparison to minimum flow discharges of the main rivers in the region.For example, the simulated average leakage for the Uruguay River was 8 m 3 s −1 while the observed absolute minimum discharge for the period 1931/2001 was 382 m 3 s −1 (its average discharge is 2300 m 3 s −1 ).This meager flow renders the verification of some modeling results very challenging.Araújo et al. (1999) postulated that the principal discharge area of the GAS was probably located between the Paraná and Uruguay rivers, although Campos (2000) raised doubts about this hypothesis.Nonetheless, this hypothesis was tested with the model and proved to be consistent with the other water budget terms.
In these four cases, stream/aquifer fluxes shown on the left of Fig. 10 were sensitive to changes in hydraulic conductivity.Streams located in heavily pumped regions, i.e.Tacuarembó River and Jacaré Papira River, switched from a gaining condition in early years to a losing condition over time for Z5 zonation.An auxiliary variable was defined for the analysis.Figure 10 shows the difference between leakage for a particular year-Flux(i)-and leakage for steady state-Flux(SS)-.Transient recharge relative to steady state recharge, represented as a solid line, is also shown on the same figure.The Tacuarembó and Jacaré Papira rivers show a similar behavior: flow from the river to the aquifer increases over time in response to increasing pumping.This situation has relevant connotations for conjunctive water resources management in localized areas of the aquifer and should be studied in more detail by combining field work to verify flow magnitudes and numerical simulations to predict system response under various scenarios.Recharge impacts leakage on the Tacuarembó River; however, it has little influence on the Jacaré Papira River.Recharge and pumping do not affect the stream/aquifer relationship in confined areas of the aquifer, represented in this case by the Paraná River.
For those rivers, stream/aquifer fluxes for steady state, year 4 (maximum recharge year, minimum pumping), year 30 (minimum recharge year, average pumping), and last simulated year (close to average recharge, maximum pumping) were evaluated, comparing exchange fluxes with the corresponding mean streamflow (Table 4).Leakage from the aquifer to the river decreases in the Jacaré Pepira Stream due to increasing pumping.By the same token, leakage from the stream to the aquifer in the Ypané River increases slightly with time; however, that change is no so drastic due to relatively low pumping rates in the area.The simulated condition for the Tacuarembó River changes from effluent to influent though the leakage magnitude is small.There is no pumping in the area to justify this behavior.Finally, the interaction between the aquifer and the simulated reach of the Paraná is negligible compared to the river mean discharge.

Hydraulic conductivity
For all zonations, K-values resulted from the automatic calibration algorithm available in TRANSIN, minimizing an objective function written in terms of heads, parameters and concentrations (if a transport problem is solved).In TRANSIN, the minimizing algorithm uses the Marquardt Method, an iterative algorithm to solve non-linear problems for parameter estimation by the least square method.
Partitioning the model domain into an increasing number of K zones was effective in improving the model fit, reducing calibration errors.However, calibrated K for all five scenarios were higher in the central region of the modeled area, with values above the conductivity range typically expected for sandstones, even considering scale effects (Fig. 11).Along a central corridor and some adjacent areas, K was always greater than 15 m d −1 , even higher than 30 m d −1 for scenarios Z2, Z4 and Z5, reaching a maximum of 114.8 m d −1 for scenario Z2.The maximum calibrated K reduced considerably for scenario Z5, reaching 72.6 m d −1 .Those values are consistent with lower hydraulic gradients in  the area, but they also would indicate the need for the model to conduct flows in that area through, for instance, a preferential flow zone or a connection with overlying/underlying geologic units.This option should be explored in the future, incorporating a more detailed geologic layering and their hydraulic interconnection to the conceptual model.Zones of low calibrated hydraulic conductivity along parts of the eastern boundary and the northwest coincide with high hydraulic gradients.
Except for the homogeneous case, the lower end of the K range was between 0.1 and 2.3 m d −1 , a value coherent with sandstone K-values reported in the literature (Freeze and Cherry, 1979) and for the GAS.
Whenever possible, the zones calibrated hydraulic conductivity was compared with both point data and K ranges reported by previous authors, seeking calibration values coherent with available data.Sixteen pumping test data (all concentrated in 150 km 2 , the average finite element size of the model mesh is 25 km 2 ), with K ranging from 0.17 to 19.92 m d −1 , were available at the Uruguayan-Brazilian border around the cities of Rivera/Santana (Gómez et al., 2010).
Seven pumping tests, with K ranging from 1.6 to 3.8 m d −1 , were available in southern Paraguay (BGR, 2008).For the best model fit scenario, K for the zones overlaying those available point data were 5-10 m d −1 (K22), and 2.5-5 m d −1 (K31), respectively, showing a good match between calibrated and observed values.
Based upon eleven pumping tests, highly concentrated in space, hydraulic conductivity values of around 3 m d −1 were reported for the Riberao Preto area in NE Brazil (Sracek and Hirata, 2002).Calibrated values for that area (K11) were in the range between 2.5 and 5 m d −1 .
Besides the point-to-zone comparison, an attempt was made to compile ranges of reported K-values attributed to various authors that could be used to further assess the model calibration.Transmissivity (and for that matter hydraulic conductivity) is one of many aquifer parameters that vary with the scale of measurement.This issue is well documented in the literature (Sánchez-Vila et al., 1996;Nilsson et al., 2001).An in-depth analysis of this topic is well beyond the scope of this work.Nevertheless, given the spatial extent of the constructed model, a general comment is merited.This is the first, though not exhaustive, attempt to compile published values of hydraulic conductivity for GAS sandstones in the four countries in order to assess the consistency of the model automatically calibrated K and identify spatial scale effects.Figure 12 shows the range of reported K-values compiled so far, attributed to various authors.Instead of the classical representation of log K vs. the scale of observation, the x-axis simply corresponds to a bibliographic reference number.
Depending on the use of packer tests, slug tests or pump tests, a clear scale effect in crystalline rocks, porous carbonate rocks and carbonate aquifers was identified by Sanchez-Vila et al. (1996), Schulze-Makuch and Cherkauer (1998) and Whitaker and Smart (2000), respectively.For GAS sandstones, there seems to be a scale dependence in hydraulic conductivity due to progressive incorporation of larger and better connected transmissive zones as the support volume increases.Reported K-values from aquifer test data show little variability.Available data sources do not provide enough information as to identify the test-type influence demonstrated by the aforementioned authors.Mean hydraulic conductivity for this data set is 2.15 m d −1 , with reported minimum and maximum values of 0.1 m d −1 and 4.56 m d −1 , respectively.
The mean K-value obtained from calibration of local scale groundwater flow models is 2.38 m d −1 , comparable to pump test results.This is presumably due to the classical approach of groundwater model calibration in which the model is used to predict recharge rates from information on water levels,   The group "Local scale aquifer test" includes aquifer tests data; the group "Local scale modeling" includes calibrated K from models ranging from 100 km 2 to 5000 km 2 in extent; the group "Regional scale modeling" includes calibrated K from models whose extent is larger than 400 000 km 2 .Reference number 16 corresponds to the model of the Paraguayan sector of the GAS; references 17-21 correspond to the five zonations discussed in this work.K is reported either as a single value or a K range.
hydraulic conductivity and other parameters that may lead to non-unique modeling results (Scanlon et al., 2002).Hydraulic conductivity and recharge rates are often highly correlated; consequently, calibration based only on water level data is limited to estimating the ratio of recharge to hydraulic conductivity.Hence, as Scanlon et al. (2002) stated, the reliability of recharge estimates depends on the accuracy of the hydraulic conductivity data.Then, in many local scale models aimed at estimating recharge rates, hydraulic conductivity values estimated from field data are considered rather representative while recharge rates are the main calibration parameter (Rabelo and Wendland, 2009;Gomez et al., 2010), leading to "calibrated" K-values very similar to pump test values.
Table 5. Model structure identification based on the following criteria implemented in TRANSIN: AIC (Akaike, 1974), BIC (Akaike, 1977), ø (Hannan, 1980) and d k (Kashyap, 1982).As the size of the region of interest increases, calibrated K also increases in both its magnitude and calibrated range.It should be pointed out that the Paraguayan groundwater flow model was based on a more detailed representation of geologic formations.GAS sandstones were simulated with two homogeneous layers, for which calibrated K was 0.05 m d −1 and 3 m d −1 , respectively; therefore, they should be considered single values and not a hydraulic conductivity range (see dotted line in Fig. 12).

Model structure
In this section, different criteria implemented in TRANSIN are used to evaluate the model structure.All of them are based on judging the models, i.e. different zonations, according to the maximum likelihood goodness of fit.Besides the model structure, the Akaike information criterion (AIC) by Akaike (1974) also includes the number of parameters.Akaike (1977) extended the previous approach adding the natural logarithm of the number of data defining BIC, a Bayasian information criterion.In 1980, Hannan (1980) introduced the ln (natural logarithm) of the ln of the amount of data while Kashyap (1982) judged the models goodness of fit based on the number of parameters, the ln of the number of data and the ln of the determinant of the Fisher information matrix, which is the expectation of the Hessian as well as the lower limit of the covariance matrix.The interested reader may resort to the original works of Akaike (1974Akaike ( , 1977)), Hannan (1980) and Kashyap (1982) or to the discussion by Heredia (1995) for a more in-depth explanation Table 5 shows the values of each criterion for the different models.The Z5 zonation was the best model, mainly supported by the large number of K zones.Therefore, it is expected to yield the best fit, even at the risk of model overparametrization.
Increasing the number of zones decreases the values of the different criteria, suggesting that the definition of new K zones is a consistent approach.Even though in the case of Z3 a better fit is achieved, the model is overparameterized.This is concluded from the values of ø and BIC, which allow seeing the relationship between the number of parameters and the number of data.In contrast, with only two more zones, the Z4 model improves its predecessor, suggesting that zoning is defined consistently.This is reasonable because Z3 and Z4 zonings were defined using different sources of information and criteria (see Sect. 5.2).

Conclusions
Fresh water has already become a limiting resource in many parts of the world.In the future, it will become even more limiting due to increased population, urbanization, and climate change.This limitation will be caused not just by increased demand for water, but also by pollution.Immersed within this context, the Guarani Aquifer System (GAS) is being increasingly exploited for freshwater supply, and for industrial and agricultural uses.In large aquifers, relevant for their considerable size, regional groundwater modeling remains challenging given geologic complexity and data scarcity in space and time.Yet, it may me conjectured that regional scale groundwater flow models can help in understanding the flow system functioning and the relative magnitude of water budget components, which are important for aquifer management.The present study on the GAS has shown that a transient-regional scale groundwater flow model can provide valuable insights regarding those two issues, given an extraction volume still very small in comparison with the aquifer volume.Even though the model was constructed on a simplified conceptual model, it constitutes the first attempt to simulate the entire aquifer including budget terms previously overlooked.
The hypothesis of a continuous sedimentary unit may be sustainable as a first approximation to construct a numerical model covering the full extent of the aquifer.Combined with an increasing number of K zones and an appropriate set of boundary conditions, that hypothesis yielded errors within the calibration target in a regional sense.Nonetheless, this approach was insufficient to improve calibration in areas known for the presence of structural controls that may influence groundwater flow patterns.
Given the amount and quality of data for calibration, model results were acceptable as measured by standard statistics.Surely, the availability of current piezometric levels, extraction volumes and stream discharges would help to produce a better model.
Calibrated K-values were coherent, with low/high hydraulic gradients.However, calibrated K in the central region of the modeled area was above the conductivity range typically expected for sandstones.This result highlights the need to analyze the possible hydraulic interconnections between the GAS and pre-GAS/post-GAS sediments.
The location and character of the southwestern boundary of the aquifer remains an open issue, though it was demonstrated that reaches along the Paraná and Uruguay rivers could be potential discharge zones, as postulated by previous authors.Simulated stream leakage along those reaches was very small, but so were other water budget components when analyzed in perspective considering the aquifer extent and its storage volume.The model was also instrumental for the understanding of the dynamics of the system along the western boundary, which resulted in an outflow condition, compatible with observed water level data and hydraulic gradients.The magnitude of the simulated discharge through that boundary matched independent estimates by other researchers.At a regional scale, the importance of the stream/aquifer interaction process was manifested by the 61.4 % contribution of this term to total outflow, its magnitude being second to recharge.In addition, model-calculated recharge was coherent with recent estimates from other studies.Even though an independent water budget for the entire aquifer is not yet available, these qualitative and quantitative analyses contribute to building confidence in model results.
On average, pumping represented 16.1 % of inflows while aquifer storage experienced a small overall increment.The model water balance indicates that the current rate of groundwater withdrawals does not exceed the rate of recharge in a regional sense.Notwithstanding, pumping is concentrated in heavily populated and industrialized areas where groundwater withdrawals are expected to continue rising in coming years; consequently, at local scale the situation may be reversed, even at present.It is worth noting that this regional model did not intend to quantify local-scale issues.Local models already developed in critical areas would serve that purpose.
The model presented in this work greatly improved its predecessors, integrating information recently generated and extending the model area.In all, the parameter sets and the water balance from the calibrated model add to the current understanding of the hydrodynamics of the GAS, highlighting the importance of contributing water balance terms.Moreover, the model was instrumental at identifying data and conceptual model weaknesses and uncertainties that can be grouped into three major themes: geology, role of structures on the flow system, and definition of discharge/recharge zones.Future works should be mainly directed to the following: -Analyze the hypothesis of a compartmentalization of the aquifer and its influence on the regional flow system, as suggested by recent studies.
-Analyze the role of pre-GAS and post-GAS formations on piezometric levels and hence, on groundwater flows.
-Evaluate the role of local geologic structures on the flow system.This would help to reproduce some piezometric (and hydrochemical) anomalies that could not be represented by the current conceptual/numerical models.
-Conduct more in-depth and model-independent studies of flows distribution especially related to recharge, pumping and river/aquifer interactions in outcropping areas, performing water balances at representative areas with sufficient field data to support hypotheses and conclusions.
-Simulate groundwater age in order to validate alternative hypotheses of the flow system functioning, supported by isotopic sampling and analysis.

Fig. 4 .
Fig. 4. Simulated boundary conditions, recharge areas and stream reaches (note: only the simulated reach of minor streams and major rivers like the Paraná and Uruguay rivers, are indicated; the full drainage network is omitted for the sake of clarity).

Fig. 5 .
Fig. 5. Time function for recharge and pumping applied to the transient simulation.

Fig. 8 .
Fig. 8. Residual error at calibration points for all conductivity zonations.Circles indicate model underestimation; triangles indicate model overestimation.Residual errors have been grouped into interval classes.

Fig. 10 .
Fig.10.Left: stream/aquifer flux (m 3 s −1 ) versus time for selected streams for all conductivity zonations.Right: ratio of transient recharge to steady state recharge (solid line -dimensionless) and change of stream/aquifer flux with respect to steady state condition (m 3 s −1 ) for all conductivity zonations (the latter is not included for the Paraná River because the magnitude of the change was negligible, i.e. less that 0,1 m 3 s −1 , for all zonations).

Fig. 12 .
Fig. 12.Comparison between calibrated hydraulic conductivity and hydraulic conductivity values compiled from bibliographic sources.The group "Local scale aquifer test" includes aquifer tests data; the group "Local scale modeling" includes calibrated K from models ranging from 100 km 2 to 5000 km 2 in extent; the group "Regional scale modeling" includes calibrated K from models whose extent is larger than 400 000 km 2 .Reference number 16 corresponds to the model of the Paraguayan sector of the GAS; references 17-21 correspond to the five zonations discussed in this work.K is reported either as a single value or a K range.
themselves and their use in the context of TRANSIN, respectively.
et al.: Conceptual and numerical modeling approach of the Guarani Aquifer System

Table 1 .
Recharge rates and volumes over each simulated recharge area.

Table 2 .
Goodness of fit estimators.

Table 3 .
Model mass balance (volumetric rates expressed in hm 3 yr −1 , percentages referred to average rates) -PEF: prescribed eastern flow; PH: prescribed head; WF: western flow; SF: southern flow.Maximum values correspond to the forth simulated year; minimum values correspond to the thirtieth simulated year.

Table 4 .
Comparison between stream-aquifer fluxes with streamflow for selected streams and simulation years (percentages express stream leakage relative to mean streamflow).