Assessment of Changing the Abstraction and Recharge Rates on the Land Subsidence in the Nile Delta, Egypt

Abstract

The majority of residential, agricultural, and industrial areas are situated on cohesive soil in the Nile Delta, Egypt. Excessive pumping from the Nile Delta aquifer to meet the increasing demands for water could lead to aquifer system compaction and land subsidence. Land subsidence endangers infrastructure such as buildings, bridges, canals, and roads, as well as deteriorating lands and agricultural resources. The objective of this research is to investigate the land subsidence and predict the future behavior of the middle Nile Delta. The study goals are met by using a numerical model (MODFLOW) to simulate groundwater flow and an analytical solution to calculate land subsidence conditions. In this study, three scenarios are considered including; decreasing aquifer recharge, increasing abstraction and combination of the two. The results reveal that decreasing recharge by 94.4%, 88.8%, and 83.2% led to 30-, 60-, and 90-mm land subsidence, respectively, while increasing abstraction by 146%, 193%, and 233% led to land subsidence by 190, 380, and 560 mm, respectively, in the Nile delta. However, the combination of the two scenarios led to 220-, 440-, and 650-mm land subsidence. According to the results the future land subsidence due to over pumping from the Nile Delta should be considered in the future development plans of the country which intend to increase the abstraction from the Nile Delta aquifer. Increasing abstraction could increase the land subsidence that may cause many damages in different properties.

1. Introduction

Groundwater (GW) is one of the most significant natural water resources and a source of survival for many living organisms. In Egypt, GW is one of the most valuable water resources after the Nile River. Most of Egypt’s residential areas are located on cohesive soil (clay) in Nile valley and Delta. The Nile Delta (ND) in Egypt, like all other deltas all over the world, is a densely inhabited area [1]. The Delta is home to 50 million people and has a population density of 2300 people per square kilometer, which is expected to double in the next 50 years [2]. The Nile Delta is the most significant strategic and economic region, with precious resources (such as oil, gas, black sand, and groundwater), fertile soil, and historic ancient sites [3]. It also includes several densely populated cities in Egypt as well as the majority of the industrial and commercial cities (Alexandria, Port Said, Damietta, Rosetta, and Suez) [4]. Due to the population explosion, economic development, and the impact of climate change, a poor future is predicted for the Nile Delta [5,6]. The Nile Delta faces serious environmental challenges such as land subsidence (LS), sea level rise (SLR), and coastline erosion.

One of the most serious threats to the ND is LS. LS is a vital geological hazard that causes the land surface elevation to gradually lower [7]. It causes serious problems such as damages to buildings and infrastructures and deterioration of land and agricultural resources. Additionally, LS increases the risk of flooding and seawater incursion in coastal areas [7]. LS in the ND is influenced by a lot of natural and anthropogenic physical processes. Holocene fluvial sediment LS which is linked to global SLR and climate change, as well as localized tectonic, downslope mass movements [8]. The Nile River discharge cycle is another physical phenomenon that influences LS rates in the ND. The recharge/discharge cycle carries sediments from the Nile’s various climate regimes and causes land-loading variations as water tables rise and fall [3]. Building dams and regulators holds sediments upstream and reduces sediment loads reaching deltas and coastal area. For example, the establishment of the Aswan High Dam in southern Egypt decreased the rate of sedimentation since the 1970s [9]. Local tectonic movements can also induce to LS in the ND in the long run [10]. GW and oil and gas withdrawal, gives enough space for sediments to re-arrange, resulting in soil compaction [11].

One of the most common causes of LS in the world is GW extraction due to increasing the demand for water resources for domestic, agricultural, and industrial purposes. For example, in the United States (USA), GW over pumping is responsible for more than 80% of current subsidence [12]. Land subsidence of up to 10 m due to GW extraction is also observed in Mexico City, Mexico [13,14]. In Taiwan, the subsidence area of 2400 km² was subjected to 3.4 m LS due to massive GW use [15]. LS due to GW extraction in sedimentary environments leads to the compaction of soil layers [16]. The volume of groundwater withdrawn and the layers from which groundwater was obtained are the factors that influence the magnitude and location of compaction. The principle of effective stress can be used to describe the mechanism of soil compaction (Terzaghi, 1925) [17]. The effective stress, which is linked to the pore water pressure, supports the overall loading that is attributed to solid material and fluid. Groundwater withdrawal generates a decrease in pore water pressure, which leads to an increase in strain and effective stress that compresses the soil skeleton and compacts the soil layers when the total stress is given [16]. Therefore, controlling LS required setting a limitation of GW withdrawn volume.

A number of researchers have used geodetic techniques, numerical models, and analytical solutions to investigate the rate and patterns of LS. Investigation of LS using geodetic techniques such as Global Positioning System (GPS) and Interferometric Synthetic Aperture Radar (InSAR). The geodetic techniques have been extensively used to assess LS in various regions all over the world, such as the Mekong Delta (Vietnam) [18]; the Nile Delta (Egypt) [1,3,19,20,21,22]; the Mexicali Valley (Mexico) [23]; Miami Beach and Norfolk (USA) [24] and Las Vegas Valley (USA) [25]. The geodetic techniques’ disadvantages are: the limited level of confidence in the physical causes of subsidence, as well as the complicated processing of SAR imaging in vegetation-covered regions and wetlands due to low coherence pixels. Furthermore, a scarcity of accurate, independent geodetic data to corroborate or InSAR result to estimate the groundwater level [3]. Investigation of LS based on numerical models is useful for simulating compaction for complicated aquifer parameters, boundary geometry and pumping activities. Many studies used numerical models to assess LS in different regions all over the world such as, the Los Banos-Kettleman City area in California [26], the Emilia-Romagna coastland of Ital [27], the Suzhou area of China [28] and the Shanghai area of China [29]. The limitation of numerical model is the difficulty to measure LS for multi-layers aquifer and near pumping wells [16]. Yilong, Yuan, et al. (2020) has developed a THM model to investigate land subsidence due to gas production in south China [30]. The results indicated that the maximum land subsidence occurred around the gas production wells. Additionally, Kiryukhin, A.V., et al. (2015) investigated land subsidence at Mutnovsk geothermal field in Russia using THM model [31]. The model results showed high rates of subsidence in 2008 reaching 6–18 mm/year. Investigation of LS based on analytical models is convenient for initial investigation, sensitivity and uncertainty analysis, and the numerical model’s verification [16]. Analytical solution has been applied in many areas such as, Shen, S. and Xu, Y (2011) [32] applied analytical solution to predict the future behavior of LS in Shanghai due to groundwater over-pumping. Po-Lung Lin et al. [16] developed a one-dimensional deformation analytical model to simulate groundwater flow and land subsidence in Tawian. Abd-Elhamid et al. [33] established an analytical model to simulate the soil compaction in multi-layer aquifer considering different scenarios of abstraction rates and pumping wells location.

Studies have investigated LS in the Nile Delta of Egypt. For example, between 1993 and 2000, Aly et al. [19] used InSAR images to investigate LS in Greater Cairo. Between 1993 and 2000, the mean rate of LS was around −7 mm/y. Aly et al. (2012) [20] used PSI images to estimate the magnitude of LS from 1993 to 2000 in two main cities in the Delta, Mansura and Greater Mahala, the rates were −9 and −5 mm/y, respectively. Wöppelmann et al. [21] investigated LS in Alexandria between 2003 to 2010 based on PSI images, the rate ranged between −0.4 to −2 mm/y Gebremichael et al. [1] used PSI images between 2004 and 2010 over the whole Delta to determine LS. The results showed a stable and minor land lowering in the southern part of the Delta of −0.4 mm/y, a local uplift in the central zone of the Delta’s plains with a rate of 2.5 mm/y and a high subsidence rate of −9.7 mm/y in the Menoufia governorate, and reclaimed desert due to GW withdrawal. Saleh and Becker [22] investigated LS across the Nile Delta, using InSAR data and the PSI images. The rates from Cairo, Tanta, Mahala, Mansoura, Damietta, and Port Said were around −6.4 ± 0.4 mm/y, −4.0 ± 0.6 mm/y, −4.8 ± 1.0 mm/y, −10.0 ± 1.2 mm/y, −10.3 ± 1.6 mm/year, and −4.9 ± 1.6 mm/y, respectively. Rateb and Abotalib [3] investigated three patterns of LS due to urban loading, natural sediment compaction, and GW over-pumping. The rates of LS due to GW withdrawal ranged from −20 to −16 mm/year and −6 to −12 mm/year in newly reclaimed lands on the west and east of the Delta’s flood plains.

In this study, LS of the ND is assessed with a new technique that is significantly different from the previous studies based on analytical solution. In this study numerical and analytical models are developed to investigate LS in the Middle Nile Delta aquifer (MNDA). MODFLOW is used to predict the groundwater levels at different scenarios of GW including over-pumping, decreasing the recharge, and combination of both scenarios. Moreover, the analytical and numerical solutions are used to predict future possible LS in MNDA according to different scenarios of GW pumping and aquifer recharge. Additionally, SURFER 13.0 is used to map both GW drawdown and LS in the MNDA for different scenarios.

2. Materials and Methods

2.1. Study Area and Used Data

The Nile Delta (ND) contains the majority of fertile lands in Egypt, as around two-thirds of its area is under agriculture. Irrigation water is derived mainly from surface water via canals branching out of the Nile River as well as groundwater from the Nile Delta aquifer (NDA). The NDA is one of the world’s deepest freshwater aquifers. The Middle Nile Delta aquifer (MNDA) is considered as the study area in this paper since it is anticipated to be the most affected part of the ND by land subsidence (LS). The MNDA is located between 30°20′ and 31°50′ N latitude and 30°10′ and 31°35′ E longitude covering an area of 9000 km². It is bordered from the north by the Mediterranean Sea, from the south by Cairo, and from the east and west by the two branches of the Nile River. The apex of the MNDA is around 20 km north of Cairo. Its Mediterranean Sea base is around 245 km long, the length of the right branch (Damietta) is about 240 km and the length of the left branch (Rosetta) is around 235 km (see Figure 1a).

2.1.1. Geology of the Study Area

The NDA system is a massive and complex groundwater system. The geological units of the aquifer system are the top unit of Holocene clay (aquitard), quaternary and late tertiary gravels and sands (aquifer), and base unit (see Figure 1b). The clay cap thickness in the southern part of the Nile Delta ranges from 5 to 25 m, however it reaches 50 m in the northern part [34]. The Quaternary aquifer’s thickness ranges from 200 m near Cairo in the south to around 1000 m towards the Mediterranean Sea in the north [35]. The major aquifer of the study area is the Quaternary aquifer [36]. The top layer of the quaternary aquifer is a clay layer, so it is considered a semi-confined aquifer. The level of groundwater in the NDA increases from north to south. It varies from 1 to 2 m above MSL in the north, 3 to 4 m above MSL in the middle, and 5 m above MSL in the south [36,37].

Figure 1. The Nile Delta aquifer, (a) the location map [38] (b) hydrogeological cross section (A-A) in the NDA [39].

Figure 1. The Nile Delta aquifer, (a) the location map [38] (b) hydrogeological cross section (A-A) in the NDA [39].

2.1.2. Meteorological Data of the Study Area

The average daily temperature in the Nile Delta varies from 17 to 20 °C along the Mediterranean to more than 25 °C in the south [40]. In the study area, the evaporation rates range from 7 mm/day in the south to 4 mm/day along the northern Mediterranean Sea coast [41]. The average rainfall in the Nile Delta is fairly small and ranges from 25 mm/y in the south and middle to 200 mm/y in the north [34].

2.1.3. Population Growth and Historical Situation of Groundwater Levels

Egypt is one of the most densely populated countries in the Middle East. The Nile Delta of Egypt is overpopulated area as it represents for only 2% of Egypt’s total land area, but it is home to 41% of the country’s people [42]. The Egyptian population was 82 million in 2010 and increased to 102 million in 2020 [43] as shown in Figure 2b. The population growth in the Nile Delta is followed by the increase in agriculture and industrial activities. Egypt’s population is predicted to double by 2078, according to current projections. The population is currently rising at a rate of 1.94%, which means that every year, about 2 million residents are added to the population [43]. This increasing population growth rate causes an increase in the water demand which resulted in extensive extraction from groundwater. Around 85% of Egypt’s total GW abstraction is provided by the NDA. Over the last 30 years, total yearly GW extraction from the NDA has rapidly increased. According to the Egyptian Research Institute of Groundwater (RIGW), the abstraction rate was 0.1 billion cubic meters (BCM)/year from the year 1980 to 2010 and the abstraction rate increased to 0.2 BCM/year from 2003 to 2010. In the next years, it is predicted that the yearly pumping rate will approach 0.20 BCM/year [44] as shown in Figure 2a.

The NDA is recharged through the Nile River flow, percolation from excess irrigation, rainfall and seepage from canals. The significant source of GW recharge in the Nile Delta is the Nile River flow as the NDA receives around 35.5 km³/y from the Nile flow [45]. Excess irrigation water recharges the NDA at a rate ranging from 0.25 to 1.1 mm/d. The NDA recharge rate from rainfall is very small—it varied from 0 to 134 mm in winter season. The groundwater situation between 1980 and 2010 in the Nile Delta was complicated as it was dependent on seepage losses from canals and GW over-pumping due to reclamation projects [36]. In both the Eastern Nile Delta and Western Nile Delta, depression cones in GW tables were created due to GW over-pumping. In the Middle Nile Delta, there was no change in GW levels, which means a steady-state situation [36]. The GW recharge is highly affected by the Nile flow change. The Nile flow will change by 0.28%/year due to climate changes and the estiblashment of the Grand Ethiopian Renaissance Dam (GERD) [46]. The construction of GERD will reduce the Nile flow to Aswan High Dam (AHD). The Nile flow to AHD was estimated to decrease by 10 to 50% in 2020 [47].

2.2. Analytical Solution of Land Subsidence

The analytical solution for simulating LS based on the consolidation due to change in effective stress (σ_e). The developed solution is based on Terzaghi 1D consolidation theory [48]. This theory based on the relationship between the effective stress (σ_e), the pore water pressure (p) and the total stress (σ_t). During the groundwater withdrawal, the hydraulic head (h) and poor water pressure (P) decrease, thus the effective stress increases (σ_e) by the same value and the total stress (σ_t) remain constant with time as following:

Δσ_e = −ΔP= −γ_w Δh =γ_w Δs

(1)

The relationship between effective stress (σ_e) and strain (ε) was revealed by Das (2006) [36] as following:

Δε = α Δσ_e

(2)

where α is the compressibility coefficient of a soil.

By combining Equations (1) and (2), the following equation is obtained

Δε = α γ_w Δs

(3)

As shown in Equations (1) and (3), groundwater withdrawal causes increase in effective stress. Therefore, consolidation occurs, then land subsidence happens when the consolidation transmits to the land surface.

The compressibility coefficient (α) can be derived from the following relationship

S_s = γ_w (α + n β)

(4)

The water compressibility coefficient (β) can be ignored as it is very small compared to the soil compressibility coefficient (α). Therefore, Ss = γ_w α.

The compaction of single soil layer (δ_i) with thickness (b_i) can be calculated based on the change in groundwater drawdown (Δs) from specific storage change (Ss) as shown:

δ_i (t) = α γ_w Δs_i (t) b_i = S_S/γ_w ∗ γ_w Δs_i (t) b_i = S_s Δs_i (t) b_i

(5)

Land subsidence for multi-layer aquifer system (δ) can be calculated by accumulating the compaction of single soil layers as following:

$$\mathsf{\delta }\left(\mathrm{t}\right)={{\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{n}}\mathsf{\delta }}}_{\mathrm{i}}\left(\mathrm{t}\right)={\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{S}}_{\mathrm{si}}}{ \Delta \mathrm{S}}_{\mathrm{i}}{\left(\mathrm{t}\right)\mathrm{b}}_{\mathrm{i}}$$

(6)

2.3. Numerical Simulation of Groundwater Flow

In this study, a numerical model (MODFLOW) is built to simulate 3D, variable density, and steady-state GW flow in MNDA porous media. A mathematical combination of the water balance equation and Darcy’s law derives the governing equation of the GW flow [49]. According to the following equation [50], the MODFLOW model represents GW flow in anisotropic and nonhomogeneous media:

$$\frac{\partial }{\partial \mathrm{t}}({\mathrm{K}}_{\mathrm{x}\mathrm{x}}\frac{\partial \mathrm{h}}{\partial \mathrm{x}})+\frac{\partial }{\partial \mathrm{y}}({\mathrm{K}}_{\mathrm{y}\mathrm{y}}\frac{\partial \mathrm{h}}{\partial \mathrm{y}})+\frac{\partial }{\partial \mathrm{z}}({\mathrm{K}}_{\mathrm{z}\mathrm{z}}\frac{\partial \mathrm{h}}{\partial \mathrm{z}})+\mathrm{W}={\mathrm{S}}_{\mathrm{s}}\frac{\partial \mathrm{h}}{\partial \mathrm{t}}$$

(7)

where K_xx, K_yy, and K_zz are values of porous media hydraulic conductivity along the x, y and z coordinate axes (LT⁻¹), h is the potentiometric head (L), W is volumetric flux per unit volume rep-resenting source and /or sink water in (T⁻¹), S_s is specific storage (L⁻¹), and t is time (T).

2.3.1. Model Boundary Geometry

The 3-D model of MNDA was established using a grid system comprising 160 columns and 172 rows with cell dimensions of 1 km × 1 km as shown in Figure 3a. The simulated area (9000 km²) was divided into 11 layers, including the clay cap layer and other layers represented quaternary aquifer with equal thickness. Figure 3b,c shows two cross sections taken in the model in x and y directions. As indicated in Figure 3, the layers numbers 1, 2–5, 6–9, and 10–11 are in colors blue, green, red and purple, respectively. According to the cross sections, the aquifer thickness ranges from 200 m at the south (near Cairo) to 1000 m at the north (along the Mediterranean Sea shoreline) with a clay layer between 20 at the south to 50 m at the north.

2.3.2. Model Boundary Conditions

The flow boundary conditions were assigned as a constant hydraulic head of zero along the shoreline of the Mediterranean Sea in the north and a specified head of 16.96 m above MSL in the south. The east and west boundaries were determined by the water levels in the Nile’s Damietta and Rosita branches. The water level in the Damietta branch declined and ranged from 13.66 m above MSL in the south to 0.5 m above MSL in the north. Moreover, the water level in the Rosita branch varied from 13.17 m above MSL in the south to 0.5 m above MSL in the north. Canals and drains packages were used to assign canals and drains in the model. The study area’s main and branch canals had average depths of 2.0 and 3.0 m, respectively. Drain package head ranged from 8 m above MSL in the south to 0.25 m above MSL in the north. Figure 4a shows the MNDA boundary conditions.

2.3.3. Model Hydraulic Parameters

Numerous field experiments were conducted to estimate the aquifer parameters. Table 3 shows the initial values of hydraulic parameters in the MNDA area. Specific storage (S_s), specific yield (S_y), hydraulic conductivity (K), and effective porosity (n_eff) are all the hydraulic parameters for each layer of the study area. These data were collected from previous studies. The hydraulic conductivity was ranged from 0.1 to 0.25 m/d for clay cap and from 5 to 100 m/d for the other layers, and the storativity ranges from 10⁻⁴ to 10⁻³ [51].The model recharge from excess irrigation water and the GW extraction rate were based on the identified value by Morsy 2009 [36]. The recharge rate to the aquifer is 0.25 to 0.8 mm/day. The abstraction wells distribution in the MND is assigned according to data collected by RIGW [35,36]. Figure 4b presents the distribution of 946 abstraction wells assigned to the study area. Production irrigation water wells and drinking water wells are among the groundwater abstraction activities in the Nile Delta aquifer. The abstraction from wells produced a total volume rate of 2.22 × 10⁶ m³/day. As shown in Figure 4b, there are few abstraction wells in the northern part of the MND to avoid saltwater intrusion due to over-pumping because the northern part is near to the Mediterranean Sea.

2.3.4. Model Calibration

The model calibration is the process that shows the difference between the calculated head and the field data obtained by RIGW 2002 [37], which have been presented by 24 observation wells throughout the study area. Figure 5 shows the distribution of 24 observation wells in the study domain. The model was calibrated for steady-state using a trial-and-error process for the aquifer parameters presented in

Table 1. Proposed scenarios of different GW abstraction and surface recharge rates.

Stage	Scenarios No.	Year	Abstraction		Recharge
			Values (BCM.yr−1)	Rates (%)	Values (BCM.yr−1)	Rates (%)
Base Case	1	2010	0.81	100	1.4	100
Decreasing Recharge	2	2030	0.81	100	1.32	94.4
	3	2050			1.24	88.8
	4	2070			1.16	83.2
Increasing Abstraction	5	2030	1.18	146	1.4	100
	6	2050	1.56	193
	7	2070	1.89	233
Combination	8	2030	1.18	146	1.32	94.4
	9	2050	1.56	193	1.24	88.8
	10	2070	1.89	233	1.16	83.2

until the GW levels in the observation wells matching the levels in 2002. After model calibration, a model verification and preliminary sensitivity analysis were performed to investigate the most significant parameters affecting the model results.

2.4. Proposed Scenarios for Simulating the Land Subsidence in the MNDA

In the MNDA, the MODFLOW calibrated results for groundwater head are used as a base case for future scenarios as shown in Table 1. Due to population growth, abstraction rates will increase to 146%, 193%, and 233% in 2030, 2050, and 2070, respectively as shown in Figure 2a. Due to climate change and the construction of GERD, the Nile flow will change by 0.28% every year, reducing recharge to 94.4%, 88.8%, and 83.2% in 2030, 2050, and 2070, respectively. The predicted heads are used to evaluate the land subsidence using the analytical solution. Free version of SURFER 13.0 model was used for mapping land subsidence in the study area.

3. Results

3.1. Model Calibration Results

The model was calibrated with 24 values of GW head for steady state condition. The calibration results indicated that the maximum and the minimum residual (the difference between the calculated and the observed heads) reached −0.893 and 0.231 m, respectively, while the absolute residual mean was 0.471, the root mean square was 0.487 m and normalized RMS was 3.249%. The model’s calibration aim was 10% difference in height between the minimum and maximum GW heads, which is around 1.5 m as presented in Figure 6a.

Table 2. Result of model calibration for calculated and observed heads.

Well	Calculated Head (m)	Observed Head (m)	Difference(m)
1/A	1.42	1	0.42
2/A	1.89	2.3	−0.41
3/A	1.56	2	−0.44
4/A	2.48	3	−0.52
5/A	3.97	3.7	0.27
6/A	4.8	4.2	0.6
7/A	4.18	4.6	−0.42
8/A	3.97	4.4	−0.43
9/A	5.96	5.4	0.56
10/A	7.13	6.6	0.53
11/A	8.53	8.3	0.23
12/A	8.69	9.1	−0.41
13/A	9.96	10.4	−0.44
14/A	10.67	10.1	0.57
15/A	10.12	9.6	0.52
16/A	11.22	10.9	0.32
17/A	11.01	11.6	−0.59
18/A	11.58	12.1	−0.52
19/A	11.94	11.5	0.44
20/A	14.16	15	−0.84
21/A	14.6	14.1	0.5
22/A	14.98	15.4	−0.42
23/A	15.53	16	−0.47
24/A	15.64	15	0.64

presents the agreement between the model results and the field measurements confirms the model calibration.

The model verification comes after the model calibration to to ensure that the model is reliable under various conditions. The model was verified using a set of independent data (14 points) as shown in Figure 6b.

Table 3. Result of calibrated hydraulic parameters.

Layer No.		Layer Type	Layer Thickness (m above MSL)	Hydraulic Conductivity (m/Day)		Storage Coefficient	Specific Yield	Effective Porosity
				Kh	Kv	SS (-)	Sy (1/m)	neff (%)
Initial parameters	1	Clay	13.2 to −42	0.10–0.25	0.01–0.025	10−3	50–60	50–60
	2–5	Fine sand with lenses of clay	1.1 to −398	5–20	0.5–2	5 × 10−3	30	30
	6–9	Coarse sand quaternary	−103.2 to −757.6	20–75	2–7.5	2.5 × 10−3	25	25
	10–11	Graded sand and gravel	−218 to −939	75–100	7.5–10	5 × 10−4	20	20
Calibrated parameters	1	Clay	13.2 to −42	0.15–0.35	0.015–0.035	0.6 × 10−3	50–60	50–60
	2–5	Fine sand with lenses of clay	1.1 to −398	25–40	2.5–4	6.5 × 10−3	30	30
	6–9	Coarse sand quaternary	−103.2 to −757.6	40–130	4–13	3 × 10−3	25	25
	10–11	Graded sand and gravel	−218 to −939	130–150	13–15	2 × 10−4	20	20

presents the model verification results. The results of verification and sensitivity analysis showed that the most sensitive parameters are the hydraulic conductivity and specific storage.

3.2. Land Subsidence Calibratiom

Calibration land subsidence is the process of comparing the LS rate obtained from the analytical solution with LS rates from previous studies in the MND. In this section, the analytical model is used to calculate land subsidence in the period from 2010 to 2016. The abstraction rate of GW is indicated in

Table 4. Land subsidence model calibration scenario.

Year	2010	2016
Population (Million)	82.2	94.4
Abstraction rate (BCM/Year)	4.9	7
% Abstraction rate	100%	142%

. The analytical model results showed that the maximum LS is 170 mm from 2010 to 2016. So, the maximum rate of LS is 10.6 mm/year. Ls rates in the major cities of the MND are also measured as shown in Figure 7. LS rates in Kafr El-Sheikh, Tanta and Shibin El-Kom are 1.7 ± 0.5 mm/ year, 4.35 ± 0.6 mm/year and 8.75 ± 0.6 mm/ year, respectively. The analytical model results have been compared with the results of LS obtained by using INSAR and GPS techniques developed by Saleh and Becker [22]. Saleh and Becker investigated LS due to GW abstraction in the whole Nile Delta. Their results showed that LS rates due to excessive GW withdrawal were −6.4 ± 0.4 mm/y, −4.0 ± 0.6 mm/y, −4.8 ± 1.0 mm/y, −10.0 ± 1.2 mm/y, −10.3 ± 1.6 mm/year, and −4.9 ± 1.6 mm/y from Cairo, Tanta, Mahala, Mansoura, Damietta, and Port Said. Comparing the results from the analytical model and the results of INSAR gave a good agreement.

3.3. Effect of Decreasing Recharge Rates on LS in the MNDA

The reduction in the Nile flow results in the depletion of the recharge to the MND aquifer. In this stage of the case study, the model was used to assess the impact of decreasing recharge rates to 94.4 %, 88.8%, and 83.2% while maintaining the GW abstraction volume constant at 0.81 BCM y⁻¹. For scenario 2, the results show that when the recharge rate decrease to 94.4%, the maximum GW drawdown is 0.36 m from the starting GW level in 2010, and the maximum LS value is 30 mm (Figure 8a,b). In scenario 3, as the recharge rate drops to 88.8%, the difference between the GW levels reaches a maximum drop of 0.70 m from the starting levels in 2010, and the LS value reached a maximum of 60 mm (Figure 8c,d). Finally, in scenario 4 the recharge rate decreases to 83.2% resulting in a maximum GW drawdown of 1.05 m from the starting level in 2010, and a maximum LS value of 90 mm (Figure 8e,f). Figure 8 shows that the groundwater drawdown is displayed as concentric circles. We believe that this is due to the dense concentration of pumping wells in core of the study area. Moreover, the results indicate that the area with the highest LS rates is the core of the MND as shown in Figure 8.

3.4. Effect of Increasing Abstraction Rates on LS in the MNDA

Population growth leads to an increase in the GW abstraction from the MND aquifer [52]. In this part of the case study, the model evaluated the impact of increasing abstraction rates to 146 %, 193%, and 233 % while maintaining the GW recharge volume constant at 1.4 BCM y⁻¹. In scenario 5, the abstraction rate increased to 146% causing a maximum GW drawdown of 2.4 m from the starting levels in 2010, and a maximum LS value of 190 mm (Figure 9a,b). When the abstraction rate increased to 193% in scenario 6, the maximum GW drawdown reached a 5.0 m from the starting level in 2010, and the maximum LS value reached 380 mm (Figure 9c,d). For scenario 7, the abstraction rate increased to 233% resulting in a maximum GW drawdown of 7.0 m from the starting level in 2010, and a maximum LS value of 560 mm (Figure 9f). The intensive concentration of pumping wells in the core of the study area is responsible for presenting the groundwater drawdown as concentric circles as indicated in Figure 9. The results also indicated that the area with the highest LS rates is the core of the MND as shown in Figure 9.

3.5. Effect of Combination between Increasing Abstraction and Decreasing Recharge Rates on LS in the MNDA

In this stage, the model was used to assess the impact of both increasing the GW abstraction rates to 146%, 193%, and 233% and decreasing recharge rates to 94.4%, 88.8%, and 83.2%, respectively. The results showed that in scenario 8, when the abstraction rate increased by 146% and the recharge rate decreased to 94.4%, the maximum GW drawdown was 2.8 m from the starting level in 2010, and the maximum LS value was 220 mm (Figure 10a,b). In scenario 9, the abstraction rate increased to 193% while the recharge rate decreased to 88.8%, resulted in a maximum GW drawdown of 5.6 m from the starting levels in 2010, and the maximum LS value reached 440 mm (Figure 10c,d). For scenario 10, when the abstraction rate increased to 233% and the recharge rate decreased to 83.2%, the maximum GW drawdown was 8.5 m from the starting GW level in 2010, and the maximum LS value was 650 mm (Figure 10e,f). As shown in Figure 10, the groundwater drawdown is represented as concentric circles because of the intensive concentration of pumping wells in middle zone of the study area. In addition, the results indicated that the area with the highest LS rates is the core of the MND as shown in Figure 10.

4. Discussion

In this study, different scenarios were applied to the case study (MNDA) in order to predict the future behavior of LS in the ND. Land subsidence was predicted for both the current and future situations. By comparing the tested scenarios, it was clear that the minimum GW lowering obtained from scenario 2 with 0.36 m from the starting GW in 2010. While, the maximum GW lowering obtained from scenario 10 with 8.5 m from the starting GW in 2010 as shown in Figure 11a. Also, Figure 11b shows that the minimum LS value obtained from scenario 2 equals to 30 mm. While the maximum LS value obtained from scenario 10 equals to 650 mm. As a result of the previous scenarios, decreasing recharge rates to 94.4%, 88.2%, and 83.2% leads to increasing LS rates by 1.5 mm/y during the period of 2010–2070. Increasing abstraction rates to 146%, 193%, and 233% leads to increasing LS rates by (9 + 0.5) mm/y during the period of 2010–2070. Finally, the combination of increasing abstraction rates to 146%, 193%, and 233% and decreasing recharge rates to 94.4%, 88.2%, and 83.2%, respectively, leads to increasing LS rates by (10 + 0.5) mm/y during the period of 2010–2070.

5. Conclusions

Land subsidence is a serious phenomenon that is associated with aquifer system compaction. In this study, the GW flow model MODFLOW was merged with the analytical solution based on Terzaghi’s theory for assessing the land subsidence in the Nile Delta. Different scenarios of decreasing recharge, increasing abstraction, and the combination of the two were applied to estimate land subsidence in the MND. For the first case, the results showed that there is an obvious effect of decreasing recharge on LS in the MNDA. The maximum values of GW drawdown and land subsidence were 1.05 m, 90 mm, respectively, after 60 years since 2010 when the recharge rate decreased to 83.2%. For the second case, the results showed that there is an obvious effect of increasing abstraction on LS in the MNDA. The maximum values of GW drawdown and land subsidence were 7 m and 560 mm, respectively, after 60 years until 2010, when the abstraction rate increased to 233%. For the third stage, the results showed that there is an obvious drawdown in GW head due to combination of decreasing recharge and increasing abstraction that increased the LS in the MNDA. The maximum values of GW drawdown and land subsidence were 8.5 m and 650 mm, respectively, after 60 years, until 2010 when the recharge rate decreased to 83.2% and the abstraction rate increased to 233%. Additionally, the results showed that this stage is the most critical case. The results demonstrate that changing groundwater abstraction and recharge rates should be controlled to mitigate land subsidence and achieve proper management of GW resource. This study gave an insight to decision makers to consider the effect of changing the rates of groundwater pumping and recharge to groundwater aquifers to avoid more land subsidence which may damage many properties.