Quantification of Recharge and Runoff from Rainfall Using New GIS Tool: Example of the Gaza Strip Aquifer

The Gaza Strip forms a transition zone between the semi-humid coastal zone in the north, the semi-arid zone in the east, and the Sinai desert in the south. Groundwater is the only water source for 1.94 million inhabitants, where the only fresh replenishment water for the aquifer comes from rainfall. This study focuses on testing a newly developed GIS tool to estimate the spatial and temporal distribution of runoff and recharge from rainfall. The estimation of surface runoff was made using the Soil Conservation Services Curve Number Method, while groundwater recharge was estimated using Thornthwaite and Mather’s Soil Moisture Balance approach. The new tool was applied to the Gaza aquifer for the year 1935 and for the period from 1973 to 2016. A comparison was made between the results obtained with the developed GIS tool and the frequently used Thiessen polygon method for rainfall distribution. Runoff and recharge were estimated for the year 1935 (prior to development) to compare with the current developed conditions. It was found that the built-up and sand dune areas stand in an inverse relationship, where the former is replacing the latter (built-up area expanded from 30.1 km2 in 1982 to 92.1 km2 in 2010). Recharge takes place in the sand dune area, whereas runoff increases in the built-up area. Due to development, runoff almost tripled from 9 million m3 in 1982 to 22.9 million m3 in 2010, while groundwater recharge was reduced from 27.3 million m3 in 1982 to 23 million m3 in 2010, even though the rainfall increased between 1982 and 2010 by 11%. Comparison between the newly developed GIS tool and the Thiessen polygon-based estimation shows that the former leads to higher values of runoff and recharge for dry years.


Introduction
The Gaza Strip is a coastal area along the eastern Mediterranean Sea and lies at latitude 31 • 25 59 N and longitude 34 • 22 34 E. It forms a transition zone between the semi-humid coastal zone in the north, the semi-arid zone in the east, and the Sinai desert in the south (Figure 1). The Gaza Strip has an area of 365 km 2 , where more than 1.94 million inhabitants are living (5315 persons per km 2 ) [1]. Groundwater is the only water source for the Gaza people, where the only fresh replenishment water for the aquifer comes from rainfall. Other recharge components are not fresh, including agricultural return flow, eastern groundwater inflow (lateral flow), water and waste water network leakage, and leakage from partially treated wastewater retention basins before dumping into the Mediterranean Sea.
The Gaza coastal aquifer consists of the Pleistocene age Kurkar Group [2] and the recent (Holocene age) sand dunes. The Kurkar Group includes marine and eolian calcareous sandstone, reddish Table 1. Summary of studies carried out to estimate recharge from rainfall for the Gaza Strip *.
The objectives of this study are (1) to calculate the recharge and runoff from rainfall for the Gaza Strip based on actual collected historical data from 1973 to 2016 using a newly developed tool in GIS based on full spatial distribution; and (2) to compare the results with calculations carried out using Thiessen polygon method for rainfall distribution for the same period (1973-2014) by Mushtaha et al. [5].
The developed GIS tool makes it possible to calculate groundwater recharge from rainfall by Thornthwaite and Mather's [16] Soil Moisture Balance Method using daily records of rainfall and evaporation, and rainfall runoff calculated according to Soil Conservation Services Curve Number (SCS-CN) method [17] using Geospatial Hydrologic Modeling System (Geo-HMS) [18] based on the land use map, hydrological soil map, digital elevation map, and SCS reference table incorporated in ArcMap extension.   [19,20]).

Materials and Methods
A new GIS tool was developed to estimate spatial and temporal distribution of runoff and recharge using the CN and SMB methods, respectively.
The SCS-Curve Number (CN) approach [17] was used for runoff estimation, while Thornthwaite and Mather's Soil Moisture Balance (SMB) method [16] was used for recharge estimation. These two methods were translated for use by the ESRI-GIS interface [21] using Model Builder and the Python language. The runoff and recharge calculations were performed from the year 1973 to the year 2016 (43 years). Since human activities changed over time, four land use maps were suggested for use by Mushtaha et al. [5]. This study used the same land use maps, in addition to the land use map for the year 1935, which can be considered representative of conditions prior to the strong development of the area.
The Inverse Distance Weighting (IDW) method was used to interpolate the point data within the GIS interface, because it gives greater weights to the points closest to the location of observation, and the weights diminish as a function of distance. IDW has been shown to be an appropriate method for daily rainfall interpolation [22].

Curve Number (CN) for Estimation of Runoff from Rainfall
The curve number (CN-II) estimation [17] was based on a combination of the land use maps, hydrologic soil map ( Figure S1 in Supplementary Materials), digital elevation model (DEM) map, and SCS reference table (Table 2). Geo-HMS Hydrologic Modeling System [18] was used to produce five CN-II GRID maps based on the five land use maps (1935, 1982, 1994, 2004 and 2010). The GRID map cell size was set to 100 m × 100 m. Each cell has a unique CN-II value. Adjustments to CN-I for  [19,20]).

Materials and Methods
A new GIS tool was developed to estimate spatial and temporal distribution of runoff and recharge using the CN and SMB methods, respectively.
The SCS-Curve Number (CN) approach [17] was used for runoff estimation, while Thornthwaite and Mather's Soil Moisture Balance (SMB) method [16] was used for recharge estimation. These two methods were translated for use by the ESRI-GIS interface [21] using Model Builder and the Python language. The runoff and recharge calculations were performed from the year 1973 to the year 2016 (43 years). Since human activities changed over time, four land use maps were suggested for use by Mushtaha et al. [5]. This study used the same land use maps, in addition to the land use map for the year 1935, which can be considered representative of conditions prior to the strong development of the area.
The Inverse Distance Weighting (IDW) method was used to interpolate the point data within the GIS interface, because it gives greater weights to the points closest to the location of observation, and the weights diminish as a function of distance. IDW has been shown to be an appropriate method for daily rainfall interpolation [22].

Curve Number (CN) for Estimation of Runoff from Rainfall
The curve number (CN-II) estimation [17] was based on a combination of the land use maps, hydrologic soil map ( Figure S1 in Supplementary Materials), digital elevation model (DEM) map, and SCS reference table (Table 2). Geo-HMS Hydrologic Modeling System [18] was used to produce five CN-II GRID maps based on the five land use maps (1935, 1982, 1994, 2004 and 2010). The GRID map cell size was set to 100 m × 100 m. Each cell has a unique CN-II value. Adjustments to CN-I for dry conditions and CN-III for wet conditions were calculated according to equations given by Chow et al. [23], based on Antecedent Moisture Conditions (Table 3). GRID maps for CN-I, CN-II and CN-III were calculated for all land use maps (3 × 5 GRID maps). The CN value was used to calculate runoff.

Soil Moisture Balance (SMB) for Estimation of Groundwater Recharge from Rainfall
The Soil Moisture Balance method was used to calculate the recharge quantity to the groundwater after the plants and crops have taken their needs from water in the root zone and after reaching soil field capacity. Food and Agriculture Organization (FAO) [24] determined the specific water content at saturation, the field capacity, and the permanent wilting point, depending on soil type. The effective precipitation was obtained based on the runoff calculations that were performed (actual precipitation minus runoff depth) [5].
Plant Available Water (PAW) is one of the main parameters that need to be found in order to estimate recharge using Thornthwaite and Mather's Soil Moisture Balance (SMB) method [16]. PAW (in mm) is a function of field capacity (FC in %), permanent wilting point (PWP in %) and plant root zone depth (Zr in mm): PAW = 10 (FC − PWP) × Zr [24]. Five PAW GRID maps were produced based on the different land use maps. FC and PWP were estimated using the FAO [24] reference table.

Rainfall and Evaporation
Rainfall for the year 1935 was taken from an average annual rainfall map , as shown in Figure 2. Twelve rainfall stations are distributed in the study area ( Figure S2 in Supplementary Materials), with daily rainfall data since 1973. Daily pan evaporation data from 1986 to 1992 are available, while average daily evaporation was estimated before 1986 and after 1992. Daily rainfall and evaporation data were prepared from 1973 to 2016 as shape files on a yearly basis. The shape file attribute data include day of the year, rain gauge station coordinates (X and Y), daily rainfall (in mm), daily evaporation (in mm), and 5-day antecedent rainfall (in mm).

Land Use
Five land use maps (1935, 1982, 1994, 2004 and 2010) were prepared to estimate CN, in addition to the DEM, with a grid size of 25 m × 25 m. Mushtaha et al. In Reference [5], the land use map for 1982 was developed according to the spatial distribution of the population based on monitoring population growth since 1973, in addition to the trend of land use activities from existing land use maps (1994, 2004 and 2010). The land use map for the year 1935 was conceived as being for a pre-developed area with 72,000 inhabitants [26]. Figure 3 shows the new 1935 land use map, and Figure S3 in the Supplementary Materials shows the four different land use maps.

Rainfall and Evaporation
Rainfall for the year 1935 was taken from an average annual rainfall map , as shown in Figure 2. Twelve rainfall stations are distributed in the study area ( Figure S2 in Supplementary Materials), with daily rainfall data since 1973. Daily pan evaporation data from 1986 to 1992 are available, while average daily evaporation was estimated before 1986 and after 1992. Daily rainfall and evaporation data were prepared from 1973 to 2016 as shape files on a yearly basis. The shape file attribute data include day of the year, rain gauge station coordinates (X and Y), daily rainfall (in mm), daily evaporation (in mm), and 5-day antecedent rainfall (in mm).

Land Use
Five land use maps (1935, 1982, 1994, 2004 and 2010) were prepared to estimate CN, in addition to the DEM, with a grid size of 25 m × 25 m. Mushtaha et al. In Reference [5], the land use map for 1982 was developed according to the spatial distribution of the population based on monitoring population growth since 1973, in addition to the trend of land use activities from existing land use maps (1994, 2004 and 2010). The land use map for the year 1935 was conceived as being for a predeveloped area with 72,000 inhabitants [26]. Figure 3 shows the new 1935 land use map, and Figure  S3 in the Supplementary Materials shows the four different land use maps.

Development of GIS Tool
Pre-defined GRID maps (CN-I, CN-II, CN-III and PAW) have been called during the calculations according to the condition statement. Runoff and recharge equations were translated into GIS as a loop to perform calculations on daily basis for each cell. Firstly, rainfall, evapotranspiration (PET) and 5-day antecedent rainfall (AMC) were interpolated using the Inverse Distance Weighting (IDW) method in GIS using the same grid size (100 m × 100 m). Potential maximum retention (S) was produced based on a conditional statement to select the proper CN value according to Equation (1) and Table 3. The runoff (Q in mm) GRID map was produced based on the daily rainfall (P in mm) in each cell using Equation (2).
Recharge calculation was performed based on three main components: (1) calculation of the Accumulated Potential Water Loss (APWL) based on the APWL of the previous time step, and augmented by the deficit of the new time step; (2) calculation of the effective rainfall (Peff in mm = P -Q); and (3) calculation of soil moisture (SM) based on APWL. Recharge is calculated as (Peff − PET) after the soil has reached field capacity. Initial pre-defined GRID maps for APWL and SM are zero and PAW value respectively, and they are changed according to the calculations at every loop based on Equations (3) and (4) in the Thornthwaite and Mather method [25]. Figure S4 in the Supplementary Materials shows a screenshot of the newly developed GIS tool for runoff-recharge estimations.
In the wet season (Peff > PET), SM increases by (Peff − PET) up to the soil FC and thus SM is known; APWL can be found using Equation (3). Water will percolate to the saturated zone producing groundwater recharge and AET will be equal to PET.

Development of GIS Tool
Pre-defined GRID maps (CN-I, CN-II, CN-III and PAW) have been called during the calculations according to the condition statement. Runoff and recharge equations were translated into GIS as a loop to perform calculations on daily basis for each cell. Firstly, rainfall, evapotranspiration (PET) and 5-day antecedent rainfall (AMC) were interpolated using the Inverse Distance Weighting (IDW) method in GIS using the same grid size (100 m × 100 m). Potential maximum retention (S) was produced based on a conditional statement to select the proper CN value according to Equation (1) and Table 3. The runoff (Q in mm) GRID map was produced based on the daily rainfall (P in mm) in each cell using Equation (2).
Recharge calculation was performed based on three main components: (1) calculation of the Accumulated Potential Water Loss (APWL) based on the APWL of the previous time step, and augmented by the deficit of the new time step; (2) calculation of the effective rainfall (P eff in mm = P − Q); and (3) calculation of soil moisture (SM) based on APWL. Recharge is calculated as (P eff − PET) after the soil has reached field capacity. Initial pre-defined GRID maps for APWL and SM are zero and PAW value respectively, and they are changed according to the calculations at every loop based on Equations (3) and (4) in the Thornthwaite and Mather method [25]. Figure S4 in the Supplementary Materials shows a screenshot of the newly developed GIS tool for runoff-recharge estimations.
In the wet season (P eff > PET), SM increases by (P eff − PET) up to the soil FC and thus SM is known; APWL can be found using Equation (3). Water will percolate to the saturated zone producing groundwater recharge and AET will be equal to PET.
In the dry season (P eff < PET), APWL increases by (PET − P eff ) and thus APWL is known; SM is calculated using Equation (4). AET = P eff + ∆SM, and no groundwater recharge occurs.
where Potential Evapotranspiration (PET) in mm: Effective rainfall (P eff ) in mm, Soil Moisture (SM) in mm, Accumulated Potential Water Loss (APWL) in mm, Actual Evapotranspiration (AET) in mm.

Results
In the year 1935, the average rainfall was 112.8 million m 3 In the dry season (Peff < PET), APWL increases by (PET − Peff) and thus APWL is known; SM is calculated using Equation (4). AET = Peff + ∆SM, and no groundwater recharge occurs.
where Potential Evapotranspiration (PET) in mm: Effective rainfall (Peff) in mm, Soil Moisture (SM) in mm, Accumulated Potential Water Loss (APWL) in mm, Actual Evapotranspiration (AET) in mm.

Results
In  The long-term average recharge is 24.5 million m 3 /year (19.2% of the rainfall). Recharge strongly varies with rainfall and intensity, but the average over the considered periods shows a tendency to decrease (Figure 4). Figure 5 shows the long-term average spatial distribution of runoff and recharge (1973-2016) as produced from the GIS tool. Meanwhile, Table 5 shows the tabulated results of rainfall, runoff and recharge results using GRID calculations per year.
Over the entire period, 24 years were below the long-term average rainfall and 19 years were above the long-term average rainfall. Annex I and Annex II in the Supplementary Materials show the 43-year GRID maps for runoff and recharge respectively.
Built-up areas increased three times from the year 1982 to 2010 (from 30.1 to 92.1 km 2 ), while the sand dune area decreased by more than 70% (114.8 to 31.5 km 2 ) [5]. This dramatic change in land use activities increased the runoff from 9 to 22.9 million m 3 , and recharge decreased from 27.3 to 23 The long-term average recharge is 24.5 million m 3 /year (19.2% of the rainfall). Recharge strongly varies with rainfall and intensity, but the average over the considered periods shows a tendency to decrease (Figure 4). Figure 5 shows the long-term average spatial distribution of runoff and recharge (1973-2016) as produced from the GIS tool. Meanwhile, Table 5 shows the tabulated results of rainfall, runoff and recharge results using GRID calculations per year. Over the entire period, 24 years were below the long-term average rainfall and 19 years were above the long-term average rainfall. Annex I and Annex II in the Supplementary Materials show the 43-year GRID maps for runoff and recharge respectively.   Built-up areas increased three times from the year 1982 to 2010 (from 30.1 to 92.1 km 2 ), while the sand dune area decreased by more than 70% (114.8 to 31.5 km 2 ) [5]. This dramatic change in land use activities increased the runoff from 9 to 22.9 million m 3 , and recharge decreased from 27.3 to 23 million m 3 during the same land use evolution period ( Figure 6).

Discussion
In

Discussion
In 1935, under pre-developed conditions in the Gaza Strip, runoff was very low, at 1.3 million m 3 compared to the runoff calculations in the period from 1973 to 2016. The results of the analysis (43 years) show that the long-term average rainfall was 119 million m 3 over the surface area of the Gaza Strip (326 mm/year), leading to long-term average runoff of 13.4 million m 3 and long-term average recharge of 24.5 million m 3 .
The changes in runoff were due to expansion in built-up area at the expense of sand dune area [5], as shown in Figure 6, with an overall trend of increasing runoff by more than 1600% (from 1.3 million m 3 to 22.9 million m 3 ), while recharge decreased by more than 38% (from 37.1 million m 3 to 23 million m 3 ) over the entire period ( Figure 4 and Table 5). Meanwhile, the runoff volume changed from one year to another within the same land use group (land use maps  . Long rainfall periods resulted in high recharge volumes ranging from 50% to 150% greater than the long-term average recharge volume (for example, the year 1982/1983, which had 52 rainy days). Annex II in the Supplementary Materials shows the yearly spatial distribution of recharge for the past 43 years (1973 to 2016). The analysis shows that when the rainfall duration was greater than 30 days and the rainfall quantity was higher than the long-term average of 326 mm, the recharge and runoff quantities were higher than the long-term average (for example, the year 1982/1983), see Figure 4 and Table 5.
On the other hand, in the year 2013/2014, more than 70% of the rainfall fell in five days (10-14 December 2013) and produced a runoff quantity of 60.2 million m 3 , exceeding the recharge quantity of 31 million m 3 . This runoff created flooding everywhere in the Gaza Strip.
The sand dune area decreased over time, as shown in the different land use maps. The change in the sand dune area from 114.8 km 2 (31.46%) of the total Gaza Strip area in 1982 to 92.9 km 2 (25.45%) in 1994 slightly increased runoff by 18% (from 9 to 10.6 million m 3 /year), while between 1994 and 2010, the sand dune areas decreased by 66% (92.9 km 2 to 31.5 km 2 ), leading to an increase in runoff by 116% (from 10.6 to 22.9 million m 3 ). Annex I shows the yearly spatial distribution maps for runoff.
A comparison of the results obtained from the developed GIS GRID method and those from the Thiessen polygon method [5] for the same period (1973-2014) is found in Table 6. The comparison shows that the overall long-term average rainfall and recharge over the area did not change, while a wide variation occurred in the long-term average of runoff (42.6% difference) due to the difference in spatial rainfall estimation and the influence of parameters such as soil texture and changes in land use (Figures 7 and 8). With the Thiessen polygon method, average values of PAW and CN were calculated for each polygon area (ranging in size from 3 km 2 to 82.5 km 2 ), while with the GRID method, the average was taken for every area of 100 m × 100 m. On the other hand, comparing the long-term average recharge with previous studies (see Table 1) shows a lower value than in previous studies, which is due to either the limited data used in previous studies and/or the use of assumptions, such as a soil coefficient for recharge. The differences between the two methods (GIS GRID vs. Thiessen polygon method) show a wide range for both yearly runoff and recharge. The largest differences occur in the dry years, where rainfall is lower than the long-term average (326 mm), as shown, for example, in the year 1975/1976. Figures S7 and S8 in the Supplementary Materials show a comparison between the two methods for selected years of lower and higher than long-term average rainfall, respectively. In general, it is indicated that the GIS GRID method provides a higher estimation than the Thiessen polygon method for both runoff and recharge when the rainfall is lower than the long-term average rainfall. In years with rainfall that was higher than the long-term average rainfall, the two methods' estimations are similar.

Conclusions
This study uses spatial daily rainfall and evapotranspiration variations for the past 43 years to estimate the runoff and recharge in the whole Gaza Strip. The long-term average (43 years) rainfall is 326 mm, average runoff is 13.4 × 10 6 m 3 /year (9.6% of rainfall), and average recharge is 24.5 × 10 6 m 3 /year (19.2% of rainfall). These results were obtained using a newly developed GIS tool based on daily records of rainfall and evapotranspiration with pre-defined GRID maps for topography, soil texture and land use.
The results obtained from the new GIS tool were compared with the Thiessen polygon method for the same period. This shows that the Thiessen polygon method is not accurate in estimating recharge or runoff because it does not provide a meaningful spatial distribution of rainfall, and because land use and soil texture, in reality, are changing within the Thiessen polygons. The GIS tool results in higher values of runoff and recharge for low rainfall years compared to the Thiessen polygon method.
The main factor of increasing runoff is the evolution of built-up area, where before 1992, the total built-up area was 30.1 km 2 of the Gaza Strip and the average runoff (1973 to 1992) was 7.9% (9 million m 3 ) of the rainfall; after development of the area, the total built-up area was 92.1 km 2 of the Gaza Strip, and the runoff had increased to 18.1% (22.9 million m 3 ) of the rainfall.
On the other hand, the decreasing size of the sand dune area affects the recharge quantities, where before 1973, the sand dune area was 114.8 km 2 of the Gaza Strip, and the average recharge (1973 to 1992) was 23.9% (27.3 million m 3 ) of the rainfall; after development of the area, the total sand dune area was 31.5 km 2 of the Gaza Strip, and the recharge had decreased to 18.2% (23 million m 3 ) of the rainfall.
Water planners and decision-makers in the Gaza Strip will benefit from the collected historical data and results of this research in securing runoff and recharge spatial distribution and volumes when making a proper plan to beneficially use the runoff to recharged the groundwater, while at the same time identifying vulnerable areas where most recharge takes place, in order to avoid changing these areas into built-up areas.