An Integrated Statistical, Geostatistical and Hydrogeological Approach for Assessing and Modelling Groundwater Salinity and Quality in Nile Delta Aquifer

1. Introduction

Coastal aquifers represent an important source of water especially in arid and semiarid regions. Stresses on coastal aquifers are related to salinization due to seawater intrusion (SWI) [1,2,3,4,5]. The main causes of SWI are the over-pumping of groundwater [6,7,8], sea level rise as a result of climate change and its subsequent global warming [9,10,11,12], as well as overpopulation along coastal areas [13,14,15]. Accordingly, SWI is seriously damaging the groundwater quality as well as soil, cultivation, and ecosystem [6,16,17,18], therefore, the management of SWI and resulting environmental problems have become a research interest [19,20,21]. SWI in coastal has been assessed through different methodological approaches, such as hydrogeochemical investigations [7,22], geophysics [23,24], isotope hydrogeology [25,26], and data analysis and numerical modeling [27,28,29,30], machine learning algorithms [31,32] and geostatistical modelling [33,34].

The Nile Delta aquifer in Egypt is one of the most important aquifers in the world; however, it suffers from SWI in the coastal parts of the Mediterranean Sea as a result of anthropogenic impacts (Figure 1) [7,35,36].

Due to the complexity of spatially and temporally variable processes that determine the salinity of subsurface aquifers, a reliable assessment of salinity-related risks is highly dependent on the ability of statistical tools to accurately capture spatial variability and to estimate the spatial relationships between salinity and other indicators of groundwater quality such as chloride and sodium content and electrical conductivity. The spatial association of the estimates of all these indicators at non-sampled locations with multivariate interpolation methods also makes it possible to delineate homogeneous areas by salinisation level to be submitted to site-specific (differential) management. For these reasons, in this study we used multivariate geostatistics as a branch of applied spatial statistics, which is based on a spatial dependence model to provide the prediction of indicators and the uncertainties of these predictions [37]. Above all, it is capable of combining different types of data, with different physical and statistical properties, taking into account their proximity both over space and time and the different support of the measurements [38].

The objective of this paper was to delineate zones with different salinization level over time and investigate the factors controlling seawater intrusion of the Nile Delta aquifer, using multivariate geostatistics compared with the classical analysis of the raw chemical data, with the perspective also to design effective groundwater management policies.

2. Materials and Methods

2.1. Hydrogeological Setting

The Nile Delta aquifer of Egypt is one of the most promising freshwater aquifers in the world. However, it suffers from saltwater intrusion from the north, where it is directly in contact with the Mediterranean coast. During the past 30 years, the aquifer has been severely pumped for domestic and irrigation purposes, which has been causing a direct impact on groundwater levels and its drawdown. This overexploitation of groundwater resources has caused significant salinization and successive well abandonment, especially in the central and northern parts of the Nile Delta.

The aquifer system is mainly represented by the sandy-gravelly Pleistocene aquifer with interbedded clay lenses (Mit Ghamr Formation) (Figure 2). This aquifer is overlaid by the upper Holocene sandy clay aquiclude layer (Bilqas Formation) and underlaid by the aquiclude Pliocene clays in some places and the Miocene limestone and marl in others [32,39](Figure 2a). The upper sandy clay cap layer is composed of Nile silty clay, silt and sandy clay with varying thickness from about 20 m in southern and central parts to about 70 m in the northern part [40]. The thickness of the Pleistocene aquifer increased towards the coastal area [7,41,42], where it usually varies between 250 (southern portion) and 900 m (northern portion). The Pleistocene aquifer varies between the unconfined condition in the south-eastern parts, from semi to unconfined condition in most eastern and northern parts and confined conditions in central and western parts of the study area. This variation in the aquifer nature is related to the lithological composition and /or missing of the upper Holocene cap cover, which serves in some areas as aquitard layer of coarser grains allowing upward leakage or downward flow [35,39].The distributed fresh water canals (Figure 2b) on the River Nile which dissect the Holocene cap layer represent the main source for the aquifer recharge in addition to the contribution of irrigation return flow [35,43].

The hydrogeological setting summarized above, combined with the fact that all the wells analysed are limited to a maximum depth of 130 m (with an average of 57 m), allows us to consider a study volume contained in a single multilayer aquifer mostly developed between Holocene sandy clays and Pleistocene clayey sands. This assumption is corroborated by the presence of intervals of clayey sands within the Holocene succession and sandy clays within the Pleistocene succession. Groundwater levels in the Pleistocene aquifer decreases towards the northern parts and ranges from 8 to 14 m above sea level (m asl) in the southern parts of the Nile Delta to less than 1 m asl in the northern parts [39,40,42,44].

2.2. Data Collection

Hydrochemical and hydrogeological data of the Nile Delta aquifer in the study area of three sampling campaigns of years: 1996, 2007 and 2018, were used which considered different monitoring networks (Table 1) and whose sampling points had variable total depths. The monitoring networks also included hand dug wells (18 – 45 m depth), irrigation wells (35–90 m wells) and drinking wells (100–120 m depth).

Before collecting groundwater samples (year 2018), each well was purged for at least 30 minutes to remove stagnant water. Total dissolved solids (TDS), pH and electrical conductivity (EC) were measured in the field prior to the sample collection by using calibrated field Orion portable meters (Thermo Scientific, Singapore). The collected samples were stored in iceboxes and sent to the central laboratories of the Productive Sufficiency Institute (PSI) in Zagazig University for different hydrochemical analyses. The data were used to construct a detailed database of the aquifer water quality over time.

Ion chromatography (Dionex ICS-1100) was used to analyze the major cations and anions (Mg²⁺ , K⁺, Ca²⁺ , Na⁺, , SO4^2- , and Cl^-), while HCO₃ ^- was determined by titration using H₂SO₄ (0.01 N).

2.3. Raw Chemical Data

In order to characterize the hydrogeochemical structure of the aquifer, and to compare geostatistical approach with the traditional one, geostatistical calculations were supported by analysis techniques typical of hydrochemistry. For this purpose, the hydrochemical data set was processed, classified and compared, for each year of survey, using the Schoeller – Berkaloff graph.

The Schoeller Berkaloff graph can be used both to classify a single water sample and to classify and compare more than one sample. It consists of a graph with the main cations and anions dissolved in water (Ca, Mg, Na+K and Cl, SO4, HCO3+CO3) positioned on the horizontal axis and the concentrations of the ions, expressed in meq/l, on the logarithmic vertical axis. Each water sample will therefore be represented by a broken line joining the concentrations, which will highlight the presence of dominant, secondary and minor ions. The analysis of two or more broken lines (groups of broken lines) will allow to perform comparisons between the represented waters, based on chemistry. In this type of graph, waters with the same concentration but different chemistry will then be easily distinguishable. The graph can also show the relationships existing between the various ions in a water or a group of waters. The raw water chemistry data were also analyzed with regard to their spatial distribution with various geospatial techniques, to provide spatial distribution maps for each ion and each survey year, which will be discussed in section 2.4. The analysis of the chemistry of the waters with the Schoeller Berkaloff graph allowed us to select the dominant ions of the waters which together with salinity (represented by Electrical Conductivity) were used in the creation of the spatial distribution maps. To verify any piezometric and groundwater flow variations, the maps of hydraulic head over the three years were also estimated and compared with the hydrochemical ones.

The evolution of water chemistry over time was analyzed besides with the Schoeller Berkaloff graphs, also with the QQ plots of hydraulic head, EC, Cl, Na and Ca for the pairs of contiguous surveys 1996 - 2007, and 2007 -2018. The QQ Plot graphs (which relate the quantiles of a physical or chemical property in two different surveys) were not created with raw data but from the geostatistical estimations, due to the non-congruence of measurement locations among different surveys. The graphs compare two data distributions and, if there is not perfect consistency, the data depart from the bisector of the first quadrant.

2.4. Data Processing

2.4.1. Exploratory Data Analysis

Exploratory analysis was aimed at describing statistically the data distributions of the variables. Conventional basic statistics, including maximum, minimum, mean, median, standard deviation (SD), skewness, and kurtosis, were calculated.

Q-Q plot was used to compare the data distribution of the same variable over time. Although the interpolation techniques of geostatistics are not parametric, in the case of strongly skewed distributions, to facilitate the variogram fitting process and produce an estimate of uncertainty [46], it was preferred to transform all variables characterized by measurement unit and magnitude different into standardized Gaussian variables with zero mean and unit standard deviation. For this purpose, Gaussian anamorphosis, by using an expansion of Hermite polynomials truncated at the first 100 terms [38,47] was performed.

2.4.2. Multivariate Geostatistical Analysis

For a detailed description of the multivariate techniques of geostatistics, reference is made to the many textbooks on the subject [37,47,48]. Only those techniques most closely specific to the case study will be briefly mentioned here. For each survey date (1996, 2007 and 2018), a co-regionalisation dataset was created as follows: at each well location together with the hydrochemical sample data three auxiliary raster variables were collocated by migration of the pixel to the nearest sampling point. The auxiliary variables were: the minimum distance from the Damietta branch, which is a branch of the Nile River located in the west-northern part of the study area as an indicator of freshwater resource; the distance from the Manzala lake located in the north-eastern part of the study area and the Suez Canal, located in the eastern part of the study area, as an indicator of saltwater resource (Table 2 and Figure 3). The distances of the three auxiliary variables were calculated using the Near tool, which is available in the ArcGIS Pro software, developed by ESRI. The tool calculates the distance and additional proximity information between the input features and the closest features in another layer or feature class using the tools available in ArcGIS Pro 2.7.0, USA.

The whole set of both experimental direct and cross-variograms was fitted to a Linear Model of Coregionalization (LMC) [49], which considers all the studied variables as the result of the same independent physical processes acting over different spatial scales. All direct and cross semi-variograms of the variables are then modelled as a linear combination of standardized semi-variograms of unit sill, each one corresponding to an identified spatial scale.

The sills of the direct and cross variograms of each scale constitute a matrix called coregionalization matrix, which must be a symmetric positive semi-definite matrix. Each coregionalization matrix can be interpreted as a variance-covariance matrix specific to a given spatial scale [50].

Multi-Collocated Cokriging

Ordinary cokriging (OCK) is one of the most basic multivariate geostatistical methods of interpolation, where the local mean is assumed to be constant but of unknown value. A way of integrating secondary finer resolution grid information into primary sparse variable modelling is multi-collocated cokriging [51]. It is a technique is very similar to ordinary cokriging but with the only difference being the neighbourhood search. Since using all the secondary exhaustive information contained in the neighbourhood can lead to an intractable solution due to too much information, the secondary variable is used only at the target location and also at all locations where the primary variable is defined. It is less accurate than full cokriging, because it does not use all the auxiliary information contained within the neighbourhood, but it is much less computationally intensive. Nevertheless, as the co-located secondary datum tends to screen out the influence of more distant secondary data, there is actually little loss of information.

Cross-validation

This involves deleting each sample in turn and then kriging it (z*) independently of all other points in the estimation neighbourhood [52].

Two statistics[51], defined as follows, were used to evaluate the performance of the LMC using cross-validation:

$$ME={\displaystyle \frac{1}{N}} \sum _{i=1}^N\left(z_i-z^{*} \right)$$

$$MSSE={\displaystyle \frac{1}{N}} {\sum _{i=1}^N\left({\displaystyle \frac{z_i-z_i^{*}}{{\sigma }_i}} \right)}^2$$

where N is the number of active observations and σ_i the cokriging standard deviation of the estimate $z_i^{*}$.

For an unbiased and accurate estimate, ME should be near zero and MSSE near 1, because the latter approximately represents the ratio of an experimental and a theoretical variance [53]. However, if MSSE is different from 1 but within the tolerance interval $(1-3\sqrt{{\displaystyle \frac{2}{N}}};1+3\sqrt{{\displaystyle \frac{2}{N}}},)$[54] the estimation may be considered accurate.

• Factor cokriging analysis (FCA)

To synthesise the complex multivariate variation of the study area in a limited number of roughly homogeneous macro-zones, which could therefore be subjected to differential management, factor cokriging analysis (FCA) was applied using the same LMCs previously estimated for each sampling date.

FCA is quite similar to the traditional Principal Component Analysis (PCA) but with the difference that each co-regionalization matrix (corresponding to a given spatial scale) is broken down into the eigenvalues and eigenvector matrices [48]. In FCA, only the eigenvectors, called regionalised factors (FC), corresponding to eigenvalues greater than one were retained, because they describe a larger portion of spatial variance than that of each variable standardized to variance 1. The regionalised factors corresponding to the nugget effect were omitted, as they are mostly affected by measurement errors.

Mapping the estimated retained regionalized factors provides a synthetic illustration of the behaviour and relationships among variables at different spatial scales. The interpolation of the FCs is performed by solving a modified co-kriging system, as described by[55].

3. Results

3.1. Exploratory Data Analysis

Table 3 shows the descriptive statistics of the whole dataset for the three different surveys. This has been done to investigate the shape of the data distributions before carrying out the geostatistical analysis. The skewness values reported in Table 3 indicate that most of the study variables for each survey are skewed either positively or negatively except for total depth in 1996 and 2007 and water level in 2018. Therefore all variables were transformed using the Gaussian anamorphosis prior to geostatistical analysis.

3.2. Raw Hydrochemical Data

The chemical analyses of all samples for each survey, processed through the Schoeller - Berkaloff diagrams (Figure 3), show the presence of two main groups of water: predominantly calcium and/or sodium bicarbonate (Figure 3 B, E, H ) and predominantly sodium chloride (Figure 3 C, F, I). The first group includes waters with chemical formula:

Na, Ca, Mg – HCO3, Cl, SO4; Ca, Na, Mg – HCO3, Cl, SO4; Mg, Ca, Na – HCO3, Cl, SO4.

The second group includes waters with chemical formula:

Na, Ca, Mg – Cl, SO4, HCO3; Na, Mg, Ca – Cl, SO4, HCO3; Na, Ca, Mg – Cl, HCO3, SO4.

The cations of the first group (bicarbonate-calcium/sodium waters) do not show the prevalence of a particular ion but, moving from the first to the last sampling campaign, they show a decrease in the samples with high Na and K contents. As regards anions, a constancy in the HCO₃ contents over time and a variability in the Cl and SO₄ contents are observed. Concentrations are between 1 and 10 meq/l.

In the second group (sodium chloride waters) the concentrations are higher and vary from 1 to over 100 meq/l and therefore more than 10 times the concentrations of the first group. This situation is clearly visible in the Scholler Berkaloff diagrams which show water types of the second group located in the upper part of the graph. The concentrations also show a slight increase in values over time and a high variation in SO₄ content, so much as to highlight the presence of some sulphate-alkaline samples.

3.3. Geostatistical Results

• Linear model of Coregionalization

Three linear models of coregionalization of the transformed Gaussian variables were applied for three periods 1996, 2007 and 2018 (Table 4, Table 5 and Table 6 respectivly).

Since it was not possible to check whether any anisotropy existed given the small number of samples, it was preferred to fit an isotropic model of spatial dependence. LMC developed for the dataset collected in 1996 involves three spatial structures: a nugget effect; a spherical model with range of 1200 m and a Bessel-K model with a range of 100000 m; the one of 2007 involves 3 spatial structures: a nugget effect; a spherical model with range of 10000 m and a Bessel-K – scale with a range of 100000 m (Table 7); the one of 2018 involves 3 spatial structures: a nugget effect; a spherical model with range of 1000 m and an isotropic Bessel-K – scale with a range of 100000 m.

As can be seen, the spatial structures apart from absolute values of the parameters remain invariant enough over time, especially that related to the long range, most likely associated with hydrogeological characteristics that do not change over time. On the contrary, it is different for the shorter-range structure, which has a variable range over time because it is more related to anthropological activity and possible movements of groundwater flows (saline water intrusion and/or freshwater adductions). The appropriateness of the selected LMC models was tested by cross-validation. These results fell within the standard range (results not shown), indicating that the model performance was accurate and precise.

3.4. Spatial Distribution

The multi-collocated co-kriging based on the LMC allowed estimating the spatial distribution of hydraulic heads, EC, Cl, Na, and Ca. Pumping-induced hydraulic head variation is usually considered the main anthropogenic cause of SWI, whereas the four physico-chemical variables can be considered as hydrogeochemical proxies of aquifer salinization, as they describe the increase in groundwater salinity and water-type [56].

In Figure 4 the hydraulic head evolution is reported from which it is possible to note the following:

- In 1996 a regular piezometric distribution with almost homogeneous gradients, flow lines perpendicular to the coastline and parallel to the Damietta Branch of the Nile.

- In 2007, a less regular distribution of the piezometric surface was observed with a lowering of the hydraulic head in the central and southern part and the presence of concentrated recall zones in the central part.

-In 2018, a further generalized reduction in hydraulic head was observed in the southern and northern portions of the study area and a greater number of water recall areas compared to 2007. In the central part, however, a recovery in the hydraulic head was observed which gave rise to a lobed shape of the piezometric surface (circled in red in Figure 4) which indicates an alternation of areas of probable attenuation of the pumping.

The thematic maps of groundwater salinity (expressed in EC) are shown in Figure 5, where a clear trend of increasing salinity towards the coastal area can be seen in the 1996 survey. In the 2007survey, the level of salinization increased even further, especially nearer the Manzala Lake. The pattern of salinity in the 2018 survey is somewhat similar to that of 2007, with salinity increasing towards the north-west part of the study area. It is worth noting that in 2018 survey two zones of ground water salinity less than 2 dS/m appeared in the eastern part of the study area corresponding to the lobate shape of the piezometric surface described above.

The spatial distribution of chloride (Figure 6) during the reference years was consistent with the EC one. The lower concentrations of chloride in the southwestern part might be attributed to the proximity of fresh water sources, such as the Damietta branch and the Ismailia canal. It is worth noting that in 2007 chloride concentrations in the northern part increased dramatically from around 1000 mg/l to more than 1500 mg/l, if compared to the chloride concentrations in the 1996 map. This evidence likely relates the groundwater overexploitation for crop irrigation in this cultivated area. This interpretation is further supported by the decrease in hydraulic head between 1996 and 2007 (Figure 4).

The sodium spatial distributions (Figure 7) show that in 2007 sodium concentrations increased from the northern part towards the southwestern part, covering about 40% of the study area and reducing the water quality of the wells.

For both chlorides and sodium in 2018 there was a slight decrease in the areas of highest concentration. Furthermore, the spatial distribution shows some zones located in the northeastern part with lower sodium concentrations of less than 100 mg/l, once again at the lobate piezometric surface.

Figure 8 shows the spatial maps of calcium, where an opposite pattern to that of EC, Cl and Na was observed. Its concentration decreased especially in the northern part in 2007. Comparing the calcium concentrations between 2007 and 2018, it can be seen an increase in calcium concentrations in the north-eastern part. New zones of high calcium concentrations appeared in the central part of the 2018 map corresponding to the lobed zone of the piezometric surface.

Considering the spatial distribution obtained from the Cokriged maps described above and to better describe the evolution of chemistry over time, the Q-Q plot graphs (Figure 9 and Figure 10) describing the evolution of the estimated EC, Cl, Na and Ca are shown below in relation to the hydraulic head.

Moreover, to analyze the evolution of SWI at a regional scale and over the considered period, we took advantage of Q-Q plots to compare the estimated proxy variables (i.e., hydraulic head, EC, Cl, Na, and Ca) and interpret the results from a hydrogeological and hydrogeochemical standpoint (Figure 9).

As for the estimated hydraulic heads between 2 and 4 m a.s.l., no significant changes were observed between 1996 and 2007. However, for values between 4 and 10 m a.s.l. there was a significant decrease in 2007 compared to 1996, indicating a decrease in the hydraulic head.

In the second period from 2007 to 2018 (Figure 10) the hydraulic heads decreased but in a less evident way. This probably means a decrease in pumping rates compared to the previous period.

For groundwater salinity the estimates up to 5 dS/m remained almost the same for both years 1996 and 2007 (Figure 9). However, there was an increase in water salinity especially for the estimates between 5 and 25 dS/m in 2007, indicating that the water salinity level increased during the period 1996-2007. When comparing 2007 and 2018, no significant changes were observed. The only change was the increase in salinity for values above 23 dS/m in 2018. A similar trend to that observed for salinity was observed for chloride, which shows that it increased between 1996 and 2007 (Figure 9), whereas remained slightly higher in 2007 than in 2018 up to the value of 8000 mg/L, beyond which the trend reversed.

There was an increase in sodium concentration in 2007 compared to 1996 (Figure 10) whereas an opposite pattern was observed in 2018 compared to 2007, indicating a decreasing Na concentration during the period 2007 - 2018, quite similar to the behaviour of Cl. Regarding the calcium content, the Q-Q plot (Figure 10) of 1996 vs 2007 shows that there was a decrease in its concentration in 2007 compared to 1996, especially for the values of 120 mg/L and above. Differently, calcium concentration increased in 2018 compared to 2007.

Ultimately, the Q-Q plot graphs graphically render what is already evident in the analysis of the spatial distributions in the three study periods. However, it is a concise though non-spatial way of characterising the entire distribution of waters at the three survey dates.

• Delineation of homogenous zones based on factorial kriging

show the structure of the regionalised factors whose eigen value is greater than one by spatial scale and survey. In 1996 (Table 8) on the first long-range factor with eigen value greater than 1, which explained 46.56 % of the spatial variance at that scale (100000 m), the distance from the freshwater source (Damietta branch), EC and Cl weigh more and positively, whereas the distances from the saltwater sources (distance from Lake Manzala and the Suez Canal (Table 4) negatively. This factor could therefore be interpreted as an indicator of the salinity level of the waters resulting from the equilibrium of lateral fresh and salt water supplies. On the second long-range factor with eigen value greater than 1, which explained 21.18% of the spatial variance at this scale, calcium, potassium and distance from Suez Canal weigh mainly and positively, whereas total depth and distance from Damietta branch negatively. This factor could therefore be interpreted as an indicator of water quality, particularly for calcium and potassium content. On the third long-range factor with eigen value greater than 1, which explained only 14.5% of the spatial variance at this scale, the distances from Damietta branch and lake Manzala, bicarbonate, and water level weigh negatively. Indeed, giving a physical interpretation to this factor is rather difficult, showing an inverse relationship with water table depth and bicarbonate content. The fourth long range factor with an eigenvalue greater than 1 explained a variance rate of less than 10%, so it was decided to neglect it. At short range there were no factors with eigen values greater than 1.

In 2007 (Table 9) the first short-range factor with an explained variance of 35.87 % was positively correlated with Ca and Mg and at a less extent with Cl and EC, and negatively correlated with pH. It could then be interpreted as an indicator of the level of solutes in the waters. The retained second short-range factor, with an explained variance percentage of 28.55, was negatively correlated with Cl, EC, Na, pH and total depth. It could then be assumed as an inverse indicator of the salinization level and then of the water quality.

In the same survey there were 2 regionalised long-range factors with eigenvalues greater than one. The first factor explained 64.7% of the variance at this scale and was positively correlated with Cl and Na and negatively correlated with distance from Lake Manzala and water level. It could then be assumed as an indicator of saline water intrusion from lake Manzala. The second factor, which explained only about 28% of the spatial variance, was positively correlated with distance from the Suez Canal (salt water) and negatively correlated with distance from the Damietta branch (fresh water). Therefore, also this factor could be interpreted as an indicator of the salinisation level of the waters (inverse indicator of water quality) resulting from the balance of saline and fresh water supply.

Finally, in 2018 survey (Table 10) the first long range factor, which explained 71.9% of the scale variance, was mostly positively correlated with Cl, EC, Na, Mg and SO4 and negatively correlated with distance from Lake Manzala. It is then a clear indicator of the salinity of the waters resulting from the proximity of the salty lake. The second factor explained 20.55 % of the scale variance and was positively correlated with distance from Suez Canal and negatively correlated with bicarbonate content, of which it could be taken as an inverse indicator. The third factor explained less the 10% of the scale variance and was then omitted. At short range there were no factors with eigen values greater than 1.

• Mapping first factors in the three surveys

It was therefore decided to take the first long-range factor for the three surveys as a scale-consistent indicator of the level of salinisation of the waters. A partition the study area into three macro-areas of equal size, each comprising approximately 33% of the total number of pixels, was made for each survey. Although the level of salinisation and its gradient change in the three surveys (the three colour scales are not the same), as it is evident also from the previous investigations, there was a consistent downward trend in salinity in the NE/SW direction during the whole period 1998 – 2018 of study. This decrease is manifested in a compact south-western area and in a more irregular central-eastern area with alternations of zones with different salinity (Figure 11Figure 12 and Figure 13).

4. Discussion

Schoeller Berkaloff diagrams (Figure 3) generally show an increase of sodium chloride over time and a related increase in total salinity (graphs C, F and I in Figure 3). Furthermore, they show a decrease in calcium bicarbonate over time (graphs B, E and H in Figure 3). The comparative analysis with the spatial distribution (Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8) and with the QQplots (Figure 9 and Figure 10) allows us to affirm that the increase in sodium chloride water occurred due to migration towards the southwest at the expense of volumes of calcium bicarbonate water. The comparative analysis with the reconstruction of the piezometric surface over time (Figure 4) also allows us to ascribe the aforementioned variations to a general decrease in the hydraulic head starting from a situation of equilibrium in 1996 up to a situation of general pumping with the presence of localized pumping areas capable of drawing brackish water from Lake Manzala.

In general, it is possible to note that over 22 years there has been an increase in sodium chloride, a general increase of total salinity of the waters and a migration of water salinization towards the southwest.

Compared to this general trend, the central band, approximately 35 km wide and extending for approximately 10-20 km north of the Ismailia Canal, presents piezometric and chemical structures between 2007 and 2018 that are peculiar and different from the general structure described above. The piezometric surface from 2007 to 2018 evolves into a lobed shape which highlights the alternation of piezometric highs and lows which testify the decrease, albeit not homogeneous, in pumping. In this area we also observe:

- in 2018 the appearance of two zones of ground water salinity less than 2 dS/m (Figure 5);

- a decrease in the concentration of Cl and Na (Figure 6 and Figure 7);

- an increase in Ca concentration in 2018 (Figure 8).

The Ca concentration shows an inverse trend compared to that of Na and Cl. In 2007 the Ca concentration decreased compared to that of 1996 because the sodium chloride waters replaced the bicarbonate ones with an increase in salinity due to the decrease in the hydraulic load which draws salt/brackish water from Manzala Lake. In 2018, at the same time as the uneven decrease in pumping, a recovery in the concentration of Ca and a slight decrease in Cl and Na were observed, probably due to the recall of bicarbonate surface water (i.e. not sodium chloride) from the Ismailia channel. The piezometric surface reconstructed in 2018, in fact, shows a water flow from the Ismailia canal towards the northern area.

The analysis of the raw chemical data and their spatial distribution estimated in three time windows show, therefore, a generalized decrease in the hydraulic head with a consequent increase in salinity and Na and Cl content over time and towards the southwest. This trend appears to be mitigated in the central part of the study area where the decrease in pumping rates (even if not in a homogeneous way) and the recall of surface water from the Ismailia Canal attenuates salinization and generates an increase in Ca content.

The above described geospatial and multi-temporal distribution of ion concentration, salinity and hydraulic head (piezometric surface) was implemented with the identification of homogeneous zones based on the kriging factorial analysis that considered the distance from the salty (Suez Canal and Manzala Lake) and fresh (Damietta Branch and Ismailia Canal) surface water bodies as auxiliary variables. The result of such analysis clearly highlights the influence of the above-mentioned surface water bodies in the evolution of the chemistry of the groundwater both over time and space. In fact, a general increase in water salinity from the northeast to the southwest was shown (Figure 11 and Figure 12), due to the lowering of the hydraulic head (Figure 4) and the consequent recall of water from Manzala Lake in the period from 1998 to 2007. Differently, a decrease in salinity in 2018 (Figure 13) was estimated, due to probable decreases in pumping rates which caused a slight decrease in salinity in the southwestern part and an irregular decrease in salinity, Cl and Na concentration in the central eastern part (Figure 5, 6, 7 referred to 2018). This was probably due to the decrease in pumping rates and the recall of fresh water from Ismailia Chanel, as highlighted by the distribution of piezometric heads (Figure 4).

5. Conclusions

The proposed multidisciplinary approach, based on the chemical analysis of raw data, their spatial interpolation using multivariate geostatistical techniques, allowed to evaluate the temporal evolution and the extent of salt water intrusion in the eastern portion of the Nile Delta.

In particular, the dependence of the groundwater chemical parameters on the distance between salt and fresh surface water bodies was introduced for the first time. The distances to salt and freshwater bodies in the Linear Model of Coregionalization, as auxiliary variables, made the spatial distribution of key variables obtained by Ordinary Cokriging more reliable and physically based.

Factorial Kriging analysis using Gaussian estimates of Ordinary Cokriging provided the multi-temporal distribution of the aquifer salinization process.

Among the limitations of the analysis carried out, which have conferred uncertainty in the evaluation of the space-time variability, it is worth highlighting: the chemical low-frequency monitoring campaigns that have not allowed a continuous recording of the phenomenon, and the non-coincidence of the monitoring locations in the three survey campaigns.

A more accurate analysis with more robust results would make it possible to define an optimal uniformly distributed network of survey points, coupled with regular and more frequent monitoring of groundwater even along the vertical dimension. Hopefully, continuous monitoring of salinity at least should be foreseen.

Such a monitoring network would allow for a physically-based management of groundwater and surface water resources with the possibility of modulating demand with availability, considering the risk of salinization but also the possibility of blunting the salinization phenomenon through irrigation with fresh water from canals. The availability of data in space and time would enable the use of density-dependent mathematical modelling aimed at rational and efficient deep and surface water management.