4.1 WASA calibration
Calibration-free WASA simulation presented monthly river flow higher than the stream gauge observations. Gradual increase of soil saturated hydraulic conductivity and porosity increased the amount of water infiltrating the soil profile, therefore decreasing runoff. The optimal fit 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 identified 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.
4.2 WASA recharge simulation
WASA provides a separate output file giving daily deep groundwater recharge for each sub-basin in m3. 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 flow time series, a high interannual variability of groundwater recharge was observed, with a maximum annual recharge of 267 mm simulated for sub-basin 186 in 2004 and a minimum annual recharge of less than 1 mm simulated in the sub-basin 172 in 1993, 2012, 2015 and 2016 (Figure 6).
When comparing the values used in previous studies with the mean annual recharge of each sub-basin for the whole calibration period, the recharge rates of 148.2 mm/y and of 196 mm/y estimated in the Golder-Pivot (2005) and Rede de Pesquisa (2007) studies, respectively, are much higher than the results obtained by WASA (186: 86.8 mm/y; 172: 40.2 mm/y; and 170: 47.6 mm/y), probably because the former studies assumed a very high constant groundwater recharge rate from the infiltrated rainfall. The long-term recharge values from WASA are more similar to those calculated by the World bank (2011): 60 mm/y in the western part of the aquifer and 16 m/y in the eastern part. Moreover, the latter study showed higher values of recharge in the western part of the Medium aquifer, which corresponds to the sub-basin 186.
4.3 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 sub-basin 172, although this sub-basin is mainly located over a crystalline basement. This lower surface runoff was compensated by a higher actual evapotranspiration rate.
4.4 Water balance of the groundwater
The order of magnitude of most of the inflows and outflows of the groundwater system could be quantified 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).
Table 1
Water balance components of the Medium aquifer for the whole period of WASA simulation (1990-2016).
Component
|
Value (mm/y)
|
Catchment rainfall infiltration
|
64.6
|
Superior aquifer percolation*
|
1.4
|
(Santana aquiclude percolation)
|
1.1
|
(Recharge from springs)
|
0.3
|
Groundwater abstraction
|
33
|
*The sum of Santana aquiclude percolation and recharge from springs. |
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 inflow to the Medium aquifer (2,166 km2) was 2.3 x106 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 x105 m³/y or 0.3 mm/y, considering a 7% infiltration 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 infiltration. The catchment rainfall infiltration 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, significantly 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 x106 m3 every year. Moreover, it is expected that during a prolonged drought period, water demand increases, resulting in increased groundwater abstraction.
4.5 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 flow in which inflow was equal to outflow 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 flow in the Medium aquifer after the steady state calibration showed a general flow direction from S-SW to N-NE. The flow 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 calibration resulted in a mean residual head of -0.62 m and a root mean squared (RMS) residual head of 9.48 m. The largest residual heads were observed in the central part of the model, with a maximum residual of almost 26 m. The lowest residual heads were mainly observed in the western part of the model, with a minimum residual of 0.03 m observed. The coefficient of correlation of 0.925 between observed and simulated head reflects a satisfactory representation of the groundwater flow at the beginning of the transient simulation.
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 flow 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 flow.
Table 2
Calibrated parameters for the steady state simulation.
Parameter
|
Optimal value (m/d)
|
Optimal value (mm/y)
|
Hydraulic conductivity - layer 1
|
11.8
|
-
|
Hydraulic conductivity - layer 2
|
0.1
|
-
|
Hydraulic conductivity - layer 3
|
0.1
|
-
|
Hydraulic conductivity of the streambed
|
5.0
|
-
|
Recharge rate - Zone 1
|
-
|
1.4
|
Recharge rate - Zone 2
|
-
|
188.7
|
Recharge rate - Zone 3
|
-
|
31.1
|
Recharge rate - Zone 4
|
-
|
152.9
|
Recharge rate - Zone 5
|
-
|
365.0
|
Recharge rate - Zone 6
|
-
|
18.7
|
4.6 MODFLOW integrated with WASA: Transient calibration
Simulation of the groundwater flow in transient mode requires imposing stresses on the system varying with time. 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 flow 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 flow 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 specific 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.
Table 3
Optimal parameters of the transient calibration
Parameter
|
Optimal value (m/d)
|
Hydraulic conductivity - Layer 1
|
0.79
|
Hydraulic conductivity - Layer 2
|
1.03
|
Hydraulic conductivity - Layer 3
|
3.10
|
Hydraulic conductivity of the streambed
|
2.21
|
Specific yield
|
0.16
|
Figure 9 shows the changes of hydraulic head over the simulation period. No large change of regional flow direction or water level could be observed, in accordance with the low recharge computed for the period 2012-2016 (Figure 6). The main change from Figure 9a to Figure 9b was a decrease of the groundwater dome presented in the central part of the model domain. This change may be attributed to the overestimation of groundwater level in the initial head distribution. A steady decrease of the hydraulic head could be observed in the model over the period of simulation in the northwestern part of the model, which could be explained by the high density of pumping wells in this area (Figure 3).
The capability of the calibrated model to reproduce the intra-annual and interannual variations of the hydraulic head at monitoring wells differed in the study area. The model accurately reproduced the interannual head variations in the northwest (corresponding to WASA sub-basin 186), such as shown in Figure 10a-b. On the other hand, performance was much poorer in the central part (equivalent to sub-basin 170) as shown in Figure 10c-d, and in the southeast (sub-basin 172) as shown in Figure 10e-f. This could be partly explained by the performance of WASA in estimating recharge in the three sub-basins (see NSE of the calibrated WASA model). WASA performed better at reproducing discharge at the outlet of sub-basin 186. The groundwater modelling uncertainties may also be related to the simplified representation of the geology (e.g. no variation of lithology and sediment thickness) and the groundwater abstraction data scarcity.
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) and the World Bank (2011). For transient simulations, our model performed much better than that from Rede de Pesquisa (2007). Also, for transient simulations, the WASA-MODFLOW modelling performance was similar to the model performance in the study of the World Bank. However, 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 .
Table 4
Summary of model performances of the MODFLOW-WASA and three previous models
Simulation
|
Residuals (m)
|
Integrated MODFLOW-WASA
|
Golder-Pivot (2005)
|
Rede de Pesquisa (2007)
|
World Bank (2011)
|
Steady-state
|
Mean residual head
|
-0.62
|
11.8
|
-
|
-2.65
|
RMS residual head
|
9.48
|
21.15
|
-
|
14.13
|
Maximum head difference
|
25.9
|
40.8
|
-
|
-34.9
|
Minimum head difference
|
0.03
|
1.1
|
-
|
-23
|
Transient
|
Mean residual head
|
-1.15
|
-
|
-
|
-0.51
|
RMS residual head
|
10.27
|
-
|
27.88 (CPRM wells); 14.4 (COGERH wells)
|
9.757
|
Maximum head difference
|
25.9
|
-
|
-
|
29.58
|
Minimum head difference
|
0.03
|
-
|
-
|
-1.35
|
4.7 Simulation of the São Francisco River transfer for the recharge of the Medium aquifer
CAC was added to the modelling as a new SFR2 segment. For small discharges, such as the CAC releases, we assumed that flows were concentrated in streambed paths that have much lower hydraulic conductivity than the total streambed area and river banks (Lange 2005; Souza and Costa 2014). Therefore, assuming the presence of clogging layers with hydraulic properties similar to silty clay soils, the hydraulic conductivity of the streambed for the five streams receiving water from CAC was lowered from 2.21 m/d (calibrated value) to 0.01 m/d, which avoided having all the water released infiltrating into the aquifer in the first stream cells. The potential of the SFR water transfer-CAC implementation was investigated by running the calibrated MODFLOW in transient mode for 5 years beyond the end of the head time series, i.e. a scenario from August 2016 to January 2022. We compared the scenario results with and without CAC. The head distribution of the calibrated transient model for August 2016 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 x104 m3, which represented an additional recharge of 6.89 x106 m³/y in the context of a prolonged drought. 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 (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.