Impact of Climate Change on Groundwater Resources in the Klela Basin, Southern Mali

Investigations of groundwater resources in order to understand aquifer system behavior are vital to the inhabitants of the Klela Basin, Mali, because groundwater is the only permanent water resource and is used for drinking water and irrigation. Due to climate change, this vital resource is being threatened. Therefore, MODFLOW was applied in this study to simulate groundwater dynamics. The aim of this study was to evaluate the impact of climate change on groundwater resources in the Klela basin using the RCP4.5 (Representative Concentration Scenario 4.5 W/m 2) climate scenario. Climatological, geological, hydrogeological, hydraulic and demographic data were collected and used as model input data. Groundwater recharge was estimated to be approximately 165.3 mm/year using the EARTH (Extended model for Aquifer Recharge and soil moisture Transport through the unsaturated Hardrock) model. Recharge was then used as groundwater model input. The sandstone aquifer in the study area was simulated in steady and transient conditions. The results showed that hydraulic conductivity values varied from 1.1 to 13.9 m/day. The model was used for scenario quantification after model calibration and verification using three different piezometer data sets. The results of the simulated MODFLOW model showed a decrease in groundwater levels over time.


Introduction
Groundwater resources provide important sources of drinking water throughout the world, especially in developing countries where most human activities (e.g., domestic water use, small-scale irrigation, etc.) depend on groundwater.The investigated basin in Mali is an agriculture zone, growing cotton and rice during the rainy season and potatoes during the dry season.The surface area used for agriculture has been extended in the basin to assure food security for the growing population.In addition, most of the rural population uses water from traditional wells that are recharged by rainfall and dry out a few months after the recharge period.Moreover, water sustainability in the town of Sikasso's supply system remains a challenge due to groundwater exploitation in the context of increased water demand and climate variability and change.Therefore, better understanding and management of this vital resource are the principal requirements of responsible water use.Because of the growing population and climate change, a scientific tool should be used to assess present and future difficulties related to water supply and demand in the study area.MODFLOW is one of the most popular hydrologic models in the world and can achieve the objectives of this study.MODFLOW is a modular, three-dimensional finite-difference groundwater model that was developed by the United States Geological Survey (USGS) [1].It uses the groundwater flow equation, which is a combination of the continuity equation and Darcy's law [2].Over the past few years, the application of MODFLOW to describe and predict aquifer system behavior has increased significantly [3].MODFLOW has been used by many authors [4][5][6][7] throughout the world to simulate groundwater dynamics; however, few studies have been conducted in West Africa, especially in Mali, where some studies [8,9] have predicted rainfall decreases.For example, Bricquet et al. [10] showed a reduction in groundwater storage in the Bani basin (which contains several subbasins, including the Klela basin) due to decreases in rainfall and runoff.Additionally, Bokar et al. [11] demonstrated that groundwater levels decreased from 2 to 15 cm per year from 1940 to 2008 in the Kolondieba catchment, a subbasin of the Bani basin.The interaction between surface water and groundwater has been the topic of many discussions [12][13][14] in Mali, which have concluded that the river system drains the aquifer in the southern part of the country.Until now, no study has been performed in the Klela basin to confirm or refute these results or to quantify groundwater storage.
The partitioning of rainfall into groundwater (recharge) and the interaction between surface water and groundwater is unknown in the Klela basin.However, knowledge of these processes is required to facilitate water resources management in the Klela basin in the context of climate change.
Therefore, the purpose of this study is to evaluate the impact of climate change on groundwater resources in the Klela basin using the RCP4.5 climate scenario.EARTH (Extended model for Aquifer Recharge and soil moisture Transport through the unsaturated Hardrock) and PMWIN 5.3 (Processing MODFLOW for Windows) developed by Chiang and Kinzelbach [15] is applied to estimate groundwater recharge and simulate groundwater behavior under steady and transient state conditions.

Study Area
The Klela basin is located in the Sikasso Region in southern Mali, with a small part located in western Burkina Faso.Situated between 5 ¥ 55 I 58.8" and 5 ¥ 16 I 12" longitude and 11 ¥ 40 I 58.8" and 10 ¥ 59 I 45.6" latitude, it covers approximately 3680 km 2 .The elevation of the Klela Basin decreases from south to north, from 748 m above sea level (masl; south) to 305 masl (north) (Figure 1).The main river of the area is the Lotio River, which has periodic flows and is a tributary of the Banifing, which is a tributary of the Bani River basin.The only assured water supply is groundwater, and all the inhabitants depend directly or indirectly on this resource.The minimum monthly temperature varies from 12.3 ¥ C to 26.8 ¥ C, the maximum temperature ranges from 28.8 ¥ C to 39.9 ¥ C and the average temperature is 27.4 ¥ C (Figure 2).The climate belongs to the Soudano-Sahelian zone and is characterized by a dry season (November-May) dominated by a dry wind from the Sahara (harmattan) and a rainy season (June-October) with a wet wind from the Guinean Gulf (monsoon) [16].The mean annual precipitation is approximately 1131 mm/year.The average annual potential evapotranspiration is approximately 2060 mm/year (Figure 2) based on the Blaney-Criddle method.The geology of the study area is mainly constituted by sandstone and dolerite formations (Figure 3).However, the Tabular Infra-Cambrian fractured sandstone formation is predominant.Groundwater can be found in the fractured bedrock at depths of up to 100 m.The recent formations that cover the sandstone bedrock are constituted by clay, sand and laterite, with a maximum covering thickness of 64 m.

Methodology
The data and methods used in this modeling study are summarized in Figure 4.

Methodology
The data and methods used in this modeling study are summarized in Figure 4.

Methodology
The data and methods used in this modeling study are summarized in Figure 4.The hydrogeological data (e.g., drilling logs, pumping tests, water levels, etc.) were provided by the DRH (Direction Regionale de l'Hydraulique de Sikasso), and climate data are from the DNM (Direction Nationale de la Meteorologie).GIS was used to draw all the maps and interpolate data using spatial analysis.PMWIN version 5.3 was used to model the aquifer behavior in the study area.

EARTH Model
EARTH (Extended model for Aquifer Recharge and soil moisture Transport through the unsaturated Hardrock) is a lumped parametric model that can be used to estimate groundwater recharge.It was developed by van der Lee et al. in 1989 in a project in Botswana and tested under different climate conditions [19].It is a one-dimensional model designed for semi-arid areas.EARTH combines direct and indirect methods.The direct method describes the recharge process from the atmosphere to the soil zone.The indirect method uses groundwater level fluctuations as indicators of the recharge process.The direct method determines recharge using physical processes above the groundwater table, whereas the indirect method calculates the groundwater level that results from recharge estimated by the direct method [19].The EARTH model has five modules: MAXIL, SOMOS, LINRES, SUST and SAFLOW (Figure 5).MAXIL, SOMOS and LINRES comprise the direct part of the model and simulate recharge in the unsaturated zone.They can be calibrated with measured time series of soil moisture data.SATFLOW is the indirect part of the model and calculates the water level from the recharge simulated by the direct part.SUST calculates the surface runoff.The hydrogeological data (e.g., drilling logs, pumping tests, water levels, etc.) were provided by the DRH (Direction Regionale de l'Hydraulique de Sikasso), and climate data are from the DNM (Direction Nationale de la Meteorologie).GIS was used to draw all the maps and interpolate data using spatial analysis.PMWIN version 5.3 was used to model the aquifer behavior in the study area.

EARTH Model
EARTH (Extended model for Aquifer Recharge and soil moisture Transport through the unsaturated Hardrock) is a lumped parametric model that can be used to estimate groundwater recharge.It was developed by van der Lee et al. in 1989 in a project in Botswana and tested under different climate conditions [19].It is a one-dimensional model designed for semi-arid areas.EARTH combines direct and indirect methods.The direct method describes the recharge process from the atmosphere to the soil zone.The indirect method uses groundwater level fluctuations as indicators of the recharge process.The direct method determines recharge using physical processes above the groundwater table, whereas the indirect method calculates the groundwater level that results from recharge estimated by the direct method [19].The EARTH model has five modules: MAXIL, SOMOS, LINRES, SUST and SAFLOW (Figure 5).MAXIL, SOMOS and LINRES comprise the direct part of the model and simulate recharge in the unsaturated zone.They can be calibrated with measured time series of soil moisture data.SATFLOW is the indirect part of the model and calculates the water level from the recharge simulated by the direct part.SUST calculates the surface runoff.Details of the five modules are given by van der Lee and Gehrels [19] and are briefly described as follows: 1. MAXIL (Maximal Interception Loss) estimates the water intercepted by the Earth`s surface.2. SOMOS (Soil Moisture Storage) describes water storage in the root zone and calculates changes in soil moisture storage by removing actual evapotranspiration, percolation, evaporation and surface runoff to determine the infiltrated precipitation.3. SUST (Surface Storage) calculates ponding and runoff.It uses SUSTmax, which is the maximum ponding volume that can be stored at the surface.If the amount of ponding is greater than SUSTmax, then runoff occurs.4. LINRES (Linear Reservoir Routing) redistributes percolation (output of SOMOS) over time in the unsaturated hard rock or soil beneath the root zone using the parametric transfer function.Details of the five modules are given by van der Lee and Gehrels [19] and are briefly described as follows: 1.
MAXIL (Maximal Interception Loss) estimates the water intercepted by the Earth's surface.

2.
SOMOS (Soil Moisture Storage) describes water storage in the root zone and calculates changes in soil moisture storage by removing actual evapotranspiration, percolation, evaporation and surface runoff to determine the infiltrated precipitation.

3.
SUST (Surface Storage) calculates ponding and runoff.It uses SUSTmax, which is the maximum ponding volume that can be stored at the surface.If the amount of ponding is greater than SUSTmax, then runoff occurs.4.
LINRES (Linear Reservoir Routing) redistributes percolation (output of SOMOS) over time in the unsaturated hard rock or soil beneath the root zone using the parametric transfer function.

SATFLOW (Saturated Flow
) is the last part of the EARTH model.It is a simple one-dimensional parametric model that uses the outputs from the previous modules.SATFLOW calculates the groundwater level using the recharge estimated by the direct part of the model.
A general mathematical explanation of the EARTH model derived from the equation of de Vries 1974 cited in [20] is given by [21], as follows: where fh/ft is the change in water level head during one month (m/month), h is groundwater level (m), DR is the drainage resistance (a lumped, site-specific parameter), R is the recharge (m/month) and S is the specific yield.
The numerical solution of Equation ( 1) (linear transfer function) yields Equation ( 2): where DR = drainage resistance (L 2 /βT), L = length of flow path, β = 2 for radial flow or β = 4 for parallel flow and T = transmissivity.Monthly observed groundwater levels from two piezometers and rainfall data were used as model inputs.Groundwater levels collected during two years (2012-2013) at wells were used for calibration.

Groundwater Model
In MODLFOW, the three-dimensional movement of groundwater with constant density through porous media is described by the following partial-differential equation [22]: where K xx , K yy , and K zz are of the hydraulic conductivities in the x-, yand z-directions (LT ¡1 ); h is potentiometric head (L); w is a volumetric flux per unit volume and represents sources and/or sinks of water (T ¡1 ); S s is the specific storage of the porous material (L ¡1 ); and t is time (T).

Spatial Discretization
In numerical groundwater modeling, an aquifer system is replaced by a discretized domain consisting of an array of nodes and associated finite difference blocks (cells).In PMWIN, the model domain is divided into a rectangular mesh comprising columns, rows and layers [15].
In this study, the model area was discretized into 239 columns and 256 rows, resulting in a total grid cell number of 61,184, with a grid size of 300 m ¢ 300 m, taking into account the aquifer extension and hydrogeological data availability.The total number of active grid cells describing the active model domain was 40,992.The aquifer of the study area was modeled under an unconfined condition represented by a single numerical layer.

Aquifer Geometry
The aquifer geometry step establishes the top and bottom of an aquifer.Determination of aquifer geometry is critical in numerical modeling because it influences model calibration results [23].Drilling log data and an SRTM-based Digital Elevation Model (DEM) from HydroSHED.orgwere used to estimate the top and bottom of the aquifer and aquifer thickness in the study area.The DEM was used as the aquifer top and the aquifer thickness was arbitrarily set to 300 m.
A scatter plot of surface elevation vs water table (Figure 6) shows that in most wells, the water table is up to 20 m below ground level.This assumes a possible strong interaction between groundwater and surface water bodies within the basin.

Aquifer Properties
The hydraulic properties of the aquifer, including transmissivity, hydraulic conductivity and storativity, determine water movement through the aquifer as well as storage in the aquifer.In this study, the horizontal hydraulic conductivity value was calibrated until a good correlation between simulated and observed heads was obtained.For this reason, the model domain was divided into six hydraulic conductivity zones.
Storativity data were not available in the study area.Therefore, specific yield data were estimated using the Cooper-Jacob method [24] and used for transient calibration.

Boundary Conditions
Boundary conditions specify how an aquifer interacts with the environment outside the model domain [25].Selection of boundary conditions is a key step in model design because they significantly control flow patterns [26,27].Because there is not enough hydrogeological information in the study area, the topographical divide was considered a no-flow boundary condition.Thus, there was no hydrogeological interaction between the aquifer inside and outside the basin.Drainage of the basin was limited to groundwater-surface water interaction.

Initial Hydraulic Head
Historical water level data from more than 100 boreholes and data interpolated using the Kriging method were assigned to the model as initial hydraulic heads in the steady-state simulation.The head distribution computed using the calibrated steady-state simulation was used as the starting head distribution in the transient simulation.

Groundwater Sources/Sinks
As previously explained, the average recharge value (0.00044 m/day) estimated by the EARTH model was assigned to the entire model domain.

Aquifer Properties
The hydraulic properties of the aquifer, including transmissivity, hydraulic conductivity and storativity, determine water movement through the aquifer as well as storage in the aquifer.In this study, the horizontal hydraulic conductivity value was calibrated until a good correlation between simulated and observed heads was obtained.For this reason, the model domain was divided into six hydraulic conductivity zones.
Storativity data were not available in the study area.Therefore, specific yield data were estimated using the Cooper-Jacob method [24] and used for transient calibration.

Boundary Conditions
Boundary conditions specify how an aquifer interacts with the environment outside the model domain [25].Selection of boundary conditions is a key step in model design because they significantly control flow patterns [26,27].Because there is not enough hydrogeological information in the study area, the topographical divide was considered a no-flow boundary condition.Thus, there was no hydrogeological interaction between the aquifer inside and outside the basin.Drainage of the basin was limited to groundwater-surface water interaction.

Initial Hydraulic Head
Historical water level data from more than 100 boreholes and data interpolated using the Kriging method were assigned to the model as initial hydraulic heads in the steady-state simulation.The head distribution computed using the calibrated steady-state simulation was used as the starting head distribution in the transient simulation.

Groundwater Sources/Sinks
As previously explained, the average recharge value (0.00044 m/day) estimated by the EARTH model was assigned to the entire model domain.
Groundwater abstraction was based on household (population and livestock) water consumption and irrigation water needs from pumping wells.

Streamflow Routing
In this study, the streamflow-routing package [28] was used to simulate the exchange of fluxes (sources and sinks) between groundwater and surface water.The Streamflow-Routing package is designed to account for the amount of flow in streams and to simulate the interactions between surface streams and groundwater [29].This package requires parameters such as streambed hydraulic conductance, elevation of the streambed top, elevation of the streambed bottom, width of the stream channel, slope of the stream channel, manning's roughness coefficient, etc.For more details, see [15].
In MODFLOW, the flow from streams to an aquifer (losing stream) or from an aquifer streams (gaining stream) is calculated depending on the groundwater and streamwater tables.When groundwater table is above the streamwater table, the volumetric flux between streams and an aquifer is calculated using Equation (4) [30]: where is the length of the river within a cell, W (L) is the width of the stream, M (L) is the thickness of the streambed, h S (L) is the head in the stream and h (L) is the groundwater head.
When the groundwater table h is below the elevation of the streambed bottom, the volumetric infiltration flux from the river to the aquifer is calculated using Equation (5).
The value of C STR was determined to be 500 m 2 /day.The hydraulic conductivity of the streambed is unknown; therefore, it was set to 0.5 m/day based on the dominant soil material in the study area.The length L within the cell is 300 m, the width is 10 m, and the streambed thickness is 3 m.

EARTH Model
The calibrated results in both piezometers are shown in Figure 7.The correlation coefficients (R 2 ) exhibit a good correlation between simulated and observed groundwater levels.A trial and error calibration process was performed as follows.First, the specific yield was fixed and the other parameters, such as drainage resistance and percentage of rainfall (%R), were adjusted to obtain the best estimation between the observed and simulated groundwater levels (Table 1).To evaluate the model quality, three efficiency criteria were used, namely, the correlation coefficient, root mean square error (RMSE) and Nash-Sutcliffe coefficient (Table 2).

Calibration
During the calibration process, model input parameters are changed so that the model output matches measured values [31].Model calibration can be performed under steady-state and transient conditions.In this study, a two-step procedure was applied.In the first step, a steady-state simulation was used to analyze the system and generate consistent initial conditions for the dynamic simulation, which are required to perform the climate change impact study.

Steady-State Simulation
Mean water table data collected by the DRH Sikasso in 1985 (Direction Regionale de l`Hydraulique de Sikasso) and data measured in the field in 12 boreholes in 2014 were used to calibrate the steady-state simulation model.This calibration was first performed via a trial-and-error procedure.Then, the output was optimized using the PEST code integrated in PMWIN.The model performance yielded R 2 = 0.98, RMSE = 5.67, and Nr = 0.86, exhibiting a good calibration.The scatter plot comparison between simulated and observed groundwater levels is shown in Figure 8.The output heads generated by the model (Figure 9) were used as initial conditions in the transient simulation.The calibration result show that hydraulic conductivity varies from 1.1 m/day to 13.9 m/day.

Calibration
During the calibration process, model input parameters are changed so that the model output matches measured values [31].Model calibration can be performed under steady-state and transient conditions.In this study, a two-step procedure was applied.In the first step, a steady-state simulation was used to analyze the system and generate consistent initial conditions for the dynamic simulation, which are required to perform the climate change impact study.

Steady-State Simulation
Mean water table data collected by the DRH Sikasso in 1985 (Direction Regionale de l'Hydraulique de Sikasso) and data measured in the field in 12 boreholes in 2014 were used to calibrate the steady-state simulation model.This calibration was first performed via a trial-and-error procedure.Then, the output was optimized using the PEST code integrated in PMWIN.The model performance yielded R 2 = 0.98, RMSE = 5.67, and N r = 0.86, exhibiting a good calibration.The scatter plot comparison between simulated and observed groundwater levels is shown in Figure 8.The output heads generated by the model (Figure 9) were used as initial conditions in the transient simulation.The calibration result show that hydraulic conductivity varies from 1.1 m/day to 13.9 m/day.

Transient Simulation
Only the trial and error method was used in the transient simulation due to technical problems with PEST.The groundwater level of the steady-state simulation was used as the initial hydraulic head in the transient simulation (Figure 9).For the transient calculations, the streamflow-routing package was used to simulate the interactions between streams and groundwater.In the same manner, sources (recharge) and sinks (wells) were introduced into the model.Specific yield was set to 0.06 for the whole model domain.
The model was run from November 2010 to November 2014 (48 months).Therefore, the simulation time was divided into 48 stress periods (with 32 dry and 16 wet periods), and each stress period represents one month.The length of a stress period was divided into days, resulting in total time steps over 48 months and a total simulation time of 1,460 days.
Mean monthly groundwater table data from three piezometers (F7, F15 and F18) were used to calibrate the model.In total, five piezometers were monitored weekly in the study area, but only three were chosen for modeling purposes.This choice was based on the data quality and data length.The quality of the calibration results is shown in Figure 10.Four efficiency criteria (Mean Absolute Error (MAE), Root Mean Square Error (RMSE), the Nash-Sutcliffe coefficient (E) and the Correlation Coefficient (R 2 )) (Table 3) were used to evaluate the model accuracy.The model efficiency results reflect an acceptable model calibration.However, F7 is not as well simulated as the others.

Transient Simulation
Only the trial and error method was used in the transient simulation due to technical problems with PEST.The groundwater level of the steady-state simulation was used as the initial hydraulic head in the transient simulation (Figure 9).For the transient calculations, the streamflow-routing package was used to simulate the interactions between streams and groundwater.In the same manner, sources (recharge) and sinks (wells) were introduced into the model.Specific yield was set to 0.06 for the whole model domain.
The model was run from November 2010 to November 2014 (48 months).Therefore, the simulation time was divided into 48 stress periods (with 32 dry and 16 wet periods), and each stress period represents one month.The length of a stress period was divided into days, resulting in total time steps over 48 months and a total simulation time of 1,460 days.
Mean monthly groundwater table data from three piezometers (F7, F15 and F18) were used to calibrate the model.In total, five piezometers were monitored weekly in the study area, but only three were chosen for modeling purposes.This choice was based on the data quality and data length.The quality of the calibration results is shown in Figure 10.Four efficiency criteria (Mean Absolute Error (MAE), Root Mean Square Error (RMSE), the Nash-Sutcliffe coefficient (E) and the Correlation Coefficient (R 2 )) (Table 3) were used to evaluate the model accuracy.The model efficiency results reflect an acceptable model calibration.However, F7 is not as well simulated as the others.

Water Budget
One of the best ways to assess model simulation quality is by analyzing the water budget.Although a good water balance cannot guarantee a good simulation, a bad water balance indicates problems in the model [32].In this study, the percentage discrepancy in all the stress periods of the model is nearly zero.Thus, the model equations have been correctly solved [15].The average annual  One of the best ways to assess model simulation quality is by analyzing the water budget.Although a good water balance cannot guarantee a good simulation, a bad water balance indicates problems in the model [32].In this study, the percentage discrepancy in all the stress periods of the model is nearly zero.Thus, the model equations have been correctly solved [15].The average annual water budget is summarized in Table 4.Because the boundary condition is a no flow boundary, the only water source entering the aquifer is recharge.The annual recharge is slightly greater than groundwater seepage (Table 4), i.e., the amount of water that is entering the aquifer is almost the same as the amount exiting the aquifer and entering surrounding streams, which are the discharge points of the study area.According to the boundary condition, groundwater can leave the catchment only via the stream interaction.Inflow from the stream is calculated from the interaction between surface water and groundwater throughout the river system.The mean annual water budget shows that groundwater storage is slightly decreasing by 39,169,133 m 3 /year (approximately 10.6 mm/year) (calculated from Table 4).
The seasonal recharge from July to October was constant and represented the important water inflow entering the aquifer.The stream leakage that is leaving the aquifer and entering streams decreases during the dry season and increases during the rainy season (Figure 11). Figure 12 shows that during the dry season, groundwater storage decreases.It then increases during the rainy season.The maximum value of storage occurs in November.
boundary condition, groundwater can leave the catchment only via the stream interaction.Inflow from the stream is calculated from the interaction between surface water and groundwater throughout the river system.The mean annual water budget shows that groundwater storage is slightly decreasing by 39,169,133 m 3 /year (approximately 10.6 mm/year) (calculated from Table 4).
The seasonal recharge from July to October was constant and represented the important water inflow entering the aquifer.The stream leakage that is leaving the aquifer and entering streams decreases during the dry season and increases during the rainy season (Figure 11). Figure 12 shows that during the dry season, groundwater storage decreases.It then increases during the rainy season.The maximum value of storage occurs in November.

Scenario Quantification
One of the most important objectives of our groundwater modeling study is to use the scenario data to predict the future behavior of the aquifer system.Thus, the model was used to quantify groundwater levels considering scenario data after model calibration and verification using three different piezometer data sets.
Monthly precipitation data from the RCP4.5 (Representative Concentration Scenario 4.5 W/m 2 ) scenario were taken from the GCM ECHAM, downscaled to a 0.44° resolution by the Swedish Meteorological Service (Sveriges Meteorologiska och Hydrologiska Institute, SMHI) and provided by the CORDEX (Coordinated Regional Climate Downscaling Experiment) initiative.They were used to calculate the groundwater recharge (Figure 13), which was used as model input to quantify future groundwater levels in the Klela basin from June 2010 to November 2050.RCP4.5 is an intermediate pathway and a stabilization scenario based on the assumption that all countries around the world implement emissions mitigation policies [33,34].
The Thornthwaite Monthly Water Balance model [35], which is based on monthly temperature and rainfall data, was used for this recharge estimation.The model was calibrated to fit the groundwater recharge previously simulated using EARTH.To calibrate this model, direct runoff was set to zero because local scale surface runoff infiltrates on the way to the river system.There is no direct runoff contribution to the discharge of the basin.The soil moisture storage capacity was calibrated to be 375 mm by comparing the groundwater recharge of the two models (EARTH and Thornthwaite).All other parameters were unchanged.More details on this model can be found in [35].After calibration, the model was applied to simulate groundwater recharge in the past (1970-2010) and the future (2010-2050).The overall result shows that recharge decreases from 1970 to 2050 and severe droughts occur from 2030 to 2037 (Figure 13).
The results of the simulated MODFLOW model reveal a decrease in groundwater level (GWL) over time (Figure 14).From June 2010 to November 2050, GWL decreases from 365 mamsl to 348 mamsl, or to approximately 17 m, in piezometer F7.During the same period, F18 decreases from 356 mamsl to 338 mamsl, or approximately 18 m.F15 decreases from 348 mamsl to 336 mamsl, or approximately 11 m.The overall mean of all the piezometer drops was approximately 15 m.Furthermore, an important decline in groundwater storage in the 2030s was modeled in all the piezometers considered in this study.Groundwater level reductions may intensify if water extraction increases over time.

Scenario Quantification
One of the most important objectives of our groundwater modeling study is to use the scenario data to predict the future behavior of the aquifer system.Thus, the model was used to quantify groundwater levels considering scenario data after model calibration and verification using three different piezometer data sets.
Monthly precipitation data from the RCP4.5 (Representative Concentration Scenario 4.5 W/m 2 ) scenario were taken from the GCM ECHAM, downscaled to a 0.44 ¥ resolution by the Swedish Meteorological Service (Sveriges Meteorologiska och Hydrologiska Institute, SMHI) and provided by the CORDEX (Coordinated Regional Climate Downscaling Experiment) initiative.They were used to calculate the groundwater recharge (Figure 13), which was used as model input to quantify future groundwater levels in the Klela basin from June 2010 to November 2050.RCP4.5 is an intermediate pathway and a stabilization scenario based on the assumption that all countries around the world implement emissions mitigation policies [33,34].
The Thornthwaite Monthly Water Balance model [35], which is based on monthly temperature and rainfall data, was used for this recharge estimation.The model was calibrated to fit the groundwater recharge previously simulated using EARTH.To calibrate this model, direct runoff was set to zero because local scale surface runoff infiltrates on the way to the river system.There is no direct runoff contribution to the discharge of the basin.The soil moisture storage capacity was calibrated to be 375 mm by comparing the groundwater recharge of the two models (EARTH and Thornthwaite).All other parameters were unchanged.More details on this model can be found in [35].After calibration, the model was applied to simulate groundwater recharge in the past (1970-2010) and the future (2010-2050).The overall result shows recharge decreases from 1970 to 2050 and severe droughts occur from 2030 to 2037 (Figure 13).Comparing the results of the EARTH model in this study to others carried out throughout the world in arid and semi-arid areas, the results are reasonably similar.For example, Usher et al. [36] applied the EARTH model in a fractured sandstone aquifer in South Africa (Kalkveld area) and obtained a recharge rate of 7.4% of the mean annual precipitation of 559 mm.In a similar study, Baalousha [37] applied the same model in Palestine in calcareous sandstone and sand aquifers and computed the recharge to be 36.95% of the average annual rainfall of 321 mm.Values between 7.2% and 26.5% of mean annual rainfall were also noted by Wang et.al.[38] in China in unsaturated sand,  Comparing the results of the EARTH model in this study to others carried out throughout the world in arid and semi-arid areas, the results are reasonably similar.For example, Usher et al. [36] applied the EARTH model in a fractured sandstone aquifer in South Africa (Kalkveld area) and obtained a recharge rate of 7.4% of the mean annual precipitation of 559 mm.In a similar study, Baalousha [37] applied the same model in Palestine in calcareous sandstone and sand aquifers and computed the recharge to be 36.95% of the average annual rainfall of 321 mm.Values between 7.2% and 26.5% of mean annual rainfall were also noted by Wang et.al.[38] in China in unsaturated sand, Comparing the results the EARTH model in this study to others carried out throughout the world in arid and semi-arid areas, the results are reasonably similar.For example, Usher et al. [36] applied the EARTH model in a fractured sandstone aquifer in South Africa (Kalkveld area) and obtained a recharge rate of 7.4% of the mean annual precipitation of 559 mm.In a similar study, Baalousha [37] applied the same model in Palestine in calcareous sandstone and sand aquifers and computed the recharge to be 36.95% of the average annual rainfall of 321 mm.Values between 7.2% and 26.5% of mean annual rainfall were also noted by Wang et.al.[38] in China in unsaturated sand, with rainfall ranging from 511 to 557 mm.No previous studies have been conducted in the Klela basin using the EARTH model; Henry [39] applied the HELP (Hydrologic Evaluation of Landfill Performance) model in southern Mali, including the Klela basin, and obtained an annual recharge of 12.6% of the mean rainfall of 1049 mm.
Regarding the MODFLOW results, despite a lack of hydrogeological data, they appear to be acceptable based on the model performance evaluation results.One of the most important uncertainties in this modeling is the establishment of boundary conditions.In this study, the boundary conditions were set as no flow boundaries.In reality, this is very difficult to prove.For example, Belay [40] applied the PMWIN model to establish a water balance for Lake Beseka in Ethiopia and concluded that the choice of boundary conditions has a larger effect on the dynamics than does precipitation.Therefore, neglecting lateral outflow in groundwater modeling may affect the water balance.

Conclusions and Recommendations
The behavior of groundwater in the Klela basin has been evaluated in the context of climate change.This evaluation was performed using the RCP4.5 scenario data and applying the groundwater model MODFLOW and the Thornthwaite model for groundwater recharge.In this scenario, groundwater storage will decrease by approximately 10.6 mm/year.Furthermore, hydraulic heads generated by the model displayed a reduced groundwater level, increasing groundwater availability stresses in areas where agriculture is the main activity.By considering the scenario analysis in which the groundwater recharge decreases significantly compared to past recharge, the reduction of groundwater storage is important, especially in the 2030s, when severe drought events are simulated due to climate change effects.These events threaten access to groundwater for irrigation, agriculture and domestic use.
Groundwater recharge, the most important parameter in groundwater modeling, calculated using the EARTH model could have been over-or underestimated in some areas of the basin because the piezometers were not spatially distributed in a way that adequately allowed determination of spatial recharge patterns.Moreover, the true boundary of the basin was not known, potentially increasing the uncertainties in the results of this study.Despite these shortcomings, the results are reliable and can serve as a basis for watershed management.
In conclusion, a focus should be placed on determination of the hydrogeological situation in future studies.More piezometers should be implemented in the Klela basin to permit permanent monitoring of groundwater resources in the context of global change.
This study also shows that MODFLOW, which is widely used to model groundwater resources, can be applied in the Klela basin to quantify the influence of climate change on groundwater resources and can support water resources management.

Figure 1 .
Figure 1.Geographical location of the Klela Basin in southern Mali.

Figure 1 .Figure 1 .
Figure 1.Geographical location of the Klela Basin in southern Mali.

Figure 4 .
Figure 4. Methodology and data used in this study.

Figure 4 .
Figure 4. Methodology and data used in this study.

Figure 4 .
Figure 4. Methodology and data used in this study.

Figure 6 .
Figure 6.Scatter plot between surface elevation and water table.The blue stars are water table field data, the green line represents surface elevation and the red line corresponds to surface elevation minus 20 m.

Figure 6 .
Figure 6.Scatter plot between surface elevation and water table.The blue stars are water table field data, the green line represents surface elevation and the red line corresponds to surface elevation minus 20 m.

Figure 7 .
Figure 7. Correlation between observed and simulated groundwater levels in meters above mean sea level (mamsl) at (a) piezometer F7 and (b) piezometer F15.

Figure 7 .
Figure 7. Correlation between observed and simulated groundwater levels in meters above mean sea level (mamsl) at (a) piezometer F7 and (b) piezometer F15.

Figure 8 .
Figure 8. Simulated and observed hydraulic heads in meters above mean sea level based on the steady-state calibration in the Klela basin.

Figure 8 .Figure 9 .
Figure 8. Simulated and observed hydraulic heads in meters above mean sea level based on the steady-state calibration in the Klela basin.Hydrology 2016, 3, 17 11 of 18

Figure 9 .
Figure 9. Calculated hydraulic heads (mamsl) in the steady-state simulation, with borehole locations used to calibrate the model.

Figure 10 .
Figure 10.Simulated and observed heads during the transient calibration in the Klela basin from November 2010-November 2014.

Figure 10 .
Figure 10.Simulated and observed heads during the transient calibration in the Klela basin from November 2010-November 2014.

Figure 11 .
Figure 11.Monthly average water budget in the Klela Basin in Mm³/month.

Figure 12 .
Figure 12.Monthly change in groundwater storage (Mm 3 /month) in the Klela Basin.

Figure 13 .
Figure 13.Long-term groundwater recharge in million cubic meters (MCM) calculated using historical and future climate data from the RCP4.5 scenario based on the Thornthwaite Monthly Water Balance and rainfall (mm) from 1970 to 2050.

Figure 14 .
Figure 14.Long-term monthly groundwater level in meters above mean sea level and groundwater recharge (MCM) in the Klela basin calculated using RCP4.5 climate data from June 2010 to November 2050.

Figure 13 . 18 Figure 13 .
Figure 13.Long-term groundwater recharge in million cubic meters (MCM) calculated using historical and future climate data from the RCP4.5 scenario based on the Thornthwaite Monthly Water Balance and rainfall (mm) from 1970 to 2050.

Figure 14 .
Figure 14.Long-term monthly groundwater level in meters above mean sea level and groundwater recharge (MCM) in the Klela basin calculated using RCP4.5 climate data from June 2010 to November 2050.

Figure 14 .
Figure 14.Long-term monthly groundwater level in meters above mean sea level and groundwater recharge (MCM) in the Klela basin calculated using RCP4.5 climate data from June 2010 to November 2050.

Table 1 .
The values used to calibrate the model (%R: Percentage of rainfall).

Table 3 .
Summary of calibration errors in the transient simulation (MAE: Mean Absolute Error; RMSE: Root Mean Square Error; E: Nash-Sutcliffe coefficient and R 2 : Correlation Coefficient).

Table 3 .
Summary of calibration errors in the transient simulation (MAE: Mean Absolute Error; RMSE: Root Mean Square Error; E: Nash-Sutcliffe coefficient and R 2 : Correlation Coefficient).

Table 4 .
Mean annual groundwater budget in m 3 /year in the Klela basin from November 2010 to November 2014.