Aquifer storage of treated wastewater for subsequent recovery as an important strategy for sustainable water security in Kuwait

This research explores the possibility of creating an underground water reservoir by injecting treated wastewater and evaluates the recovery ef ﬁ ciency and water quality under different injection/recovery scenarios. An injection/recovery well was established in the Dammam Formation aquifer in Kuwait. Six observation wells were also drilled nearby to monitor the water levels and quality. A three-dimensional (3D) numerical hydraulic and transport model was set up using the coupled version of MODFLOW and MT3DMS. Four scenarios of cyclic and continuous injection with different injection and pumping rates stretching over 5,400 days were investigated. Two dispersivity values, low (10 m) and high (200 m), were used in each scenario in order to establish the minimum and maximum recovery limits of useful water that can be expected from the site. The combination of injection rate of 650 m 3 /d and subsequent recovery at 1,000 m 3 /d proved to be the optimum option for the storage of water, resulting in good recovery ef ﬁ ciency and acceptable water quality. Based on these results, it was concluded that aquifer storage is a feasible strategy for mitigating growing water scarcity in the State of Kuwait and other countries in the larger Middle East.

INTRODUCTION destructive, and economically cost-effective method for treating the treated wastewater effluents and naturally purifying them through filtration in soil layers (Akber et al. 2008).
ASR is applicable in a wide range of hydrogeological settings and climatic conditions (Sheng & Zhao 2015). From the mitigation of water scarcity perspective, which was the primary policy objective of this study, ASR provides long-term storage capacity without the problems associated with surface storage reservoirs such as contamination and evaporative losses, with minimum natural environment disturbances, and lower costs (Maliva et al. 2006;Khan et al. 2008). Furthermore, ASR can function as an effective clean technology in removing various forms of pollutants, such as heavy metals and microorganisms (Page et al. 2015). The rationale for the creation of strategic water storage through aquifer recharge is two-folds: (1) maintaining the potentiometric head of groundwater and (2) reuse when the need arises, paving the way to sustainable management of national and regional water resources.
Modeling tools to predict the quantity and quality of water supply sources can facilitate efforts to identify which combination among different scenarios will optimize sustainable water resources supply (Hägg et al. 2021). Previous studies conducted on an adjacent area have not included the influence of a fracture zone in the model on the transmission of the potentiometric effects of pumping and on the movement of the injected water, and did not monitor changes in the water level and the quality of recovered water simultaneously. Objectives of this study are to (i) develop a multi-aquifer flow and transport model for the recharge site, to predict the recovery efficiency under different injection-recovery scenarios and (ii) examine the feasibility of creating a strategic reserve of freshwater for use during high-demand periods at an acceptable water quality level. To realize these objectives, a numerical model was set up using a coupled version of MODFLOW and MT3DMS to investigate the changes in both the water level and the water quality during recharge and subsequent recovery at a pilot recharge site in Kuwait adjacent to the area selected by the Ministry of Electricity and Water (MEW), Kuwait, to create an underground storage of water in the Dammam Formation aquifer through artificial recharge of wells using reverse osmosis (RO) membrane treated municipal wastewater.

MATERIALS AND METHODS
The study involved numerical modeling of injection of reverse osmosis (RO) treated wastewater in the Dammam Formation aquifer through an ASR well located in Kabd, Kuwait, followed by the pumping back of water from the same well. The location of this pilot-scale recharge site is shown in Figure 1. The Visual MODFLOW (VMOD) software was used for this purpose, and system version 2011.1 was used for numeric modeling. This coupled version of the flow and transport simulation system has various options for choosing the numerical engines as per the requirements of an investigation. In the present case, MODFLOW 2005 software, released by the United States Geological Survey in 2005 (Harbaugh 2005) and based on the MODFLOW software originally developed by McDonald & Harbaugh (1984, 1988, was chosen for the flow simulation. This advanced version of MODFLOW is built on a process-based approach and has built-in parameter estimation and local grid refinement capabilities. For the transport simulation, the MT3DMS software engine (Zheng & Wang 1999), the advanced version of the original software MT3D developed in 1990 (Zheng 1990), was selected for the study. MT3DMS has a better equation solver and accommodates non-equilibrium sorption, dual-domain advection-diffusion mass transport, and add-on reaction packages. The combination is currently widely used for in-depth investigations of aquifer flow and transport problems.
A multilayer three-dimensional (3D) numerical hydraulic and solute transport model was developed using the coupled version of MODFLOW 2005 and MT3DMS. The MT3DMS module of VMODFLOW was used for simulating the water quality. The initial calibration of the model was done using the water level data collected during the pumping test. The groundwater flow models are based on the partial differential equation of McDonald & Harbaugh (2003): where x, y, and z are Cartesian coordinates aligned along the major axes of hydraulic conductivity K xx , K yy , K zz ; h is the 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).
In general, S s , K xx , K yy and K zz are functions of space (S s ¼ S s (x,y,z), K xx ¼ K xx (x,y,z), K yy ¼ K yy (x,y,z) and K zz ¼ K zz (x,y,z)), h and w are functions of space and time (h ¼ h(x,y,z,t) and w ¼ w(x,y,z,t)) so that Equation (1) describes groundwater flow in three-dimensions under nonequilibrium conditions in a heterogeneous, anisotropic aquifer. Equation (1) assumes a constant density of water. In the present case, as the difference in salinity, and hence the density, of the injected water and the groundwater is not significant, this assumption remains valid for all practical purposes.
The partial differential equation describing the fate and transport of the contaminants of the species k in 3D, transient groundwater flow systems can be written as follows (McDonald & Harbaugh 2003): where θ ¼ porosity of the subsurface medium (dimensionless); C k ¼ dissolved concentration of species k (ML À3 ); t ¼ time (T); x i,j ¼ distance along the respective Cartesian coordinate axis (L); D ij ¼ hydrodynamic dispersion coefficient tensor (L 2 T À1 ); Water Supply Vol 00 No 0, 4 corrected Proof v i ¼ seepage or linear pore water velocity (LT À1 ); it is related to the specific discharge or Darcy flux through the relationship, v i ¼ q i /θ; q s ¼ volumetric flow rate per unit volume of aquifer representing fluid sources (positive) and sinks (negative) (T À1 ); C k s ¼ concentration of the source or sink flux for species k (ML À3 ); P R n ¼ chemical reaction term (ML À3 T À1 ).
The left-hand side of Equation (2) can be expanded into two terms, i.e., where q 0 s ¼ @u @t is the rate of change in transient groundwater storage (unit, T À1 ). As no or very little chemical reaction is expected under the artificial recharge conditions, the chemical reaction term, ∑R n , may be set to 0.

Model description
A rectangular area of 18 km Â 18 km with the injection and the monitoring wells located approximately at its center was chosen as the modeled area. This has allowed the modeled area boundaries to be at enough distance to remain unaffected by the injection-recovery operation to be conducted at the site. Based on the lithological and geophysical logs of the ASR well and the six monitoring wells and the conventional cores covering the Dammam Formation recovered from well ASR-M-06, the vertical thickness covered by the model was divided into five layers as follows: Layer 1  In the study area, there is no or very little natural recharge from rainfall due to the great depth of the groundwater table (90-100 m), scanty rainfall (average 100-110 mm/year), flat, gently sloping topography and high evaporation rate (more than 2,500 mm/year) (Mukhopadhyay et al. 1996). Natural groundwater movement through the study area under the prevailing hydraulic gradient and transmissivity conditions has been taken into account by the numerical model through the use of appropriate boundary and initial conditions.

Model calibration using pumping test data
The spatial distribution of the injection/recovery well and the monitoring wells is depicted in Figure 2. The water levels in the ASR well and the monitoring wells recorded during the injection and pumping test carried out earlier were used to calibrate the model (Mukhopadhyay 2019). Hydraulic conductivity was used as a calibration parameter because model results were most sensitive to this parameter. The calibration process suggested the presence of a fracture zone with much higher hydraulic conductivity than the aquifer matrix, passing through the ASR well and some of the nearby monitoring wells. The calibrated values of the hydraulic parameters are presented in Table 1. The calibrated values of the model were used for investigation in all four scenarios. The first and the second scenarios assume cyclic injection and recovery for creating the reserve of water. The possibility of creating a long-term storage of water through continuous injection for 1,800 days (approximately five years) was investigated in the third and the fourth scenarios. Two alternative runs with dispersivity values of 10 and 200 m for each scenario were made. Mukhopadhyay (2019) suggested that this is one of the most important parameters that affect the recovery efficiency. The assumption of these two dispersivity values was based on the results of the field experiment on artificial recharge of the Damman Formation aquifer that suggested that the longitudinal dispersivity in this aquifer can vary from 10 to 500 m depending on the transmissivity of the aquifer at the recharge location, recharge duration and the distance traveled by the recharged water (Streetly 1998). Longitudinal dispersivity for specific type of geological medium governs the spreading of contaminants (Mukhopadhyay et al. 1994). Based on previous field experiments, the longitudinal dispersivity value of 10 m was appropriate for the recharge experiment conducted in the Sulaibiya water field, where the hydraulic conditions were very similar to those at the current recharge site (Mukhopadhyay 2019). Therefore, this value  corrected Proof was selected as the most plausible dispersivity value to be used for the scenario runs. As the maximum dispersivity value of 500 m was observed in an area with very high transmissivity (4,000 m 2 /d or higher) by Mukhopadhyay et al. (1994) and the prevailing transmissivity of the Dammam Formation at the selected recharge site was much lower (only about 40 m 2 /d), 200 m was taken as the maximum possible value of the dispersivity at the study location.

RESULTS AND DISCUSSION
The initial and final heads and associated total dissolved solid (TDS) concentrations during the injection and the recovery periods, respectively, for the four scenarios, are presented in Tables 2-5. The volumes of water injected and those recovered with TDS levels of 500 mg/L, 1,000 mg/L and 1,500 mg/L are tabulated in Table 6. The three levels of TDS were selected  corrected Proof based on the following three criteria: (1) the United States Environmental Protection Agency (USEPA 2016) has set a value of 500 mg/L as the maximum level for drinking water (Mukhopadhyay et al. 1994); (2) Berkey (2017) proposed 1,000 mg/L for the limit of acceptable TDS for drinking water; (3) the Canadian guideline has adopted a value of 1,500 mg/L as the maximum limit of TDS in drinking water (SDWF 2017). The evolution of potentiometric heads and TDS contents in the ASR well and the spatial distribution of the head and water quality (TDS) at different time points in the recharge and recovery operation for the four scenarios are presented in Figures 3-5. The injected water moved further down the hydraulic gradient along the fracture zone, causing elongation of the injected water body in the direction of the fracture (Figure 6). The recovery efficiency, defined as the percentage of the injected volume of water that is recovered with a certain specified quality level, for each of the scenarios was computed for the     corrected Proof three TDS content levels ( 500 mg/L, 1,000 mg/L and 1,500 mg/L) of the water recovered by pumping (Table 7). For all scenario runs with a longitudinal dispersivity value of 10 m, the recovery efficiencies for all three TDS levels were 70% or above. For the longitudinal dispersivity value of 200 m, the outcomes may be explained as follows: (1) Scenario 1: no recovery of water with TDS 500 mg/L in any of the injection/recovery cycles; negligible (0.6%) recovery for water with TDS 1,000 mg/L in the third cycle and .100% in the fourth and the fifth cycles; and .100% recovery for water with TDS 1,500 mg/L in all cycles except the first cycle (0% recovery).

) of injected and recovered water of prescribed quality in each cycle of injection and pumping for different scenario runs
(2) Scenario 2: no recovery of water with TDS 500 mg/L in the first four injection and recovery cycles, but .100% recovery in the fifth cycle; no (0%) recovery for water with TDS 1,000 mg/L in the first cycle, but .100% recovery in the rest of the cycles; and .100% recovery for water with (3) Scenario 3: 0% recovery of water with 500 mg/L and 1,000 mg/L TDS and 141% recovery for water with 1,500 mg/L TDS. (4) Scenario 4: The aquifer runs dry on the 864th day after the start of pumping and the water that is recovered has TDS . 500 mg/L resulting in a recovery efficiency of about 74% for TDS levels of 1,000 mg/L and 1,500 mg/L by the time the well runs dry. The above observations suggest that aquifer dispersivity has a big role to play in determining the overall recovery efficiency of an ASR project.
The results further indicate that cyclic injection and recovery of water from an ASR well completed in a brackish aquifer results in the creation of storage with better quality water and higher overall recovery efficiency for a pre-defined quality criterion for the recovered water when compared with those attained under the continuous injection scenario where injection continues over an extended period without any intermittent recovery of the injected water. Though both scenarios 1 and 2 indicate good recovery efficiencies as long as the acceptable TDS can be higher than 500 mg/L (Tables 2 and 3), the water level in the ASR well is expected to rise to a height of 115 m above the mean sea level (amsl), which is only 25-30 m below the ground surface at that location in the case of Scenario 2 (Table 2). In the case that the water level rises higher due to clogging of the injection face by suspended solids or dissolved air in the recharge water, it will create operational problems for sustained injection. Therefore, for the cyclic injection option, Scenario 1 will be more acceptable where the maximum rise of water level in the injection well without clogging will be to a level of about 81 m amsl.
Similarly, for the continuous injection options (Scenarios 3 and 4), the modeling not only suggests better recovery efficiencies for Scenario 3 (Table 6), it also indicates more manageable water levels (maximum rise to 89 m amsl during injection and maximum decline to 18 m below mean sea level (bmsl) during recovery) (Table 3). In the case of Scenario 4 with higher  corrected Proof injection and recovery rates, the water level in the ASR well will rise to about 132 m amsl (only 8-13 m below the ground level) even if there is no clogging and during recovery, after 864 days of pumping, the water level will reach below the base of the aquifer (86 m bmsl by which time only 74% recovery will take place. Thus, even for the continuous injection option, lower injection and recovery rates (Scenario 3) will be preferable. Consequently, to avoid operational problems during the execution of an ASR project, the injection rate and the pumping rate for recovering the injected water should be judiciously selected, taking into consideration the aquifer parameters, the prevailing potentiometric head of the aquifer to be injected, and its distance from the ground surface at the site of the ASR well. A recovery efficiency between 80% and 100% was achieved in both the cyclic and the continuous injection scenarios 1 and 3, respectively. This finding is not significantly different from that of Mukhopadhyay & Al-Otaibi (2002) who reported 70%-80% recovery efficiency for a similar (655 m 3 d À1 ) injection rate. The slight difference may be due to different aquifer properties at the recharge site studied by Mukhopadhyay & Al-Otaibi (2002), the assumed injected water salinity of 500 mg L À1 , and multi-well injection and recovery simulated by these authors. The inclusion of a fracture zone in the model, the presence of which was indicated from the analysis of the pumping test data (Mukhopadhyay 2019) in the current model and assessment of its influence on the transmission of the potentiometric effects of injection and pumping and on the movement of the injected water ( Figure 6) are novelties of this study.
The comparison of the rise and fall in water levels and the recovery efficiencies achieved in scenarios 1 and 2 with those achieved by scenarios 3 and 4 is also consistent with the finding that the cyclic injection and recovery would be much more efficient in controlling the rise and fall in the potentiometric head during injection and recovery and in attaining better recovery efficiency when compared with those parameters obtained for the continuous injection option. The improvement in the recovery efficiency with the number of cycles in the scenario runs 1 and 2 (Table 6) is also in line with the similar observation made by Pyne (1998). The reduction in the recovery efficiency as the longitudinal dispersivity value increases from 10 to 200 m in all scenarios agrees with the findings reported by Shulze-Makuch (2005), Anderson & Lowry (2004) and Merritt (1985).
As treated wastewater travels from and to the ASR well during injection and recovery, residence time plays an important role in determining the suitability of the water recovered for different intended uses. Though the effects of the length of the storage period between the end of injection and beginning of recovery have not been investigated in the present study, it may be noted here that other researchers like Ross (2018), Al-Damkhi et al. (2009) and Anderson & Lowry (2004) have observed a decrease in recovery efficiency as the length of the storage period increases. Overall, our findings suggest that water extracted during the recovery period may not necessarily be injected water alone, but a blend of injected and ambient water in varying proportions, depending on the duration of storage, the duration and rate of injection and extraction, and other factors. This observation is in line with that made by other workers in this area from different parts of the world (Anderson & Lowry 2004;Eini et al. 2020).
It may be noted here that Mukhopadhyay & Al-Otaibi (2002) carried out a study for a possible recharge site 10-12 km further to the east of the current study area, which may appear to overlap the current study. In that study however, the model used included 37 ASR wells instead of one. These wells were arranged in a hexagonal pattern within a 10 km Â 10 km area. The transmissivity of the Dammam Formation aquifer within the area varied from 40 to 200 m 2 /d. All the wells were injected and pumped simultaneously. In this study, however, for practical reasons, MEW decided to select a different area for the creation of a strategic reserve of water through artificial recharge, based on various operational and hydrogeological considerations. As there were apprehensions about the applicability of the results of this study to the area finally selected for artificial recharge due to the differences in hydrogeological conditions and groundwater quality, it was decided to carry out a fresh modeling study at the selected site that was adjacent to the area finally selected by MEW for artificial recharge of the Dammam Formation aquifer. A pilot-scale one-well recharge-recovery operation, consisting of four cycles of recharge and recovery, extended over a year at the same site was planned (Figure 1). Therefore, the outcomes of this study are of more practical relevance to the implementation of artificial recharge in Kuwait.
Moving forward, while the safe reuse of recycled water is technically feasible and that helps meet the shortage of fresh water in all parts of the world, acceptability (safety perception) for society of wastewater reuse for various purposes (potable or non-potable) is a critical factor for the successful implementation of wastewater reuse to meet growing demands. For religious reasons, there has been an intense debate in the Muslim world (where water scarcity is most acute) on the validity of the direct use of recycled water for certain social activities such as ablution (Al-Damkhi et al. 2009). Ablution is a requirement for cleaning certain parts of the body before performing prayers, a process that consumes significantly large volumes of water.
A major advantage of water extracted from ASR is that the practice of injecting treated wastewater into the ground makes it appear as groundwater in people's minds and consequently removes religious and/or psychological obstacles, leading to an enhanced public acceptability of the reclaimed wastewater. On the other hand, from a practical perspective, clogging of the injection wells and/or aquifer basins due to accumulation of suspended solids can be a major challenge (Mohamed et al. 2016). This can be a particularly serious problem if wastewater treatment is inadequate. One final point to be communicated is that there is no single technology that will solve all of our challenges and more integrated research work remains to be done regarding ASR.
A possible limitation of the study that may be argued is that over the simulation period of five years or more, changes in groundwater movement and quality due to natural recharge of the aquifers and extraction of groundwater in the adjacent areas have not been considered. In this context, it might be mentioned that it was shown by Mukhopadhyay et al. (1996) that except for a few topographic depressions, mostly located in the northern parts of the country, natural recharge from rainfall was negligible in Kuwait. Its exclusion from the model should, therefore, be acceptable. Due to the limited time available, it was not possible to collect all relevant data on the production history and the future projections for the production of the brackish water fields contiguous to the study area from MEW and the Kuwait Oil Company (KOC), the organizations operating these fields. It was, however, apparent from the latest brackish water production statistics, published by MEW (MEW 2019), that brackish water production from these fields had been drastically reduced in recent years, and the effects of the operation of these fields on the artificial aquifer recharge at the study location should not be significant. The initial and boundary conditions of the model and the model extent should have taken care of the lateral recharge and discharge across the model boundaries and the changes in potentiometry, water quality and water movement in the vicinity of the recharge site due to the injection and pumping operations. Thus, major conclusions derived from the study remain valid. A more detailed modeling of the planned development of a strategic reserve of water through artificial recharge needs to be undertaken before the implementation of the project in the future that should take into account the prevailing status of aquifer exploitation in the vicinity of the project site.

CONCLUSIONS
This study explored the potential of the ASR technique for water security in Kuwait. Numerical 3-D flow and transport modeling of the recharge of the Dammam Formation aquifer through an ASR well in the Kabd area in Kuwait has been carried out to investigate the feasibility of implementing this technique in Kuwait. The results suggest that cyclic injection and recovery should give better recovery efficiency than that obtained from uninterrupted injection over an extended period followed by recovery. However, as long as the dispersivity of the aquifer is reasonably low (about 10 m), the injecting of freshwater (TDS 70 mg/L) at the rate of 650 m 3 /d and subsequent recovery at the rate of 1,000 m 3 /d will be ideal for both the cyclic and the continuous injection options with high recovery efficiency and a reasonable rise in the potentiometric head during the injection.

DATA AVAILABILITY STATEMENT
All relevant data are included in the paper or its Supplementary Information.