Simulation of seawater intrusion in the Nile Delta aquifer under the conditions of climate change

Coastal aquifers represent an important source of freshwater, especially in arid and semi-arid regions. In such regions, groundwater resources are over-exploited to meet the development and urbanization of coastal zones. The excessive pumping of groundwater in coastal aquifers has reduced the freshwater flux to the sea and allowed the inland migration of the seawater. Several investigations have reported the abandonment of well fields due to seawater intrusion in coastal zones (Abd-Elhamid 2010).

A number of studies have been conducted to simulate the seawater intrusion in the Nile Delta aquifer using different numerical techniques. Examples of these studies include, among others, El-Fayoumi (1968), Wilson et al. (1979), Kashef (1983), Farid (1985), Sherif et al. (1988, 1990a, 1990b, 2012), Darwish (1994), Amer & Sherif (1995, 1996), Amer & Farid (1981) and Sherif & Singh (1997). Few studies addressed the possible impacts of climatic change and sea level rise on the seawater intrusion in the Nile Delta aquifer (e.g. Sherif & Singh 1999; Sefelnasr & Sherif 2014).

Mohamed (2004) studied the global and local factors affecting the rise of sea level and its impact on coastal regions and specifically deltaic zones. He assessed the impact of such factors on the Northern coasts of Egypt extending from Alexandria to El-Arish for the years 2020 and 2050. He used an artificial neural network to estimate the change in sea water level from available tide gauge and weather station data in the vicinity of the Northern Egyptian coasts. Nasar (2004) conducted a numerical and experimental study to assess the environmental impacts of extracting brackish water for desalination and recharging the brine water through deep wells to control the seawater intrusion into the Nile Delta aquifer. A number of design charts for extraction and recharge wells, including discharge and recharge rates and locations and depths of pumping and injection wells, were presented. Saleh (2009) developed a numerical model using MODFLOW and MT3D to quantitatively and qualitatively assess the possible impacts of the policies related to groundwater extraction and development activities in the Nile Delta region.

Recently, Pandian et al. (2016) presented finite element modelling of a heavily exploited coastal aquifer using FEFLOW software to assess the response of groundwater level to the changes in pumping and rainfall variation due to climate change. The study predictions indicated that 10% increase of recharge and 10% decrease of pumping causes a 3 m and 6 m increase in groundwater head in the studied aquifers, respectively. Surinaidu (2016) used different hydro-chemical mixing models and multivariate statistical techniques, including factor and cluster analysis to identify the salinity sources in the central Godavari delta. The results of the study revealed that very high salinity in pumping wells is due to up-coning of saltwater. Shapouri et al. (2016) examined the variation of stygofauna composition collected in wells, along a gradient in groundwater salinity/conductivity in a coastal aquifer from southern Portugal. The studied coastal aquifer is considered vulnerable to salinization due to seawater intrusion, caused by overexploitation of the aquifer. The study suggested that stygofauna composition and structure can be a useful complementing tool for monitoring seawater intrusion into coastal aquifers, where reduction or deterioration in groundwater resources is predicted.

In this paper, a coupled finite element model for simulating seawater intrusion in coastal aquifers is employed. The model accounts for transient density-dependent flow, heat transfer and solute transport in saturated and unsaturated soils. The model is developed and applied to simulate seawater intrusion in the Nile Delta aquifer along three different vertical cross sections. The results of the developed model were compared with those of Sherif et al. (1999a). Also, the results of the current model were compared to results of SEAWAT code for model verification. The model was also applied to investigate the position and movement of the transition zone in response to sea level rise and extensive groundwater pumping. Such pumping might be necessary in the face of climate change with diminished rainfall and surface water resources. Different scenarios of sea level rise and groundwater pumping have been considered. The selected scenarios accounted for (a) sea level rise, (b) decline of the piezometric head on the land side due to additional groundwater pumping, and (c) combined effect of sea level rise and additional groundwater pumping.

The governing differential equations describing transient fluid flow, heat transfer and solute transport in saturated/unsaturated soils are presented in terms of four primary variables; pore-water pressure ⁠, pore air pressure ⁠, absolute temperature and solute concentration (c). The governing differential equations for water flow, air flow, heat transfer and solute transport can be expressed as:

Detailed information related to governing differential equations for water flow, air flow, heat transfer and solute transport and the numerical solution of coupled fluid flow and solute transport can be found in Abd-Elhamid & Javadi (2011). The validation of the model can be found in Abd-Elhamid (2010).

Arid and semi-arid regions are more vulnerable to water shortage problems due to the scarcity of rainfall and limited natural recharge events. The tremendous increase in the population of the coastal regions and the associated increase in water consumption by different sectors have led to over-exploitation of the groundwater resources in coastal aquifers over the last few decades. This practice has accelerated seawater intrusion and caused significant deterioration of the groundwater quality in the Nile Delta aquifer. A comprehensive study is, therefore, needed to assess the current situation of seawater intrusion and predict any possible additional intrusion considering different scenarios of pumping and sea level rise. In the current study, the 2D-FEST is employed to simulate the current status and assess the future conditions of groundwater in the Nile Delta aquifer under different scenarios. For validation purposes, the results are compared with those of Sherif et al. (1999a). The developed model has then been used to identify the seawater intrusion and the dispersion zone in other vertical cross sections that have not been investigated before.

Egypt lies between latitudes 22° and 32° North and longitudes 25° and 35° East. It is bounded by the Mediterranean Sea to the North and the Red sea to the East. The Southern and Western boundaries are political boundaries with Sudan and Libya, respectively. The total area of Egypt is about one million km², of which 94% is desert. The population of the country is currently estimated as 85 million. There are three interdependent sources of water in Egypt: Nile water, groundwater, and drainage water. The average annual rainfall is around 150 mm (Abdelaty et al. 2014a, 2014b). Due to scarcity and randomness, rainfall does not contribute to the water resources budget of the country. The water of the Nile River represents the only renewable water resource in the country. The Egyptian share from the Nile water, according to the international agreements, is limited to 55.5 billion m³ per year since 1959. The per capita share from the renewable water is estimated at 650 m³/year. Currently, about 5.5 billion m³/year of groundwater is pumped from the different aquifers, 85% of which originates from the Nile water and the rest is fossil water. About 5 billion m³/year of agriculture drainage water is being reused, while about 12 billion m³ of low quality drainage water is expelled to the Mediterranean Sea and the northern lakes (Abdelaty et al. 2014a, 2014b). Expansion in the reuse of drainage water is not recommended because of its low quality and high content of pesticides and other chemicals.

The Nile Delta represents the most fertile land in Egypt. However, less than two-thirds of it is under cultivation. The uncultivated lands are found along the sea to the north and in the vicinity of Delta fringes close to the desert. Irrigation mostly exploits surface water through an extensive network of canals branching out from the Nile River. The growth of population in the Nile Delta region along with the increase in agricultural and industrial activities has imposed an increasing demand for freshwater. This increase in demand in the Delta is met by intensive pumping of fresh groundwater, causing subsequent lowering of the piezometric head and upsetting the dynamic balance between freshwater and saline water bodies in the aquifer. Like any coastal aquifer, an extensive saltwater body has intruded the Nile Delta aquifer forming the major constraint against aquifer exploitation. The Nile Delta region lies within the temperature zone which is a part of the great Desert belt. The average rainfall along the shore of the Mediterranean Sea is about 180 mm/year. This amount decreases very rapidly as one proceeds inland to about 26 mm at Cairo.

The groundwater in the Nile Delta aquifer fills a vast underground bowl situated between Cairo and the Mediterranean Sea. If it were not for the presence of a saline body of seawater along the bottom of this bowl, the Nile Delta aquifer could be easily exploited to the best interests of Egypt. Seepage from the Nile River, canals and the surplus of irrigation water percolate downward, representing the main source of aquifer recharge. It may also be recharged, nominally, by northward flow from the Nile Valley aquifer. On the other hand, the aquifer loses some of its fresh water to the Mediterranean Sea and the drainage system in the northern part of the Delta.

The Nile Delta aquifer system forms an immense and complex groundwater system. It is a leaky Pleistocene aquifer system overlain by a semi pervious Holocene aquitard (clay cap) and underlain by an impermeable Miocene aquiclude. The clay cap (aquitard), being the main source through which aquifer recharge occurs, is of great importance. Based on the thickness and vertical permeability of the upper confining clay layer, the aquifer system is identified as confined, semi-confined or phreatic. On the contrary, the aquiclude layer has no hydrologic importance, except that it forms an impermeable bottom and defines the aquifer geometry. The bulk of the Nile Delta aquifer consists of deltaic deposits 200–300 m thick on average. It is dominated by unconsolidated coarse sands and gravels with occasional clay lenses. The top boundary of the deltaic deposits, which acts as a cap for the aquifer, is a semi-pervious clay and silt layer. The thickness of the upper clay cap increases uniformly toward the Mediterranean Sea from about 5 m near Cairo and reaching about 70 m near the shore line. The clay cap is inter-fingered with the aquifer near the shore.

The Nile Delta with its fringing area covers an area of about 22,000 km². It lies between latitudes 30° 25′ and 31° 30′ North, and longitudes 29° 50′ and 30° 15′ East. The Nile Delta has its apex about 20 km north of Cairo. The length of its base at the Mediterranean Sea is about 245 km. The length of its right branch (Damietta) is about 240 km and that of the left branch (Rosetta) is about 235 km (see Figure 1). Water flowing through the Damietta and Rosetta branches accounts for less than half of the original total northward flow across the Delta plain, that is now captured by artificial waterways comprising over 10,000 km of canals, usually only 2–3 km apart. This huge surface water network interacts with the groundwater system of the Nile Delta.

Horizontal and vertical dimensions of the Nile Delta aquifer were obtained using deep oil borings and data from test bore holes (Farid 1980). The boundary between the Pleistocene aquifer and the Pliocene and/or Miocene aquiclude was defined. The thickness of the Pleistocene aquifer increases toward the Mediterranean Sea, and decreases southward, almost vanishing near El Manawat and partly separating the Nile Delta aquifer from the Nile Valley aquifer. South of Tanta in the transverse direction, the bottom boundary forms a concave shape. The aquifer thickness decreases eastward and southward with a maximum depth in the middle of the Delta. North of Tanta, the maximum thickness of the aquifer shifts toward the east.

The clay cap takes different profiles along the shore of the Mediterranean Sea. Generally, it is divided into two layers. The first layer is the upper clay cap which acts as an aquitard, about 20 m thick, with a low permeability. The second layer is the lower clayey sand layer, a few meters thick, with a higher permeability. The thickness of the upper clay cap (aquitard) and the clayey sand layer are well defined at many locations in the Delta area. Except in the vicinity of the shore line, their thicknesses are mostly less than 20 m and 15 m, respectively. No clay cap exits at the Delta fringes.

The vertical water movement affects, to a great extent, the water balance of the system. Downward movement of water through the clay cap occurs in two different stages. The first is the downward infiltration of irrigation water from the ground surface to the subsoil water table through the unsaturated zone of the clay cap. The velocity of this movement is defined by the downward infiltration velocity. The second is the movement of the subsoil water from the water table to the groundwater in the aquifer through the saturated zone of the clay cap. The velocity of the flow in this zone is governed by Darcy's law and is defined by seepage velocity. Vertical hydraulic conductivity, K_cv, of the clay cap has an average value of 2.5 mm/day (Abdelaty et al. 2014a, 2014b).

Field experiments were carried out to determine the aquifer parameters. A hydraulic conductivity of 100.0 m/d and a storativity of about 10^–4 to 10^–3 were found to be representative of the regional values for the Nile Delta aquifer (Abdelaty et al. 2014a, 2014b). Farid (1980) reported different values for hydraulic conductivity and storativity at various locations. The hydraulic conductivity of the aquifer decreases toward the south and west. An effective porosity of 0.3 is considered representative for the aquifer medium. The dispersivity of the Nile Delta aquifer has not been measured. The dispersivity accounts for the lack of information about the pore velocity fluctuation when passing from the microscopic to the macroscopic configuration of the solid liquid interface. Values between 0.1 and 500 m can be found in the literature. Based on Sherif et al. (1988), the longitudinal and lateral dispersivities for the Nile Delta aquifer are set to 100 m and 10 m, respectively. The free water table (if not measured) is generally assumed to be 1.0 m below the ground surface of the Delta. The piezometric head is monitored periodically by the Groundwater Research Institute (GWRI), National Water Research Center, Egypt (Abdelaty et al. 2014a, 2014b).

On the land side (left hand side) boundary, the concentration is constant and equal to the freshwater concentration, C_f. The freshwater flux through the landside boundary is defined and a hydrostatic pressure distribution is assumed in the vertical direction. The bottom of the aquifer is impermeable. For the upper leaky boundary, where the aquifer is recharged by freshwater from irrigation water and leakage from the Nile River and irrigation network, the concentration is equal to that of freshwater, C_f, and the concentration gradient is set equal to zero. At the seaward boundary, where the aquifer is exposed to the Mediterranean Sea, the concentration is equal to the seawater concentration, C_s, where the flux across the boundary is directed inward, and as zero dispersive flux boundary where the flow direction is outward toward the sea. The pressure is hydrostatic at this boundary with an atmospheric pressure at sea level.

The 2D-FEST model was employed to identify the dispersion zone and the flow pattern in the Nile Delta aquifer using the variable density approach. All the simulations were conducted in the vertical view. Three vertical sections in the eastern, middle and western parts of the Nile Delta aquifer were considered, as shown in Figure 2. The upper boundaries of the three sections were assumed horizontal while the lower boundaries were inclined with an increased depth toward the seaside. The geometric characteristics of the three sections are as follows:

• Section I: length is 120 km, depths are 200 m at the land side and 900 m at the sea side.

• Section II: length is 150 km, depths are 240 m at the land side and 680 m at the sea side.

• Section III: length is 70 km, depths are 200 m at the land side and 600 m at the sea side.

In scenario 1 the water level in the sea is raised by 100 cm, the free water table is kept unchanged and all other parameters are kept to their original values. In this scenario the equiconcentration line 35.0 advanced by a distance of 9 km, measured along the bottom boundary and the equiconcentration line 1.0 advanced inland by a distance of 10 km, as shown in Figure 9(a) compared against the current situation presented in Figure 3.

In scenario 2 climate change and population growth may create a shortage of surface water, resulting in additional demands for groundwater resources. To investigate the effect of additional pumping from the Nile Delta aquifer the piezometric head on the land side was lowered by 100 cm (as a result of additional pumping) and other parameters were kept unchanged. Under this scenario, the equiconcentration line 35.0 advanced inland a distance of 6 km (compared with current conditions) measured along the bottom boundary. The equiconcentration line 1.0 advanced inland by a distance of 8 km. The width of the dispersion zone increased considerably, as shown in Figure 9(b).

In scenario 3 the water level in the sea is raised by 100 cm, the piezometric head on the land side is lowered by 100 cm and all other parameters are kept to their original values. The equiconcentration line 35.0 advanced inland a distance of 14 km (compared with current conditions) measured along the bottom boundary. The equiconcentration line 1.0 advanced inland by a distance of 15 km, as shown in Figure 9(c).

The results of the simulation model indicate that the change of water level in the sea side has a significant effect on the position of the transition zone, especially if the effect of sea level rise is combined with the effect of increased abstraction from the aquifer, as is the case in scenario 3. In this case the transition zone is shifted further inland.

The Nile Delta aquifer is severely susceptible to seawater intrusion from the Mediterranean Sea causing serious environmental impacts. Groundwater resources should thus be developed and managed carefully to avoid any further deterioration in the water quality. Seawater intrusion in the Nile Delta aquifer was simulated using the 2D-FEST model and the results were compared with the results of SEAWAT code. A good agreement was observed between the results of the two models. Different scenarios of climate change and sea level rise and their impacts on seawater intrusion were simulated using the 2D-FEST model. The simulation of the seawater intrusion problem in the Delta revealed that the equiconcentration line 35 (representing seawater) has intruded inland to a distance varying from 42 to 64 km measured along the bottom boundary from the seaside. However, the equiconcentration line 1 has intruded inland to a distance varying from 68 to 112 km measured along the bottom boundary from the seaside. The results show that the rise in sea level has a significant effect on the position of the transition zone. On the other hand, the combination of sea level rise and over-pumping results in movement of the transition zone further inland. The rise in sea water levels will impose additional saline water heads at the sea side and therefore more seawater intrusion is anticipated. A 100 cm rise in the Mediterranean Sea water level will cause additional intrusion of about 10 km in the Nile Delta aquifer. Reduction in piezometric head water at the land side due to over-pumping will cause an additional intrusion of 8 km. The combination of scenarios 1 and 2 will cause an additional intrusion of 15 km.