Groundwater quality distribution by geostatistical investigation ( GIS ) , Nile Delta , Northern Egypt

The Nile Delta represents two-thirds of agricultural output with low precipitation and population growth that depletes the groundwater. The Quaternary aquifer is consisted of Holocene (aquitard) underlain by Pleistocene (aquifer), and groundwater flows from south to north, northwest, and northeast directions. The total dissolved solids (TDS) concentration increases toward the western and northeastern part in 1997 by seawater intrusion and irrigation practices. The TDS concentration anomalies were in the north and east, and are attributed to pumping within 2003, 2010, and 2013 periods. There is a good match between hydrogeochemistry (TDS content) and hydrogeology (groundwater flow) with the TDS concentration change found to decline mainly in the southwestern and northeastern part within 1997, 2003, and 2013 periods. The TDS concentration change declines toward the northern part of the eastern Nile Delta by increase in potentiometric surface from the River Nile in 2010. The geostatistical application (GIS) is used to determine the groundwater quality and hydrogeochemical parameters distribution, which is accomplished by ordinary kriging. The first produced map is the default options, while the next map (2nd) incorporated more of the spatial relationships constructed. When the latter is estimated, the exploratory spatial data investigation (ESDA) techniques check parameters extent. Trends were deleted and spatial autocorrelation was modeled. The ESDA and geostatistical techniques were used and the surfaces of the major ions were found more accurate. The third surface showed critical probability that TDS concentration threshold point increase drinking and irrigation purposes in the northeastern part and decrease in the southwestern parts of the Nile Delta. The default kriging is the best for mapping the hydrogeochemical parameters.


INTRODUCTION
Agriculture is the dominant industries in the Nile Delta area (~ 63% of the total agricultural area) and is caused by soil type (Dawoud, 2004) and watering canals included.The quality of the groundwater was affected by sea level rise, River Nile flows, overexploitation, septic tanks and agricultural practices.Assessment of spatial correlation in hydrogeochemistry is useful for investigating groundwater quality.In geostatistical application, kriging tool is the most used which estimates the piezometric surfaces and hydrogeochemical parameters distribution.Many authors investigate spatial distribution of groundwater quality by Gographic Information System (GIS) software such as Ayazi et al. (2010), Manap et al. (2012), El Afly (2012), Machiwal et al. (2012), Neshat et al. (2013), and Manap et al. (2013).The aims of the dissertation are 1) to determine the anthropogenic sources impact on groundwater quality, and 2) to investigate the application of various spatial models and data transformation to interpret the spatial distribution of the groundwater quality.In this study, kriging techniques in the framework of GIS software (ArcGIS Geostatistical Analyst) are used.

Hydrogeology
Active natural delta has stopped currently (Stanley and Warne, 1996).
The area and coastline length of Nile Delta were 22 × 10 3 km2 and 225 km, respectively (Stanley and Warne, 1996).The topography varies from 18 m (Cairo) to < 1 m (coast).A third of the Nile water was evapotranspirated and leaked into the Quaternary aquifer.Two-third of Nile water was distributed slowly in irrigation systems and drains (Sestini, 1992).The Quaternary deposits consisted of upper Holocene (aquitard) and lower Pleistocene (main aquifer) (Figure 1).The maximum thickness of Bilqas Formation (Holocene) is 71 m in the north and 77 m in the east, and decrease gradually towards the south (Zhaghloul et al., 1989).The leaky aquifer was recharged from the surface waters distributed in the Nile Delta.The mean recharge in the southern and middle parts was 0.8 mm/day, whereas at northern sector, the leakage watering canals is cut by wastewaters networks (El Arabi et al., 2001).The Quaternary deposits in the Nile Delta form a huge aquifer while thickness ranges from 200 m at southern part to 1000 m in the northern (Figure 1).
Transmissivity changes from 5000 to 25000 m -3 in the southern part to 0.0009 in the northern part (Dahab, 1993).

METHODS OF ANALYSIS
The hydrogeochemical parameters of major ions and piezometric surfaces were collected within different periods (1997-2013) (Geriesh, 2003;Geriesh et al., 2015;Saad and Mohamed, 2010;RIGW/IWACO, 1999).The historical review of TDS concentration and piezometric heads within different periods was determined by GIS contour.The impact of anthropogenic sources on groundwater quality was estimated and the correlation between hydrogeology (groundwater flow) and hydrogeochemisrty (TDS content) was achieved.The TDS concentration and piezometric surfaces changes within different periods were determined by GIS software.Spatial distribution of groundwater chemistry parameters is analyzed using a geostatistical analyst tools (GISsoftware).Geostatistic is considered as a set of numerical methods that focus on behaviors of spatial issues mainly represented by random models in such a way that the time interval investigation identifies temporal data (Olea, 2000).It focuses on spatially autocorrelated information with basic spatial forms, whereas the latter was indicated in semi-variogram investigation.Semi-variogram classified spatial correlation of the elements.The semivariogram indicates the relationship between the lag distance on X-axis and the semivariogram value on the Y-axis.Lag distance represents the length among the estimation of a special property.From the semivariogram, the spatial correlation of a spatially changing property could be clarified.The semivariogram point trends from low to high values indicating higher spatial autocorrelation at the small lag length (Nayanaka et al., 2010).Kriging is a technique concerned with making optimal, unbiased estimates of regionalized variables at unsampled locations using the structural properties of the semivariogram and the initial set of data values (David, 1977).The procedure was illustrated in Figure 2, and the following steps (Johnston et al., 2001;Mehrjardi et al., 2008) were performed: a) Exploratory spatial data analysis (ESDA) was performed using ArcGIS software, for the groundwater chemistry, to study the spatial distribution of the data.b) Spatial interpolation for groundwater chemistry data was performed using ArcGIS software, and the ordinary kriging is applied by involving the following steps, semivariogram modelling, model validation using cross-validation, and generation of the groundwater chemistry maps.

Historical review of the TDS concentration and potentiometric surface
The TDS concentration increases in the western and northeastern part during 1997 (Figure 3a), and is caused by seawater intrusion and irrigation practices.The TDS concentration anomalies appeared toward the north and east by pumping that enhanced the seawater intrusion within 2003, 2010, and 2013 periods (Figure 3b-d) periods (Figure 3b-d).The highest TDS concentration (anomalies) was in the north and east due to the subdivision of the Pleistocene aquifer into sub aquifers by aquitard Holocene sediments (Figure 1).The seawater that intrude into the sub aquifers, react with aquitard, and the TDS concentration sometimes go higher than seawater.The groundwater flow was towards the north, northwest, and northeast (Figure 4a-d).There is a good match between hydrogeochemistry (TDS concentration in Figure 3) and hydrogeology (groundwater flow in Figure 4).The TDS concentration contributed to rock water interaction during the groundwater flow from the recharge (south) to discharge areas (TDS anomalies).The water level of the Holocene aquitard is generally higher than potentiometric surface of the Pleistocene aquifer (Geriesh et al., 2015).Moreover, the Holocene aquitard pressure is higher than Pleistocene aquifer that enhances the downward leakage from agricultural activity into aquitard and finally in the Pleistocene aquifer, thus, helping to raise the TDS concentration throughout the Nile Delta area.
The TDS concentration change decline mainly in the southwestern and northeastern part (Figure 5a-c

Creating a surface using default parameters
By default, ordinary kriging and prediction surface was chosen for TDS concentration (Figure 8a).The semivariogram/covariance models (Figure 8a) checks spatial relationships among measured points where    1980, 1992, 1999, 2003, and 2010).The expected point in crosshairs is used alongside the measured locations.TDS concentration map was created and called default kriging (Figure 9a).The highest TDS concentration is more or less in the same areas in default kriging, log kriging, and trend removed maps (Figure 9ac), whereas the lowest TDS content is towards southern part, but spotted in different directions as shown in Figure 9a-c.

Exploring the data
This involves use of ESDA tools to explore data and gather information to build good interpolation models.Examine the data distribution, identify the data trends, and understand the spatial autocorrelation and directional influences.

Histogram tool:
The interpolation techniques can produce a map, which are better outputs if the data is normally distributed.The TDS and major ions concentrations skewed (lopsided) and both the mean and the median greatly differs (Figures 10 and 11), indicating that they are not normally distributed.The interpolation method of the information always focuses on their normality.The information indicated high skewness, and accordingly, they normalized by the suitable change in order to give a good accuracy.In the ESDA of ArcGIS, the logarithmic transformations applied to all hydrogeochemical elements that were attributed to positive numbers.Figures 10a-d    QQ figure (Figures 10c and 12), they transform to logarithmic (Figure 10d and 12) for normal distribution.Figure 12 is near to being a straight line in case of HCO 3 , while the rest ions depart from the line before transformation.2nd level polynomial) produces a U form (Figure 13a).The information shows positive trend from the center directed around, that is, towards the sea (north and east) and agricultural activity (south and west).Because the trend is U formed, a 2 nd level polynomial is a best select to apply as a global trend model.This trend is possibly attributed to the fact that the contamination is high at northern and eastern parts and low due southern and western parts.

Explore spatial autocorrelation and directional influences:
The semivariogram/covariance cloud (Figure 13b) checks the spatial autocorrelation among the estimated locations whereas its supposed points nearer to each other are dissimilar.Semivariogram point is the squared change among the values of each pair of points, and is figured on the y-axis.X-axis is the length inbetween each pair of determination (Figure 13b).The semivariogram/covariance cloud indicates a pair of points.Points nearer to each other are dissimilar.In the semivariogram figure, the nearest points (on the far left on the x-axis) have low semivariogram (low values on the y-axis).As the length among the pairs of points increases (moving right on the x-axis), the semivariogram increase (move up on the y-axis).Some points near to each other (near zero on the x-axis) have a higher semivariogram (high on the y-axis) than expected (Figure 13b), reflecting inaccurate chemical analysis.

Mapping with trend removed
This involves produce reform model by ordinary kriging, trend delete, and checks for anisotropy.TDS concentration exhibited a trend and a directional influence.Once again, it applies the ordinary kriging interpolation technique, but this time includes trend and anisotropy to produce good expectation.Ordinary kriging is the simplest geostatistical model because the number of assumptions behind it is the lowest.Once the trend deleted, the statistical investigation will be performed on the residuals or the short-interval variation of the map.The trend is automatically added back before the last map was determined (Figures 8c and 9c).Accordingly, the expectations are the logic outputs.
Semivariogram/Covariance modeling: The map represents the fastest difference in the northern-southern direction and a more gradual variation in the east-west direction (result in ellipse form).The northern and southern trend attributed to seawater intrusion and agricultural impact, respectively.The eastern and western trends were influenced by seawater intrusion (Suez Canal) and agricultural activity, respectively.The aim of semivariogram/covariance modeling (Figure 14) is to estimate the best fit for a TDS model in the semivariogram.By deleting the trend, the semivariogram will model the spatial autocorrelation with no application of data trend.The latter will be automatically added back to the estimations before the last map was created (Figure 9c).Semi variograms determined after normalizing the data by ArcGIS Geostatistical Analyst in Universal Transverse Mercator (UTM) coordinates system.Ordinary kriging used semivariogram modeling performed on the experimental semivariogram obtained from the data by checking the angle direction, angle tolerance, and lag size, until the model fits the data.
Directional semivariograms: Directional impacts on the locations of the semivariogram and the model will be fit.In specified directions, points nearer to each other are dissimilar than in rest directions.Geostatistical Analyst applies in the directional influences, or anisotropy in the semivariogram model.Anisotropy can be caused by lithogenic and anthropogenic sources.The directional influence (Figure 15a-c) can be statistically quantified and accounted for TDS concentration.Further, the semivariogram model increases more gradually and then flattens out (Figure 15a-c).

Cross-validation:
The aim of cross-validation is to ascertain the model that gives the best expectation.It is compared by the following: i) The expectations are unbiased as shown an average expectation error near zero.
ii) The standard errors are accurate as shown by a root- iii) The expectations do not deviate much from the estimated points as shown by root-mean-square error and mean standard error that are as small as possible.
A comparison of the two interpolation models (Figure 16) based on the previous rules show that the default kriging (Figure 9a) is the best.After fitting the semivariogram model, the cross-validation tool was applied to determine how the aquifer chemical model is expected to be at unsampled points.The cross-validation procedure deletes each data point, one at a time, expects the aquifer chemical contents, and then compares it with the estimated.In the cross-validation, the determined statistics was applied to diagnose whether the model and/or its associated elements concentration are logic.
The expected ratio between the predicted and the estimated points should be at slope line 1:1.Besides, the expectation error clarifies the change between the expectation and the actual determination.The accuracy in the expectation outputs estimated by the average error that should be near to zero; the root mean square error and the average standard error should be as small as possible; and the root mean square-standardized error should be near to one (Johnston et al., 2001).The best model identifies the aquifer chemical concentration maps (Figure 17).
The Ca and Mg concentrations increase in the northern and southeastern parts (Figure 17a-b) due to seawater intrusion and agricultural activity, respectively.It decreases towards the southwestern part attributed to inflows of River Nile.The lowest Ca and Mg concentrations scattered in eastwest line (Figure 17ab).The Na and Cl contents enrich the northern and northeastern parts (Figure 17c-d) because of seawater intrusion.The SO 4 concentration increases in the northern and northeastern parts (Figure 17e) due to its being affected by seawater intrusion.While the central decline in SO 4 concentration is caused by both agricultural activity (fertilizers) and septic tanks, the lowest Na, Cl, and SO 4 contents in groundwater was mainly scattered around River Nile (Figure 18c-e).The HCO 3 content increases in the eastern and central parts (Figure 17f) due to inflows of River Nile.

Mapping the probability of TDS exceeding a critical threshold (drinking and irrigation)
The probability of TDS concentration that exceed a drinking (1000 ppm) and irrigation (2000 ppm) WHO (2004) guidelines is created (Figure 18a-b).Any value below threshold is represented by zero while the values above threshold are indicated by one.The probability, for drinking application, increased due northeastern by seawater intrusion and anthropogenic sources and decline in the southwestern parts (Figure 18a) by inflows of the River Nile.In irrigation purpose, it increased in the northern and eastern parts and decreased in the western part (Figure 18b).

CONCLUSION AND RECOMMENDATION
The Nile Delta is influenced by anthropogenic sources.The TDS concentration change decline mainly in the southwestern and northeastern part within 1997, 2003, and 2013 periods.It is attributed to recharge from River Nile.The TDS concentration change increases due the northern and eastern parts as a result of seawater intrusion, irrigation practices, waterlogging, increase in pumping rate and rock water interaction.ArcMap with geostatistical analyst gives many methods for producing maps, which surface, analyze, and understand spatial phenomena like seawater intrusion and agricultural impact.Kriging technique is a useful tool in mapping distribution of the hydrogeochemical elements.A log transformation is applied in order to reduce the skewness of the data.On comparing the two interpolation models,      the default kriging is the best.The probability of the TDS content for drinking and irrigation purposes decreases in the southwestern parts.

CONFLICT OF INTERESTS
The author has not declared any conflict of interests.

Figure 2 .
Figure 2. Flow chart of the steps followed for the geostatistical analysis.(modified after Marko et al., 2013).

Figure 3 .
Figure 3. Spatial distribution of TDS concentration within different periods.
) within 1997, 2003, and 2013 periods and is attributed to recharge from River Nile.The TDS concentration change increases due the northern and eastern parts (Figure 5ac) due to seawater intrusion, irrigation practices, waterlogging, increase in pumping rate (Figure 6), and rock water interaction.The potentiometric surface change increases in the eastern part (Figure 7a-c) by seawater intrusion from Suez Canal, River Nile, and irrigation practices; but showed patchy decrease due western part (Figure 7a-c).The TDS concentration change increases due south, northwest, and northeast and depleted in the central and western part (Figure 5b) within 2003 and 2013 periods.The potentiometric surface change between 2003 and 2013 increases in the middle and eastern part of Nile Delta (Figure 7b), leading to increase in TDS concentration in most of the irrigated areas.The potentiometric surface change decreases generally in the western part (Figure 7b) by low recharge from River Nile.

Figure 3 .
Figure 3. Spatial distribution of TDS concentration within different periods.

Figure 4 .
Figure 4. Spatial distribution of potentiometric surfaces within different periods.

Figure 4 .
Figure 4. Spatial distribution of potentiometric surfaces within different periods.

Figure 5 .
Figure 5. Spatial distribution of TDS concentration change within different periods.
closer things are dissimilar than apart.The semivariogram model fitting is applied in order to touch the spatial correlations of the information called variography.The crosshairs show a location that has no measured value.

Figure 5a .
Figure 5a.Spatial distribution of TDS concentration change within different periods.

Figure 7 .
Figure 7. Spatial distribution of piezometric surfaces change within different periods.

Figure 8 .
Figure 8.The semivariogram/covariance model of the TDS concentration.
and 11 clarifies the difference in the skewness coefficient values pre-and post-change.Hydrogeochemical elements such as calcium, magnesium, sodium, chloride, suphate, and bicarbonate have low skewness post change, thus, they need

Figure 9 .
Figure 9. Kriging types of the TDS concentration.

Figure 10 .
Figure 10.Histogram and QQ plots of the TDS concentration before transformation (a & c) and after logarithmic transformation (b & d).

Figure 11 .
Figure 11.Histograms of the major ions before transformation and after logarithmic transformation.

Figure 13 .
Figure 13.Global trend analysis of the TDS concentration.

Figure 14 .
Figure 14.Semivarigrams of the TDS and major ion concentrations with trend removed.

Figure 15 .
Figure 15.Semivarigram with trend removed and directional influence of TDS concentration.

Figure 16 .
Figure 16.Cross validation between default and trend removed kriging of the TDS concentration.

Figure 17 .
Figure 17.Best models of the major ions.

Figure 16 .
Figure 16.Best models of the major ions.

Figure 18 .
Figure 18.Probability map of the TDS concentration based WHO guidelines (2004).