Observations and Prediction of Recovered Quality of Desalinated Seawater in the Strategic ASR Project in Liwa, Abu Dhabi

To be able to overcome water shortages, Abu Dhabi Emirate started an Aquifer Storage and Recovery (ASR) project with desalinated seawater (DSW) as source water near Liwa. It is the largest DSW-ASR project in the world (stored volume ~10 Mm3/year), and should recover potable water for direct use. DSW is infiltrated into a desert dune sand aquifer using “sand-covered gravel-bed” recharge basins. In this study, we evaluate the hydrogeological and hydrogeochemical stratification of the (sub)oxic target aquifer, and water quality changes of DSW during trial infiltration runs. We predict water quality changes of DSW after 824 d of infiltration, during 90 d of intensive recovery (67% recovered) without storage (scenario A), as well as after 10 years of storage (scenario B, with significant bubble drift). Monitoring of preceding trials revealed a lack of redox reactions; little carbonate dissolution and Ca/Na exchange; much SiO2 dissolution; a strong mobilization of natural AsO4, B, Ba, F, CrO4, Mo, Sr and V from the (sub)oxic aquifer; and immobilization of PO4, Al, Cu, Fe and Ni from DSW. The Easy-Leacher model was applied in forward and reverse mode including lateral bubble drift, to predict water quality of the recovered water. We show that hydrogeochemical modeling of a complex ASR-system can be relatively easy and straightforward, if aquifer reactivity is low and redox reactions can be ignored. The pilot observations and modeling results demonstrate that in scenario A recovered water quality still complies with Abu Dhabi’s drinking water standards (even up to 85% recovery). For scenario B, however, the recovery efficiency declines to 60% after which various drinking water standards are exceeded, especially the one for chromium.


Introduction
Water scarcity has driven many countries in arid zones, such as the Middle East and Abu Dhabi in particular, to desalinate large volumes of seawater for fresh water supply [1]. Episodic problems with seawater quality due to, e.g., harmful algae blooms [2,3] and oil spills, energy supply and fear of The ASR plant is thus not a normal ASR system, which is exclusively composed of wells performing two tasks: infiltrate and recover [20].

The Pilot
The first operation run of the pilot started on 1 October 2003 and lasted until December 2004. DSW was infiltrated via a series of 5 ASR wells (not discussed further) and a recharge basin recovery scheme consisting of a covered gravel-body with reverse drains in it, 4 recovery wells and tens of observation wells around. DSW was infiltrated via the basin for 250 d at ~250 m 3 /h, and after 48 d of storage part of the stored volume was recovered in 70 d at ~250 m 3 /h. Intensive monitoring of water quality yielded important insights in the ambient hydrochemical stratification and water quality changes during infiltration, stand-still and recovery [8].
A second infiltration run via the basin took place in 2008, but water quality monitoring was not undertaken and DSW was not pumped out. This offered the possibility to sample, in August 2014, 6-year-old DSW from the aquifer.   The ASR plant is thus not a normal ASR system, which is exclusively composed of wells performing two tasks: infiltrate and recover [20].

The Pilot
The first operation run of the pilot started on 1 October 2003 and lasted until December 2004. DSW was infiltrated via a series of 5 ASR wells (not discussed further) and a recharge basin recovery scheme consisting of a covered gravel-body with reverse drains in it, 4 recovery wells and tens of observation wells around. DSW was infiltrated via the basin for 250 d at ~250 m 3 /h, and after 48 d of storage part of the stored volume was recovered in 70 d at ~250 m 3 /h. Intensive monitoring of water quality yielded important insights in the ambient hydrochemical stratification and water quality changes during infiltration, stand-still and recovery [8].
A second infiltration run via the basin took place in 2008, but water quality monitoring was not undertaken and DSW was not pumped out. This offered the possibility to sample, in August 2014, 6-year-old DSW from the aquifer.

Quantitative Description of the Break-through Curve
The first infiltration run of the pilot yielded valuable insight into the break-through curve (BTC) of nearly all main constituents and trace elements. These observed BTCs are characterized by 3 parameters: pore volume, retardation or leach factor and (semi)permanent concentration change ( Figure 4).
The dimensionless parameter called "pore volume" (PV) forms the time axis of water quality observations and model predictions: where tINF = total infiltration period [d]; and t50 = the observed 50% breakthrough time of conservative tracer or the calculated travel time via Equation (4). One PV means that the whole aquifer, from infiltration point to the monitoring or recovery well, has been flushed with the infiltration water exactly one time. Retardation factors R or leach factors L can then be simply deduced from concentration plots against PVs (Figure 4).
Sorbing and oxidizing solutes as well as desorbing and dissolving solutes are retarded during aquifer passage compared with conservative solutes. In the latter case, raised concentrations drop to the influent level long after passage of the conservative chloride front. These delays are quantified for solute i by, respectively, the well-known retardation factor Ri, and the less well known leach factor Li [25]:

Quantitative Description of the Break-through Curve
The first infiltration run of the pilot yielded valuable insight into the break-through curve (BTC) of nearly all main constituents and trace elements. These observed BTCs are characterized by 3 parameters: pore volume, retardation or leach factor and (semi)permanent concentration change ( Figure 4).
The dimensionless parameter called "pore volume" (PV) forms the time axis of water quality observations and model predictions: where t INF = total infiltration period [d]; and t 50 = the observed 50% breakthrough time of conservative tracer or the calculated travel time via Equation (4). One PV means that the whole aquifer, from infiltration point to the monitoring or recovery well, has been flushed with the infiltration water exactly one time. Retardation factors R or leach factors L can then be simply deduced from concentration plots against PVs (Figure 4).
Sorbing and oxidizing solutes as well as desorbing and dissolving solutes are retarded during aquifer passage compared with conservative solutes. In the latter case, raised concentrations drop to the influent level long after passage of the conservative chloride front. These delays are quantified for solute i by, respectively, the well-known retardation factor R i , and the less well known leach factor L i [25]:  [-]; and r P = reaction coefficient related to (prod). For practical reasons, t i is set at ≥90% not at 100%. If for some reason, the BTC shows a partial breakthrough due to a prolonged phase of continued partial immobilization or mobilization, so the additional parameter ∆C is needed to describe the BTC (Figure 4). Equations (2) and (3) hold for a stationary retardation or leaching process, with a homogeneously distributed reactive phase in the aquifer. Of course, (reac) or (prod) should have no other sinks or sources, unless these can be properly accounted for.  [-]; and rP = reaction coefficient related to (prod). For practical reasons, ti is set at ≥90% not at 100%. If for some reason, the BTC shows a partial breakthrough due to a prolonged phase of continued partial immobilization or mobilization, so the additional parameter ΔC is needed to describe the BTC (Figure 4). Equations (2) and (3) hold for a stationary retardation or leaching process, with a homogeneously distributed reactive phase in the aquifer. Of course, (reac) or (prod) should have no other sinks or sources, unless these can be properly accounted for.

Lithological and Geochemical Stratification
The aquifer system on well fields A, B and C could be schematized into a succession of 14 (sub) horizontal layers in between ground surface and 40 m below sea level [23], based on all drilling logs (>372), geophysical logs (including the eccentered wireline NMR logging of permeability and porosity), pumping tests, infiltration pilot tests and fluid flow logging, as presented by [9][10][11][12]. For modeling purposes, the 14 layers were aggregated into 6 main aquifer layers (a-f in Figure 5), of which the very low permeability aquitard f needs less consideration. , which displaces fluid A (concentration C A ), in terms of travel time (t 50 ) and 90% break-through or leach times (t i ; t i 3, t i 4 = ti for curve 3 and 4, respectively), the dimensionless parameter PV (pore volumes flushed) and the (semi)permanent concentration change (∆C), the maximum concentration (C max ) or minimum concentration (C min ) after full break-through or leaching. Line 1 = conservative tracer (R = L = 1) without dispersion; Curve 2 = as 1, however, with dispersion; Curve 3 = compound retarded by sorption or leaching (R = L = 4); Curve 4 = compound retarded by ad-or desorption (R = L = 5), with continued removal c.q. addition. (A) (C B > C A ); and (B) (C B < C A ).

Lithological and Geochemical Stratification
The aquifer system on well fields A, B and C could be schematized into a succession of 14 (sub) horizontal layers in between ground surface and 40 m below sea level [23], based on all drilling logs (>372), geophysical logs (including the eccentered wireline NMR logging of permeability and porosity), pumping tests, infiltration pilot tests and fluid flow logging, as presented by [9][10][11][12]. For modeling purposes, the 14 layers were aggregated into 6 main aquifer layers (a-f in Figure 5), of which the very low permeability aquitard f needs less consideration. The mean geochemical composition of layers a-f was calculated from the geochemical data of 4 deep core drillings, one in each well field and one in the middle of well fields A, B and C. The data were derived from [12], containing the following information: Sample description with photographs, petrographic analysis and mineralogical counting by petrological photomicroscopic examination of a thin section impregnated with fluorescent resin, XRD analysis, porosimetric analysis, chemical analysis of main constituents (including LOI and Acid Solubility), chemical analysis of heavy metals (probably in nitric acid; no details given), and particle size distribution (by sieving).
The digital data of all 41 samples were used to calculate mean values for aquifer layers a-f, and to quantify the content of selected minerals by petrochemical calculations [23]. Petrochemical calculations were needed to correct specific data for water losses (loss on ignition), to calculate the cation exchange capacity (CEC) and to derive the mineral content from data on elements that are present in more than one mineral.

Hydrochemical Analyses in Period 2003-2013
Samples were taken from practically all recovery and monitoring wells, in both the pilot area and well fields A, B and C. The analyses include data on gases (O2, residual Cl2), turbidity, color, temperature, pH, EC, ORP (Oxidation Reduction Potential, Eh), main anions (Cl, SO4, S, HCO3, CO3, NO3, NO2, PO4, and F), main cations (Na, K, Ca, Mg, NH4, Fe, and Mn), SiO2, TOC (Total Organic Carbon), CN, selected trace elements (Ag, Al, As, B, Ba, Br, Cd, Cr, Cu, Hg, Mo, Ni, Pb, Sb, Se, Sr, V and Zn), and the stable isotopes 2 H and 18 O. Most samples were filtered in the field, and all concentrations (TOC excluded) refer to total dissolved concentrations, thus without (further) speciation. Microbial parameters and organic micropollutants were analyzed but showed negligible concentration levels.
Samples of the ambient groundwater were taken in the pilot area in 2003, and within and around well fields A, B and C in 2011-2013. Analytical data of DSW samples prior to infiltration, during and after aquifer passage were exclusively available from the pilot. All data (obtained from [9][10][11][12]   The mean geochemical composition of layers a-f was calculated from the geochemical data of 4 deep core drillings, one in each well field and one in the middle of well fields A, B and C. The data were derived from [12], containing the following information: Sample description with photographs, petrographic analysis and mineralogical counting by petrological photomicroscopic examination of a thin section impregnated with fluorescent resin, XRD analysis, porosimetric analysis, chemical analysis of main constituents (including LOI and Acid Solubility), chemical analysis of heavy metals (probably in nitric acid; no details given), and particle size distribution (by sieving).
The digital data of all 41 samples were used to calculate mean values for aquifer layers a-f, and to quantify the content of selected minerals by petrochemical calculations [23]. Petrochemical calculations were needed to correct specific data for water losses (loss on ignition), to calculate the cation exchange capacity (CEC) and to derive the mineral content from data on elements that are present in more than one mineral.

Hydrochemical Analyses in Period 2003-2013
Samples were taken from practically all recovery and monitoring wells, in both the pilot area and well fields A, B and C. The analyses include data on gases (O 2 , residual Cl 2 ), turbidity, color, temperature, pH, EC, ORP (Oxidation Reduction Potential, Eh), main anions (Cl, SO 4 , S, HCO 3 , CO 3 , NO 3 , NO 2 , PO 4 , and F), main cations (Na, K, Ca, Mg, NH 4 , Fe, and Mn), SiO 2 , TOC (Total Organic Carbon), CN, selected trace elements (Ag, Al, As, B, Ba, Br, Cd, Cr, Cu, Hg, Mo, Ni, Pb, Sb, Se, Sr, V and Zn), and the stable isotopes 2 H and 18 O. Most samples were filtered in the field, and all concentrations (TOC excluded) refer to total dissolved concentrations, thus without (further) speciation. Microbial parameters and organic micropollutants were analyzed but showed negligible concentration levels.
Samples of the ambient groundwater were taken in the pilot area in 2003, and within and around well fields A, B and C in 2011-2013. Analytical data of DSW samples prior to infiltration, during and after aquifer passage were exclusively available from the pilot. All data (obtained from [9][10][11][12] and from data files supplied by Dr. G. Koziorowski (GTZ International Services) were controlled, elaborated and stored with Hydrogeochemcal.xlsx [26].

Monitoring Campaign in August 2014
In the period of 3-7 August 2014 a sampling campaign was conducted to take 31 samples from divergent observation wells that had been sampled earlier and from DSW, in order to: (i) check potential bias in the existing hydrochemical data set of well fields A, B and C; and (ii) obtain data on chromium speciation: Cr(III) versus Cr(VI).
Important aspects that were tackled with great care, are: sufficient well purging based on purged volume and stable field parameters, sampling without applying vacuum and excluding atmospheric exposition, reducing exposure to sunlight and wind, flow regulation of the pump, field measurements (EC, pH, temp, O 2 ), filtration of water over a 0.45 µm membrane, dedicated sample preservation for specific parameter groups, cooling, nearly daily shipment to the Netherlands, and rapid analysis in the certified Vitens Laboratory (Leeuwarden). The quality of the analysis was validated using HGC 2.0 to exclude potential impact of errors.

Predictions by EL Modeling
Two models were constructed, an Excel based Easy-Leacher (EL) model [25] and a PHREEQC-2 [27] flowtube model. PHREEQC-2 was applied to model more in detail the behavior of chromate and arsenate along a small number of flowlines during infiltration phase. PHREEQC-2 and EL produced nearly identical results for chromate and arsenate behavior during infiltration, justifying the application of the simpler EL model. Further details about the PHREEQC-2 modeling and its results are given by [24] and not considered here further.
EL simplifies 3D groundwater flow into a 2D set of maximum 50 flow tubes through a maximum of 10 horizontal aquifer layers. The travel times are either derived from a hydrological model, or calculated analytically. Chemical transport is calculated on the basis of pore volumes (dimensionless time scale), retardation and leach factors (superimposed on the pore volumes) based on either mass balances or field observations, CaCO 3 equilibrium (if relevant), redox reactions (if relevant), and expert rules on among others redox reaction kinetics. EL was given the task to do the all-round water quality modeling during all ASR phases (infiltration, storage with bubble drift, and recovery), and to combine the output of a relatively high number of flowlines into a mixed output as generated by a well field.
In EL, the whole ASR system was schematized by one representative recharge basin (the average of basins A-C), 5 aquifer layers (of which the upper layers a, b and c are most important) and 7 flowlines within each aquifer layer departing from the basin towards one recovery well in each of the 7 well rings at 75-225 m radius ( Figure 5).
The expansion of the DSW bubble in each aquifer layer and the travel times along each flowline were calculated by assuming first vertical flow down to each aquifer layer and then horizontal radial flow, so that:  [8]. This simplification creates some distortion during the first 30 d, but these are of minor importance on the long run of 824 d of infiltration. With Equation (4), the travel time was calculated from recharge basin to its 7 surrounding rings of recovery wells at the distances specified in Figure 6. It is deduced that in theory all wells, also those in the outermost ring will pump DSW after 824 d of infiltration. In layers D and E, probably little DSW will be present. During 10 years of storage, the infiltrated DSW bubble is predicted to move laterally down the regional hydraulic gradient, with the following velocity (vB,N), assuming an equal gradient in each layer: where ΔH/ΔX = mean regional hydraulic gradient in the aquifer at well fields A, B and C during storage phase [m/m]. Vertical bubble drift by buoyancy has been ignored in accordance with FEFLOW model predictions [9]. The water quality evolution during injection phase was calculated for each flowline where it crossed its destination well (node point), and also for the "fictive", mixed sample taken from all 35 node points, in proportion to: (i) the preset pumping regime (the inner wells pumping more than the outer wells); and (ii) the transmissivity of each main contributing aquifer layer (a-c). This mixed sample thus represents the output from the whole well field, when pumping out a negligible amount of water without disturbing the continuous DSW bubble expansion.
EL in forward mode assumes the following: (1) input quality (DSW) is constant; (2) retardation factors, leach factors and the (semi)permanent concentration changes are derived from observations during the pilot; (3) redox reactions are absent as observed in the pilot; (4) reactive minerals such as calcite, dolomite and silicates are not depleted; and (5) reaction kinetics do not play an important role.
The recovery phase was modeled by moving backward in the time series that displays the quality evolution of this "imaginary", mixed (averaged) water sample. Contrary to the infiltration phase, this "imaginary", mixed sample now becomes the true output of the well field during recovery phase, showing after some time an increasing instead of decreasing contribution of native groundwater. This is in harmony with theory and the predictions by [9,11].
As during the 90 d of recovery 67% of the infiltrated water volume will be pumped out, the way back in the expansion time series needs to be as long as 0.67 × 824 = 552 d (= ΔtBACK-1). This way, water quality at the start of recovery, without storage phase, will be the water that flushed each of the 7 well rings prior to pumping (on day 824), and this water only needs to be mixed in proportion to the pumping rates of each well ring ( Figure 6). At the end of pumping (after 90 d) we obtain about the same water as predicted to surround the wells on day 824 − 552 = 272. Water compositions in between can be calculated by just following the predicted water quality from 824 to 272 d back in time. In order to plot this water quality evolution during 90 d of recovery, the expanded time scale (from 90 to 552 d) needs to be compressed back again to 90 d, by multiplying it with 90/552, and to mirror it from backward into forward mode ( Figure 7). During 10 years of storage, the infiltrated DSW bubble is predicted to move laterally down the regional hydraulic gradient, with the following velocity (v B,N ), assuming an equal gradient in each layer: where ∆H/∆X = mean regional hydraulic gradient in the aquifer at well fields A, B and C during storage phase [m/m]. Vertical bubble drift by buoyancy has been ignored in accordance with FEFLOW model predictions [9]. The water quality evolution during injection phase was calculated for each flowline where it crossed its destination well (node point), and also for the "fictive", mixed sample taken from all 35 node points, in proportion to: (i) the preset pumping regime (the inner wells pumping more than the outer wells); and (ii) the transmissivity of each main contributing aquifer layer (a-c). This mixed sample thus represents the output from the whole well field, when pumping out a negligible amount of water without disturbing the continuous DSW bubble expansion.
EL in forward mode assumes the following: (1) input quality (DSW) is constant; (2) retardation factors, leach factors and the (semi)permanent concentration changes are derived from observations during the pilot; (3) redox reactions are absent as observed in the pilot; (4) reactive minerals such as calcite, dolomite and silicates are not depleted; and (5) reaction kinetics do not play an important role.
The recovery phase was modeled by moving backward in the time series that displays the quality evolution of this "imaginary", mixed (averaged) water sample. Contrary to the infiltration phase, this "imaginary", mixed sample now becomes the true output of the well field during recovery phase, showing after some time an increasing instead of decreasing contribution of native groundwater. This is in harmony with theory and the predictions by [9,11].
As during the 90 d of recovery 67% of the infiltrated water volume will be pumped out, the way back in the expansion time series needs to be as long as 0.67 × 824 = 552 d (= ∆t BACK-1 ). This way, water quality at the start of recovery, without storage phase, will be the water that flushed each of the 7 well rings prior to pumping (on day 824), and this water only needs to be mixed in proportion to the pumping rates of each well ring (Figure 6). At the end of pumping (after 90 d) we obtain about the same water as predicted to surround the wells on day 824 − 552 = 272. Water compositions in Water 2017, 9, 177 9 of 25 between can be calculated by just following the predicted water quality from 824 to 272 d back in time. In order to plot this water quality evolution during 90 d of recovery, the expanded time scale (from 90 to 552 d) needs to be compressed back again to 90 d, by multiplying it with 90/552, and to mirror it from backward into forward mode (Figure 7).
In the case of a 10-year storage phase, lateral bubble drift can be taken into account by first extending the backward period of 552 d (∆t BACK-1 ) with the period (=∆t BACK-2 ) that would be needed to get the retrograded bubble face back on its position at 552 d without bubble drift, and then resetting the time scale to 90 d by multiplying it with 90/(552 + ∆t BACK-2 ), and mirror it from backward into forward mode (Figures 7 and 8).
The calculation of ∆t BACK-2 is then as indicated in the realistic example elaborated in Figure 8. In order to simplify the calculations, the weighted average value of ∆t BACK-2 is taken, being 174 d in the example of Figure 6. Addition of ∆t BACK-2 = 552 yields a total set-back time of 726 d. This means that the quality of the water recovered after 90 d of pumping, is to be looked up in the quality output list on day 824 − 726 = 98, in case of bubble drift during 10 years of storage.
How this example with resulting time shifts works out in the %DSW and TDS concentration of the water recovered, is shown in Figure 7. The underlying calculations were performed in EL, and match the predictions by [9] quite well.
EL in backward mode assumes that: (a) the forward evolution is reversed at 6 times higher speed; (b) the 6 times higher recovery rate does not provoke serious upconings (as shown by 3D modeling results [9]; and (c) during storage and recovery no further reactions with the aquifer are taking place.
Fluxes in the different aquifer layers, from the recharge basin towards the recovery wells and beyond them, were set equal to their contribution to the aquifer's transmissivity.
Model calibration was done on: (1) available data from the pilot study in 2003-2004 when DSW was infiltrated via a recharge basin; and (2) groundwater quality as observed in the same pilot area in August 2014, after about 6 years of storage in the local aquifer system (since a second recharge run in 2008).
Water 2017, 9, 177 9 of 25 In the case of a 10-year storage phase, lateral bubble drift can be taken into account by first extending the backward period of 552 d (ΔtBACK-1) with the period (=ΔtBACK-2) that would be needed to get the retrograded bubble face back on its position at 552 d without bubble drift, and then resetting the time scale to 90 d by multiplying it with 90/(552 + ΔtBACK-2), and mirror it from backward into forward mode (Figures 7 and 8).
The calculation of ΔtBACK-2 is then as indicated in the realistic example elaborated in Figure 8. In order to simplify the calculations, the weighted average value of ΔtBACK-2 is taken, being 174 d in the example of Figure 6. Addition of ΔtBACK-2 = 552 yields a total set-back time of 726 d. This means that the quality of the water recovered after 90 d of pumping, is to be looked up in the quality output list on day 824 − 726 = 98, in case of bubble drift during 10 years of storage.
How this example with resulting time shifts works out in the %DSW and TDS concentration of the water recovered, is shown in Figure 7. The underlying calculations were performed in EL, and match the predictions by [9] quite well.
EL in backward mode assumes that: (a) the forward evolution is reversed at 6 times higher speed; (b) the 6 times higher recovery rate does not provoke serious upconings (as shown by 3D modeling results [9]; and (c) during storage and recovery no further reactions with the aquifer are taking place.
Fluxes in the different aquifer layers, from the recharge basin towards the recovery wells and beyond them, were set equal to their contribution to the aquifer's transmissivity.
Model calibration was done on: (1)

Hydrogeology
The hydrogeological schematization is presented in Figure 5. It shows an upper aquifer with transmissivity 768 m 2 /d, formed by layers a-c, which consist of reddish to yellowish brown, eolian sand(stone) of Quaternary age. It is unsaturated in its upper 30 m, and consists there mostly of uncemented or weakly cemented sand. From the groundwater table downwards to the aquifer base at 62 m ASL, an alternation is observed of semi-consolidated and weakly to moderately cemented sands and sandstone.
Aquitard 1 and Aquifer 2, between 62 m ASL and 14 m BSL, likely correspond to the Pleistocene Medinat Zayed formation, consisting of mainly eolian and some fluvial and lacustrine deposits [11]. It is composed of yellowish brown to red brown to gray calcarenaceous sandstone with intercalations of siltstone, mudstone, marl and thin sand lenses.
Aquitard 2, at the base of the considered aquifer complex, is composed of a Neogene formation [11], consisting of light gray to white calcarenite, limestone and dolomite with interbeds of white to pinkish gypsum, marl, and chalk.

Geochemistry
The mean geochemical stratification (Table 1) follows more or less the hydrogeological subdivision ( Figure 5). The interpreted data reveal the presence of the following reactive phases in the aquifer system, with their main potential interaction with water within brackets: Bulk Organic Material (denitrification, sorption), clay minerals (sorption), iron (hydr)oxide coatings of sand grains (source of Fe and oxyanions like chromate, vanadate and arsenate), calcite and dolomite (or dolomitic limestone; source of Ca, Mg, Sr and HCO3), feldspars (source of Na, Ca and SiO2), gypsum (source of Ca, Sr and SO4), pyroxene (source of oxyanions like chromate, vanadate and arsenate), and (fluoro)apatite (source of Ca, PO4 and F). Gypsum and dolomite were concentrated in aquitard 2 (at the base), feldspar and iron (hydr)oxides in aquifer 1 (the top).
The eolian sand is typically coated with iron (hydr)oxide, and more or less cemented mainly by calcite. The iron hydroxide coatings are considered to be the main (but genetically the secondary) source of oxyanions like chromate [Cr(VI)], vanadate, molybdate, selenate, arsenate and phosphate. These negatively charged ions are (chemi)sorbed to these positively charged coatings, and may desorb from it under specific conditions.

Hydrogeology
The hydrogeological schematization is presented in Figure 5. It shows an upper aquifer with transmissivity 768 m 2 /d, formed by layers a-c, which consist of reddish to yellowish brown, eolian sand(stone) of Quaternary age. It is unsaturated in its upper 30 m, and consists there mostly of uncemented or weakly cemented sand. From the groundwater table downwards to the aquifer base at 62 m ASL, an alternation is observed of semi-consolidated and weakly to moderately cemented sands and sandstone.
Aquitard 1 and Aquifer 2, between 62 m ASL and 14 m BSL, likely correspond to the Pleistocene Medinat Zayed formation, consisting of mainly eolian and some fluvial and lacustrine deposits [11]. It is composed of yellowish brown to red brown to gray calcarenaceous sandstone with intercalations of siltstone, mudstone, marl and thin sand lenses.
Aquitard 2, at the base of the considered aquifer complex, is composed of a Neogene formation [11], consisting of light gray to white calcarenite, limestone and dolomite with interbeds of white to pinkish gypsum, marl, and chalk.

Geochemistry
The mean geochemical stratification (Table 1) follows more or less the hydrogeological subdivision ( Figure 5). The interpreted data reveal the presence of the following reactive phases in the aquifer system, with their main potential interaction with water within brackets: Bulk Organic Material (denitrification, sorption), clay minerals (sorption), iron (hydr)oxide coatings of sand grains (source of The eolian sand is typically coated with iron (hydr)oxide, and more or less cemented mainly by calcite. The iron hydroxide coatings are considered to be the main (but genetically the secondary) source of oxyanions like chromate [Cr(VI)], vanadate, molybdate, selenate, arsenate and phosphate.
These negatively charged ions are (chemi)sorbed to these positively charged coatings, and may desorb from it under specific conditions. Another mineral frequently mentioned by [12], in concentrations of ≤1%, is pyroxene. Pyroxenes are of particular interest due to their potential contribution to Cr. Their composition is as follows: XY(Si,Al) 2 O 6 with X representing Ca, Na, Fe(II), Mg and more rarely Zn, Mn and lithium, and Y representing Cr, Al, Fe(III), Mg, Mn, Sc, Ti and V. There is indeed a good positive correlation between Cr, Al, Fe, Ni, Ti, V and Zn [23], indicating that pyroxenes could be a significant primary source of Cr. No mention was made of the occurrence of olivine, another mineral typical of ophiolites, which [28] suspect as being the primary source of Cr in Abu Dhabi's groundwater. This mineral has a high weathering thus low resistance potential, and has possibly therefore not been found.

Chemistry of Native Groundwater
The strong vertical zonation of water qualities in the aquifer system is mainly linked to a rising salinity (TDS) with depth ( Figure 5). The latter is dictated by the presence of shallow "fresh" groundwater on top of deep seated brackish to saline paleowater [11]. The "fresh" groundwater, with a low Cl/Br ratio and high stable isotope content, is considered to be derived from local precipitation during relatively warm and humid climatic conditions with high evaporation losses, without significant contributions from evaporite dissolution [11,29]. The brackish to saline paleowater, with a high Cl/Br ratio and low stable isotope content, is considered to be derived from rainwater that dissolved evaporitic rock during relatively cold climatic conditions with less evaporation losses [11,29].
The following parameters follow the salinity (TDS) increase with depth: EC, Cl, SO 4 , NO 3 , Na, Fe, Mn, Al, B, Br, Cr, Cr(VI), Mo, Ni, Sr and Cl/Br ratio, whereas HCO 3 , Ba and Ca/Mg ratio decline with depth ( Table 2). All groundwater is (sub)oxic (containing little O 2 and much NO 3 ), more or less in equilibrium with calcite and dolomite (as expected on the basis of geochemistry), slightly undersaturated with respect to barite, and (strongly) undersaturated with respect to fluorite and gypsum (base of layer e excluded where near equilibrium with gypsum).
The current mean-annual rainfall of 40-60 mm/year [30] would be sufficient to sustain natural recharge by local precipitation at a rate of~7.6 [11] or 11 mm/year [8]. Tritium data on samples taken close to the groundwater table, stable isotope data, fluoro-chlorinated hydrocarbons and reaction patterns of the groundwater table to recharge events confirm that recharge of groundwater is actually taking place [8].
Between well fields A, B and C some differences exist [23]. Well field B is the most saline at all depth levels. Well field C is the least saline and most homogeneous of all. Well field A forms an intermediate, and is therefore taken as a representative for the whole plant.  Table 2. Simplified chemical stratification of the native groundwater, averaged over both pilot and ABC plant areas. Aquifer 1 = layers a-c; Aquitard 1 = layer d; Aquifer 2 = layer e. Abu Dhabi drinking water standards according to [13], WHO guidelines according to [31], both with red numbers for parameters considered at risk during ASR recovery. All samples showed: Ag < 1, Be < 0.1, Pb < 5 µg/L.

Results of Sampling Campaign in August 2014
The sampling campaign of 3-7 August 7 2014 revealed [23] that there was only limited bias in the 2003-2013 hydrochemical data set of well fields A, B and C, and that practically all (>95%) of dissolved Cr consists indeed of Cr(VI).
The following bias was deduced for the large data set 2003-2013: O 2 concentrations for all monitoring wells were too high due to oxygenation, and suspended fines raised the concentrations of PO 4 , Fe, Mn, NH 4 , Al, and possibly also Zn in many samples. In addition, the minimum detection limits (MDLs) could be lowered for Cd, Hg, Ni, and Sb. All hydrochemical data presented here have been corrected for the demonstrated bias, where possible.

Results of Pilot
The available monitoring data at the pilot in 2003-2004 and 2014 (some of which are in Table 3) reveal among others the following quality changes of DSW in the aquifer: These changes were subsequently plotted against dimensionless time parameter PV (Equation (1)), as shown for a selection of constituents in Figure 9. The resulting patterns were translated into leaching and retardation factors and semi-permanent concentration changes (Table 4). These observations assisted in fine-tuning the EL model, which is partly based on both calculated and observed values of L, R and DC.   The plots in Figure 9 show the observed concentration trends for selected constituents as function of the dimensionless time parameter "pore volume" (PV). The following classification of elements with similar behavior can be made on the basis of their concentration/PV-plot.

•
Ca, Mg and Sr (Ca and Mg shown in Figure 9) concentrations start low and with increasing PVs rise to an asymptotic value. All 2014 data plot below those of 2003-2004, which is partly explained by the lower DSW input in 2008. Only Ca and Mg show an initial concentration below DSW input, which is probably linked to Na exchange, with Ca developing an asymptotic value at or below DSW input. Contrary, Mg and Sr approach an asymptotic value little (Mg) or far (Sr) above DSW input. Obviously, Mg and even more so Sr are dissolved from the aquifer, and Ca is not. Sources of Mg and Sr could be dolomite and/or silicate minerals. • SiO 2 , K and V concentrations (SiO 2 and K shown in Figure 9) decline in an exponential way, but relatively slowly compared to the next group, and the low concentration asymptote is also relatively high and in most cases far above DSW input. This suggests that these elements are initially desorbed and that later on they could be dissolving from the aquifer, possibly from quartz (SiO 2 ), K-feldspar (K and SiO 2 ) or for instance pyroxenes (SiO 2 and V). A long storage time yields much higher concentrations for SiO 2 , pointing at slow dissolution kinetics. Vanadium is likely present as VO 4 3− , which could desorb from iron (hydr)oxide coatings.
• Na, As, B, Cr, F and Mo (Na and Cr shown in Figure 9) concentrations also decline in an exponential way, but much more rapidly than the previous group. The low concentration asymptote is relatively high and in any case significantly above DSW input, especially for As, B and Cr. The asymptote approaches DSW input better for Na, F and Mo. This suggests that all these elements are very rapidly desorbed and that later on especially As, B and Cr could be steadily dissolving from the aquifer matrix. Except for Na and B (as H 3 BO 3 ), this group is composed of anions, As as AsO 4 3− , Cr as CrO 4 2− , F as F − and Mo as MoO 4 2− . The source mineral of chromate and vanadate could be pyroxenes [23]. • Al, Fe and Ni concentrations (Al shown in Figure 9) decline in an exponential way, approximately as rapidly as the previous group. Their concentration seems to be dictated by colloidal particles that become initially mobilized. On the one hand, lack of filtration may have contributed to this. On the other hand, it is well known [32], that clay minerals in aquifers tend to be mobilized by deflocculation when the sodium adsorption ratio (SAR) of the infiltration water and of the native groundwater is high, their salinity low, the clay mineral content high, and the dominant type of clay minerals unfavorable (smectite > illite > kaolinite). The SAR value of the native groundwater seems to be more important than the one of the infiltration water. With SAR values of~20, the risk of mechanical clogging of the aquifer when the particles strand in the pore necks is relatively high. This risk is estimated higher in aquifer layer C than in A and B because SAR of native groundwater is higher and permeability is lower.
In Table 4, the estimated retardation (R) or leach factors (L) as derived from the PV-plots ( Figure 9) are listed, together with the estimated (semi)permanent concentration increase (DC). Table 4. Interpretation of water quality observations at the pilot, in terms of retardation factor (R), Leach factor (L) and (semi)permanent concentration change after break-through (∆C).

Parameter
It can be concluded that those parameters that could cause problems with drinking water standards (Table 2), such as EC, TDS, Na, Cl, SO 4 , F, As and Cr, are swiftly leached out from the aquifer, namely most of them within 1-3 PVs. The leach factor for Cr corresponds well with the chromate retardation factor as observed in a column dosage experiment [33], and with observations elsewhere [34,35].
The (semi)permanent concentration increases can, in general, be attributed to mineral dissolution or prolonged desorption from low permeability but high concentration pockets within the aquifer. The significant Ca decrease and Mg increase during the six-year storage period could point at dedolomitization, in which possibly also Sr 2+ is involved: where X + Y + Z = 1, possibly with Mg = 0.38 and Sr = 0.02. The relatively strong increase for SiO 2 during the six-year storage period could at least partly be explained by quartz dissolution at the relatively high water temperatures (35 • C). The following relation by [36] can be used to predict quartz solubility at temperatures 0-200 • C, yielding 14.3 mg/L at 35 • C: SiO 2 = 10 4.83− 1132 temp+273.15 (7)

Model Settings
The SWSR ASR system was modeled for just one of three recharge basins, assuming this to be representative for all three basins. The applied hydrogeological and geochemical stratification of the aquifer system into five layers follows the information given in Section 3. The low TDS DSW was selected as constant input to the recharge basin. Its quality is shown in Table 3.
The travel times to all recovery wells (condensed into one well per ring) and flux contributions from each aquifer layer to the well discharge, were calculated as indicated in Section 2.7. The combination of travel times leads to the so-called response curve of the whole well field. This is defined as the cumulative frequency distribution of underground travel times from the recharge basin to the wells during infiltration. The resulting response curve is shown in Figure 10.
Water quality changes in the recharge basin were completely ignored, because of exclosure from both sunlight and direct atmospheric contact, the high quality of DSW and the expected ultralow concentration of suspended solids. The formation of bottom muds in the basin was excluded. The complete lack of redox reactions necessitated the setting of the reactivity of bulk organic material at zero.

Quality Evolution during Nonstop Infiltration
If the infiltration via one of the recharge basins of the SWSR facility would last for 100 years without interruption, and the recovery wells would periodically pump out a water sample during this period varying between 0.1 and 100 years, then the abstracted water quality would be as predicted in Table 5. This table also shows the quality of DSW prior to infiltration and of native groundwater as abstracted by all ABC recovery wells.
Many things can be concluded from this table, because the recovered water quality also covers those moments that are representative for the 90 d recovery with or without drift. That aspect follows in Section 6.3.
Here, it is concluded that continued flushing leads to an initially fast and then slow elution of native groundwater, a significant leaching of the trace elements As, B, Ba, Cr, F and Sr, and a slow breakthrough of Cu, Ni and Zn which show higher concentration levels in DSW than in native groundwater. For selected elements, the breakthrough or leaching curves have been plotted in Figure 11. The bulk of changes clearly happens during the first year when the native groundwater is displaced.
The predicted water quality after 0.27 year is also representative of water quality after 90 d recovery following 10 years of lateral drift after 824 d of DSW infiltration. The quality after 0.74 year is as the one after 0. 27 year, but without lateral drift (thus, with direct recovery). The quality after 2.26 years is the water that is recovered on the first day after 824 d of DSW infiltration.

Quality Evolution during Nonstop Infiltration
If the infiltration via one of the recharge basins of the SWSR facility would last for 100 years without interruption, and the recovery wells would periodically pump out a water sample during this period varying between 0.1 and 100 years, then the abstracted water quality would be as predicted in Table 5. This table also shows the quality of DSW prior to infiltration and of native groundwater as abstracted by all ABC recovery wells.
Many things can be concluded from this table, because the recovered water quality also covers those moments that are representative for the 90 d recovery with or without drift. That aspect follows in Section 6.3.
Here, it is concluded that continued flushing leads to an initially fast and then slow elution of native groundwater, a significant leaching of the trace elements As, B, Ba, Cr, F and Sr, and a slow breakthrough of Cu, Ni and Zn which show higher concentration levels in DSW than in native groundwater. For selected elements, the breakthrough or leaching curves have been plotted in Figure 11. The bulk of changes clearly happens during the first year when the native groundwater is displaced.
The predicted water quality after 0.27 year is also representative of water quality after 90 d recovery following 10 years of lateral drift after 824 d of DSW infiltration. The quality after 0.74 year is as the one after 0.27 year, but without lateral drift (thus, with direct recovery). The quality after 2.26 years is the water that is recovered on the first day after 824 d of DSW infiltration.

Quality Evolution during Nonstop Infiltration
If the infiltration via one of the recharge basins of the SWSR facility would last for 100 years without interruption, and the recovery wells would periodically pump out a water sample during this period varying between 0.1 and 100 years, then the abstracted water quality would be as predicted in Table 5. This table also shows the quality of DSW prior to infiltration and of native groundwater as abstracted by all ABC recovery wells.
Many things can be concluded from this table, because the recovered water quality also covers those moments that are representative for the 90 d recovery with or without drift. That aspect follows in Section 6.3.
Here, it is concluded that continued flushing leads to an initially fast and then slow elution of native groundwater, a significant leaching of the trace elements As, B, Ba, Cr, F and Sr, and a slow breakthrough of Cu, Ni and Zn which show higher concentration levels in DSW than in native groundwater. For selected elements, the breakthrough or leaching curves have been plotted in Figure 11. The bulk of changes clearly happens during the first year when the native groundwater is displaced.
The predicted water quality after 0.27 year is also representative of water quality after 90 d recovery following 10 years of lateral drift after 824 d of DSW infiltration. The quality after 0.74 year is as the one after 0.27 year, but without lateral drift (thus, with direct recovery). The quality after 2.26 years is the water that is recovered on the first day after 824 d of DSW infiltration.

Quality Evolution during Recovery
The standard EL output list had to be elaborated, as indicated in Figure 7, in order to depict the water quality evolution during the 90 d of recovery (with or without 10 years standstill before). The result is shown for selected elements (Na, Cl, NO 3 , As, Cr and F), TDS and the contribution of DSW in Figure 12, and for many more water quality parameters in Table 5.
It is concluded that TDS and the concentration of most dissolved constituents slowly rise, while the percentage of DSW slowly decreases during progressive recovery. The main reason of this worsening water quality is the admixing of native groundwater. Without a preceding 10-year storage phase and thus practically without bubble drift, recovered water quality complies with the UAE drinking water standards (Table 2), with only As approximating the maximum permissible concentration near the end.
Ten years of storage with bubble drift have a more negative impact on water quality, as expected. This becomes manifest only halfway during the 90 d recovery phase (Figure 12), and starts to become critical in between 80 and 90 d, especially regarding Na, As, Cr and F. These constituents are predicted to exceed the drinking water MPCs: Na 170 (MPC 150 mg/L), As 11.2 (MPC 10 ug/L), Cr 70.7 (MPC 50 ug/L) and F 1.67 (MPC 1.5 mg/L). drinking water standards (Table 2), with only As approximating the maximum permissible concentration near the end.
Ten years of storage with bubble drift have a more negative impact on water quality, as expected. This becomes manifest only halfway during the 90 d recovery phase (Figure 12), and starts to become critical in between 80 and 90 d, especially regarding Na, As, Cr and F. These constituents are predicted to exceed the drinking water MPCs: Na 170 (MPC 150 mg/L), As 11.2 (MPC 10 ug/L), Cr 70.7 (MPC 50 ug/L) and F 1.67 (MPC 1.5 mg/L).

Discussion
The relatively simple modeling approach as presented in Section 2.7 and applied in Section 6, strongly relies on: (i) observations done during trials in the pilot area; and (ii) reversibility of the predicted water quality evolution during forward (infiltration) mode. This means that complications due to shorter or longer detention times in the aquifer (during storage and recovery) or due to processes like biofouling of the basin proximal zone should not occur.

Discussion
The relatively simple modeling approach as presented in Section 2.7 and applied in Section 6, strongly relies on: (i) observations done during trials in the pilot area; and (ii) reversibility of the predicted water quality evolution during forward (infiltration) mode. This means that complications due to shorter or longer detention times in the aquifer (during storage and recovery) or due to processes like biofouling of the basin proximal zone should not occur.
Problematic deviations from our predictions could arise, when the relatively high recovery rate (six times higher than the infiltration rate) would provoke serious upconing of brackish to saline groundwater from Aquifer 2, or when pockets of lower permeability in aquifer 1 would deliver more of their interstitial, low quality groundwater. A more dedicated spatial model, in which density effects and reactive transport (e.g., [19]) are combined, is needed to assess the potential impact of these processes.
There are other important questions still to answer, probably necessitating column studies and prolonged field monitoring during operation. The first question is: Which mineral sources could be responsible for the remarkable, semipermanent background mobilization of especially Cr, V, As and Sr, and how can this mobilization be mitigated? The second is whether clay minerals will deflocculate upon flushing of the target aquifer, and thereby clog it. There are weak signs that this might have happened to a minor degree during the pilot. The risk of clay mobilization will be higher during the operational phase, because still lower TDS infiltration water (DSW) will be used, and the ABC area contains some higher salinity subareas than the pilot area.

Conclusions
In this paper, we present the peculiar characteristics of an eolian-fluvial, (sub)oxic, calcareous sand(stone) aquifer system in a desert environment, the unique water quality changes during an ASR pilot including a six-year storage phase, and a relatively simple modeling approach to rapidly predict water quality of the recovered water, which is to be distributed without treatment.
In the applied Easy-Leacher model, pilot observations on retardation, leaching and (semi)permanent concentration changes are combined with a strong schematization of the complex basin-ASR system, groundwater flow (including bubble drift) and hydrogeochemical processes.
The hydrogeochemical simplifications are justified by the low reactivity of the aquifer, which is explained by the very low solute content of desalinated seawater (DSW), the infiltration of (sub)oxic DSW into an already (sub)oxic aquifer (thus, preventing redox reactions), a rapid leaching and, at high pH (7.8-8.5), high mobility of the problematic trace elements Cr (>95% as CrO 4 2− ), B (as H 3 BO 3 ), F, Mo (as MoO 4 2− ), and V (as VO 4 3− ).
The pilot observations and modeling results demonstrate that recovered water is most likely complying with Abu Dhabi's drinking water standards in both studied ASR scenarios. Without preceding 10-year storage (and thus hardly any drift, but with replacement of native groundwater), recovered water quality is expected to comply with the drinking water standards after 824 d of infiltration, even beyond the 67% recovered in 90 d (up to~85% recovery). The same conclusions also hold for the modeled scenario with 10-year storage (and thus significant bubble drift), however not for 90 but for 80 d. This means that the recovery efficiency will drop from~85 to 60%. After 80 d (60% recovered), MPC exceedances are then expected for: Na (20 above 150 mg/L), As (1 above 10 µg/L), Cr (21 above 50 µg/L) and F (0.2 above 1.5 mg/L). The main reason for their exceedances is the admixing of native groundwater.
This prediction has contributed to the implementation of some adaptations to the current ASR scheme operational since May 2015. The injection rate was raised by 20%, and during storage 568 m 3 /d will be infiltrated continuously during any storage period. By adjusting pumping schemes based on real time monitoring, the recovery efficiency can probably be maintained at about 80%.
If drinking water standards would become more strict, as is currently discussed in the Netherlands regarding As (from 10 to 1 µg/L) and CrO 4 2− (from 50 to 0.3 µg/L; [37]), then Abu Dhabi's drinking water standards will be exceeded in all modeled cases. This may be acceptable, however, during a short emergency situation after a calamity. Potential risks of upconing of brackish/saline deep, native groundwater, aquifer clogging by clay mobilization and basin clogging were not addressed, and may necessitate further research.
Acknowledgments: This study was funded by the Environment Agency of Abu Dhabi (EAD). Figures 1-3 were kindly provided by GTZ-Dornier. George Koziorowski (GTZ) is acknowledged for sharing some of his water quality data and for valuable discussions. Anonymous reviewers helped to improve the text.
Author Contributions: Pieter J. Stuyfzand evaluated all existing data, developed and applied the ASR model for predicting future recovered water quality, and wrote the paper. Koen G. Zuurbier conducted the field campaign in 2014. Niels Hartog performed the PHREEQC-2 modeling. Ebel Smidt and Mohammed A. Dawoud formulated the framework for this research and guided all steps in conducting this research, contributed by supplying information and edited texts.

Conflicts of Interest:
The authors declare no conflict of interest.