Potential Use of Treated Wastewater as Groundwater Recharge Using GIS Techniques and Modeling Tools in Dhuleil-Halabat Well-Field/Jordan

: Due to limited rainfall and precipitations, different developing countries depend on groundwater (G.W.) resources to challenge water scarcity. This practice of continuous and excessive G.W. pumping has led to severe water shortages and deteriorated water quality in different countries. Recharging of treated wastewater (TWW) into G.W. provides a critical solution for solving water scarcity, extending the well’s service life, and maintaining the G.W. supply. However, effective injection practice requires accurate tools and methods to determine the best location for groundwater recharge (GWRC). This work offers a new tool based on GIS–Multi-Criteria Analysis to identify the potential site and locations for GWRC with TWW. The developed methodology was applied to one of the most used well-ﬁeld areas in Jordan (Dhuleil-Halabat). The G.W. ﬂow for the B-B2/A7 formation system in the area of study was simulated using Processing Modﬂow (version 8.0). The analysis combined six thematic maps produced following the environmental, technical, and economic criteria to draw conclusions and recommendations. Both steady and transient conditions were used to predict the future changes that might occur under different stresses and after continuous GWR. The study evaluated three possible scenarios of artiﬁcial GWRC to evaluate the process efﬁciency and determine the effect on the water table level. The results revealed that only 0.05% (0.14 Km 2 ) of the total surface area of 450 Km 2 is suitable for GWRC. A GWRC with TWW at a rate of 3.65 Mm 3 /year (MCMY) would provide a good G.W. table recovery to 39.68 m in the year 2025, maintain a steady-state water table ≥ of 50.77 m for up to six years, and secure water supply for future generations. The proposed methodology can be used as a useful tool that can be applied to regulate the GWRC practice worldwide.


Introduction
Groundwater (G.W.) has been used as the main water supply in arid and semi-arid areas of highly populated countries. However, the continuous pumping at a rate greater than the natural recharge significantly reduced the water table and decreased the availability of water. Continuous monitoring showed that the low annual precipitation slows down the natural recharge of these aquifers. Artificial groundwater recharge (GWRC) was highlighted as unique and sustainable G.W. management to preserve the water table and extend the service life of the aquifers [1,2]. In this context, the managed aquifer recharge (MAR) is a planned water recharge for subsequent recovery and maintaining G.W. supplies [3][4][5][6][7]. Several techniques were used in MAR, including injection wells, aquifer storage and recovery (ASR), and infiltration basins. Zhang and Kanyerere [7] presented a review of the most current MAR practices. The review evaluates the importance of MAR in maintaining the service life of aquifers. It was highlighted in this review that the concept of MAR might have stemmed over two thousand years. This review also identified the common problems and challenges that render MARs from systematic clogging mechanism and prevention to theoretical seepage and infiltration calculations as well as purification mechanism and the use of big data and artificial intelligence. Bouwer [8] presented a work that focuses on the design of a system for artificial GWRC considering infiltration rates and the unsaturated zone between the land surface and the aquifer to establish adequate permeability and fast treatment. Asano and Cotruvo [9] presented a review that discusses the GWRC recharge and its management with a focus on the health and regulatory aspects of this practice using reclaimed municipal wastewater. The review shows some uncertainties concerning health risk considerations when reclaimed municipal wastewater is used for GWRC, especially when a large portion of the wastewater is recharged into the G.W. In the present work, the term MAR is referred to as artificial GWRC through recharging the basins with TWW.
GWRC only occurs when a sufficient amount of water is available to flow underneath the ground and infiltrate into the saturated zone [8,10,11]. The selection of suitable sites for MAR is dependent upon several factors, including occurrence and movement of G.W., hydrogeological factors of the recharge area (structure, slope, drainage, pattern, landform, land use/land cover), and climate [12][13][14][15]. Thus, it is important to consider these factors while exploring the potential for GWRC.
The selection of the most appropriate locations for GWRC with TWW depends on different criteria (environmental, technical, and economic) as outlined in Table 1. The environmental criteria include the residential areas as well as the water supply resources. In contrast, the technical criteria deal with soil characteristics, aquifer depth, land use, and slope, as well as the soil type and texture. The economic criteria, on the other hand, include the cost of pipes and pumping devices. Spatial analysis has been introduced as a powerful tool for identifying the suitable and ideal areas for the development and application of MAR. The site selection, in general, can be determined by integrating and weighing the geospatial data characterizing the study area according to the aims of the research. A geographic information system-multicriteria decision analysis (GIS-MCDA) has been developed for different settings and scales around the world to support the selection of MAR sites [15][16][17][18][19][20][21][22][23]. The geophysical imaging approach was widely used for assessing, improving, and identifying the highest permeability areas that are appropriate for aquifer artificial recharge in MAR sites [24,25]. Table 1. Environmental, technical, and economic criteria used for the selection of the appropriate locations for GWRC (adapted from [26]).

Environmental
The limitations set to prevent the pollution by infiltration of treated wastewater and direct interaction with human settlements: the safeguard distance was set at: 200 m away from urban residential areas; 50 m away from irrigation water resources; 100 m away from human drinking water supply sources.

Technical
The slope of the recharging basins should be between 0 and 12 percent. Higher slopes would increase runoff, soil erosion, soil instability, and replenishment costs [27]. Soil texture: the recharging basins should be situated in permeable soils to allow for high infiltration rates and suitable wastewater treatment [27][28][29][30]. The unsaturated zone should be greater than 5 m to avoid G.W. contamination [29]. Low permeability is highly recommended in this zone to restrict the downward movement of water and clog the soil. Soil Types: Rocky soil should be avoided except for Anthrosol soils because at depths ≥ 1 m it permits infiltration at a rate of 1 m/s. The remaining contaminants are largely eliminated mainly within the top one meter [31]. Aquifer depth: A G.W. level of at least five meters below the ground surface is recommended to provide an adequate distance for TWW purifications [28,[30][31][32]. Land Use: All uncultivated zones (e.g., bush area) have been considered for GWRC.

Economic
The cost includes the piping system with a maximum distance of 8 km as well as transporting of the TWW to the recharging site [28,30,33].
Artificial GWRC with TWW is becoming increasingly important as a G.W. management alternative. Various suitability criteria have been presented throughout the research literature [26,27,30,[32][33][34][35][36][37][38][39]. For example, high-quality TWW is a steadily and widely available natural resource that can be used in GWRC as a vital alternative to sustain G.W. potential, resource conservation, and management of large volumes of wastewater. During this process, the TWW is infiltrated through the unsaturated zone, where it undergoes a natural purification process and maintains the supply of high-quality G.W. Previous works have highlighted different challenges for GWRC practice, including the effect of recharged volume of TWW as well as water quality parameters on the quality of the produced G.W. [9,[40][41][42]. Although different research works have shown that obtaining a direct measurement of the precise amount of the TWW to be recharged is nearly impossible, other research works have suggested the use of indirect methods or the development of sophisticated models/tools that can be used for such objective [21,31].
The G.W. modeling techniques have become appropriate tools to simulate the behavior of G.W. systems and to achieve optimum management and uses. Besides this, these modeling techniques play an important role in assessing the MAR scenarios and screening the potential GWRC locations. Ringleb et al. [5] had collected and evaluated 216 studies published over the past 30 years addressing the importance of G.W. modeling in MAR practical applications.
The aforementioned literature indicated that GWRC could be amongst the possible solution to resolves water scarcity and secure clean, safe, and steady water resources for the coming generation. There is limited research on the use of indirect methods in determining the volume and location of TWW as a potential source for artificial GWRC. As such, there is a knowledge gap on the topics related to (1) the use of artificial GWRC with TWW for the management of G.W. aquifers in water scarcity countries, (2) the development of sophisticated models and tools to determine the required volume of GWRC, and (3) utilizing restricted water resources efficiently and allocating the G.W. resources properly.
Given that, the present study aims to investigate the use of GIS-based multi-criteria analysis to identify potential site locations for GWRC with TWW. This case study was carried in the Dhuleil-Halabat well-field in Jordan (30.5852 • N, 36.2384 • E), which sits on the eastern highland of the country inside the Amman Zarqa Basin (AZB). Jordan is one of the water-poorest countries in the world due to the limited annual rainfall amount, which has shown high fluctuation from year to year. On the other hand, the high temperatures during the summer have increased the evaporation rate from surface water, creating big stress on the regularly available water. Presently, the available per capita freshwater in Jordan is around 85 L/day, which is significantly low compared with the international standards and other countries. According to the Water Stress Index (WSI)-defined as the value of annual rainfall divided by the total population (m 3 /capita/year)-a value ≥1800 reflects the norm, while a value ≤1700 reflects existing stress, a value <1000 reflects scarcity, and a value <500 is described absolute scarcity. On the WSI, Jordan has a value of 450, which is classified as absolute scarcity [43]. It is predicted that the difference between water supply and demand will increase in the coming 15 years, which will increase the harsh daily realities and make it difficult to access safe water resources. This will certainly impact human needs, the country's improvement, and could hinder the development in the industrial, agricultural, and economic sectors. In the mid-1960s, the G.W. from the Dhuleil-Halabat well-field was extracted from six wells and increased to 113 wells in 2019. Previous research work related to Dhuleil-Halabat well-field has focused on the use of G.W. modeling, GWRC, sources of pollution, aquifer characterization, and determining the water quality parameters of the [44][45][46][47][48][49][50][51][52][53]. Due to over-pumping, the aquifer equilibrium was disrupted, and the volume of G.W. decreased, which also led to a decrease in the quality of the G.W. The model developed in this study can be used as a useful tool to control and manage the volumes and location of the TWW that could be used in the artificial GWRC of this area for economic and resource preservation purposes.  Figure 1 presents the study area, Dhuleil-Halabat well-field (260 to 290 East and 165 to 180 North) according to the Palestine Grid, which is situated to the northeast of the capital Amman and covers an area of approximately 450 km 2 . This area is mainly located in the AZB and stretches from Khirbet As-Samra in the west to Qaa' Khanna in the east [54]. The topography of the study area is irregular with elevated ranges between 580 m and 700 m above the mean sea level (AMSL) as indicated in Figure 2a.

Study Area and Climate Conditions
leil-Halabat well-field has focused on the use of G.W. modeling, GWRC, sources of pollution, aquifer characterization, and determining the water quality parameters of the [44][45][46][47][48][49][50][51][52][53]. Due to over-pumping, the aquifer equilibrium was disrupted, and the volume of G.W. decreased, which also led to a decrease in the quality of the G.W. The model developed in this study can be used as a useful tool to control and manage the volumes and location of the TWW that could be used in the artificial GWRC of this area for economic and resource preservation purposes. Figure 1 presents the study area, Dhuleil-Halabat well-field (260 to 290 East and 165 to 180 North) according to the Palestine Grid, which is situated to the northeast of the capital Amman and covers an area of approximately 450 km 2 . This area is mainly located in the AZB and stretches from Khirbet As-Samra in the west to Qaa' Khanna in the east [54]. The topography of the study area is irregular with elevated ranges between 580 m and 700 m above the mean sea level (AMSL) as indicated in Figure 2a.  The climate of this area is characterized by warm summers and cold winters (i.e., Sahara type), with average precipitation in the range of 120 mm/year in the southeastern part to 160 mm/year in the northwest area. The monthly rainfall distribution has the highest values during the rainy season in December (an average of 29 mm) and January (an average of 33 mm), and a sharp decrease occurs between March and May (4-20 mm). The monthly variation in temperature peak up in July (min 18 • C, max 44 • C, average 28 • C) and decrease in January (min 5 • C, max 23 • C, average 10 • C). The annual surface evaporation rate, based on a class-A-pan test, ranges from 2450 to 3850 mm/year with an average value of 2949 mm/year, varying between. The monthly evaporation rate ranges from 116 mm in January to 395 mm in July. Most of the rainfall occurs during very short intervals and forms flash flooding events. The volume of runoff in this area was estimated using the runoff coefficient for wet, dry, and average conditions to be 2.2 MCMY [53]. The climate of this area is characterized by warm summers and cold winters (i.e., Sahara type), with average precipitation in the range of 120 mm/year in the southeastern part to 160 mm/year in the northwest area. The monthly rainfall distribution has the highest values during the rainy season in December (an average of 29 mm) and January (an average of 33 mm), and a sharp decrease occurs between March and May (4-20 mm). The monthly variation in temperature peak up in July (min 18 °C, max 44 °C, average 28 °C) and decrease in January (min 5 °C, max 23 °C, average 10 °C). The annual surface

Properties of the Aquifer System
The main aquifers in this study are the Basalt aquifer and the underlying Amman-Wadi Es Sir aquifer (B2/A7). Both aquifers are part of a combined system (B-B2/A7 formation) as they are directly connected. The basalt aquifer is restricted to the eastern area and a channel crossing the study area is the central area in the west-east direction. Typically, basaltic aquifers are characterized by hydraulic anisotropy and discontinuous heterogeneity. The large variations consist in the hydraulic conductivity of basalt aquifer systems. However, the B2/A7 aquifer is the most important in Jordan because of its vast reaches and favorable aquifer properties (high karstification) as described by Margane et al. [55].
The G.W. flow to the area comes from Jabal Al Arab in Syria, which is diverted to the G.W. divide. This divide was previously recognized within the framework of the G.W. assessment of Northern Jordan [55]. The divide was explained as a local recharge mound. In 1965, the G.W. gradient in the area was likely 0.0015. An inflow component comes from the East Amman area, which reaches the study area south of the Halabat well-field. The G.W. flow in the area of study in 1965 is presented in Figure 2b. The current G.W. gradient in this area is 0.004. In 1965 the G.W. flow came from the south/southeast and had a gradient of 0.002. In the western portion of the study area, the current G.W. gradient is shallow, whereas, in 1965, the G.W. was flowing westward with a gradient of 0.0024. The previous observations show that, due to the heavy exploitation, the Dhuleil-Halabat G.W. gradients have vastly changed over the past four decades.
The basalt aquifer overlies the B2/A7 aquifer in the Wadi Dhuleil area where it is hydraulically connected to it. The thickness of the basalt aquifer is in the range of a few meters up to more than 200 m. Until roughly 40 years ago, the basalt flows reached the Wadi Zarqa where the Sukhneh springs discharged water from the basalt aquifer. The transmissivity of the basalt aquifer was determined by several wells [53]. The transmissivity ranges between 10 and 100 m 2 /d, while the hydraulic conductivity is between 0.2 and 31.1 m/d.
The B2/A7 aquifer extends over much of Jordan and has a high permeability, large storage capacity, and receives GWRC from precipitation. The thickness of the B2/A7 formations can reach up to more than 300 m in the Dhuleil-Halabat region. Numerous exploratory wells were drilled in the Dhuleil-Halabat area and were used to determine the hydraulic characteristics of the B2/A7 aquifer [53]. The transmissivity of the B2/A7 aquifer was determined by performing the pumping test for several wells, and the results indicate that the transmissivity ranges between 476 and 3400 m 2 /d, while the hydraulic conductivity is between 14.7 and 47.6 m/d [53].
According to the Water Authority, the amount of G.W. extraction from this well-field was around 4.5 MCM in 1965 and reached up to 27.5 MCM in 1994 [53]. From 1994 to 2007, no further wells were put into operation, and the extraction decreased to 15.5 MCM yearly. Thereafter, the extraction remained constant until 2019. The reduction of extraction in the subsequent years results from the high decline in water levels and less productivity of the wells. The soil types in the Dhuleil-Halabat area are mostly uniform, which can be divided based on soil texture into four groups: Sand (6.97%), Silty Clay (36.73%), Silty Loam (45.15%), and Clay (11.15%), as exemplified in Figure 2c. The land use of the study area is composed of 78.07% bare soil, including urbanized areas, 9.33% field crops, tree crops, and vegetables, and 5.17% of the study area is pastureland ( Figure 2d).  Figure 2a exhibits the slope map that was produced for this region as well as how the overall slope in the study area is changing from west to east, showing where hills cover a wide area in the west and around the study area.

Data Analysis
Each criterion thematic map was transformed into raster form with 30 m resolution using spatial analysis (Table 1). Each pixel (30 m × 30 m) in the selected raster thematic map was assigned a value (0 or 1) according to the suitability criteria. Those coded with 1 are suitable for recharging, and those coded with 0 are unsuitable for recharging. The raster calculator tool in ArcGIS 10.3 was used to perform a Boolean operating analysis between the grid cells for each thematic map to produce a final suitability map for the recharging basin.

Groundwater Flow Modeling
The processing Modflow (PM) version 8.0 was used to investigate the behavior and response of groundwater flow following GWRC with TWW [56]. The model was developed based on simplification or extraction of the complex physical reality and its processes. As the validity of the forecast data depends upon how close the model is to the actual field conditions, high-quality field data was used to validate the used model. The two aquifers in the Dhuleil-Halabat region (B-B2/A7 and Kurnub) were divided by the Lower Ajlun Aquitard (A1/6) layer, as shown in Figure 3. All the withdraws were assigned to B-B2/A7 formation, while the Kurnub aquifer in most parts of Dhuleil-Halabat was set off operation. Under such circumstances, the analysis is only conducted on the B-B2/A7 formation.    [57]. These values were compared to well-logs obtained from previously drilled wells throughout the study area. The B-B2/A7 formations have a thickness of approximately 516 m and 200 m, respectively.
Based on the net flow of the G.W. and the structural maps of the B-B2/A7 formation, the constant head boundary was assumed in the north-eastern, south-eastern, and western areas of the model domain with 530 m, 530 m, and 505 m water levels, respectively. Dry spots in the north and south of the study area were considered a no-flow boundary (Figure 4c). The initial condition for the current model was identified based on the water level contour lines in the equilibrium phase before the aquifer disturbance due to pumping.
Based on the available data, the water level in 1965 was chosen as the initial condition to calibrate the model in the steady-state condition (Figure 2b).
Based on previous research literature and pumping tests, the horizontal hydraulic conductivity (HHK) values in the B-B2/A7 formation was ranged between 0.043 and 250 m/day [53,58]. The HHK values were modified during the calibration of the model. The isohyetal map for the precipitation distribution in the study area for normal conditions was used for estimating the natural recharge, as shown in Figure 4d [56]. The normal year conditions were assumed and applied for all of the years of the calibration period. The average recharge percent was calculated to be 3.7% of the rainfall based on the standardized precipitation index approach (SPI) [59].    Based on the net flow of the G.W. and the structural maps of the B-B2/A7 formation, the constant head boundary was assumed in the north-eastern, south-eastern, and western areas of the model domain with 530 m, 530 m, and 505 m water levels, respectively. Dry spots in the north and south of the study area were considered a no-flow boundary (Figure 4c). The initial condition for the current model was identified based on the water level contour lines in the equilibrium phase before the aquifer disturbance due to pumping. Based on the available data, the water level in 1965 was chosen as the initial condition to calibrate the model in the steady-state condition (Figure 2b).
Based on previous research literature and pumping tests, the horizontal hydraulic conductivity (HHK) values in the B-B2/A7 formation was ranged between 0.043 and 250 m/day [53,58]. The HHK values were modified during the calibration of the model. The isohyetal map for the precipitation distribution in the study area for normal conditions was used for estimating the natural recharge, as shown in Figure 4d [56]. The normal year conditions were assumed and applied for all of the years of the calibration period. The average recharge percent was calculated to be 3.7% of the rainfall based on the standardized precipitation index approach (SPI) [59].

A. Steady-State Calibration
The model calibration was achieved through the trial and error method. The procedure includes adjusting the hydraulic parameters during the sequential runs and carrying out calculations until the computed water levels match the actual water level contours [60][61][62].

A. Steady-State Calibration
The model calibration was achieved through the trial and error method. The procedure includes adjusting the hydraulic parameters during the sequential runs and carrying out calculations until the computed water levels match the actual water level contours [60][61][62].

B. Transient Calibration and Validation
In the transient calibration, the simulated head and the HHK values that were obtained from steady-state calibration were used. The transient calibration was achieved by changing the specific yield and storage values. The transient calibration of the area was performed using the existing water-level data from two observation wells (AL1043 and AL2698) for the period 1989-2017. The available data for the period 1988-2008 was used for the model calibration, while the available data for the period 2009-2019 was used for the model validation.

C. Model Prediction
The transient model was used to perform the predictive simulations and the sustainability of the G.W. system. After the transient calibration, three different future prediction scenarios were developed up to the year 2035. The scenarios considered the water level change in the aquifer system as follows:

WWTP Characterization
As-Samra wastewater treatment plant (WWTP), which opened in 2005 with activated sludge technology, is the largest treatment processed in Jordan, serving 60% of the Jordanian population [62]. The design capacity of this process is 360,000 m 3 /day [62]. Table 2 presents the effluent water quality parameters from As-Samra WWTP. It can be seen that most of the water quality parameters comply with Jordanian standards for GWRC except for the Chemical Oxygen Demand (COD). With such good effluents characteristics, clogging in the infiltration basin is anticipated. However, this form of clogging is relatively simple to remedy with routine basin maintenance [63][64][65]. The excess COD can be effectively and reliably removed from treated water during the long-term operation of MAR systems (Infiltration Basin). Different researchers have proven that more than 50% of COD could be removed at the top 30 cm of the soil layers [65][66][67][68]. Table 2. Effluent water quality parameters from As-Samra WWTP.

Parameters
As  Figure 5 presents the six thematic maps categorized as water points (W.P.), urban areas (U.A.), land use (L.U.), soil type (T.S.), slopes (S), and WWTP distance (DWWTP) generated using ArcGIS 10.3 (spatial analysis tool). For each thematic map, the acceptable areas were coded as 1 and the rejected areas were coded as 0. Figure 5a,b present the environmental criteria which include a buffer area of 100 m radius for water consumption wells and a 200 m radius for the residential area, both of which were coded as 0. Figure 5c presents the slope thematic map for the study area, which was reclassified into two categories. The first category was related to the slopes with a value less than or equal to 12% (coded as 1), while the second category was related to the slopes with values >12% (coded as 0). The soil texture was identified as 1 for and 0 for the rest of the areas, as shown in Figure 5d. The suitable water table depth for GWRC must have a static G.W. level of at least 5 m below the ground to have a sufficient unsaturated zone for W.W. purification. Results showed all parts of the study area were considered suitable and coded as 1.

GIS Analysis
The land use thematic map of the studied area was produced considering all of the bush areas (i.e., uncultivated areas) as a suitable site for recharging and coded as 1 as presented in Figure 5e. The economic criteria include an 8 km buffer zone around the four WWTPs, which are coded as 1 as shown in Figure 5f. The raster GIS maps were generated using the raster calculator to create the suitability map for the GWRC basins ( Figure 6). Table 3 outlines the suitable area for GWRC and the percentage from total area based on each criterion. It was observed that water point is suitable for about 97.86% of the total land area. Up to 99.5% of the studied area,~267.29 km 2 has the required slope, 97.97% of the studied area~263.03 km 2 has the required urban contribution, and 97.86% of the area (262.76 km 2 ) has the right water points. However, only 55.56% of this land (≈149.19 km 2 ) has the required soil texture, and just 7.96% of the land (18.68 km 2 ) has a good opportunity in terms of close WWTP. wells and a 200 m radius for the residential area, both of which were coded as 0. Figure 5c presents the slope thematic map for the study area, which was reclassified into two categories. The first category was related to the slopes with a value less than or equal to 12% (coded as 1), while the second category was related to the slopes with values >12% (coded as 0). The soil texture was identified as 1 for and 0 for the rest of the areas, as shown in Figure 5d. The suitable water table depth for GWRC must have a static G.W. level of at least 5 m below the ground to have a sufficient unsaturated zone for W.W. purification. Results showed all parts of the study area were considered suitable and coded as 1. The land use thematic map of the studied area was produced considering all of the bush areas (i.e., uncultivated areas) as a suitable site for recharging and coded as 1 as presented in Figure 5e. The economic criteria include an 8 km buffer zone around the four WWTPs, which are coded as 1 as shown in Figure 5f. The raster GIS maps were generated using the raster calculator to create the suitability map for the GWRC basins  A GIS-based MCDA, combined with sensitivity analysis, informs decision-makers about which criteria are genuinely required and how the results are affected by weight changes. Consequently, it can maximize the criterion combination, adjust the conceptual framework, and minimize the deviation, leading to enhanced and more persuasive results, which will be considered in future research. Figure 7 presents the comparison between the actual and the computed water level contour map for the B-B2/A7 formation in the Dhuleil-Halabat region produced by Modflow. Previously reported data on the HHK [51,53,55,58] and pumping tests [69] were used as initial values in the steady-state simulation. After that, the HHK values were continuously updated during the calibration process of the model. Results from the calibration phase revealed that the HHK values for the study area range between 5 and 35 m/day. Low conductivity values of 9 and 5 m/day were noticed on the south-central and south-west sides of the area (see Figure 8).  A GIS-based MCDA, combined with sensitivity analysis, informs decision-makers about which criteria are genuinely required and how the results are affected by weight changes. Consequently, it can maximize the criterion combination, adjust the conceptual framework, and minimize the deviation, leading to enhanced and more persuasive results, which will be considered in future research. Figure 7 presents the comparison between the actual and the computed water level contour map for the B-B2/A7 formation in the Dhuleil-Halabat region produced by Modflow. Previously reported data on the HHK [51,53,55,58] and pumping tests [69] were used as initial values in the steady-state simulation. After that, the HHK values were continuously updated during the calibration process of the model. Results from the calibration phase revealed that the HHK values for the study area range between 5 and 35 m/day. Low conductivity values of 9 and 5 m/day were noticed on the south-central and south-west sides of the area (see Figure 8).     Conversely, high conductivity values~35 and 25 m/day were observed in the East and West regions, respectively (Figure 8). The water balance for the B-B2/A7 formation at the steady-state conditions (Figure 9) indicates that the inflow is greater than the outflow by 0.01 MCMY with a discrepancy of 0.19%. The calculated total recharge amount for the study area was determined to be 10.32 MCMY (9.11 MCMY) with inflow from Jebal Al Arab and Azraq basin (i.e., southern constant head boundary (Figure 4c), and 1.21 MCM from the areal recharge. In contrast, the outflow from the eastern boundary towards Wadi Zarqa was determined to be 10.31 MCMY. In 1994, when the pumping reached 27.5 MCMY, a large amount of water flowing to pumping wells comes from G.W. storage. This causes a rapid decrease in the G.W. levels and creates a great hydraulic gradient around pumping wells. So, the hydraulic gradient forces G.W. to flow towards the pumping wells and reduces G.W. discharge to the western boundary, and creates recharge from the eastern and southern constant head boundary. With extra pumping and increasing the hydraulic gradient, the capture (low discharge and high recharge) increased and eventually equaled the pumping rate of 15.5 MCM in around 35 years as determined from Modflow calculations. In this situation, the capture sustains the pumping rate, and no more G.W. storage is removed. A new equilibrium state is achieved at a sustainable yield of 15.5 MCMY instead of 10.31 MCMY as calculated by steady-state calibration processes. Conversely, high conductivity values ~35 and 25 m/day were observed in the East and West regions, respectively (Figure 8). The water balance for the B-B2/A7 formation at the steady-state conditions (Figure 9) indicates that the inflow is greater than the outflow by 0.01 MCMY with a discrepancy of 0.19%. The calculated total recharge amount for the study area was determined to be 10.32 MCMY (9.11 MCMY) with inflow from Jebal Al Arab and Azraq basin (i.e., southern constant head boundary (Figure 4c), and 1.21 MCM from the areal recharge. In contrast, the outflow from the eastern boundary towards Wadi Zarqa was determined to be 10.31 MCMY. In 1994, when the pumping reached 27.5 MCMY, a large amount of water flowing to pumping wells comes from G.W. storage. This causes a rapid decrease in the G.W. levels and creates a great hydraulic gradient around pumping wells. So, the hydraulic gradient forces G.W. to flow towards the pumping wells and reduces G.W. discharge to the western boundary, and creates recharge from the eastern and southern constant head boundary. With extra pumping and increasing the hydraulic gradient, the capture (low discharge and high recharge) increased and eventually equaled the pumping rate of 15.5 MCM in around 35 years as determined from Modflow calculations. In this situation, the capture sustains the pumping rate, and no more G.W. storage is removed. A new equilibrium state is achieved at a sustainable yield of 15.5 MCMY instead of 10.31 MCMY as calculated by steady-state calibration processes.

Transient Calibration and Validation
The transient calibration was achieved by changing the specific yield and storage values and perform several computer runs until the observed and simulated heads matched. Calculations determined that the average specific storage and average specific yields were 1.36 × 10 −4 m −1 and 0.025, respectively. Figures 10 and 11 compare the observed versus simulated fluctuations in G.W. heads in two observation well stations (AL1043 and AL2698). Results showed a good agreement between the observed and simulated fluctuations, confirming the accuracy of the developed model. There are very low differences between the measured and simulated G.W. heads in both the calibration and validation stages because of withdrawals in the model domain. It is clear that the water levels in the two observation wells are strongly influenced by the G.W. extraction from the surrounding pumping wells.

Model Evaluation
The calibration processes (steady-state and transient) between the measured and simulated heads were quantitatively assessed and are outlined in Table 4. The results in table 4 show that the developed model successfully predicts the measured values for

Model Evaluation
The calibration processes (steady-state and transient) between the measured and simulated heads were quantitatively assessed and are outlined in Table 4. The results in table 4 show that the developed model successfully predicts the measured values for

Model Evaluation
The calibration processes (steady-state and transient) between the measured and simulated heads were quantitatively assessed and are outlined in Table 4. The results in Table 4 show that the developed model successfully predicts the measured values for both wells. The statistical analysis showed that the maximum and minimum residual of the three measurements ranged between 0.765 to 4.799 and −0.815 to −4.683, respectively. The RMSE showed very low values (0.463 to 0.736) and R 2 values close to one (0.931 to 0.996). The quantitative measurements provide a different way of assessing the calibration results beyond the visualized results previously provided. Typically, these measurements indicate that the match between the measured and simulated heads for the calibration process is more accurate than the validation process. This result was anticipated as the calibration stage for the hydraulic conductivity specific yield values hanged until the best fit was reached. However, for the validation stage, the calibrated values from the calibration stage were used. The grid size of 500 m by 500 m was used for steady-state calibration as the spatial resolution of data, and the monthly temporal resolution of data was used for the calculation of performance metrics (Transient State).  Figure 12 presents the model prediction for the three studied scenarios at the observation well (AL1043). The model was used to predict and assess the sustainability of the future GWRC as per the proposed scenario up to the year 2035. For the first scenario, which assumed constant G.W. withdrawal for 2025, 2030, and 2035 without any artificial GWRC recharge, the results showed that the maximum drawdown would be concentrated in the Southern region, where many of the pumping wells exist. These wells would reach roughly 70.12, 75.17, and 78.28 m in the years 2025, 2030, and 2035, respectively. The model also predicted that in 2030 there would be some dry areas in the South-West region with no-flow (dry areas), which are the areas that had saturation thickness ≤ of 35 m in 1965. For the second scenario assume artificial GWRC at a rate of 1.825 Mm 3 /year while keeping the current extraction rates constant up to the year 2035. The results indicated that the maximum drawdown would replicate scenario 1, as the average distance between the recharging sites and the location of the maximum drawdowns is 20 km (see Figure 6). The results also showed that with such practice the G.W. levels at observation well AL1043 would increase by roughly 12.82, 25.63, and 28.07 m in 2025, 2030, and 2035, respectively. The obtained results suggest an interesting recharging rate that will preserve the G.W. supply and sustain water resources. The results of the third scenario (Artificial GWRC at a rate of 3.65 MCMY and constant extraction rate up to the year 2035) show that the maximum drawdown would replicate scenario 1, as the average distance between the recharging sites and the location of maximum drawdowns is 20 km. With this scenario, the  The obtained results (see Figure 12) suggest a good strategy for sustainable maintenance of the G.W. via recharge with TWW. This strategy has several advantages, including the management of excess volumes of TWW that cannot be used for irrigation or landscaping and sustainable GWRC to maintain resources for the future generation. Although scenarios 2 and 3 show promising results in maintaining the G.W. up to 2035, the effect of water quality parameters on the produced G.W. and the cost associated with this practice required further investigation. Studying the fate and transport of different constituents in the TWW as a function of time is of great interest and would be the focus of our research work in the future. Besides this, developing a well-defined procedure for the injection and distribution of the TWW within the recharge area still requires further investigation. It is also worth mentioning that, herein, the proposed model work focused on simulation in the upper aquifer. The lower layer and the unsaturated zone were not included in this study. It will be highly encouraged to develop simulation work to cover these areas and verify the results. The effect of the unsaturated zone on hydrologic problems is often ignored or treated with crude approximations due to the scarcity of data in this zone. During the GWRC process, the recharging water typically reaches an aquifer by flowing through an unsaturated zone overlying it. Due to localized recharge processes, analytical and numerical models for coupled unsaturated-saturated flow models have received little attention. It was assumed during the use of the GWRC package from MODFLOW software that all the infiltrated water is transmitted to the water table as recharge. With such an assumption, the recharged flow can locally increase the water table more than the case which considers the flow through the unsaturated zone.

Conclusions
A case study was developed to identify the potential site location for G.W. recharge in Dhuleil-Halabat in Jordan. Six thematic maps were produced following environmen- The obtained results (see Figure 12) suggest a good strategy for sustainable maintenance of the G.W. via recharge with TWW. This strategy has several advantages, including the management of excess volumes of TWW that cannot be used for irrigation or landscaping and sustainable GWRC to maintain resources for the future generation. Although scenarios 2 and 3 show promising results in maintaining the G.W. up to 2035, the effect of water quality parameters on the produced G.W. and the cost associated with this practice required further investigation. Studying the fate and transport of different constituents in the TWW as a function of time is of great interest and would be the focus of our research work in the future. Besides this, developing a well-defined procedure for the injection and distribution of the TWW within the recharge area still requires further investigation. It is also worth mentioning that, herein, the proposed model work focused on simulation in the upper aquifer. The lower layer and the unsaturated zone were not included in this study. It will be highly encouraged to develop simulation work to cover these areas and verify the results. The effect of the unsaturated zone on hydrologic problems is often ignored or treated with crude approximations due to the scarcity of data in this zone. During the GWRC process, the recharging water typically reaches an aquifer by flowing through an unsaturated zone overlying it. Due to localized recharge processes, analytical and numerical models for coupled unsaturated-saturated flow models have received little attention. It was assumed during the use of the GWRC package from MODFLOW software that all the infiltrated water is transmitted to the water table as recharge. With such an assumption, the recharged flow can locally increase the water table more than the case which considers the flow through the unsaturated zone.

Conclusions
A case study was developed to identify the potential site location for G.W. recharge in Dhuleil-Halabat in Jordan. Six thematic maps were produced following environmental, technical, and economic criteria. The G.W. recharge potential procedure is demonstrated using the GIS-based multicriteria, which include different G.W. recharge potential factors. Three-dimensional G.W. flow models were used to calculate the G.W. budget, aquifer characteristics, and predict the aquifer response to water levels for the coming 15 years (up to 2035). Out of the 450 Km 2 of the study area, the potential site for the G.W. recharge area was determined to be 0.14 km 2 of the Eastern region.
The horizontal hydraulic conductivities of the aquifer (B-B2/A7) ranged between 5 and 35 m/day. The distribution of horizontal conductivity is irregular due to the complexity of the aquifer. The total GWRC, inflow from the eastern boundary through the basalt of the B-B2/A7 of the aquifer, and infiltrations due to rainfall were 10.32 MCMY, 9.11 MCMY, and 1.21 MCMY, respectively. The outflow across the western boundary reaches up to 10.31 MCMY towards Wadi Zarqa. According to the transient state calibration, the specific storage and the specific yield of the B-B2/A7 aquifer were found to be 1.36 × 10 −4 m −1 and 2.5%, respectively. The optimal use of G.W. resources should not exceed 15.5 MCMY.
Lastly, three scenarios were performed to predict the aquifer system response to GWRC with TWW. The results exemplified that Artificial GWRC at a rate of 3.65 MCMY and constant extraction rate up to 2035 has the best results in terms of aquifer recharging, maintaining a steady-state G.W. level of 39.68 m extends the well-field service life, and manages high volumes of excess TWW. In contrast, withdrawal without a recharge is the least effective due to high G.W.