Anthropogenic wetlands due to over-irrigation of desert areas: a challenging hydrogeological investigation with extensive geophysical input from TEM and MRS measurements

.


Introduction
The rise of groundwater due to farmland over-irrigation, i.e., waterlogging, is a widespread occurrence with a large global impact on the available water resources.This comes with often contrasting effects on the water quantity, which increases, and the water quality, which deteriorates (Scanlon et al., 2007).Waterlogging occurs due to a number of causes, including large percolation losses from the irrigation water applied to the fields and seepage losses from channels providing significant localized sources of recharge water.One example of the groundwater level rising in semiarid regions due to man-made activities can be found in parts of the Xinjiang Province in western China.The diversion of river water for irrigation raised the groundwater levels from more than 7 m below the ground surface in the 1950s to about 1 m below the ground surface in the late 1980s (Sheng and Xiuling, 2001).Another example is in the Thar Desert of India, running along the western border with Pakistan.A mean water table rise of 1.1 m year −1 was recorded following the construction of the Indira Gandhi Nahar Pariyojana (IGNP) irrigation system.This system was started in 1960, and an area of nearly 3960 km 2 became susceptible to waterlogging at the end of 2001 (Choubey, 1997;Sharma, 2001).A third example is the farmland along the Lower Arkansas River Basin in Colorado that has been continuously irrigated since the 1870s.It began to develop shallow water tables by the early part of the 20th century, with an average water table depth of less than 2 m below the surface in 1999 (Gates et al., 2002).Waterlogging is usually accompanied by the salinization of the surface and subsurface environments (Leaney et al., 2003).
When over-irrigation is combined with low-permeable soils relatively close to the ground surface, low relief or depressions, and the absence of natural drainage either through surface water or groundwater systems, waterlogging can create anthropogenic perennial inland wetlands even in arid or semiarid regions.These "arid wetlands", i.e., humid zones in an arid or semiarid climate, which may seem to be a contradiction (Lemly et al., 1993), have been the focus of hydrological research over the last few years.This is because they represent regions of high conservation value, are crucial refuges for native wildlife (fish, amphibians, snails, and plants), and provide a habitat for migratory birds.Worldwide arid wetlands are threatened by increasing anthropogenic pressure, leading to water contamination through the development of agriculture in the surrounding area and/or water shortages due to surface water and groundwater use caused by population growth (Ashley et al., 2004;Li et al., 2015;Minckley et al., 2013).
An unequivocal contradiction arises when wetlands are generated and grow in desert areas due to the overuse of water for irrigation.In such cases, the wetlands reflect the consequences of extremely inefficient water management with the potential waste of a precious resource.Egypt is one of the countries where this occurrence is becoming widespread.Since the 1952 revolution, Egypt has tried to increase its agricultural area through the reclamation of desert land.Land reclamation in the Egyptian context means converting desert areas into agricultural land and rural settlements primarily by "adding water".The canals fed by the Nile River are extended through existing agricultural areas to supplement the new reclamation zones (Adriansen, 2009).Hundreds of deep wells have been drilled to tap the Nubian Sandstone Aquifer System in order to support agricultural megaprojects developed within the depressions in the middle of the Western Desert.Excessive pumping, accompanied by a lack of catchment hydrogeological and geomorphological studies, has caused the formation of large surficial ponds.An example of this occurs in the Baharia depression where El Bastawesy et al. (2013) have shown by remote-sensing analyses that an increase in cultivated land from 40.9 km 2 in 1984 to 95.5 km 2 in 2011 corresponded with a wetland enlargement from 0.3 to 4.7 km 2 .A similar situation has occurred in the El Fayoum depression, located about 100 km south of Cairo, which is one of the most important agricultural areas in Egypt.The development of the new reclaimed areas located in the desert land on the boundaries of older cultivated zones has been associated with excessive irrigation, seepage losses from canals carrying Nile water to the farmland, ineffective subsurface drainage, and lack of proper land development.The movement of seepage water into local depressions has resulted in the accumulation of surface ponds, particularly in the cultivated low-lying areas, resulting in damage to historical places such as the Hawara and Lahun pyramids (El Abd and El Osta, 2014).Using remotely sensed multitemporal data, similar occurrences have been highlighted by Arnous and Green (2015) in newly reclaimed lands to the east of the Nile Delta, a few tens of kilometers from the Suez Canal.Wetlands in low-lying areas are being created by the seepage of water from adjacent irrigated uplands.
In this paper, attention is focused on the Nubariya depression (30 • 41 N, 29 • 52 E) located at the edge of the Western Desert in the proximity of the western border of the Nile Delta (Fig. 1).The hydrological evolution of this zone is exemplary of the other reclaimed desert areas that have experienced over-irrigation.A large portion of the Nubariya zone has been converted from desert into irrigated farmland since the 1980s.An extensive network of irrigation canals has been constructed, transporting water from the Nile River.Overirrigation and seepage from the bottom of the reclamation canals has significantly raised the groundwater level, resulting in a surface water accumulation in the most low-lying zone since the early 2000s.Sharaky et al. (2008) noted that worsening groundwater quality has accompanied the growth of wetlands in this area.
In the framework of the EU SWIM IMPROWARE project (Innovative means to protect water resources in Mediterranean coastal areas through re-injection of treated water; http://www.improware.eu/),an extensive hydrological and hydrogeophysical study was performed in 2013-2014 in this area.With the goals of explaining the groundwater evolution in the depression and proposing a few scenarios for the reuse of treated wastewater (TWW) to reduce the aquifer salinization, a hydrogeological investigation with extensive geophysical input and numerical modeling was undertaken.The effectiveness of coupling hydrogeophysical techniques with numerical models in hydrological applications has been demonstrated in recent years over multiple scales, from lab to catchment (Binley et al., 2015).For instance, controlled infiltration experiments were monitored using time-lapse (TL) electrical resistivity tomography (ERT), crosshole TL-ERT, and TL ground-penetrating radar (GPR) with their outcomes properly assimilated in partially saturated groundwater flow models (Busch et al., 2013;Rossi et al., 2015).TL-ERT and airborne electromagnetic (EM) surveys were coupled to variable density groundwater flow models to investigate the freshwater-saltwater mixing in flat coastal areas (Comte et al., 2010;Teatini et al., 2010).Airborne EM data were used to distribute the hydraulic conductivity of catchment-scale groundwater flow models (Dickson et al., 2014;Marker et al., 2015).
Here, based on our previous experience in the use of hydrogeophysical methods for groundwater characterization (Behroozmand et al., 2013;Vilhelmsen et al., 2014), a joint magnetic resonance sounding (MRS) and ground-based transient EM (TEM) survey was performed on a ∼ 9 km 2 zone at the southwestern side of the Nubariya pond.The TEM method measures the electrical resistivity of the subsurface, which can be linked to water quality and soil type, in a threedimensional (3-D) setting (Foged et al., 2014).The MRS method is used to measure free water content (i.e., the total volume of water that can freely move and is not bound to the grain surface) directly and provides information on pore size distribution (Legchenko et al., 2004).The free water content is measured in percentage, and hence, in the case of a saturated porous media, the free water content measurement indicates the effective porosity.By combining both methods, a detailed hydrogeological model of the investigated area (see Sect. 2 for details) in terms of hydrogeological structure and water quality was built down to a depth of about 100 m, which is the maximum depth of the investigation.
With the main goal of reconstructing the past multiyear hydrological evolution of the area driven by the need for reclamation and its possible future behavior in relation to artificial aquifer recharge using TWW, two numerical models were developed: a 3-D groundwater flow model for the entire depression area (25 × 35 km) and a 3-D density-dependent groundwater flow and transport model for a smaller zone (2.7 × 2.9 km) at the southwestern side of the pond area caused by over-irrigation (see Sect. 2 for details).The models were constrained by using the results of the geophysical survey, an interpretation of pump tests, a digital elevation model of the depression surroundings, and information on the evolution of the pond and the irrigated area derived from satellite images acquired since 1984.The local model is nested into the depression-scale model in the sense that the results from the latter in terms of the piezometric head are used as boundary conditions in the former.
In developing countries and many other parts of the world, these reclaimed desert environments are generally quite challenging to access and investigate.Moreover, the development of anthropogenic wetlands as a byproduct of over-irrigation is a process that is generally not well investigated in the hydrological community.Here, we present a unique case study where the evidence of the occurrence is profound.Finally, from a more technical point of view, the advantages of integrating MRS and TEM methods to discriminate between different deposits and groundwater quality are clearly highlighted for the first time, to our knowledge, in wetland environments.
The paper begins with a section that describes the environmental setting and hydrological evolution of the study area.Second, the MRS and TEM methods are briefly introduced, the results from the surveys carried out in Nubariya are summarized, and a hydrogeological model obtained by the joint inversion method is shown.The setup of the depression-scale and local-scale numerical models is presented and the simulation outcomes are shown.The final sections discuss the value of the proposed approach and draw the conclusions.

The Nubariya depression: evolution and challenges
The study area is located in a complex sedimentary setting along the transition between the Quaternary deposits of the Nile River Delta to the west and the Pliocene and Miocene deposits of the Western Desert to the east (Fig. 1).The area is characterized by low relief with an elevation from 0 to 100 m above the mean sea level (a.m.s.l.).The main landforms developed through the interaction of the geological structures, the processes of wind and surface water erosion, and the climatic conditions, such as temperature and humidity (Dawoud et al., 2005;Hassan et al., 2012;RIGW and IWACO, 1991).Hydrogeological investigations at the regional scale (Sharaky et al., 2008) highlight the presence of loose coarse Miocene sands with clay lenses overlain by Pliocene deposits in the upper 200 m of subsoil.The Pilocene units consist mainly of estuarine clayey facies at the base and fluviomarine and shallow marine limestones at the top.The uppermost units are exposed in the lowest parts of the landscape and their vicinities.Quaternary deltaic sediments from sands to silt dominate in the eastern part of the study area.The Quaternary layer gradually decreases in thickness from the Rosetta branch of the Nile River to almost zero along the linear depression from Nubariya to the north and Wadi El-Natrun southward, where they interdigitate more or less sharply with the Pliocene deposits.
The Nubariya depression represents the northernmost tip of the Wadi El-Natrun elongated depression situated approximately 50 km to the south.Unlike the latter, which has been extensively investigated from the hydrological, hydrogeological, and environmental perspectives over the last few decades due to the presence of alkaline lakes (Atwia et al., 2012;Khalil and Santos, 2013), very little hydrogeological information is available on Nubariya, even though it has been subject to more recent reclamation and development.
With a minimum elevation of about 7 m a.m.s.l., the Nubariya depression is surrounded by sandy hills with a mean elevation ranging between 30 and 100 m a.m.s.l.(Fig. 2).A digital elevation model (DEM) of the area has been obtained from the ASTER GDEM, which was calibrated by a kinematic DGPS in situ survey carried out in 2014.Over the last three decades, large portions of this desert area have been converted into productive farmland.A remote-sensing investigation with a series of Landsat images acquired between 1984 and 2014 clearly shows that a pond started to develop in 1999 and expanded significantly until the 2014 extent reached about 2.5 km 2 (Fig. 3).The pond enlargement corresponds to the expansion of the cultivated areas surrounding the depression, providing evidence of a connection between the two (Fig. 3).The reclamation efforts in the Nubariya area relied on Nile water diverted through a set of main (the Nubariya Canal) and secondary canals, which are usually unlined.The surface water systems in Nubariya, in most cases, cut through sands; therefore, a direct connection between the surface and the groundwater exists in most parts of the irrigation system.Seepage from the groundwater to the surface water system and vice versa is very common depending on the water levels.Overall, the Nubariya canal and its main branch (the Nasr Canal) recharge the groundwater system in the study area.Because of an inadequate network of drainage ditches, excessive irrigation (via surface irrigation techniques), and seepage from the bottom of the reclamation canals, the groundwater level rose and formed a large water accumulation in the depression.The wetland enlargement was followed by a general deterioration of the surface and subsurface environment.Farmhouses and roads were permanently flooded (Fig. 4).The water quality was degraded by evapotranspiration and the accumulation of salts in the soils associated with shallow water tables and fertilizer leaching into the underlying aquifers that discharge to the pond (Table 1).Despite the nature of the pond, coliforms are practically absent (Fadlelmawla, 2014) due to the high salinity (Dawe and Penrose, 1978).
The artificial pond has grown despite the climatic conditions.The typical climate that prevails in the area is characterized by high temperatures that span from a high of 30 • C in the summer to a few degrees above 0 on winter nights.Rainfall occurs from December to February with an average of 100 mm year −1 recorded in the northwest part of the area.Usually, the relative humidity ranges from 88 % in winter to 40 % in summer.The mean monthly value of potential evaporation ranges between 15.4 mm day −1 in June and 1.6 mm day −1 in January (Dawoud et al., 2005).

Hydrogeological survey
Hydrogeological and chemical investigations were carried out to collect information about the link between the subsurface lithostratigraphy and the water quality in the pond, thus providing specific ancillary data for the calibration of geophysical measurements.With this aim, five boreholes (Fig. 5) were drilled with continuous core sampling and pumping and recovery tests.The depth of the drilling ranges between 25 and 55 m from the ground surface.Moreover, a number of slug tests were carried out to characterize the infiltration rate, and hence the soil type, of the uppermost deposits.
The lithological analyses of the recovered core revealed the presence of a top 10 to 15 m thick sandy unit underlain by a 15 to 20 m thick clay layer overlying the limestone (fine sand) at the base borehole.The interpretation of the well tests provides estimates of the limestone hydraulic conductivity K and the elastic storage S at 1 × 10 −4 -6 × 10 −4 m s −1 and ∼ 8 × 10 −3 , respectively.The slug tests and the data in the literature (Dawoud et al., 2005) suggest K values of approximately 1 × 10 −4 and 1 × 10 −7 m s −1 for the top sand and the clayey deposits, respectively.
The results of the chemical analyses of the groundwater and pond water samples are summarized in Table 1.The measurements highlight a deterioration of the groundwater quality, mainly in terms of NaCl concentration.Moreover, the similar proportional concentrations of various anions and cations in the water samples taken in the pond and in the boreholes suggest significant mixing between the two water bodies.

Transient electromagnetic survey
The transient electromagnetic (TEM) method has been widely used in groundwater studies (Auken et al., 2003;Danielsen et al., 2003;Siemon et al., 2009) and offers a relatively fast and cost-effective way to obtain information about the ground electrical resistivity down to depths of a few hundred meters.The ground-based TEM method employs a transmitter loop deployed on the ground surface, which generates a primary magnetic field.After a short duration (∼ 10 ms), the current is switched off abruptly (within a few µs), and the secondary rate of change with the time of the magnetic fields from the induced eddy currents in the ground is measured using an induction receiver coil on the surface.Hence, the observed data are voltage data as a function of time.The early-time data, typically measured from a few µs after the current is shut off, primarily reflect the resistivity of the top layers.The late-time data, typically measured until a few ms after the current is shut off, provide the resistivity information of the deeper layers.These data, measured at different times after the current is shut off, are then interpreted in terms of the subsurface resistivity as a function of depth.For basics of the TEM method, see Christiansen et al. (2006).

Field survey and measurement setup
The TEM survey was carried out in December 2013 using the WalkTEM system (Nyboe et al., 2010).The full survey was comprised of 127 soundings, 110 of which were located on a 200 m spacing grid within the 2.7 km × 2.85 km groundwater model area.The remaining 17 soundings are located outside the grid to investigate regional geological changes (Fig. 5).The WalkTEM system was configured in a central loop with a dual-moment setup (low moment and high moment) in order to obtain soundings containing high-quality early-time measurements using a low magnetic moment and late-time measurements with a high magnetic moment.The two moments cover the entire depth range from the surface down to approximately 150 m.Furthermore, two different setups were used for the survey, one using a 40 m × 40 m transmitter loop and another using a 100 m × 100 m transmitter loop, resulting in a larger moment and hence a larger depth of investigation.Table 2 documents the used WalkTEM system setup.The instrument was calibrated prior to the survey following the approach of Foged et al. (2013).

Results
The quality of the collected data is high and appropriate for the planned analyses.For the majority of the measurements, there are usable data points from 5 µs to 10 ms for the 40 m × 40 m loop and up to to 20 ms for the 100 m × 100 m loop, respectively.Consequently, for most of the measurements, the high signal-to-noise ratio required only a limited amount of processing.However, in some parts of the survey area, the signal level for the late-time gates fell below the noise level and these gates were removed.In the following processing, the data were inverted using the AarhusInv inversion code (Auken et al., 2015) with a multilayer model setup (here 20 layers) in which the layer boundaries are fixed and the resistivities of the neighboring layers are tied together via smoothing constraints.In addition to the characterization of the overall structural model of the subsurface, multilayer models help assess the homogeneity of each unit.
2 * E c is electrical conductivity and TDS is total dissolved solids.NA: not available.
Figure 6 shows a west-to-east cross section through the resistivity structure obtained from the TEM data.The location of the profile is shown as a black line in Fig. 5.The color scale is faded white below the depth of investigation, i.e., the depth to which the model can be trusted (Christiansen and Auken, 2012).The results of the TEM survey suggest four distinct layers in the upper 100-150 m.The top layer is not continuous, vanishing in the easternmost part of the profile.It is characterized by resistivity values greater than 10 m and a maximum thickness of approximately 20 m.Fig. 6 represent the MRS results in the vicinity of the profile, as discussed in the following section.

Magnetic resonance sounding survey
Magnetic resonance sounding (MRS; also called surface NMR) is an electromagnetic geophysical method that noninvasively measures water content and pore structure.The method works based on the physical principle of nuclear magnetic resonance (NMR) and is directly sensitive to water molecules and their interaction with the pore space.In short, using Earth's magnetic field, resonance excitation is induced by passing a tuned AC current into a large wire loop deployed on the ground surface.The corresponding energizing magnetic fields propagate in the subsurface; at any position in the subsurface, a component of the energizing field excites the water molecules.The measurement continues by turning off the energizing pulse and recording the responses from all subsurface excited protons as they gradually return to their initial (equilibrium) state.This experiment is called a free induction decay (FID) and is mostly used for MRS data acquisition.The experiment is repeated for a number of energizing pulses at different pulse moments (defined as the product of the current amplitude and the pulse duration) through which different volumes of the subsurface are sampled, and therefore depth information is provided.For more information about the principles and application of the MRS method, see Behroozmand et al. (2015).
In near-surface geophysics, MRS is commonly used to estimate free water content (or porosity in the case of saturated porous media) and hydrogeological properties such as pore size and hydraulic conductivity.The initial amplitude of the MRS decaying signal is proportional to the water content, while its relaxation rate provides information about the pore structure.A small relaxation time typically indicates finegrained material, whereas a high relaxation time indicates coarse materials.The method is well established for nearsurface characterization (Chalikakis et al., 2008;Günther and Müller-Petke, 2012;Knight et al., 2012).Different studies have also dealt with estimating hydrogeological parameters from MRS and have found good correlation between MRSderived parameters and those from borehole aquifer tests (Boucher et al., 2009;Herckenrath et al., 2012;Plata and Rubio, 2008;Vilhelmsen et al., 2014Vilhelmsen et al., , 2016;;Vouillamoz et al., 2012Vouillamoz et al., , 2015)).

Field survey and measurement setup
The field campaign consisting of six MRS soundings was performed in January 2014 under the same hydrogeological conditions as the TEM survey.The locations of the MRS soundings are shown in Fig. 5. Prior to the survey, a noise scouting study was carried out over the entire area to investigate the ambient electromagnetic noise condition and to propose potential locations for acquiring MRS data.The results of our noise scouting presented low noise levels in the area, which suggests that it is a good site for measuring MRS data using a relatively low number of repeated measurements (free induction decays; FIDs).We used the NU-MISPoly system (IRIS Instruments, Orleans, France) for the data acquisition, and the system was configured in coincident loop geometry with a 100 m side square loop as both a transmitter and a receiver.In addition, up to two reference loops were deployed for recording noise.measurement configuration used in the study.The measurement sequence consisted of a noise record before excitation, an energizing pulse, and a measurement dead time (to switch from the transmit to receive phase) followed by a recorded FID (see Table 3 for details).We measured FIDs for 16 pulse moment values, and a stack size of 30 was more than enough to obtain high-quality data.

Results
The data were processed following the approach described by Dalgaard et al. (2012) and Larsen et al. (2013) and inverted jointly with the TEM data following Behroozmand et al. (2012a, b).We used the AarhusInv code (Auken et al., 2015) for the inversion of the data.As an example, the bars on the top of the TEM plots in Fig. 6  simulated data (column two), and the weighted data residuals (column three) as well as the TEM data fit (column four).
The estimated model fits both datasets very well, as shown in columns three and four.Row two shows the inversion results (black lines) in terms of resistivity (column one), free water content (column two), and relaxation time (column three).
The gray error bars denote the model parameter uncertainties, shown as the 68 % confidence intervals.The model is well determined and is in good agreement with the lithological information obtained from a nearby borehole (column four).It is noteworthy that the MRS results vary along the profile in Fig. 6, suggesting that the clay content of the second layer may vary.
Overall, the outcomes of the hydrogeophysical surveys integrate satisfactorily with the lithological and hydrological characterization available from the wellbore information.
The derived hydrogeophysical model was used as an input in the hydrological model of the area, as will be discussed in the following section.

Hydrogeological modeling
Numerical modeling was developed to gain some fundamental insight into the processes that explain the observed pond development and provide some preliminary information on the possibility of using TWW for managed aquifer recharge (MAR) in the area.Based on the hydrogeophysical characterization of the area, we set up a simplified groundwater flow depression-scale model to verify the possible connection between the desert reclamation and the wetland growth and a density-dependent flow and transport local-scale model to test the potential and the effect of recharging the saline aquifer by fresh TWW.
The models, particularly the one developed at the local scale to investigate possible MAR scenarios, rely strongly on the outcome of the hydrogeophysical investigations.Both the geometry of the hydrogeological layers (maps of the top and bottom of the various units) and the distribution versus depth of the groundwater quality have been derived through the integration of the TEM and MRS results.

Model setup
The objective was to investigate the possibility of a surplus of irrigation water being responsible for the formation and growth of the pond at the boundary of the desert area by modeling the hydrological response of the zone to the reclamation carried out from 1984-2014.By simulating the decadal rise of the water table, the model also aimed to show the inefficient use of the available water resources.Furthermore, the results from the regional-scale hydrological model can be used to provide boundary conditions for the local-scale MAR model.
For this purpose, we used a 3-D Richards' equation solver that is part of the more general code CATHY (Camporese et al., 2010).The model domain is showed in Fig. 2 by the black box.The area size amounts to 800 km 2 , i.e., approximately 32 and 25 km in the southwest-northeast and northwestsoutheast directions, respectively, and is centered on the depression.Such a large area has been selected in order to limit the effects of the large uncertainties related to the boundary conditions on the zone of main interest, i.e., the pond surroundings.According to the available piezometric information (Dawoud et al., 2005) and the digital elevation model (DEM), the domain orientation has been selected in order to use a no-flow condition along the southwest-northeast (A-B and C-D, Fig. 2) lateral boundaries of the model.The DEM represents the top surface of the model with a no-flow bottom bound placed at −34 m a.m.s.l., i.e., approximately in the middle of the first confined aquifer.The boundary conditions defined on the current land surface and along the A-D and B-C (Fig. 2) lateral sides are quite complex.Regarding the ground surface and keeping in mind the multidecade reference period spanned by the simulations, the following conditions were used: a seepage face condition on the depression zone, i.e., for the surface nodes around the pond with an elevation of less than 15 m a.m.s.l.; a net specific recharge of 400 mm year −1 is uniformly distributed on the reclaimed zone.Following the development of the farmland as retrieved from remotesensing data (Fig. 3), the recharge zone was assumed to increase from 1984 to 2014 as shown in Fig. 2. The recharge amount was obtained by combining the average estimates of the aquifer recharge by downward percolation of the irrigation surplus water (in desert areas, relatively high 365-550 mm year −1 leakage rates have been quantified for basin, furrow, and sprinkler irrigation with much lower 40-190 mm year −1 rates for drip and central pivot irrigation (Dawoud et al., 2005), rainfall infiltration, and exfiltration due to evaporation); a null infiltration/exfiltration in the desert zone, i.e., the portion of the domain complementary to the reclamation area.Due to the relatively large depth to the water table and the lack of specific data, a null balance between rainfall and evaporation has been assumed for simplicity.
Fixed and time-dependent Dirichlet conditions were prescribed along the northwest-southeast bounds toward the desert (A-D, Fig. 2) and the Nile Delta (B-C, Fig. 2), respectively.The water table was set to 0 m a.m.s.l.along A-D with a sensitivity analysis that has demonstrated a negligible influence of the selected value on the piezometric evolution in the area of interest.The information at the regional scale reveals that the water  region (El-Sayed et al., 2012) .However, due to the lack of specific data, we have used the rise of the water table along B-C as a calibration parameter in order to match the date of the pond formation and the pond enlargement versus time.This represents the best information available to calibrate the model at the scale of interest.The approach suggested a linear rise from +10 m a.m.s.l. in 1984 to +17 m a.m.s.l. in 2014, i.e., to only a few meters below the land surface.These values seem plausible on account of the available regional information (El Molla et al., 2005;Mabrouk et al., 2013) and a few surveys carried out in the area.
The DEM was used to derive the grid discretization in the horizontal direction, with a characteristic element dimension that reduces from 1000 m on the external boundaries to 50 m around the pond to allow for a more accurate representation of the infiltration/exfiltration processes and changes in the time and space of the groundwater table and saturation degree.The 3-D finite element (FE) grids consists of 146 275 nodes and 830 808 elements with 24 layers of tetrahedra that are used to represent the three main lithostratigraphic units defined from the hydrogeological and hydrogeophysical investigations (Fig. 8).
The hydrological properties (K and S) of the three main hydrostratigraphic units were assigned according to the values described above, with the parameters of the van Genuchten retention curves for the top sandy unit derived using the Rosetta code (Schaap et al., 2001); ψ s (capillary pressure) = −0.37 m, n = 4, and S wr (residual water saturation) = 0.15.

Results
A steady-state simulation was initially performed to fix the initial condition in the whole domain for the transient run.The solution was obtained by prescribing a water table elevation of 10 and 0 m a.m.s.l. on the B-C and D-A (Fig. 2) boundaries, respectively.The initial distribution of pressure and saturation degree S w , thus obtained, can be considered representative of the hydrological condition in 1984 when the area was still totally desert with negligible reclamation (Fig. 3a).With the prescribed setting, the depth to the water table is approximately 10 m below the most depressed portion of the Nubariya area (Fig. 9a).
Then, the model was applied over the period between 1984 and 2014.Figure 9 shows the evolution of S w along the vertical cross section E-E traced in Fig. 2, as computed using the CATHY code.The distribution of the saturation clearly describes the changes in the water table that evolved due to water leakage from the land surface into the aquifer.The recharge raises the water table in the eastern part of the domain with a groundwater level that reached the bottom of the depression around 2007-2008.From 2007 to 2014, the water table rose further, completely saturating the lowest por-tion of the depression.This caused the waterlogging of the most depressed zone and the pond development and growth.A comparison between the model outcome and the pond's extent acquired from satellite images is presented in Fig. 10 for the years 2007 (a), 2010 (b), and 2014 (c).The contour line S w = 0.9 can reasonably be assumed as the boundary of the waterlogged zone.The modeled evolution of the pond extent satisfactorily matches the remotely sensed information.It is noteworthy that some westward portions of the pond located outside the simulated zone with S w > 0.9 are waterlogged due to deep excavation for construction material.

Expected evolution linked to MAR scenarios
MAR using treated wastewater is one of the main challenges in Egypt with the reduction of freshwater due to climate change, the growth in water demand due to population increase, and the enlargement of desert areas reclaimed for agricultural purposes.The reuse of TWW (municipal and to some extent industrial wastewater) is considered an effective water saving measure in areas where this water would otherwise be lost from the Nile system.The primary use of TWW is for the irrigation of green areas (landscape development) and non-food agriculture or the improvement of aquifer quality by mitigating or countering saltwater intrusion and contamination.The planned extent based on treated wastewater is some 250 000 feddan (approximately 1000 km 2 ).Most of the planned area is located in the western and eastern delta with greater Cairo and Alexandria as the main suppliers of treated wastewater (United States Agency for International Development; USAID, 2010).
Within this general context, one of the aims of the IM-PROWARE project was to perform a preliminary evaluation of the effectiveness of MAR in the Nubariya zone using TWW provided by the Nubariya wastewater treatment plant (WWTP).An existing WWTP located in the area surrounding the pond has been updated with a constructed wetland (CW) tertiary system in order to improve the quality of the plant effluent to a level compatible with MAR.The CW treatment capability in the present condition amounts to 160 m 3 day −1 with a potential increase of up to 6000 m 3 day −1 , representing the working WWTP flow rate of domestic wastewater.

Model setup
A number of preliminary simulations were performed at a local scale to understand the effect of recharging the aquifer system in the Nubariya.This local scale coincides with the area where hydrogeophysical surveys have led to an in-depth understanding of the geologic setting.The US Geological Survey SUTRA code (Voss and Provost, 2002), which can handle density-dependent flows under saturated-unsaturated conditions, was used to simulate a number of scenarios of artificial aquifer recharge.Based on the hydrogeological set- ting of the study area, the limestone unit identified in the depth range between 0 and −100 m a.m.s.l.(Fig. 8) was selected as the target aquifer for water injection because of its uniform and large thickness, moderate salt contamination, and relatively high hydraulic conductivity.At this preliminary stage of the study, the following major assumptions were used: (i) the injected water is fresh compared to the formation water; (ii) well clogging is neglected; and (iii) the hydrological regime at the model boundary is constant in time and derives from the outcome of the depression-scale model as of 2014.
The model domain extends approximately 2.7 km in the southwest-northeast direction and 2.85 km in the northwestsoutheast direction (Fig. 2), roughly coinciding with the area characterized by the MRS and TEM surveys at the southwestern tip of the pond.The model spans the depth range from the ground surface to the bottom of the fine-sand confined aquifer, approximately −100 m a.m.s.l.The latter was mapped by interpolating the models derived from the few deep TEM soundings (Fig. 11).The domain was characterized by the three main hydrostratigraphic units, i.e., sandy, clayey, and limestone layers, consistent with the available data and the depression-scale model.The hydrogeological parameters were also assigned accordingly, in agreement with the values used in the regional model.The maps of the depth of the unit interfaces (top and bottom surfaces) were obtained from the TEM and MRS measurements (Fig. 11).
The domain was discretized into prismatic elements.A 2-D grid was initially developed composed of 7077 nodes and 6965 quad elements with the dimension ranging between 50 m on the external boundaries and 20 m in the central area where injection is planned.Then, the 2-D FE grid was "projected" vertically to generate a 3-D FE mesh for a total of 219 387 nodes and 208 950 elements.A vertical discretization into 30 layers 0.5 to 15 m thick allows for an accurate reconstruction of the geological formations detected in the area (Fig. 12).
The distribution of the salt concentration c has initially been assumed constant in each formation, as suggested by the layering provided by the geophysical investigations.The c values have been derived from the analyses of the collected water samples (Table 1).The solute concentrations, expressed in terms of the dissolved mass fraction (M s / M, where M s is the mass of the dissolved components and M is the fluid mass), are the following: c top sand = 0.013 and c limestone = 0.016 kg salt / kg water , i.e., 16 g L −1 .An intermediate value has been assigned to the clay layer, and we have supposed c = 0 in the unsaturated zone.The following boundary conditions were prescribed: -Null groundwater and concentration flux along the southwest-northeast boundaries (side A -B and C -D in Fig. 2).These bounds have been properly selected along a direction parallel to the main flow direction as emerged from the depression-scale model; -Dirichlet constant conditions on the northwestsoutheast bounds (side A -D and B -C , in Fig. 2).The pressure distribution provided by the depression-scale

Results
A steady-state simulation was initially performed to acquire an equilibrated condition in terms of pressure and concentration in the whole domain to be used as the initial state for the transient runs.Due to the lack of specific tracer tests, a sensitivity analysis on the hydrodynamic dispersivity (longitudinal dispersivity α L and transverse dispersivity α T ) was also performed, spanning the range 10 < α L < 100 m as suggested by Gelhar et al. (1992) for a characteristic problem dimension on the order of 100 m.The outcomes provided in the following were obtained using α L = 20 m and α L /α T = 10 with the extension of the mixing zone m L between fresh and in situ brackish water in the longitudinal (horizontal) direction that ranges between 280 m (α L = 10 m) and 570 m (α L = 100 m); m L = 340 m with α L = 20 m.Three major scenarios were used to address the variability on the injection layout: S 1 ) a single well; S 2 ) three wells along an alignment spaced between d = 60 m (S 2A ) and 100 m (S 2B ) with the central well in the same position as the well in S 1 ; and S 3 ) three wells disposed in a radial configuration with the circle radius varying between r d = 60 m (S 3A ) and 100 m (S 3B ) and the center in the same position as the well in S 1 .A 20 m long well intake is positioned in the lowest portion of the limestone aquifer, approximately between −90 and −70 m a.m.s.l., in order to keep the freshwater volume stored in the subsurface as compact as possible by taking advantage of the lighter injected freshwater with respect to the denser brackish groundwater.The injection rate of Q inj = 6000 m 3 day −1 is simulated for a period of 2 years.Q inj is evenly distributed between the three wells in scenarios S 2 and S 3 .
The results in terms of c distribution at the end of the simulation period are presented in Figs. 13 and 14. Figure 13 shows the salt concentration at the top of the injected aquifer for the various scenarios listed above.The evolution of c in time and space on a vertical cross section through the barycenter of the well layout is presented in Fig. 14 for scenario S 3A .The injected water moves upward to the bottom of the clay layer, which stops the vertical migration, and then spreads horizontally, forming a significant volume of freshwater.The freshwater gently flows southwestward due to the natural hydrodynamic regime.The value d = r d = 60 m represents the maximum distance between the wells that allows for a compact volume of freshwater after 2 years of injection.The use of three wells instead of one allows for a reduction in the maximum overpressure in the injection cells from 2.6 m to about 1.1 m.
The fate of the injected wastewater can play an important role in water reuse because it may result in a further quality improvement due to the well-known phenomenon of soil aquifer treatment (SAT) (Idelovitch et al., 2003).Therefore, we have performed a last simulation in order to investigate the possibility of withdrawing the water previously injected into the limestone aquifer.With reference to the well layout implemented in scenario S 3A , a pumping well located in correspondence to the barycenter of the injection boreholes (indicated by the X in Fig. 13) is intended to extract groundwater from the top portion of the aquifer starting after 2 years of water injection, i.e., from the final condition reached by scenario S 3A .Injection and withdrawal simultaneously occur for a 2-year period.The intake of the pumping well is 10 m long.The results are summarized in Fig. 15, showing the evolution versus time of the salt concentration in the grid cell where the pumping well is located.After the first 2 years, when c reduces from the initial groundwater value to approximately zero, the groundwater withdrawal at a rate equal to the cumulative injection rate (Q withdrawn = Q inj = 6000 m 3 day −1 ) yields a certain deterioration of the groundwater quality in the well surroundings, and hence of the extracted water, with a salt concentration that rises to approximately 4 g L −1 in steady-state conditions.Therefore, Q withdrawn must be reduced to make the groundwater suitable for agricultural purposes.The numerical simulation points out that almost fresh groundwater (c < 1 g L −1 ) can be permanently obtained by the system if the pumping rate is reduced to 2000 m 3 day −1 .Notice that the injected TWW takes about 80 days to flow from the injection to the withdrawal intakes, which are about 100 m apart.Both the time period and the distance are significant in terms of suspended solid filtration and organic matter biodegradation using secondary or tertiary wastewater effluents (Idelovitch et al., 2003;Kopchynski et al., 1996).

Model evaluation
Generally, the numerical models are simplifications of real aquifer systems.The model results are affected by (1) the numerical approximations used to solve the groundwater flow and transport equations, (2) the discretization of the modeled area, and (3) the availability and accuracy of the hydrogeological data used to define the spatial distribution of the Hydrol.Earth Syst.Sci., 21, 1527Sci., 21, -1545Sci., 21, , 2017 www.hydrol-earth-syst-sci.net/21/1527/2017/ physical parameters, the boundary conditions, and the factors forcing the system (Idelovitch et al., 2003;Tsang, 2005).Another important limitation of calibrated numerical models is the non-uniqueness of their solutions.
The models presented here rely on a number of assumptions, starting from the hydrogeological structure of the system and the distribution of the hydrological properties (K, S, α L , and α T ).The majority of the available stratigraphic information is located in the area surrounding the contemporary pond with little lithological data collected in the Western Desert zone.However, the results from the hydrogeophysical surveys carried out at the study site and the hydrogeological characterization carried out in other low-lying zones to the south, along the Wadi El-Natrun depression (Ammar, 2010;Khalil and Santos, 2013;Monteiro Santos and Sultan, 2008), have pointed out a remarkable lateral continuity and relative homogeneity of the main hydrological formations addressed by the modeling study.
Regarding the issue related to lithological and hydraulic heterogeneity, the two models did not address the spatial variation of the hydrological parameters, namely the hydraulic conductivity, elastic storage, and longitudinal and transversal dispersivity.The construction of the models in this manner would have required detailed hydrogeological information from aquifer tests and decadal piezometric records, which are not available.However, an accurate calibration of the depression-scale model is beyond the scope of the study.Indeed, its main objective was to evaluate whether a groundwater flow model is able to reproduce the formation and growth of the pond as observed by satellite data by im-plementing reasonable values of K and S (derived from the few pumping and recovery tests performed in the area) and estimating the factors forcing the system.The pond evolution over time allows for a global evaluation of the model representativeness of the hydrological processes in the area.
As reported above, the main assumption in the modeling of MAR is related to the disregard of aquifer clogging in the surrounding area of the injection wells.The deterioration of the aquifer properties due to biological, chemical, and physical clogging strongly depends on the quality of the water used in MAR (Bouwer, 2002).A number of methodologies have been proposed in the literature to account for the decrease in porosity and thus the conductivity of the medium due to clogging (Pérez-Paricio and Carrera, 2000).However, at the present stage of this study, the lack of specific information on the quality of the effluent from the Nubariya WWTP and the chemical and physical characteristics of the limestone composing the aquifer do not allow for a specific numerical investigation on the clogging effects.
Concerning the dispersivity values, a sensitivity analysis on α L and α T was carried out using the local model.The effect of the dispersivity values on the shape of the injected freshwater volume was investigated to overcome the lack of dispersion and tracer tests.Due to the relatively small dimension of the local model (on the order of 2-3 km) and the homogeneity of the limestone aquifer suggested by the geophysical results, the assumption of the parameter homogeneity seems more reliable here than in the depression-scale model.Finally, we highlight that the hydrogeophysical data acquired in this study are of high quality with good spatial coverage of the layered aquifer model.In a more complex environment, e.g., with a larger lateral heterogeneity, a higher density of hydrogeophysical measurements might be required to provide accurate hydrogeological models to be used in hydrological modeling.

Conclusions
In this study, we evaluated the possible reuse of treated wastewater in the Nubariya depression in Egypt using extensive hydrogeophysical input and numerical modeling.The geophysical surveys successfully characterized the aquifer system in the study area where distinct layers of unsaturated sand, saturated sand, and sandy clay were found.In addition, the combined use of the MRS and TEM methods provided information about the spatial distribution of the layers as well as the water quality and clay content.Within the sandy layers, we found that the water quality is inappropriate for use as drinking water or for agricultural purposes.In addition, a relatively shallow confining clay layer was found as one of the reasons for the evolution of the artificial pond in the area.Below the clay layer, a 100 m thick aquifer with a high free water content was nominated as the potential aquifer for recharge.
The estimated hydrogeophysical model was used as an input for building the hydrological model of the area.The simulations carried out at the depression scale clearly pointed out the inefficient use of the freshwater resources in the study area.Although not very well known, the problem of waterlogging and pond development within desert zones caused by over-irrigation has been reported for several locations worldwide.This study investigates a phenomenon that represents a contradiction in areas where long-term programs are underway to effectively reuse treated wastewater: large volumes of good quality freshwater are wasted and large economic investments are made to improve and reuse TWW.
The same low-permeability layer responsible for the formation of the artificial pond in Nubariya confines a relatively thick limestone aquifer and bounds the upward migration of the treated water injected into the system.The continuity of the clay unit, an important requirement in view of the recharge, was clearly highlighted by the geophysical surveys.Although preliminary, a local-scale density-dependent groundwater flow and transport model allowed for the development of an optimized MAR scheme.
This research highlights the hydrological challenges related to the effective management of water resources in reclaimed desert areas and the effectiveness of using advanced geophysical and modeling methodologies to characterize the subsurface environment, investigate naturally and/or anthropogenically driven exchanges between groundwater and surface water, and plan appropriate interventions aimed at the efficient use of available water.

Figure 1 .
Figure 1.(a) The location map of the Nubariya depression between the Western Desert and the western Nile Delta in Egypt and a Landsat image available from the US Geological Survey.(b) A map showing the main hydrogeological and aquifer units (modified after Sharaky et al., 2008).

Figure 2 .
Figure 2. A digital elevation model of the Nubariya depression and its surroundings derived from the ASTER GDEM (Advanced Spaceborne Thermal Emission and Reflection Radiometer; Global Digital Elevation Model) (Tachikawa et al., 2011) and calibrated by a kinematic DGPS in situ survey (red alignment).The boundaries of the regional and local models developed in the framework of the IMPROWARE project are shown by the black and white boxes, respectively.The blue lines and associated dates, which separate the desert (to the southwest) from the irrigated areas (to the northeast), indicate how the latter has encroached southwestward over time.This information was obtained form an analysis of Landsat images.The results from the 3-D FE pond-scale groundwater flow model are shown along the black alignment E-E.

Figure 3 .
Figure 3. Landsat images showing the time evolution of the farmland (blue lines) and the pond in Nubariya between 1984 and 2014.The extent of desert reclamation is illustrated with reference to the red line, which marks the contemporary (2014) extent of the irrigated area.The Landsat images are available from the US Geological Survey.

Figure 4 .
Figure 4.A photo of the Nubariya pond dated November 2013 showing houses and a road abandoned due to the pond expansion.The photo location and direction are shown in the inset.

Figure 5 .
Figure 5.The location map of the boreholes (red triangles) and the TEM (blue dots) and MRS (green squares) measurements carried out in the surrounding area of the Nubariya pond.The black line represents the trace of the cross section through the resistivity model shown in Fig. 6.The location of MRS03 is highlighted.The domain of the local-scale model is shown by the yellow box.The background is a Landsat image acquired in 2014 available from the US Geological Survey.

Figure 6 .
Figure 6.A west-to-east section of the resistivity model obtained from the TEM data along the black line in Fig. 5.The color is faded white below the depth of investigation (gray line).The bars on the plots represent three MRS results in the vicinity of the profile in terms of the free water content (top) and relaxation times (bottom).Figure 5 provides the location of the MRS soundings.

Figure 7 .
Figure 7.The inversion results for sounding MRS03 (see Fig. 5 for the location of the sounding).The MRS and TEM data were inverted jointly.Row one shows the MRS observed (column one) and simulated (column two) data and the weighted data residuals (column three).Column four shows the TEM data fit.Row two shows the inversion results (black lines) in terms of resistivity, free water content, and relaxation time (columns one to three, respectively) together with the model parameter uncertainties (gray error bars).Column four shows the lithological information obtained from a nearby borehole.

Figure 8 .
Figure 8.The 3-D mesh of the regional hydrogeological model.A 2-D triangular grid was projected vertically to generate the 3-D FE tetrahedral mesh used in the CATHY simulations.

Figure 9 .
Figure 9.The evolution of the saturation degree along the vertical cross section E-E traced in Fig. 2, i.e., along the main groundwater flow direction as obtained by the 3-D depression-scale model.

Figure 10 .
Figure 10.The contour lines of the saturation degree at the land surface around the depression in (a) 2007, (b) 2010, and (c) 2014 as obtained by the 3-D pond-scale model.The backgrounds are Landsat images acquired in 2007, 2010, and 2014, respectively, which are available from the US Geological Survey.

Figure 11 .
Figure 11.Elevation maps (m a.m.s.l.) of (a) the top of the clay unit, (b) the bottom of the clay unit, i.e., the top of the limestone confined aquifer, and (c) the bottom of the limestone aquifer as provided by the hydrogeophysical (deep TEM) surveys and implemented in the local FE model.The background is a Landsat image acquired in 2014 that is available from the US Geological Survey.

Figure 12 .
Figure 12.The 3-D FE grid used in SUTRA to simulate possible scenarios of artificial aquifer recharge at Nubariya: (a) the horizontal view, (b) perspective view, and (c) vertical cross section along the A-A alignment traced in (a).

Figure 13 .
Figure 13.The salt concentration (g L −1 ) on a horizontal plane at −30 m a.m.s.l.(i.e., at the top of the injected aquifer) after 2 years of aquifer recharge for the various scenarios addressed by the study.The location of the injection wells are shown by the white dots.The production well in scenario S 3A is marked by a yellow X.

Figure 14 .
Figure 14.Scenario S 3A showing the evolution of the salt concentration (g L −1 ) on a vertical section through the barycenter of the injection configuration.The vertical exaggeration is 10.

Figure 15 .
Figure 15.The salt concentration versus time at the pumping well intake for the two selected values of Q withdrawn .

Table 1 .
Chemical analyses of groundwater samples collected in 2003 (after Sharaky et al., 2008) and 2011 (after Masoud, 2014) and groundwater and pond water samples collected in 2014 within the IMPROWARE project * .
Table 3 summarizes the

Table 3 .
The MRS measurement configuration used in this study.