Integrated surface-groundwater modelling for assessment of aquifer recharge and groundwater availability in a semi-arid region in the context of inter-basin water transfer

The São Francisco River (SFR) ow has been partially transferred to dryland catchments in Northeastern Brazil (NEB), to help deal with recurrent regional water shortages. However, the inuence of this water transfer on overexploited aquifer systems had not been investigated. Our goal was to assess the groundwater recharge and the potential of the SFR water transfer to increase groundwater availability by channel transmission losses in the Medium aquifer (2,166 km²) in the state of Ceará (NEB). We developed integrated surface-groundwater modelling to improve estimations of the groundwater resources (GR), considering the spatio-temporal variability of inltrated rainfall for aquifer recharge. Surface hydrology and aquifer recharge were simulated by a process-oriented hydrological model (WASA), and groundwater ow by MODFLOW. The annual recharge (66 mm/y) was approximately twice the amount of water abstracted (1990-2016). This abstraction rate indicates a sustainable use of groundwater. However, the annual recharge dropped to 13.9 mm/y from 2012 to 2016, a drought period. MODFLOW accurately reproduced the decrease of groundwater level, driven by both low recharge and high pumping rates during the drought. We found that SFR transfer could decrease the stress on GR during a prolonged drought. However, this additional recharge (6.89x10 6 m³/y) did not compensate for the relevant decrease of the hydraulic heads in areas that do not surround the receiving streams. Therefore, the water transfer would represent an additional recharge source, but it would not alone solve the issue of the unsustainable management of GR.


Introduction
Water management is a challenge as it has to meet the growing needs of the current generations for social and economic growth, without compromising the ability for future generations to meet their needs. The challenge is even greater in water-limited environments such as semiarid northeastern Brazil (NEB), where annual precipitation levels range from 400 to 800 mm, and evaporation levels reach over 2000 or 3000 mm in some areas (Campos 2015). This large dryland (982,566 km 2 ) is characterized by its dry and rainy seasons, and recurrent droughts, i.e. consecutive years of low precipitation. Because of these natural dry conditions and increasing water demand, NEB is regularly confronted with water insecurity issues (Gutiérrez et al. 2014).
To help distribute water resources more equally within NEB, an inter-basin water transfer scheme has been designed by the Brazilian Federal Government. This project aims to transport water from the São Francisco River (SFR) to several intermittent rivers and reservoirs in the states of Ceará, Rio Grande do Norte, Paraíba and Pernambuco (Ministério da Integração Nacional 2004; Lima 2005; Tortjada 2006). SFR is a perennial allogenic river that begins its course in a water-rich region of Southeast Brazil before owing through the southern part of NEB. SFR discharge represents 3.7 times the surface runoff of the receiving region.
In the state of Ceará, the SFR's water rst reaches the rivers in the Salgado watershed, which hosts the major groundwater resources of the state within the Araripe sedimentary basin, where groundwater is abstracted by pumping wells for drinking water, irrigation, industry and leisure. However, environmental changes, such as the drying up of creeks and permanent vegetation changes, have been reported due to unsustainable groundwater abstraction rates (Mendonça et al. 2005). The Water Agency of the state of Ceará (COGERH) estimated that the amount of groundwater extracted is much larger than the mean annual recharge (Governo do Estado do Ceará 2011). Therefore, the receiving rivers from the SFR water transfer ow through an overexploited groundwater region in Ceará. While studies have been focused on the potential environmental and social impacts of such water transfer projects in Brazil (MI 2004;Caúla and Moura 2006;Henkes 2014), none of them has yet studied their potential in uence on groundwater availability due to aquifer recharge by channel transmission losses (Costa et al. 2012;2013a). These losses can play an important role in the recharge of the groundwater resources in the Salgado watershed, as losing streams often occur in such areas in northeast Brazil.
The main objective of this study was to assess the groundwater recharge and the potential of the SFR water transfer to increase groundwater availability in the Salgado watershed by channel transmission losses. In order to accomplish this goal, we developed an integrated modelling approach (WASA-MODFLOW) to improve estimations of the groundwater resources in this watershed, taking into account the spatio-temporal variability of in ltrated rainfall for aquifer recharge.
Groundwater models, such as MODFLOW (McDonald and Harbaugh 1984), can be useful to assess scenarios of surface water and groundwater management. For instance, aquifer recharge by channel transmission losses occurring at streambeds and banks can be represented in MODFLOW using a built-in model package that calculates the exchange ux between streams and groundwater. Three studies have already applied MODFLOW to quantify the groundwater resources in the Araripe sedimentary basin. They showed, however, poor performance in reproducing the observed groundwater levels, with up to 40 m difference between simulated and observed groundwater head (Golder-Pivot 2005; Rede Cooperativa de Pesquisa 2007; the World Bank 2011).
One reason for that poor model performance may be the oversimpli cation of the groundwater recharge process.
Those studies assumed a constant groundwater recharge rate from the in ltrated rainfall, which neglects both the rainfall temporal variability and the heterogeneity of the soil and land cover. Thus, it could be argued that the choice of using constant recharge rates results in an incomplete representation of groundwater recharge. Therefore, to model the regional groundwater ow system in the Salgado watershed, we included the temporal and spatial variation of the in ltrated rainfall within the watershed landscape. This process is modelled by WASA, which is a semi-distributed rainfall-runoff model that has previously been successfully applied in Ceará (e.g. Güntner  The Salgado River watershed ( Figure 1) covers an area of 12,900 km 2 , with a population of about 870,000 inhabitants distributed over 23 municipalities (IBGE 2012). It contains Ceará's major groundwater reserves within the Araripe sedimentary basin, which is formed of aquicludes and aquifer systems ( Figure 1). The focus of this study is the Medium aquifer system, because it is the main groundwater resource of the Araripe basin with estimated permanent reserves of 21.024x10 9 m 3 and annual average recharge of 27.7x10 6 m 3  The study area is characterized by high spatial and temporal variability of rainfall. Intra-annual rainfall variability is distributed in two distinct periods: a wet season from January to May with over 70% of mean annual rainfall, and a dry season from June to December. Mean annual precipitation is 920 mm, with higher values of around 1,030 mm in the western part due to the orographic effect of the Araripe plateau (Golder-Pivot 2005). The area is characterized by high potential evapotranspiration up to 1,566 mm/y as measured in the municipality of Barbalha (Costa et al. 1998).

Geology
The Salgado watershed comprises the Araripe sedimentary basin in the south, which covers 37% of the watershed area (see the aquicludes and aquifer systems in Figure 1); crystalline basement covers the remaining area. The former consists of a superposition of sedimentary formations, enclosed by the crystalline basement. The sedimentary formations result from several geological events that occurred during the rifting of Gondwana, leading to the opening of the southern Atlantic, with ages varying from Silurian to Cretaceous (Morais Neto et al. 2006). According to Ponte and Ponte-Filho (1996), the sedimentary sequence can be described as follows: The Mauriti formation of Siluro-Devonian age forms the base of the sequence in direct contact with the crystalline basement and is composed mainly of coarse quartz uvial sandstone, conglomerate and ne eolian sandstones.
The pre-rift and rift sequence are composed of deposits from upper Jurassic to lower Cretaceous, called the Cariri Valley Group with the formations of Brejo Santo, Abaiara and Missão Velha. The Brejo Santo formation is characterized by ne sandstones, siltstones and claystones. The Missão Velha is mainly composed of coarse sandstones and conglomerates. Abaiara formation consists of ne to medium sandstones and claystones, corresponding to a rifting uvio-lacustrine depositional environment.
The post-rift sequence is composed of sub-horizontal deposits of Medium Cretaceous age, covering the previous formations affected by extensive tectonic events. This sequence, called Araripe Group, comprises the Riacho das Batateiras, Santana, Arajara and Exu formations. The Riacho das Batateiras formation is mainly characterized by medium to coarse uvial sandstones. The Santana formation is mainly composed of ne deposits such as claystones, gypsum and lacustrine clays-silts. The Arajara formation is represented by ne sandstones, siltstones and claystones and composes the cliff of the Araripe plateau along with the Santana formation. The sedimentary sequence is completed by the Exu formation, covering the plateau, composed of red sandstones, claystones and conglomerates.
The sedimentary formations of the Araripe basin have been strongly subjected to the tectonic movements of the pre-rift and rift stages, with a northeast-southwest extension direction. This results in a structural succession of horsts and grabens, which affects the Cariri Valley Group, causing heterogeneities in the thickness of each formation (Golder-Pivot 2005).

Surface Water and Groundwater Resources
In the Araripe sedimentary basin, most springs, which feed the tributaries of the Salgado River, originate from the slope of the Araripe plateau at the contact of the Santana and Arajara formations. Fifteen strategic surface water reservoirs, managed by COGERH, are located in the Salgado watershed. However, none is located directly within the Araripe basin because of the high permeability of the sedimentary formations. Nonetheless, part of the surface water ow into the Araripe basin is regulated by several strategic reservoirs, namely Atalho, Quixabinha, Gomes, Manoel Balbino and Thomas Osterne ( Figure 1).
The Araripe basin contains a system of aquifers created by alternating sandstone and mudstone formations; there are three aquifer systems separated by aquicludes ( Figure 1). Recharge of the Superior aquifer system occurs through direct in ltration of rainfall. The natural discharge of this aquifer occurs at springs emerging at the contact with the underlying Santana aquiclude on the plateau slope. The Santana aquiclude is a very low permeability formation composed of clayey materials with an estimated hydraulic conductivity of 8.6x10 −6 m/d (Mendonça 2001). Machado et al. (2007) estimated that about 20% of recharge percolates through the fractures of the Santana aquiclude to the Medium aquifer system. The Medium aquifer system is composed of the Riacho das Batateiras, Abaiara and Missão Velha formations. Mendonça (2001) estimated that the three formations show a similar hydrogeological behavior with an average hydraulic conductivity of 4.3 m/d and a total thickness of about 500 m. The recharge to this system is the sum of direct in ltration of rainfall, channel transmission losses from rivers, and percolation through the Santana aquiclude. Arti cial discharge occurs through abstraction from numerous boreholes in the area. It is expected that channel transmissions losses are much larger than groundwater discharge into the rivers due to aquifer overexploitation (Governo do Estado do Ceará 2011; Souza and Costa 2014). The Brejo Santo aquiclude is a low permeability formation mostly composed of ne and clayey materials, which separates the Medium aquifer from the largely con ned Inferior aquifer.

The São Francisco River transfer in Ceará
In the state of Ceará, the SFR water transfer was extended with the Ceará's Water Belt project (Cinturão das Águas do Ceará or CAC), consisting of channels to supply water to twelve hydrographic regions. The rst part Jatí Reservoir -Cariús River ( Figure 1) has a total length of 145.3 km in the Araripe sedimentary basin, with a design discharge of 30 m 3 /s. As of early 2022, its construction is close to completion, but its in uence on the groundwater resources in the Salgado watershed has not yet been studied.

Material And Methods
The general approach was rst to calibrate WASA using observations from stream ow gauges to estimate groundwater recharge. The results from WASA were used as input for the MODFLOW simulations, which were then calibrated using head time series from monitoring wells. Finally, the integrated modelling was used to assess the potential of the SFR water transfer project to recharge groundwater by channel transmission losses in the Salgado watershed.

Groundwater model parameterization
There are several open source codes available to model groundwater ow, such as FEFLOW, OpenGeoSys or AquiMod. For this study, the choice of MODFLOW was made because of its capability to simulate all considered subsurface ow processes, and to be consistent with existing studies in the area thus facilitating comparison with these studies. We used the Groundwater Modelling System (GMS) software package to set up MODFLOW.
The model domain was designed to represent the Medium aquifer system. Its spatial discretization, a rectangular structured grid, which was required to apply the MODFLOW's nite volume technique, was selected to simulate the regional groundwater ow without imposing a computational burden. Thus, we set a xed grid size of 500 x 500 m. The resulting 3D grid (XYZ) had 208 cells in the X direction, 159 cells in the Y direction and 3 vertical layers, which was based on the main hydrogeological formations in the Medium aquifer system ( Figure 2).
The area of 2,166 km², in which the Medium aquifer is uncon ned, was used to de ne the ow domain. The cells inside this polygon were active for calculation of groundwater head, while the cells outside this area were inactivated. This resulted in 26,277 active cells in the grid. The aquifer depth was set to a mean value of 500 m, based on Mendonça (2001). Previous studies assumed homogeneous hydrogeological properties throughout the aquifer. However, we divided the total depth in three layers of 100, 100 and 300 m from surface to bottom to account for a possible decrease of hydraulic conductivity caused by sediment compaction (based on Medonça 2001). The digital elevation model (DEM) was obtained from Shuttle Radar Topography Mission (SRTM) data with a resolution of 30 m (http://earthexplorer.usgs.gov/). Its reliability was con rmed by Governo do Estado do Ceará (2009) that performed a validation of the DEM topography using elevation points collected by a Differential Global Positioning System (DGPS), which reported a good correspondence between satellite and eld data.
The lateral boundaries of the Medium aquifer were characterized by the Santana aquiclude in the east, the Brejo Santo aquiclude in the west and the crystalline basement at the northern and southern side. Based on these properties, the lateral boundary conditions were set as no-ow boundaries. However, the Superior aquifer interacts with the Medium aquifer system by percolation through the Santana aquiclude and in ltration of water emerging from springs at the plateau slope at the contact of the Arajara and Santana formations. This interaction was represented in the model by the MODFLOW Recharge (RCH) package.
The contribution of natural springs located at the cliff of the Araripe plateau to the recharge of the Medium aquifer was estimated using a database from COGERH and the National Department of Mineral Production (DNPM). The database lists 320 natural springs with data of annual water production and annual authorized water withdrawal.
Information of pumping wells were retrieved from the database of the Brazilian Geological Survey (CPRM), called SIAGAS (http://siagasweb.cprm.gov.br/layout/pesquisa_complexa.php). Wells with no data concerning depth or pumping rate were excluded. In total, 532 wells were registered within the model domain ( Figure 3). They were implemented in the model by the MODFLOW Well (WEL) package, using elevation, depth and pumping rate (m 3 /d). Because screen depths were not available for these wells, either a default screen of 10 m was chosen or the wells were assigned automatically to the rst layer only ( rst 100 m) without using screens.
The SIAGAS database also provided information about the type of groundwater use (public, domestic, irrigation, leisure or multiple uses), as well as the allowed hourly pumping rate. The daily pumping rate implemented in MODFLOW was based on Golder-Pivot (2005), Rede Cooperativa de Pesquisa (2007) and the World Bank (2011).
For each well with a pumping rate inferior or equal to 50 m 3 /h, a daily pumping regime of 8 hours was assumed, corresponding to private wells for domestic water use. For wells with a pumping rate superior to 50 m 3 /h, a daily pumping regime of 20 hours was assumed, corresponding to public and irrigation wells, likely to be pumping water more intensively. The sum of these 532 pumping rates was the arti cial discharge of the Medium aquifer.
Assuming that the yearly water demand remained constant during the MODFLOW simulation period (2010-2016), the wells represented an abstraction capacity of 7.14 x10 7 m 3 of water per year.
Streams were added to the groundwater model using the MODFLOW Stream ow Routing (SFR2) package. We adopted the simpli ed Manning equation, which required values of hydraulic conductivity of the streambed, width, sinuosity, roughness and value of incoming ow for each stream segment. The package automatically computes ow between the stream and the aquifer, including the channel transmission losses. Values for hydraulic conductivity of the streambed and stream geometry for each segment were retrieved from eld data (FUNCEME 2018). For stream reaches, where no eld data was available, the values were inferred from the closest measured stream. We considered the reservoir water release as upstream input discharge for some stream reaches.
The recharge to the Medium aquifer works as the exchange variable between the surface water model WASA and the groundwater model MODFLOW, i.e. the in ltrated rainfall. Recharge was implemented in the groundwater model as a recharge rate assigned to a polygon, which was provided by the GMS package. Then, the recharge rate was applied to the most upper active cell. The interaction with the Superior aquifer was also represented using the RCH package at the western boundary of the model as aforementioned. The estimated contribution from the Superior aquifer using RCH was added to the recharge calculated by WASA. Machado et al. (2007) estimated that 4% of mean annual rainfall in ltrates into the Superior Aquifer System in forested areas, and only one fth of this value in the case of deforested areas. They also estimated that about 20% of the total recharge percolates through the Santana Aquiclude and effectively contributes to the recharge of the Medium aquifer.

Hydrological model parameterization
The use of WASA allows representation of the temporal variability of in ltrated rainfall by calculating rainfall in ltration and soil moisture redistribution at a daily time scale. This model also permits sub-dividing the study area into spatially-distributed modelling units, including different soil-land cover patterns, to account for the spatial variability of the hydrological behavior within a watershed.
WASA was parameterized by Costa et al. (2013b) for the whole state of Ceará. The model parameterization was organized through ve different scales in cascade, as follows: sub-basin, landscape unit, terrain component, soilland cover component and soil pro le . It also considered the major reservoirs located at sub-basin outlets and the impact of small reservoirs on the simulated runoff within the sub-basins ). Costa et al. (2013b) delineated 215 sub-basins for the state of Ceará based on the location of 155 major reservoirs, 13 stream ow gauges and major river con uences. We adopted their parameterization in this study. The area of the Medium aquifer is covered by three WASA subbasins (sub-basin IDs: 170, 172 and 186, Figure 3). There is a stream ow gauge at the outlet of these WASA subbasins. There are predominantly eight land cover and eight soil types in these sub-basins (Costa et al. 2013b), which produce a complex mosaic of soil-land cover patterns. Note that there are upstream sub-basins that can ow into the sub-basins 186 and 172 ( Figure 3). The inter-basin hydrological connectivity happens if the reservoirs at the outlet of these upstream sub-basins release water downstream. Reservoir releases occur during the dry seasons in order to allocate water downstream or during very moist rainy seasons due to uncontrolled spill over ow. Both processes were simulated by WASA in this study.

Hydrological model calibration
For the three aforementioned sub-basins covering the area of the Medium aquifer, WASA was run from January 1990 to December 2016, using the rst ve years as model warm-up period. As the aim of the study was to model the regional groundwater ow over time periods of several years, the calibration of WASA primarily focused on the monthly time scale. Data available were daily discharge time series, measured at stream ow gauges located at the outlet of the sub-basins. Thus, the average river ow simulated by WASA was compared to the average river discharge at the stream ow gauge on a monthly basis.
The model was then calibrated to t the simulations to the observations, by adjusting the saturated hydraulic conductivity and the porosity of the soil pro les in a manual trial-and-error process. These parameters are the most sensitive for WASA (Güntner 2002). It was expected that WASA would produce more runoff than that observed in the study area, since the state of Ceará is mainly composed of low permeability crystalline rocks, unlike the sedimentary formations comprising the Medium aquifer. The soil parameters could therefore be adjusted to t the hydrological conditions in the study area by gradually increasing the soil saturated conductivity and the porosity, generating more in ltration and less runoff.
The performance of WASA in simulating the monthly surface water processes was evaluated using the Nash-Sutcliffe e ciency (NSE), calculated for each sub-basin (170, 172 and 186). The best t is the simulation resulting in the highest NSE (close to 1) for at least the two largest sub-basins, which cover the majority of the Medium aquifer. After calibration, the simulated water balance was assessed to ensure that the model results were hydrologically realistic. Given the fact that the study area is located in the Brazilian semi-arid, it was expected that about 70 to 80% of rainfall would be lost by evapotranspiration, 10 to 15% by runoff and about 5 to 10% lost as deep percolation into the deeper groundwater system (Ponce 1995).

Integrated surface water and groundwater modelling
WASA was integrated with MODFLOW, using recharge as the exchange variable. For each sub-basin, WASA yielded a monthly time series of the mean daily recharge to groundwater, which was obtained after stream ow calibration. These values were implemented in MODFLOW using the RCH package, by de ning polygons following the subbasins contours. Three additional polygons were created at the western boundary, where the estimated contribution from the Superior aquifer was added to the WASA recharge rates, resulting in 6 recharge zones ( Figure  4).
Then, MODFLOW simulated the groundwater ow in the Medium aquifer at a monthly scale (su cient to observe seasonal variations) with a monthly varying recharge rate in each recharge polygon. However, the scarcity of eld measurements of hydraulic conductivity and storage coe cients added uncertainty to these model parameters. 3.6 Potential of the São Francisco River transfer for the recharge of the Medium aquifer Once the integrated model was calibrated, CAC could be represented in the model, as a new channel with the SFR2 package of MODFLOW, releasing water to the main streams in the area. This new input of water to the streams might generate recharge into the Medium aquifer by channel transmission losses. By carrying out a scenario analysis, i.e. comparing simulations with and without CAC, the effects of this new water management system on the groundwater ow could be assessed. For the scenario analysis, we consider a context of a prolonged drought period, i.e. low groundwater recharge. Scenario analysis approach has been adopted to study the impact of human interventions on surface/groundwater resources in other drylands (Bauer et al. 2006;Izady et al. 2015).
The CAC route was retrieved from a COGERH shape le and used to create a new SFR2 segment. The channel elevations were estimated from the SRTM raster. The channel width was uniformly set to 20 m, based on design documents. The channel is concreted to preserve the ow of water, thus hydraulic conductivity of the streambed was set to 0. The design in ow of 30 m 3 /s was set at the entrance of the segment. An evaporation value of 1.5 mm/d was also set to account for the evaporation at the water surface. No data regarding the potential amount of water released from the canal to the sub-basin streams were available. Thus, an arbitrary value of 0.29 m³/s (5,000 m³/d) was set as in ow to each of the ve streams crossed by CAC, in order to preserve the ow to the next canal segments. Since the CAC segment considered here is only a part of this water transfer project (145.3 out of 1,250 km), it was important to ensure that the water release simulation from CAC in the Salgado River watershed does not result in a reduction of ow to the next watersheds crossed by this project.

WASA calibration
Calibration-free WASA simulation presented monthly river ow higher than the stream gauge observations. Gradual increase of soil saturated hydraulic conductivity and porosity increased the amount of water in ltrating the soil pro le, therefore decreasing runoff. The optimal t was obtained after increasing the saturated hydraulic conductivity and porosity by 65%. Figure 5 shows the interannual variability of the monthly discharge. One distinctive period from 2003 to 2009 was characterized by higher peak discharge during the wet season (January to May). On the other hand, a much drier period could be identi ed from 2010 until 2016, which was characterized by low peak discharge during the wet season. These two periods can be explained by the interannual variability of precipitation, with the presence of a prolonged meteorological drought in the latter period (not shown).
The simulated river discharge time series ( Figure 5) resulted from the calibration with the highest NSE. Because of the availability of the measured discharge time series, the chosen calibration period for the sub-basins 186 and 172 was 1995-2016, and for the sub-basin 170 was 2000-2016. Non-calibrated simulations from Güntner (2002) resulted in a NSE of -0.54 for the sub-basin 186, which was improved to 0.73 after calibration, and a NSE of 0.34 for the sub-basin 172, which was improved to 0.61. However, we found that the performance for the sub-basin 170 did not improve after calibration, with a NSE of 0.70 without calibration and 0.19 with calibration. This can be explained by the fact that most of this sub-basin area ( Figure 3) is located in the crystalline domain, i.e. outside the Araripe sedimentary basin, which corresponds more closely to the parameterization in the non-calibrated version of WASA. As a result, the recharge to groundwater that is input to MODFLOW might be an overestimation in the case of the sub-basin 170. However, this sub-basin forms only a small proportion of the total area, which leads to still a low value of groundwater recharge, when compared to the recharge calculated from the other sub-basins.

WASA recharge simulation
WASA provides a separate output le giving daily deep groundwater recharge for each sub-basin in m 3 . These values were converted to mm by dividing by the sub-basin area to match the input requirements of the MODFLOW RCH package. As expected from the precipitation and river ow time series, a high interannual variability of groundwater recharge was observed, with a maximum annual recharge of 267 mm simulated for sub-

Water balance of the sub-basins
In order to ensure the reliability of the WASA calibration results, we assessed whether the contribution of each dominant process to the catchment water balance made hydrological sense in the context of the study area (e.g. semi-arid climate, sandy soils and intense rainy seasons). To do so, we summarized the proportional values of the WASA outputs: actual evapotranspiration, surface runoff and recharge regarding rainfall in each sub-basin for the whole period of simulation (Figure 7).
We found that 75-85% of rainfall was lost by evapotranspiration. The remaining part of rainfall was then distributed over the system mainly in surface runoff contributing to river discharge (5 to 18%) and in recharge to groundwater (5 to 10%). These values of actual evapotranspiration, surface runoff and recharge are typical for semi-arid regions like northeast Brazil (Ponce 1995). Other processes involved in the surface water balance were e.g. subsurface runoff, storage contents of the soil layer and interception by plants, representing altogether less than 4% of the total contribution.
Although the overall water balance was similar for the three sub-basins, some differences in the relative contribution of each process could be related to the input rainfall data combined with the characteristics of the study area. For instance, the higher proportion of surface runoff in sub-basins 186 and 170 compared to that from 172 could be explained by higher precipitation in the former. Based on the mean annual precipitation in the three sub-basins (sub-basin 186: 975 mm; sub-basin 172: 737 mm; and sub-basin 170: 869 mm), differences in runoff responses are expected. This explained partly the lower contribution of surface runoff to the water balance in subbasin 172, although this sub-basin is mainly located over a crystalline basement. This lower surface runoff was compensated by a higher actual evapotranspiration rate.

Water balance of the groundwater
The order of magnitude of most of the in ows and out ows of the groundwater system could be quanti ed before MODFLOW application. We present here the overall water balance of the Medium aquifer for the whole period of WASA simulation (1990-2016), based on WASA outputs, available regional data and previous studies (Table 1). The amount of water percolating through the Santana aquiclude to the Medium aquifer was estimated using the average rainfall time series over the Araripe plateau. For the period 1990-2016, this resulted in a mean annual percolation of 6.5 mm in forested areas and 1.3 mm in deforested areas. When applied to an area of 600 km², corresponding to the overlying Santana formation, the estimated in ow to the Medium aquifer (2,166 km 2 ) was 2.3 x10 6 m³/y or 1.1 mm/y. The potential amount of water available from the springs located at the cliff of the Araripe plateau was estimated by calculating the difference between water production and withdrawal. The contribution of these springs to recharge was only 6 x10 5 m³/y or 0.3 mm/y, considering a 7% in ltration rate, which was based on the recharge rates obtained with WASA. The sum of the contribution of the Santana aquiclude percolation and the recharge from springs is equal to 1.4 mm/y, which is the total contribution of the Superior aquifer percolation.
We found that the total contribution of the Superior aquifer percolation could be considered negligible in comparison to the catchment rainfall in ltration. The catchment rainfall in ltration was the most relevant component of the Medium aquifer recharge, being circa 47 times the Superior aquifer percolation on average. The annual recharge was approximately twice the amount of water abstracted. This abstraction rate would indicate a sustainable use of groundwater. However, the results from WASA showed a high interannual variability of groundwater recharge. While periods with high rainfall resulted in recharge reaching 184.7 mm/y like in 2008, annual rainfall was much lower from 2012 to 2016, resulting in a mean annual recharge dropping to 13.9 mm. Therefore, when considering this drought period, the groundwater abstraction becomes unsustainable, as abstraction is now circa twice the amount of recharge.
In fact, groundwater abstraction might be much higher and, consequently, signi cantly unsustainable. The study of Golder-Pivot (2005) estimated an annual increase of water demand of 1.77%, which would represent an increase of the abstraction capacity by 1.3 x10 6 m 3 every year. Moreover, it is expected that during a prolonged drought period, water demand increases, resulting in increased groundwater abstraction.

MODFLOW integrated with WASA: Steady state simulation
First simulations were run in steady state to reproduce the initial head distribution in the study area, i.e. a stationary situation with a limited storage change. The abstraction wells were not included to simulate groundwater ow in which in ow was equal to out ow over one hydrological year. We calibrated the steady state model using PEST to reproduce the initial head distribution observed at monitoring wells on August 1st 2012. The calibration of hydraulic conductivity, recharge rate and hydraulic conductivity of the streambed resulted in the head distribution presented in Figure 8.
The groundwater ow in the Medium aquifer after the steady state calibration showed a general ow direction from S-SW to N-NE. The ow direction was mainly driven by topography, following the downhill slope from the Araripe plateau and converging to the main river network, which are the natural outlets of the aquifer. Channel transmission losses occurred mainly in upstream areas.
The model capability to reproduce the head distribution observed could be analysed from the model residuals. The values of the parameters obtained after the steady state calibration ( Table 2) were quite different from the initial values set for the Medium aquifer from previous studies. For instance, the calibrated hydraulic conductivity for layer 2 and layer 3 (0.1 m/d) seemed to be very low to represent sandy sedimentary formations, which were estimated at about 5 m/d by Mendonça (2001), using pumping test data. In addition, the calibrated recharge was quite different between the recharge zones, whereas the results from WASA were more homogeneous between the sub-basins. In fact, it seems unrealistic to observe such drastic changes. This could be explained by a high contribution of groundwater ow to storage in the initial head distribution. However, besides the uncertainties involved, the low mean residual head obtained after steady state calibration meant this simulation provided satisfactory initial head conditions for the transient calibration of the regional groundwater ow. Here, the recharge was applied as a monthly varying stress entering the system to account for the intra-annual variability caused by distinct seasonality and the interannual variability caused by the highly variable annual precipitation. Groundwater ow was then simulated using a monthly time step from September 2012 to September 2016. As stated earlier, the dataset resulting from the steady state calibration in September 2011 was exported and used as the initial heads for the transient simulation of the groundwater ow in the Medium aquifer.
The main issue was compensating the error introduced by the initial discrepancy in storage, which led to an enormous amount of water leaving the aquifer at the streams in the lowest elevation part of the model domain (NE part, downstream end of the main river network), where the groundwater level was simulated as being several meters above the streambed elevation. Ultimately, we performed a transient calibration of the model, adjusting the hydraulic conductivities of the 3 layers, the hydraulic conductivity of the streambed, and the speci c yield. The model calibration showed satisfactory results with a mean residual head of -1.15 m and a RMS residual head of 10.27 m for the total simulation period, with optimal parameters presented in Table 3.  The comparison between the performance of the WASA-MODFLOW modelling and three previous models is summarized in Table 4, which shows that the model performance was an improvement in steady-state simulation compared to the studies of Golder-Pivot (2005)  this comparison of performances should be considered with caution since the aim of study, size of modelled area, modelling approach and data used for calibration were quite different between the considered studies . was set as the initial condition. The recharge time series for the 2011-2016 simulation, which was a period of severe drought, was replicated for the 2016-2022 scenario, to assess the potential of the water transfer in a context of a prolonged drought. Figure 11 shows that the implementation of CAC led to an increase of the hydraulic head by up to 3 m at the outlet of the Medium aquifer after 5 years (2016-2022 scenario). We also found an average increase of the hydraulic heads over the study area of 0.5 to 1 m. The daily amount of groundwater recharge by channel transmission losses resulting from water release of CAC was estimated as 1.89 x10 4 m 3 , which represented an additional recharge of 6.89 x10 6 m³/y in the context of a prolonged drought. However, it is important to note that the CACdriven recharge input did not compensate for the relevant decrease of the hydraulic heads in areas that do not surround the receiving streams (Figure 11), resulting from the joint effect of low natural recharge (drought) and high groundwater pumping. Therefore, the water transfer would represent an additional source of recharge to the Medium aquifer, but it would not solve alone the issue of unsustainable management of groundwater resources.

Conclusion
We investigated the groundwater recharge of the Medium aquifer (2,166 km²) in the state of Ceará, NE Brazil, and the potential in uence of the São Francisco River inter-basin water transfer on these groundwater resources, using an integrated surface-groundwater modelling approach.
We found that the average groundwater recharge was 66 mm/y for the period of 1990-2016 (27 years), which was 7.7% of the average precipitation (860.3 mm/y). Although the presence of springs in the catchments and the direct Santana aquiclude percolation to the Medium aquifer, the catchment rainfall in ltration represented 97.9% of the total groundwater recharge. The annual recharge was approximately twice the amount of water abstracted (33 mm/y). This abstraction rate would indicate a sustainable use of groundwater. However, the mean annual recharge dropped to 13.9 mm from 2012 to 2016. Therefore, when considering this drought period, the groundwater abstraction becomes unsustainable, as abstraction is now circa twice the amount of recharge, which drove the observed decrease of the groundwater level in the study area. This groundwater level trend was accurately reproduced by MODFLOW in the northwest of the domain, where there is the highest density of pumping wells.
Part of the water transfer system's extension called the Ceará's Water Belt (CAC) was added to the integrated WASA-MODFLOW modelling, using the MODFLOW SFR2 package. After 5 years of scenario simulation, the CACdriven recharge input due to channel transmission losses led to an increase of the hydraulic head of 0.5 to 1 m over the study area or 6.89x10 6 m³/y in a context of a prolonged drought. This result shows that SFR water transfer could decrease the stress on groundwater resources in the Medium aquifer. However, it is important to note that the CAC-driven recharge input did not compensate for the relevant decrease of the hydraulic heads in areas that do not surround the receiving streams, resulting from the joint effect of low natural recharge and high groundwater pumping. Therefore, the water transfer would represent an additional source of recharge to the Medium aquifer, but it would not solve alone the issue of the unsustainable management of groundwater resources.
The CAC implementation in the integrated WASA-MODFLOW modelling yielded encouraging results for the use of the developed modelling as a support tool in decision making. Moreover, an advantage of the integrated WASA-MODFLOW modelling developed here is that it can be easily modi ed to investigate future scenarios, beyond the inter-basin water transfer. For instance, different scenarios of climate change could be investigated by using projected meteorological time series as input in the calibrated version of WASA. Future scenarios of increase in water demand could also be investigated. However, the poorer groundwater model performance in the central and southeast parts of the domain may limit reliable conclusions over the whole study area. The modelling uncertainties may be related to the large area being modelled, poor representation of hydrological uxes in the central part of the model domain, the simpli ed representation of the geology (e.g. no variation of lithology and sediment thickness) and the groundwater abstraction data scarcity. Further work should focus on the reduction of modelling uncertainties.  Location of pumping wells with data and monitoring wells used for MODFLOW calibration. Also shows WASA subbasins and stream gauges used for WASA calibration. Recharge zones in the groundwater model. Zone 1: sub-basin 170 with contribution of Superior aquifer; zone 2: sub-basin 170; zone 3: sub-basin 172; zone 4: sub-basin 172 with contribution of Superior aquifer; zone 5: subbasin 186 with contribution of Superior aquifer; zone 6: sub-basin 186.

Figure 5
Simulated vs observed stream ow discharge (m 3 /s) after model calibration.

Figure 6
Simulated annual groundwater recharge (mm) for each WASA sub-basin and annual average precipitation over the sub-basins.   Evolution of the head distribution (m) for the calibrated transient simulation.

Figure 10
Plots of simulated and observed head (m) at 6 monitoring wells (a. and b. are located in the western part of the model domain where the model reproduces the observations most accurately; c. and d. are located in the central part of the domain where the model performance to reproduce observed heads are poor; e. and f. are located in the east where there is a large discrepancy between model and the heads at the monitoring wells).

Figure 11
Evolution of the head difference (m) between the simulation with CAC active and without CAC at different times of simulation.