An Integrated Approach for Deciphering Hydrogeochemical Processes during Seawater Intrusion in Coastal Aquifers

Abstract

For managing the freshwater in the worldwide coastal aquifers, it is imperative to understand the hydrogeochemical processes and flow patterns in the mixing freshwater/saltwater zone. The Egyptian Nile Delta aquifer is a typical example. The management of seawater intrusion (SWI) requires detailed investigations of the intrusion wedge and the dynamic processes in the mixing zone. Thus, a multidisciplinary approach was applied based on holistic hydrogeochemical, statistical analysis, and DC resistivity measurements to investigate the lateral and vertical changes in groundwater characteristics undergoing salinization stressor. The results of cross plots and ionic deviations of major ions, hydrochemical facies evolution diagram (HFE-D), and seawater mixing index (SMI) were integrated with the resistivity results to show the status of the SWI where the intrusion phase predominates in ~2/3 of the study are (~70 km radius) and the compositional thresholds of Na, Mg, Cl, and SO₄ are 600, 145, 1200, and 600 mg/L, respectively, indicating that the wells with higher concentrations than these thresholds are affected by SWI. Moreover, the results demonstrate the efficiency of combining hydrogeochemical facies from heatmap and resistivity investigations to provide a large-scale characterization of natural and anthropogenic activities controlling aquifer salinization to support decision-makers for the long-term management of coastal groundwater.

1. Introduction

Coastal aquifers undergo several stressors due to global climatic changes, sea-level rise (SLR), the world’s growing population, and increasing developmental activities [1,2,3,4]. The most common stressor is groundwater salinization that is driven by water level decline and SLR [5,6,7,8]. The sources of salinity in coastal aquifers are manifold, including invasion of seawater to the fresh groundwater [9] and mixing with highly saline fossil water [10] or with wastewater [11,12]. Accordingly, coastal aquifers require costly effective management plans to limit freshwater deterioration.

In arid and semi-arid regions, seawater intrusion (SWI) is considered a global issue that restricts groundwater pumping in coastal aquifers. Therefore, understanding complex hydrogeochemical processes and factors that control the evolution of brackish groundwater in coastal aquifers as well as the sources of salinization is crucial. However, the great complexity of coastal aquifers in terms of aquifer heterogeneity, spatial and temporal variations in the flow regime, and surface water–groundwater interactions complicate the interpretation of the controlling factors in fresh–saline water distribution [13]. Numerous studies were implemented to assess SWI in coastal aquifers using different methods—e.g., hydrogeochemical investigations [14,15], geophysical methods [16,17,18,19,20], isotope tracers [21,22,23], and statistical analysis [24,25]. SINTACS, DRASTIC, AVI, and GALDIT models are popular methods for mapping coastal aquifer vulnerability to natural and anthropogenic factors based on hydrogeological and hydrogeochemical parameters (e.g., aquifer type, hydraulic conductivity, water level, electrical conductivity, chloride concentration, Revelle index (Cl⁻/HCO₃⁻ + CO₃⁻²)) [26,27,28,29]. SINTACS, DRASTIC, and AVI are mostly used to determine aquifer vulnerability to surface contamination as wastewater, while the GALDIT is mainly related to seawater intrusion [27,28]. The Nile Delta aquifer is among the foremost important aquifers in the world; however, it suffers from SWI from the northern parts at the Mediterranean Sea because of severe pumping, especially in the central parts of the Nile Delta area [30,31,32]. Sherif et al. [31] demonstrated that seawater has intruded into the Nile Delta aquifer about 100 km from the coast. The seawater interface toe location was at ~23 km in the eastern Nile Delta before the 1950s when pumping from the Nile Delta aquifer was minimal; then, in the 1990s with pumping at 1.92 billion m³/yr [33], the interface surface moved 40 km landward. In the year 2000, with the pumping of 2.4 billion m³/yr, the interface moved 116 km landward from the coastal line in the eastern Nile Delta [31,32].

Previous studies on the Nile Delta aquifer system have involved hydrogeologic and hydrogeochemical studies dealing with SWI [15,30,34,35]. These studies described the extent of SWI into the aquifer, in addition to using vulnerability mapping and mathematical simulation for seawater encroachment to the aquifer. However, there is still a gap in understanding the hydrochemical processes during SWI into the aquifer. From this point of view, the present work aimed to understand the hydrogeochemical processes occurring during the groundwater–saline water mixing process as well as the origin of solutes and sources of salinity. To fill this knowledge gap, a holistic study using hydrogeochemical indicators, multivariate statistical analysis, and geophysical mapping was proposed to delineate and explain the chemical changes of the groundwater due to SWI.

The assessment of aquifer vulnerability to saltwater intrusion or pollutants using water sampling from wells is a reliable tool for monitoring water quality due to the seawater intrusion based on the total depth of the collected samples. However, unplanned wells with different depths cannot configure the continuity of the freshwater/saltwater full interface and may be difficult when attempting to understand the intrusion mechanism in the coastal aquifers. In such cases, surface geophysical methods such as the vertical electrical resistivity (VES) offer a potentially valuable, non-destructive, and cost-effective tool for mapping the seawater interface, where the electrical properties of the freshwater and seawater show large contrast in resistivity values.

The proposed approaches include hydrogeochemical analysis using ionic deviations of major ions, hydrochemical facies evolution diagram (HFE-D), groundwater quality index for seawater intrusion (GQI_SWI), seawater mixing index (SMI), and multivariate statistical analysis. Additionally, geophysical measurements were employed to validate the hydrogeochemical results and to delineate the vertical boundary of the fresh/saltwater interface. Finally, the geospatial distribution of the different parameters was addressed for visual investigation. The integration between these approaches will add valuable information and understanding of SWI in the Nile Delta aquifer system and can be applied in similar circumstances elsewhere as a good qualitative practice for coastal aquifers management.

1.1. Study Area

The study area is represented by the Sharkia Governorate in the Eastern part of the Nile Delta of Egypt. It is located between longitudes 31°15′30.7″–32°10′28.71″ E and latitudes 30°9′17.84″–31°10′2.91″ N (Figure 1) and is characterized by arid to semi-arid climate and has low relief surface sloping gently toward the north direction [15]. The ground surface is covered mainly by deltaic sediments that are used in different agricultural activities and dissected by a large network of freshwater canals and drainage systems [36]. It is bounded from the north by El-Manzala Lake, which represents the main source of SWI to the Nile Delta aquifer, especially with severe pumping of the groundwater [34,36].

1.2. Geological and Hydrogeological Setting

The study area is covered mostly by the Quaternary sediments with varying thickness. It is composed of sand, silt, gravel, clay, silty, and sandy clay deposits with lateral discontinuity [37]. The Quaternary sediments are underlain by the older Pliocene clays in most parts of the study area, while it is underlain by the Miocene and Oligocene deposits in the east and south parts [37]. These deposits form the Nile Delta aquifer that represents the main groundwater-bearing unit in the study area (Figure 2).

The Nile Delta aquifer is classified into an Upper Holocene clay cap unit and a Lower Pleistocene aquifer unit, which represents the main groundwater storage unit in the study area. The Pleistocene unit is classified based on lithology into Late and Early Pleistocene units that vary in thickness from 200, 500, and 800 m in the south, central, and northern parts of the study area, respectively [39] (Figure 2a). The Late Pleistocene aquifer is overlain directly by the Holocene clay cap unit and formed of medium sand with intercalations of clay lenses with an average thickness of 50 m. It is underlain by the Early Pleistocene aquifer unit, which is composed of coarse sand and gravels with clay intercalations. These two units are hydraulically connected and not easily differentiated, especially when the Holocene clay cap is missing in some parts of the Nile Delta [37]. The Holocene layer is composed of silt, clayey sand, and sandy clay, with fine sand intercalations. Its thickness varies from about 5 m and 70 m in the southern and northern parts, respectively (Figure 2b,c). It acts as an aquitard layer in the southern and central parts, while it acts as an aquiclude layer in the northern parts. The aquifer is recharged through the seepage from the distributed network of irrigation canals cutting the Holocene layer, as well as the Ismailia Canal that runs in the Pleistocene sediments. The groundwater levels in this aquifer decreased from southern (water level: 8–14 m asl) to northern parts (water level: <1 m asl) of the study area, with general flow direction from the south to the north and northeast directions [40] (Figure 2b).

2. Materials and Methods

Groundwater samples were collected through the year of 2018 from 60 drilled wells with variable total depths (18–120 m below the ground surface). Each well was pumped for about 20 to 30 min before sampling. On-field measurements were made after stabilization of the meter readings to determine the total dissolved solids (TDS), electrical conductivity (EC), and pH, using standard calibrated Orion portable meters. Each groundwater sample was added to polyethylene bottle after cleaning well with water sample at the sampling location. The concentration of HCO₃⁻ was measured by titration method using H₂SO₄ (0.01 normality). The concentrations of the other major anions and cations were detected using ion chromatography (Dionex, ICS-1100) (Table S1).

The hydrogeochemical processes prevailing during SWI were investigated using the groundwater quality index for seawater intrusion (GQI_SWI), ionic deviations, seawater mixing index (SMI), Piper diagram, hydrochemical facies evolution diagram (HFE-D), binary diagrams, and statistical analysis. The obtained results were further validated using geophysical tools.

2.1. Groundwater Quality Index for Seawater Intrusion (GQI_SWI)

The GQI_SWI is an index developed by [41] to facilitate the geospatial analysis of groundwater quality (in terms of seawater salinization) rather than traditional pattern diagrams. It is calculated according to [42] as in Equations (1)–(3):

$$GQ{I}_{SWI}=\frac{GQ{I}_{Piper\left(mix\right)}+GQ{I}_{fsea}}{2}$$

(1)

GQI_Piper(mix) = (Ca²⁺ + Mg²⁺ + HCO₃⁻) × 50 (in % meq/L)

(2)

GQI_fsea = (1 − f_sea) ×100

(3)

where GQI_{Piper (mix)} is the freshwater seawater mixing index of the piper diagram; GQI_fsea is seawater fraction index.

The seawater fraction (f_sea) is calculated using the Cl⁻ concentration of the sample, where Cl⁻ is considered as a conservative parameter. Accordingly, the fraction of seawater is calculated as follows:

$$f_{sea}=\frac{m_{Cl}{,}_{sample}-m_{Cl}{,}_{fresh}}{m_{Cl}{,}_{sea}-m_{Cl}{,}_{fresh}}$$

(4)

If Cl⁻ is only derived from seawater, then m_{Cl, fresh} should equal to zero so that

$$f_{sea}=\frac{m_{Cl}{,}_{sample}}{566}$$

(5)

where the denominator is equal to the concentration of Cl⁻ in mmol/L of 35‰ (grams of salt per kilogram) seawater [41,43].

2.2. Ionic Deviations (m_{i, react})

The chemical reactions taking place during fresh/seawater interaction could be deduced by calculating the expected ions’ compositions based on conservative mixing of saltwater and freshwater, and then comparing the result with the measured concentrations in the water sample [43]. Thus, the concentration of an ion (i) by conservative mixing of seawater and freshwater is

$$m_{i,mix }=f_{sea}\cdot m_{i,sea}+\left(1-f_{sea}\right)\cdot m_{i,fresh}$$

(6)

where m_i is the concentration of i (mmol/L), f_sea is the fraction of seawater in the mixed water, and subscripts mix, sea, and fresh indicate a conservative mixture, seawater, and freshwater, respectively. Any change in the concentration m_{i, react} due to reactions then becomes simply

$$m_{i_{,}react }=m_{i,sample}-m_{i,mix}$$

(7)

Thus, the calculated difference or ionic deviations (m_{i, react}) between the mixed (m_{i, mix}) and the measured compositions indicates the possible reactions [43] where the Cl is used in these calculations due to its conservative nature. The use of ionic deviations in hydrogeochemical modeling better explains the possible reactions rather than the original concentrations.

The ionic deviations of the major ions were used considering sample 57 as the fresh groundwater endmember where it is characterized by Ca-HCO₃ water type, lowest Cl, and lower TDS concentrations and is representative of the fresh surface water of the River Nile and the Ismailia freshwater canal (Table S2).

2.3. Hydrochemical Facies Evolution Diagram (HFE-D)

This diagram was used to study the hydrogeochemical changes in groundwater experiencing salinization–freshening phases in coastal aquifers [44,45]. It is used as an alternative to other hydrochemical diagrams, yet it has more benefit as it helps in the identification of the hydrochemical facies without amalgamations (parts of the diagram necessarily representing two or more different ions, e.g., Piper and Durov diagrams). In addition, it permits visualization of possible groupings of samples and trends of evolution occurring in the aquifer [46,47,48,49]. In this diagram, the abscissa reproduces cation percentages and their evolution, considering the two most important cations (%Na and %Ca (or %Mg)) that interpret base-exchange reactions while the ordinates denote the mixing process, where %HCO₃ (for recharge waters with a bicarbonate facies) or %SO₄ (for recharge water with a sulfate facies) indicate the dominant anion in the recharge water. Details about the different hydrochemical facies in the HFE diagram are given in [44], where four heteropic facies—NaCl, CaCl₂, NaHCO₃, and CaHCO₃—are defined and each heteropic facies is represented by four other subfacies that are similar to them. Accordingly, sixteen isopic facies could be extracted from the diagram. The main hydrochemical processes occurring during seawater intrusion could be recognized by the excess or deficit compared with the line of conservative mixing between freshwater and seawater (conservative mixing line (CML)) in the freshening and intrusion stages, respectively.

The isopic facies that are related to freshening substages are f1 to f4, while those of the intrusion substages are i1 to i4. The i1 and f1 substages characterize the distal facies that define the boundary between the two invasive flows where this boundary must not be mistaken with the conservative mixing line represented in the HFE diagram. On the contrary, the proximal facies are f4 and i4 and are closely linked to the end members FW and SW, while the substages (i2, i3, f2, and f3) represent the intermediate levels where cation exchange is usually more evident. The samples sited in the freshwater field are considered as belonging to the freshening phase, although they are in the intrusion field. Similarly, samples located in the saltwater field represent the intrusion stage, even if they are in the freshening field [44].

The endmembers of seawater (SW) and freshwater (FW) can be theoretical or actual [50]. The FW endmember should be critically selected as all processes are evaluated from the starting point of the seawater–freshwater mixing line, and the only variable is the recharge water because seawater is considered stable [51].

2.4. Seawater Mixing Index (SMI)

The seawater mixing index (SMI) estimates the degree of seawater mixing [52], where the major ionic constituents of seawater (Cl⁻, SO₄²⁻, Na⁺, and Mg²⁺) were used to calculate this index (Equation (8)).

$$SMI=a\times \frac{C_{Na}}{T_{Na}}+b \times \frac{C_{Mg}}{T_{Mg}}+c \times \frac{C_{Cl}}{T_{Cl}}+d\times \frac{C_{{SO}_4}}{T_{{SO}_4}}$$

(8)

where a, b, c, and d are constants representing the relative concentration proportion of Na⁺, Cl⁻, Mg²⁺, and SO₄²⁻ in seawater where their values are 0.31, 0.04, 0.57, and 0.08, respectively. C represents the measured concentration of these ions in mg/L, where T signifies the values of regional threshold of the measured ions, which can be calculated by using the cumulative probability curves and where the inflection points in the curves were used for the concentrations of Cl⁻, SO₄^2−, Mg²⁺ and Na⁺ in groundwater samples [53]. Samples of SMI value greater than 1 characterize groundwater affected by seawater–freshwater mixing process, where SMI values lower than 1 signify freshwater.

2.5. Base Exchange Index (BEX)

The BEX was used to assess a dolomite-free aquifer system experiencing freshening or salinization; Stuyfzand [54] defined this index as BEX = Na⁺ + K⁺ + Mg²⁺ − 1.0716Cl⁻ (meq/L). A positive BEX represents freshening, a negative BEX represents salinization, and a BEX of zero indicates no base exchange.

2.6. Sodium Adsorption Ratio (SAR)

This parameter is used as an indicator of the salinity hazard. High values imply a hazard of the water to the soil. It represents the amount of sodium relative to calcium and magnesium in the water extract from saturated soil (Equation (9)).

SAR = Na/(Ca + Mg/2)^ 0.5

(9)

2.7. GIS Spatial Mapping

The measured data were used to build a geodatabase inside ArcGIS 10.3 software (Esri, Redlands, CA, USA) [55]. The locations of the collected samples were geo-referenced, converted to shapefile format, and interpolated to produce spatial distribution maps of the studied parameters in a raster format. These maps were produced using the Kriging interpolation method inside GIS Spatial Analyst Tool [55,56,57]. Kriging is a multistep procedure that comprises exploratory statistical data analysis, variogram modeling, and surface creation (Figure S1). Kriging is an advanced and accurate interpolation geostatistical procedure, and it is commonly used in soil, water resources, and geology sciences [55].

2.8. Geophysical Measurements

Electrical resistivity can be used as an indicator of SWI, where electrical conductivity (the inverse of electrical resistivity) is based mainly on the amount of TDS in groundwater. In the coastal aquifers as in the Nile Delta, the low resistivity anomalies of less than 2 Ωm indicate a seawater intrusion zone [16]. The authors of [27,34] carried out a DC resistivity survey to study the effect of pumping on SWI in the eastern and central Nile Delta, respectively. Their results categorized the groundwater into three classes based on the TDS concentration and resistivity values: 1) freshwater (35–100 Ωm), brackish water (5–35 Ωm), and saltwater (1–8 Ωm). Accordingly, the surficial DC resistivity surveys, accompanied with limited hydrochemical data are widely used for mapping subsurface problems relevant to groundwater quality to describe the spatial distribution of freshwater and saltwater in coastal regions [27,58,59,60,61]. Reasonably, the vertical electrical sounding (1D-VES) measurements were collected at 49 sites using Schlumberger array to validate the existing freshwater boundary identified through the hydrochemical facies modeling and calculated SMI at different depths. The VES technique outlines the vertical variation of the electrical resistivity and can easily detect the freshwater–saltwater interface compared with the traditional water sampling [62]. Basically, during the field work, an electrical current was released into the ground between a pair of steel electrodes (A and B), and voltage was measured between another two electrodes to estimate the subsurface electrical resistivity distributions. The subsurface electrical resistivity has generally been deployed to identify subsurface layer distributions and groundwater potentialities at depths that vary between a few meters up to hundreds of meters. The maximum current spacing (AB) ranges from 300 m near El Manzala Lake in the north to 800 m in the southern part of the surveyed area. According to the available geological information, the applied spreading is long enough to penetrate the expected fresh/seawater interface. Some of the measured soundings were conducted close to some existing shallow wells to define the parametric resistivities of the subsurface layers and the depth to the water table. The collected apparent resistivities and AB/2 measurements were inverted into layer resistivities and thicknesses using the IPI2WIN 3.01software [63]. The inversion results show the vertical variations in the soil and water quality; therefore, it can map the interfaces of freshwater and saltwater more precisely, as seen in Figure 3.

3. Results and Discussions

3.1. Groundwater Classification and Salinity

The studied groundwater samples were classified on the basis of GQI_SWI values (Table S1) into highly fresh ‘FGH’ (GQI_SWI > 80; 4 samples), slightly fresh ‘FGS’ (GQI_SWI from 75 to <80; 11 samples), slightly mixed ‘MGS’ (GQI_SWI from 70 to <75; 5 samples), fairly mixed ‘MGF’ (GQI_SWI from 65 to <70; 11 samples), highly mixed ‘MGH’ (GQI_SWI from 60 to <65; 8 samples), extremely mixed ‘MGE’ (GQI_SWI from >50 to <60; 12 samples), and saline groundwater ‘SG’ (GQI_SWI from 10–50; 9 samples) (Table S1). For simplicity, we used acronyms to describe the different groundwater classes where FG, MG, and SG are for fresh groundwater, mixed groundwater, and saline groundwater, respectively. The third letter of the acronym denotes the first letter of the water class. The groundwater salinity distribution is spatially comparable with the GQI_SWI distribution (Figure 4). Additionally, the composition of the groundwater samples was investigated using a Piper diagram [42], where it evolves from Ca-HCO₃ and Na-HCO₃ (in the south) to mixed Mg-Na-Ca-Cl-HCO₃ and mixed Na-Ca-Mg-HCO₃-Cl, to Na-Cl and Na-Cl-SO₄ (in the central parts), and finally to Mg-Na-Cl water types (in the north) (see zones A, B, and C in Figure 5). The dominant groups, Na-Cl and mixed Mg-Na-Ca-Cl-HCO₃, dominate the northern and central parts of the study area, where the Ca-HCO₃ and Na-HCO₃ water types characterize samples of the southern parts, where it has very low concentrations of TDS and Cl. In general, Na-Cl and Mg-Na-Cl-type waters reflect the higher influence of seawater in intrusion in coastal aquifers [43,64]. Additionally, the conservative mixing line between sample 57 (freshwater endmember) and seawater (the Mediterranean Sea water of [65,66]) was shown on the Piper diagram (Figure 5), where the plotting of samples above this line indicates the effect of salinization, while those below this line indicate freshening processes.

3.2. Descriptive Statistics and Coefficient of Variation

The studied variables show very wide ranges and high standard deviations (

Table 1. Descriptive statistics and coefficient of variation (CV) of hydrogeochemical data (concentrations unit is ppm; EC is in μS/cm).

	Minimum	Maximum	Mean	Std. Deviation	CV%
pH	6.9	8.9	7.6	0.4	5
TDS	283	19550	3595	5397	150
EC	389	30030	5429	8135	150
K+	4.5	146	24.3	26.5	109
Na+	16	4350	734	1131	154
Mg2+	13	960	146	242	166
Ca2+	15	1554	227	366	161
Cl−	30	10120	1459	2557	175
SO42−	15	4560	654	1147	175
HCO3−	72	810	227	129	57

). In particular, the TDS and EC values have a wide range, with means 3595 ppm and 5429 μS/cm and standard deviations of 5397 and 8135, respectively (Table 1). Similarly, the ranges and standard deviations of K, Na, Mg, Ca, SO₄, and HCO₃ ions are also wide, indicating that complex hydrogeochemical processes and/or multiple sources could be responsible for the composition of these waters (Table 1).

To confirm the variable sources of ions, the coefficient of variation (CV) was calculated for the studied variables (Table 1). Generally, it is applied to express the precision and repeatability of a geochemical analysis [67]. The CV represents the ratio of the standard deviation (SD) of an element to its average concentration, giving a standardized measure of the dispersion of a probability or frequency distribution [68]. Generally, a CV near zero indicates that the population is homogeneous (Table 1). Accordingly, the TDS, EC, and the major ion concentrations have higher CV values representing the heterogeneity of these variables and confirming multiple sources of mineralization. Additionally, the log-normal distribution of K, Ca, Mg, Cl shows a tail at the high concentration ranges, while most ions show a polymodal distribution, where the population of higher concentrations is possibly attributed to seawater mixing (Figure 6).

3.3. Pearson’s Correlation and Hydrogeochemical Reactions

Pearson’s correlation analysis [69] assesses the correlation between different variables in a system. Several studies focused on the ion exchange of major cations (Na⁺, K⁺ Ca²⁺, Mg²⁺) in coastal aquifers [43,70,71], where the major hydrogeochemical processes could be identified. Additionally, cross plots of major ions were used to decipher the hydrogeochemical reactions operating in aquifers [72,73,74,75]. In this study, integration between the results of Pearson’s correlation and cross plots were made to study the possible sources of ions in the studied groundwater samples (Table S3; Figure 7). Generally, when seawater intrudes fresh groundwater, Na⁺ is adsorbed on clay surfaces in the aquifer matrix and Ca²⁺ is released into the water leading to change of Na-Cl water type to the Ca-Cl₂ water type [64]. On the contrary, aquifer freshening is indicated by Na-HCO₃ water type, where ion exchange between Ca in the water and Na on the clay matrix occurs.

The correlation was applied to selected variables (total depth, water level, pH, EC, distance from shoreline, and Cl), saturation indices (SI_Hal, SI_Gyp, SI_Dol, SI_Cal), SAR, BEX, and the calculated ionic deviations (m HCO₃-react, m SO₄-react, m Ca-react, m Mg-react, m Na-react, m K-react) (Table S3). A correlation factor r ≥ 0.5 indicates a significant level of correlation and is used in this work to preliminarily evaluate the main source of salinity based on the used variables [76,77].

The Cl contribution to groundwater salinity is evidenced by the strong linear relationship between Cl and TDS (Table S3, Figure 7a). However, the samples of the MGE and SG lie slightly above the theoretical mixing line between freshwater (sample 57) and seawater (Figure 7a), indicating that the salinity of groundwater is derived not only from marine sources but also from other sources (e.g., anthropogenic sources from drains and return flow of evaporated surface water (Figure 7a). Na⁺, K⁺, and Mg²⁺ correlate well with Cl⁻ for most groundwater types at Cl⁻ < 150 epm (Figure 7b–d). At higher Cl⁻ levels (the SG water samples), Na⁺ and Mg²⁺ show deficit and surplus, respectively, suggesting reverse ion exchange prevailing during the seawater intrusion. This is also confirmed by the strong negative correlation between EC and Cl with reacting sodium and potassium (Table S3). Few samples plot along the 1:1 line, indicating possible halite dissolution (Figure 7b), while most samples lie slightly above and below the 1:1 line, indicating the possibility of silicate weathering and/or anthropogenic effects, respectively [78], (Figure 7b). Halite dissolution could be due to the occurrence of several clay lenses in the Quaternary aquifer [37]. Additionally, the variable sources of Na⁺ and Cl⁻ could also be deciphered from the relation between Na⁺ + K⁺ and Cl⁻+ SO₄²⁻ (Figure 7k). In this relation, the continuous increase in these ions along the 1:1 suggests a geogenic origin [79]; however, the concentrations of these elements in the studied samples do not fall along the 1:1 line (Figure 7k), indicating other sources. Mg²⁺ and Ca²⁺ show excessive concentrations for the MGF, MGH, MGE, and SG, where they deviate from the theoretical mixing line between seawater and freshwater (Figure 7e), indicating reverse ion exchange reactions. However, reverse ion exchange is not the only process that drives the hydrochemistry of the mixed groundwater, where the majority of the mixed type groundwater lies below 1:2 line in the relation between HCO₃⁻ and Ca²⁺, indicating possible CO₂ dissolution according to reaction 10 or organic matter oxidation that derives the production of CO₂ or due to refreshing recharge of surface water (Figure 7h).The increased Ca²⁺ and HCO₃⁻ concentrations in the mixed- and saline groundwater derive the saturation index of calcite and dolomite to positive values that could lead to their precipitation (see the scatter plot in Figure 7f and Figure 8a,b; SMI in the ordinate of Figure 8 will be discussed later). The SO₄²⁻ concentrations for most samples are exceeding those expected from mixing of seawater and freshwater (Figure 7g), indicating multi-sources.

$$C{O}_2+H_{2}O \leftrightarrow H_{2}C{O}_3\leftrightarrow H^{+}+HC{O}_3{}^{-}$$

(10)

When SO₄²⁻ concentrations in groundwater exceed 100 mg L⁻¹, the possible sources are attributed to anthropogenic inputs and/or lithogenic sources (sulfate evaporites, pyrite or organic matter oxidation) [80]. However, gypsum dissolution and iron sulfide oxidation are precluded for the MGE and SG and the possible sources could be related to anthropogenic impacts (Figure 7i,j). Additionally, the ionic deviations of SO₄²⁻ increase in the MGE and SG as the distance from the shoreline decreases as well as the well depths decrease proposing mixing with seawater (Figure S5c,k). Conversely, samples from distant and deep wells (more than 50 m depth) have low reacting SO₄²⁻ contents in the mixed groundwater, indicating that the well depth seems not to be influencing the SO₄²⁻ content in the mixed and fresh groundwaters (Figure S5c,k). The FG and the MGS seem not to be potentially affected by anthropogenic sources of sulfate where the samples plot above the 1:1 line (Figure 7i). Some samples of the MGS and MGF lie along the 1:1 line, suggesting possible evaporite dissolution that is associated with clay lenses in the Quaternary aquifer [37] (Figure 7i). The increase in Ca and SO₄²⁻ during aquifer salinization drive the saturation index of gypsum to positive values (Figure 8c, Table S2).

The water level shows negative correlation with SI_Hal (−0.71), SI_Gyps (−0.67), SAR (−0.57), Cl (−0.56), and EC (−0.56), which explains the dynamics during seawater mixing, as during high water levels, all of these parameters decrease and vice versa (Figure S6e,f,g). Generally, the chemistry of the SG (represented by Na⁺, Mg²⁺, and SO₄²⁻) is significantly affected by the distance from shoreline, water levels, and well depths. On the other hand, the chemistry of the rest of groundwater samples (represented by Na⁺, Mg²⁺, and SO₄²⁻) show no definite trend with shoreline, water levels, and well depths (Figure S2).

Another interesting result is the negative correlation between the BEX and Cl (−0.63), BEX and reacting Ca (−0.65), and the positive correlation between BEX with Na and K (0.93, 0.67, respectively), where a positive BEX represents freshening and a negative BEX represents salinization and confirming the cation exchange processes prevailing during aquifer salinization (Supplementary Table S3).

3.4. Factor Analysis

Factor analysis (FA) is a multivariate statistical technique used to quantify the interrelationships among a large number of variables into a smaller set of dimensions with a minimum loss of information [80]. It is used to identify important factors contributing to the groundwater hydrochemical processes [81,82,83,84]. The correlation matrix was created based on standardized data sets. Then, the eigenvalues and factor loadings of the correlation matrix were defined. Factors were extracted when eigenvalues were greater than 1 [85,86].

The sampling adequacy and the feasibility of applying the FA were tested using Bartlett’s sphericity tests and the Kaiser–Meyer–Olkin (KMO) [87] (

Table 2. The Kaiser–Meyer–Olkin (KMO), Bartlett’s sphericity tests, and the factor analysis for the analyzed groundwater samples.

Kaiser–Meyer–Olkin Measure of Sampling Adequacy				0.79
Bartlett’s Test of Sphericity		Approx. Chi-Square		2604.09
		df		153
		Sig.		0.000
Rotation Sums of Squared Loadings	Total	8.96	3.41	3.11
	% of Variance	49.78	18.95	17.29
	Cumulative %	49.78	68.72	86.02
		1	2	3
	EC	0.97	0.22	−0.01
	TDS	0.97	0.22	−0.02
	Mg	0.97	0.16	−0.01
	Cl	0.95	0.21	−0.04
	Na	0.95	0.26	−0.02
	Ca	0.95	0.15	−0.02
	SO4	0.89	0.21	−0.05
	SI Gypsum	0.75	0.54	0.12
	SI Halite	0.73	0.62	0.05
	K	0.73	0.19	0.19
	SAR	0.69	0.51	−0.02
	Distance	−0.40	−0.85	−0.08
	Total Depth	−0.05	−0.83	0.26
	Water level	−0.35	−0.83	−0.11
	SI Dolomite	0.22	0.05	0.96
	SI Calcite	0.19	0.00	0.96
	HCO3	−0.08	0.02	0.81
	pH	−0.40	−0.16	0.69

). The best adequacy results were obtained using eighteen variables: distance from the shoreline, total depth, water level, pH, EC, TDS, SAR, K⁺, Na⁺, Mg²⁺, Ca²⁺, Cl⁻, SO₄²⁻, HCO₃⁻, and saturation indices of halite, gypsum, dolomite, and calcite. Accordingly, the KMO test result was 0.79, confirming the adequacy of PCA in reducing the dimensionality of the dataset due to inter-correlation between variables. Additionally, the results from Bartlett’s sphericity test (chi-square = 2604.09; degrees of freedom ‘df’ = 153; and p < 0.001 “=zero”) confirm that there is common variance shared among the studied variables. Based on eigenvalues (more than one) and varimax rotation, three factors were discerned to explain most of the variability, with a total cumulative variance of about 86% for the groundwater samples (Table 2 and Figure S3).

Factor 1 accounts for about 50% of the total variance. It is represented by a strong positive loading on TDS, EC, Mg²⁺, Cl⁻, Na⁺, Ca²⁺, SO₄²⁻, K⁺, SAR, and saturation indices of halite and gypsum (Table 2). This indicates that the salinity of groundwater is attributed to manifold sources (seawater mixing and possible anthropogenic sources as evidenced by high SO₄²⁻ and Cl⁻).

Factor 2 explains about 19% of the total variance, in which a strong negative loading is found on distance (−0.85), total depth (−0.83), and water level (−0.83). Moderate positive loading on SI_Hal (0.62) and SI_Gyp (0.54) exist (Table 2). This factor signifies that under conditions of proximity to the shoreline and the existence of shallow wells and lower water levels, the composition of the groundwater changes, being enriched in some elements that derive the saturation indices of gypsum and halite to positive values, representing the effect of seawater mixing on the groundwater quality.

Factor 3 accounts for about 17% of the total variance, where a strong positive loading on SI_Dol (0.96) and SI_Cal (0.96), HCO₃⁻ (0.81), and pH (0.69) exist. This factor represents the freshening and water–rock interaction factor that characterizes sample populations close to the Ismailia freshwater canal due south of the study area.

3.5. Spatial Distribution of Groundwater Facies

The seawater contribution in the studied samples ranged from 0% in sample 57 to ~48% in sample 33 (Table S2). The plot of the calculated ionic deviations for Na and Ca (and Mg) against chloride concentrations shows the onset of cation exchange reactions in the aquifer at Cl concentrations (more than 500 ppm, Figure 9). The degree of seawater contribution in the groundwater samples was studied using an HFE diagram and SMI. In this diagram, sample 57 was used to represent the FW endmember, while the Mediterranean water was used as the SW endmember [65,66].

The distribution of hydrogeochemical facies of the studied samples in the HFE-D shows that water samples influenced by SWI are represented by red color, while samples representing freshwater are in green color (Figure 10a). Ten hydrochemical substages are recognized representing the freshening stage (24 samples representing 40.7% of the total samples) and intrusion stage (35 samples representing 59.3%). Details on the percentages and facies of each substage are given in (Table S4). The freshening phase of the studied samples starts from an initial situation in which the aquifer is totally salinized by seawater ingress. Going from inland to the coast, the recharge water from surface canals begins to rinse the aquifer that is in equilibrium with the Na-Cl facies and is represented by the distal freshening facies ‘f1’ (MixNa-Cl, representing MGF, MGH, and MGE; samples 58, 2, 4, 32, 13, 18, 26, 30, 31; see Table S1 for classification, Figure 10a). After that, the water evolves to substages ‘f2’ (MixNa-MixCl, representing MGF, MGH, and MGE and encompassing samples 25, 54, 55, 17, 56) and ‘f3’ (MixCa-MixHCO₃, Ca-MixHCO₃, representing FGS that is represented by samples 9, 59) and finally to substage ‘f4’ that represents the proximal freshening facies ‘f4’ (Ca-HCO₃, representing FGH and FGS of samples 38, 10, 21, 60, 7, 8) (Tables S1 and S4, Figure 10a). Ultimately, the aquifer total recovery is recognized by the heteropic Ca-HCO₃ facies representing FGH (samples 36 and 40) (Tables S1 and S4).

Conversely, in the intrusion phase when the aquifer is filled with fresh Ca-HCO₃ water, the water increasingly progresses through substage ‘i1’, representing Ca-MixHCO₃ facies of the distal intrusion (represented by the FGS, sample 37, (Table S1 and Table S4; Figure 10a)). After that, the water evolves to substage ‘i2’ (MixCa-MixCl, Ca-MixCl, represented by the FGS and MGS, samples 14, 20, 52, 53, 24). Later, the water evolves to substages ‘i3’ (MixNa-Cl, MixCa-Cl, Ca-Cl₂, representing MGF, MGH, and MGE, samples 6, 12, 19, 22, 39, 3, 28, 50, 27). At this time, some samples of the ‘i3’ substage lie away from the conservative mixing line representing the active cation exchange reactions occurring during aquifer salinization. These samples are represented by the MGF and MGS (samples 51, 11, 15, 23). The final substage is ‘i4’ that represents the proximal intrusion facies MixNa-Cl and is shown by sample 1 of MGF. Finally, in this phase, the main heteropic Na-Cl facies dominates the aquifer and is represented by MGE and SG (samples 44, 45, 46, 16, 33, 34, 35, 42, 48, 49, (Tables S1 and S4, Figure 10a)). A heatmap was constructed to show the spatial distribution of the different substages for every phase, showing the predominance of the intrusion phase in ~ 2/3 of the study area (Figure 10b).

3.6. Seawater Mixing Index (SMI)

The cumulative probability percentages of selected hydrochemical parameters were plotted against their logarithmic concentration value (Figure 11). The intersection point of each plot characterizes the threshold values (T) of the equivalent ion [52,53]. Accordingly, the determined T values of Cl, SO₄, Na, and Mg for the studied samples are 600, 145, 1200, and 600 mg/L, respectively. The threshold values and the relative concentration proportions for the ions in seawater (a = 0.31, b = 0.04, c = 0.57, d = 0.08) were substituted in Equation 8 to calculate the SMI values for each groundwater sample. If the SMI value was more than 1, the water may have been influenced by seawater mixing [52].

The SMI values show that all groundwater samples of the FG, MGS, MGF, MGH, and part of the MGE water types have SMI values smaller than 1, while part of the MGE and SG types have values larger than 1 (Figure 10c). Accordingly, 13 water samples (21.6%) are significantly affected by the mixing with seawater components, while the rest of the 47 samples (78.3%) are not (Table S2). The spatial distribution of SMI values for coastal groundwater in the study area is relatively comparable with the heatmap, where the northern parts of the study area are affected to variable degrees by SWI, while the central area in the heat map represents the dynamic zone that is affected by variable degrees of mixing with seawater, and the southern zone represents the fresh groundwater zone (Figure 10c).

In addition, a bivariate relation between SMI values and the distance from the shoreline, water level, and depth of each well is shown in Figure 12. A negative correlation exists between the SMI and distance from the shoreline, where the fresh groundwater shows nearly constant SMI with distance (the green trend line in Figure 12). On the other hand, the SMI of the mixed groundwater increases with decreasing distance (the orange trend line, Figure 12a). Samples with SMI values greater than 1 tend to be preferentially located within about 45 and 12 km from the shoreline in the west and east parts of the study area, respectively (see the blue vertical line in Figure 12a and the spatial distribution in Figure 10c)). This variation between the west and eastern parts could be attributed to the presence of scattered zones of perched water with lower salinities, as in samples 18, 28, and 55. However, the trend of SMI with water level shows that the FG has low ranges of SMI at higher water levels (green trend line), while at water levels <~5 m, the SMI increases (orange trend line in Figure 12b). The relation between SMI and well depth indicates that well depths seem not to be influencing the degree of mixing at shallow wells, as both high and low SMI exists (Figure 12c).

3.7. Ionic Deviations and SMI

The relation between the ionic deviations and SMI shows that Mg²⁺, Ca²⁺, and SO₄²⁻ increase gradually with increasing SMI for the FG and the MG classes (Figure 13a,b,e). This is attributed to the mixing processes in the central parts of the study area. The SG and some of the MGE show SMI above 1, indicating seawater intrusion in the northern parts close to El Manzala Lake (Figure 13a,b,e). Na⁺ and K⁺, of the SG and some of the MGE, are depleted dramatically at SMI more than 1, indicating the reverse ion exchange that operates during SWI (Figure 13c,d). No distinct trend of HCO₃⁻ is shown for the groundwater samples (Figure 13f).

3.8. Spatial Resistivity Distribution and Hydrogeochemical Facies

To validate the spatial distribution of the obtained hydrochemical facies and the calculated SMI, the spatial distribution of the interpreted resistivity from geophysical surveys was used in comparison with both maps, as illustrated in Figure 10d. This comparison shows that the resistivity distribution can be read in terms of the pattern of SWI and that it displays a sufficient similarity with HFE and SMI maps. Accordingly, the area could be differentiated into three zones with sharp boundaries based on resistivity and hydrochemical facies distributions.

The inspection of the resistivity map shows that the northern part of the aquifer of less resistivity (<2 Ωm) is highly impacted by SWI that hydraulically connected with El Manzala Lake. Reasonably, this zone is characterized by high Cl concentration (>1000 ppm), and the dominant hydrogeochemical facies is related to i4 (>66% Cl and >50% Na, Figure 10a) due to SWI. By contrast, the southern part of the area displays a maximum resistivity range (35–75 Ωm) with minimum Cl content, and the dominant facies is freshwater f4 (Ca-HCO₃), where the SWI is totally diminished to the south. Notably, the aquifer in the southern part is recharged from the Damietta Nile branch in the northwest and the Ismailia canal to the south (Figure 1). From the comparison of resistivity and hydrochemical maps (Figure 10), the mixing zone was clearly observed in the middle part with resistivity ranges from 8 to 30 Ωm where the area was moderately influenced by SWI to confirm the obtained mixed geochemical facies of FW and SW (i2, i3, f2, and f3; Figure 10a).

As mentioned before, the mixing zone represents the dynamic response of the aquifer to seawater/freshwater, and the resistivity values can change in lateral and vertical directions in a short distance, reflecting the rapid variations in water facies. In the resistivity map (Figure 10d), the local low resistivity anomalies in the middle parts may have been influenced by SWI as a result of intensive pumping that, accordingly, creates vertical up-coning of seawater from the deeper zone. By contrast, the southeastern area in this zone that showed high resistivity values has been influenced by freshwater from the Ismailia canal and irrigation networks. To understand the seawater intrusion mechanism with depth, vertical resistivity section A-A’ was constructed in a north–south direction perpendicular to the coastline (Figure 14a). For better deployment of the geophysical data, the SMI and HEF-D domains from seven wells with different depths were demonstrated in the upper part of the resistivity section, respectively (Figure 14a). In addition, the integrated results of geophysical and hydrogeochemical data were compared with the pervious results of saltwater simulation in the east and middle regions of the Nile Delta aquifer system [27,31,32,33]. Accordingly, the very low resistivity zone in the resistivity section (<2 Ωm) on the northern side extended in an approximately 75 km radius at a maximum depth of 100 m, representing the impact of SWI and forming a wedge-shaped bottom layer (Figure 14a). In this very high conductive zone, the main hydrogeochemical facies are related to i4 and i3, and the SMI is more than 1, indicating seawater saturation zone at that depth. The total depth of drilled wells in this zone varies between 20 and 33 m according to the distance from the shoreline; therefore, the TDS varies between 19,950 ppm at well 33 and 5136 ppm at well 16 towards the middle parts, with a very high concentration of Cl (Figure 9c). Here, the superiority of vertical resistivity measurements compared with the traditional water sampling could be noticed for mapping the brackish and seawater interface at a great depth of 120 m in the middle part of the area where the available wells do not exceed 35 m depth. In contrast and without the borehole information, the DC resistivity method failed to distinguish between the upper clay cap and the lower saltwater saturated sand in the area close to the coastal zone (>5 SMI value) due to the seawater intrusion at shallow depths. This challenge gradually decreases towards the middle parts where the hydrochemical facies changes from i4 to i3 and SMI values decrease from 5 to 1 (Figure 14a). The most interesting element along the cross section is the mixing zone between the extreme high resistivity values (>35 Ωm) in the south and very low values near the coastal line, forming a typical wedge shape of the seawater–freshwater interface. According to similar studies (i.e., 27, 32, 33, 34, and 59) for simulating seawater intrusion, the depth to wedge-shaped interface is influenced by the groundwater pumping (rate and locations) and freshwater recharge, in addition to the flow patterns towards the shoreline under the impact of SLR. Accordingly, the vulnerability of the aquifer to the seawater intrusion in the area is different from the west to east based on irrigation networks and groundwater abstraction, as illustrated in the SMI and resistivity maps (Figure 10c,d, respectively). In the mixing zone, the local resistivity anomalies with a limited thickness (10–20 m) indicate perched water bodies with low TDS (800–2000 ppm) and mixing hydrochemical facies due to the freshening of these lenses from the local irrigation system. These results are compatible with the SWI model as adapted by [31,32] in the eastern Nile Delta aquifer (Figure 14b). In this model, the sharp interface toe is positioned at 116 km (1000 ppm TDS) from the coast with a depth of 400 m [31]. The toe point was moved due to the intensive pumping where the current pumping rate in the mixing zone will not be sustainable and increases the vulnerability to SWI.

4. Conclusions

To manage coastal aquifers salinization due to the worlds’ growing population, global climatic changes, and increasing developmental activities, it is important to study the controls governing salinization. Comprehensive hydrogeochemical, statistical, and geophysical studies were performed in a coastal shallow aquifer ever known by its high vulnerability to seawater salinization (the Nile Delta aquifer, Egypt). The study included 60 groundwater samples collected from shallow and deep wells and analyzed for major elements. The groundwater was classified based on GQI_SWI into fresh, mixed, and saline groundwater—FG, MG, SG, respectively. Ionic compositions of these samples showed heterogeneity, confirming the presence of multiple sources of mineralization. The sources of solutes were studied using Pearson correlation and cross plots. The relation between EC and Cl versus Ca, SO₄, and Mg indicated that the increased salinity is associated with increasing these ions’ concentrations as a result of mixing with seawater. The ionic deviations of SO₄²⁻, Ca²⁺, and Mg²⁺ also increase in the extremely mixed and saline groundwater as the distance from the shoreline decreases as well as when the well depths decrease, confirming the mixing with seawater. Additionally, the results of factor analysis indicate that there are three factors controlling the elements’ distribution in the groundwater. These factors are: (a) seawater/anthropogenic, (b) salinization, and (c) freshening processes. A heatmap was constructed to show the spatial distribution of the compositional substages for freshwater, mixed water, and seawater phases, showing the predominance of the intrusion phase in ~2/3 (~70 km radius) of the study area. Additionally, the compositional thresholds of Na, Mg, Cl, and SO₄ for the studied samples are 600, 145, 1200, and 600 mg/L, respectively, indicating that wells of concentrations higher than these thresholds are affected by seawater intrusion. The relation between the seawater mixing index (SMI) and distance from the shoreline indicates that SMI values >1 (indicating salinization) tend to be preferentially located within about 45 and 12 km from the shoreline in the west and east parts of the study area, respectively. This variation between the western and eastern parts of the study area could be attributed to the presence of scattered zones of perched water with lower salinities. Additionally, the trend of SMI with water level shows that the fresh groundwater has nearly low ranges of SMI at higher water levels, while at water levels ~<5 m asl, the SMI increases. The resistivity distribution at different depths, interpreted from the surficial measurements integrated with available geological and hydrogeological information, successfully exposed the interaction of seawater and freshwater in the coastal aquifer. The geophysical results indicate that seawater mainly intruded inland up to 75 Km at 100 m depth from the shoreline. However, these results were restricted in the near-coast environment due to high seawater intrusion located close to the coastline. Therefore, the assessment of aquifer vulnerability to seawater intrusion using hydrochemical data based on traditional water sampling needs the combination of geophysical and hydrogeological data to evaluate seawater intrusion more precisely with depth.

The current results based on geophysical and hydrogeochemical modeling are consistent with previous simulation studies in the eastern and middle parts of Nile Delta to indicate that the vulnerability to seawater intrusion currently seriously threatens the Nile Delta aquifer. Reasonably, a suggested plan for controlling the invasion of seawater is to increase water levels >5 m asl as well as to maintain this level through controlled pumping to avoid aquifer deterioration. This could be attained through artificially recharging selected moderately saline wells in the mixing zone in the middle part of the study area in a radius of about 12–75 km (of resistivity 8–35 Ωm). Sources of recharge could be desalinized groundwater or treated wastewater. Additionally, an automated Internet of Things (IoT)-based monitoring system would be useful for selected parameters (e.g., EC, Ca, Cl, Br) to build a prediction model using continuous readings. Therefore, the results of this study could represent a base for building management systems for preventing the future advance of seawater wedge and may contribute to decision-making for the long-term management of coastal groundwater resources.