Numerical groundwater modelling as an effective tool for management of water resources in arid areas

Abstract

Groundwater, possibly of fossil origin, is used for water supply in some arid regions where the replenishment of groundwater by precipitation is low. Numerical modelling is a helpful tool in the assessment of groundwater resources and analysis of future exploitation scenarios. To quantify the groundwater resources of the East Owienat area in the southwest of the Western Desert, Egypt, the present study assesses the groundwater resources management of the Nubian aquifer. Groundwater withdrawals have increased in this area, resulting in a disturbance of the aquifer’s natural equilibrium, and the large-scale and ongoing depletion of this critical water reserve. Negative impacts, such as a decline in water levels and increase in salinity, have been experienced. The methodology includes application of numerical groundwater modelling in steady and transient states under different measured and abstraction scenarios. The numerical simulation model developed was applied to assess the responses of the Nubian aquifer water level under different pumping scenarios during the next 30 years. Groundwater management scenarios are evaluated to find an optimal management solution to satisfy future needs. Based on analysis of three different development schemes that were formulated to predict the future response of the aquifer under long-term water stress, a gradual increase in groundwater pumping to 150% of present levels should be adopted for protection and better management of the aquifer. Similar techniques could be used to improve groundwater management in other parts of the country, as well as other similar arid regions.

Editor D. Koutsoyiannis; Associate editor X. Chen

Citation El Alfy, M., 2014. Numerical groundwater modelling as an effective tool for management of water resources in arid areas. Hydrolological Sciences Journal, 59 (6), 1259–1274. http://dx.doi.org/10.1080/02626667.2013.836278

Résumé

Les eaux souterraines, quelquefois très anciennes, sont parfois utilisées pour l’alimentation en eau dans les régions arides, où la recharge des eaux souterraines par les précipitations est faible. La modélisation numérique est un outil précieux pour l’évaluation des ressources en eaux souterraines et l’analyse des scénarios d’exploitation futurs. Pour quantifier les ressources en eaux souterraines de la région de l’Est Owienat dans le Sud-Ouest du désert occidental égyptien, la présente étude évalue la gestion des ressources en eaux souterraines de l’aquifère nubien. Les pompages d’eau souterraine ont augmenté dans cette région, entraînant une perturbation de l’équilibre naturel de l’aquifère et la dégradation continue et à grande échelle de cette réserve en eau cruciale. Des impacts négatifs, comme la baisse du niveau piézométrique et l’augmentation de la salinité, sont connus. La méthodologie consiste en l’application de la modélisation numérique des eaux souterraines en permanent et en transitoire pour différents scénarios de prélèvement. Le modèle de simulation numérique développé a été appliqué pour évaluer les réponses piézométriques de l’aquifère nubien à différents scénarios de pompage pour les 30 prochaines années. Les scénarios de gestion des eaux souterraines retenus sont évalués afin de trouver une solution de gestion optimale pour la satisfaction des besoins futurs. Basée sur l’analyse de trois projets de développement qui ont été formulés pour prévoir la réponse future de l’aquifère au stress hydrique à long terme, une augmentation progressive de pompage des eaux souterraines à 150% des niveaux actuels devrait être adoptée pour la protection et une meilleure gestion de l’aquifère. Des techniques similaires peuvent être utilisées pour améliorer la gestion des eaux souterraines dans d’autres régions du pays, ainsi que dans d’autres régions arides similaires.

Key words:

• groundwater management
• numerical modelling
• Nubian aquifer
• East Owienat
• Egypt

Mots clefs:

• gestion des eaux souterraines
• modélisation numérique
• aquifère nubien
• Est Owienat
• Egypte

INTRODUCTION

Approximately one-third of the world’s land area is arid or semi-arid (Rogers 1981). These regions face major challenges in the management of scarce freshwater resources under the pressures of population and climate changes. Groundwater maybe the only source of water supply in such regions, and so needs quantification to manage its limited resources. This study presents a model of groundwater exploitation management that develops an optimal pumping rate under an expected increase in irrigation water requirements. The East Owienat area, in the extreme southwest corner of the Western Desert of Egypt, was selected in this study of groundwater resource quantification due to the fast growing need for utilization of non-renewable groundwater resources. Groundwater is mainly used for irrigation, domestic, and small industrial purposes. Climate conditions have changed over the last 10 000 years from more humid conditions after the last Ice Age to arid conditions during the last 5000 years (Edmunds and Wright 1979, Edmunds 1994, Issar 2003). Under present climatic conditions, the rate of groundwater recharge is very low, so the groundwater resources will be exhausted in the near future if abstraction continues with present trends. Due to the increase in population and expansion of the irrigated area, groundwater management is under environmental constraints. Determining the availability of groundwater resources in the Owienat area to meet present and future water demands is done using a simulation model to understand the physical behaviour of the Nubian aquifer system in southwest Egypt. It recognizes the need for sustainable use of water, balancing the long-term use with the long-term water availability, and adopts the best practice from traditional and new methods of water management. This requires social recognition of water scarcity and decision-makers will confront difficult societal choices and need good scientific understanding to develop and support sustainable management solutions.

The Oweinat area lies within the Great Desert belt in southwest Egypt (Fig. 1) which is characterized by high temperatures and a hyper-arid climate zone; the aridity index is <0.05 (UNEP 1997). This work describes the history of groundwater pumping from the Nubian aquifer, the hydrogeological setting, the resulting negative impacts on groundwater levels and the various stages involved in the development of a three-dimensional (3D) simulation model of the interactive aquifer system followed by a prediction stage to evaluate the effects of long-term pumping on groundwater in the study area. An improved groundwater pumping scheme is suggested to improve the long-term management of the aquifer system in the area. Due to growing population, there have been many development projects in this area within the last two decades, including land reclamation and new urban communities. Due to the mismanagement of aquifer systems by over pumping from a large number of wells clustered in small areas, severe environmental problems have grown, with a continuous decline of the piezometric surface by >8 m during the studied 16 years.

Fig. 1 Location map of the East Owienat area and of wells tapping the Nubian aquifer. A–A′ and B–B′ and C–C′ are the locations of hydrogeological cross sections.

The objective of this study is to develop the optimal pumping rate under the expected increase in irrigation water requirements covering the newly reclaimed areas. This will be achieved through numerical groundwater modelling. The model results document a method to quantify the sustainable yield of a groundwater flow system. The transient calibration was used to simulate potential future water use scenarios.

GEOLOGICAL BACKGROUND

The Nubian sandstone in the Western Desert attains a thickness of more than 1000 m. It is subdivided into six formations (the Six Hills, Abu Ballas, Sabaya, Maghrabi, Taref and Mut formations) ranging from Jurassic to Upper Cretaceous (Klitzsch et al. 1979, Klitzsch and Lejal-Nicol 1984). This thick sandstone unconformably overlies the basement rocks (Fig. 2) that are crosscut by an extensive E–W fault system in southern Egypt (Issawi 1971, 1978). Lithological logs of water wells tapping the Nubian aquifer show that the aquifer is composed of intercalations of sandstone, varicoloured and medium to very coarse grained, with interbeds of grey to green, sticky shale. The Nubian aquifer is unconformably overlain by Quaternary deposits that are composed of aeolian sand dunes and sand sheets, alluvial deposits, playa deposits and inland sabkha (Mansour et al. 1982, Diab 1984, Heinl and Thorweihe 1993, Ebraheem et al. 2004). The E–W fault system can be identified in cross sections, but is difficult to trace on the ground surface due to the aeolian sand dunes and sand sheets. There are local basins within the Nubian sandstone, which are characterized by the presence of fault planes exhibiting no-flow boundary conditions. Some discrete local structural lows are bounded by no-flow boundary conditions, as revealed by long-duration pumping tests (GARPAD 1994).

Fig. 2 (a) Key map of cross sections; (b) hydrogeological cross section A–A′; (c) hydrogeological cross section B–B′ (El Nahry et al. 2010); and (d) hydrogeological cross section along the line C–C′ in East Owienat.

METHODOLOGY

The approach combines information and procedures from various backgrounds to simulate the development of the groundwater resources and account for the complex boundaries and spatial variability. Understanding the Nubian aquifer system required (a) delineation of the aquifer boundary limits, (b) determination and correlation of the stratigraphic horizons, (c) characterization of the structural setting, and (d) visualization of the aquifer geometry. Well logs from 20 boreholes have been used to delineate the horizontal and vertical stratigraphic unit offsets. Several geological cross sections were drawn in different directions (Figs 1 and 2).

The water table data and land-use information were collected during two field trips between 2007 and 2011. Archived and collected data were analysed by means of statistical methods to obtain a good understanding of the prevailing hydrogeological conditions. These involved collection of pervious data (GPC 1985, GARPAD 1994, MALR 1997, 2007, NMA 2005), which were used in constructing water-level contour maps. To mathematically simulate the behaviour of groundwater in the East Owienat area, detailed conceptual understanding of the hydrogeology was gained through constructing different hydrogeological cross sections, using water-level contour maps of different dates and reviewing the aquifer hydraulic parameters. The modelling process began with data integration and description of the hydrogeological situation of the study area, such as the boundary conditions, then definition of the conceptual framework of the area for analysis and finally the numerical simulation. To establish an interpretive model, the groundwater flow model was calibrated for steady-state conditions using the geological and hydrogeological characterization. Then the calibrated transient condition model was established using the calibrated steady-state model.

Table 1 Changes in hydraulic head, drawdown and depth to groundwater for the different scenarios.

Scenario	Well abstraction OUT (m^3 year^-1)	Constant head IN (m^3 year^-1)	Storage IN (m^3 year^-1)	Hydraulic head (m a.s.l.)	Drawdown (m)	Depth to groundwater (m b.g.l.)
I	6.9 × 10^8	4.9 × 10^8	2.0 × 10^8	234.0	15.2	42.7
II	1.1 × 10^9	7.3 × 10^8	3.2 × 10^8	225.7	23.3	50.8
III	1.4 × 10^9	9.7 × 10^8	4.3 × 10^8	219.0	30.1	57.6

Hydrogeological setting

The Owienat area is located in the hyper-arid region of Western Desert, Egypt, with minimal precipitation. Nour (1996) and Sadek (1996) investigated the desert climate in this area showing that: (1) average temperature ranges from less than 6°C in winter to more than 40°C in summer, (2) the precipitation is scarce and does not exceed 6 mm year^-1, natural evaporation rate ranges from 10 mm d^-1 during January to 31 mm d^-1 during July, with an average of 22.2 mm d^-1, and (3) the relative humidity ranges from 11% in May to 41% in December. The groundwater of the study area is fossil water and is non-renewable. In the Dakhla and Kharga oases, 400 km downstream to the north of the study area, isotope studies have shown that climate conditions changed and the range of water ages derived by krypton-81 and chlorine-36 is 0.2–1.0 × 10⁶ years (Sturchio et al. 2004). In the study area, the Nubian aquifer water is dated as 14 000 years old (Schneider and Sonntag 1985). This contradiction can be explained because the study area is located 400 km upstream of the Dakhla and Kharaga oases, so nearer to the recharge areas in Chad and Sudan. There is no evidence of aquifer recharge from the Nile which has a younger isotopic profile; this is because the basement rocks (exposed on ground surface) to the east of the study area act as a solid barrier preventing recharge from the Nile.

More than 720 production wells were drilled within the area up to 2001. According to its stratigraphic position, permeability and storage capacity, the Six Hills Formation is considered to be the main water-bearing strata. The horizontal and vertical extension and boundaries of the aquifer were defined as a prerequisite to representing the geometric framework of the groundwater flow modelling. Data pertaining to the subsurface geology in the study area were derived from well logs available in the archive. The Nubian aquifer is represented by a thick sandy section, composed of sandstone with a subordinate kaolinitic sandstone. The maximum thickness of the sandy section, reported in GPC well no. VII, is 720 m. The basement rocks outcrop at Bir Safsaf and Qaret El Mayet, and are encountered at various depths ranging from 200 to 720 m at well nos 32 and VII, respectively. The saturated thickness of the aquifer is highly variable ranging between 137.25 and 659 m in wells 7/1/12 (El Awkaf) and VII (GPC), respectively. In general, where covered by Quaternary sands, the Nubian aquifer is unconfined but semi-confined conditions occur at the extreme eastern and northwestern parts where it is overlain by a considerable thickness of clayey beds (Fig. 2). The lower boundary of the aquifer system is represented by the impervious basement rocks. The depth to the basement rocks, and corresponding thickness of the Nubian aquifer, increases northwestward. The northeast boundary is represented by the basement rock outcrops (Qaret El Mayet, Nusab El Balgum and Atmour El Kebyish). The western boundary of the aquifer is unclear, because the basement rocks of G. Kamel have a relatively small peak shape and local extent, and the basin is extends westward. Consequently, the western part of the East Owienat area may have a considerable reserve of groundwater. The interpolated depth to water map shows that the depth increases from the northeast part of the study area to the southwest (Fig. 3), while the ground elevation increases toward the southwest. The depth to water ranges from 2.3 to 58.69 m at Bir Safsaf and well no. VII (GPC), respectively.

Fig. 3 Depth to water contour map of East Owienat area (1984).

Water-level measurements from 1984 to 2011 were used to determine groundwater flow directions and the location of recharge and discharge areas. Three regional iso-potentiometric maps referenced to sea level were prepared for the aquifers for 1984, 2003 and 2011. The water-level maps for 1984, 2003 and 2011 were interpolated using fitted variogram models (Fig. 4(a), (b) and (c), respectively). Owing to low abstraction rates prior to 1984, the static water-level contour is considered to represent the steady-state conditions of the aquifer (Fig. 4(a)). The piezometric level slopes from southwest to northwest, from 264 to 244 m a.s.l., with non-uniform gradients of 0.02% to 0.04% from southeast to north, respectively. The gradients imply that the aquifer is recharged from the southwest, and naturally discharges to the northeast. A dramatic change in the groundwater flow regime since 1984 has been reported showing a regular northward flow. In 2003, a general northward flow persists but with a discrete water trough at the center, which might be attributed to abnormal increases in the rate of groundwater extraction at well no. 2 (Fig. 4(b)). The hydraulic gradient varies between 0.02% and 0.07%. The higher gradient south of this water trough is perhaps due to a higher rate of groundwater abstraction. In 2011, the water trough progressively migrated towards the northwest at well no. 6 and was bigger than that developed in 2003 (Fig. 4(c)). The groundwater flow changed from a general northward flow to flow into the trough. The water trough becomes more linear and oriented from the southeast to the northwest. This linearity may suggest that a fault constrains the flow direction, as it might be considered as a path of accelerating water flow with advancing groundwater extraction at well no. 6. Thus, it is inferred from Fig. 4(a), (b) and (c) that the aquifer system is recharged from the southwest part of the study area, and naturally discharges to the northeast.

Fig. 4 Iso-potentiometric level maps of Nubian sandstone aquifer (a) 1984, (b) 2003, and (c) 2011.

Water-level monitoring from 1995 until 2011 was carried out in several observation boreholes (Fig. 5). Measurements in most of the boreholes record a decline of water level with time; evidently the heavy pumping of the aquifer has caused significant drawdown. The abstraction rate was relatively low until 1984; since then it has increased tremendously, which seemingly corresponds to the progressive population increase and advancing artificial irrigation of farmland (Fig. 4(a), (b) and (c)).

Fig. 5 Depth to water (m) in wells 32, 7, 9 and 11 of the Nubian aquifer.

The average increase of depth to water, relative to ground elevation, is about 8 m over the last 16 years (Fig. 5), which corresponds with Ghoubachi’s (2004) conclusion that the water level declined by 0.5 m year^-1. Evidently, the former remarkable decline of water level is due to the drilling of more than 720 production wells within the study area, and the consequent increase in extraction rate of Nubian aquifer groundwater from 3.6 × 10⁴ m³ to 1.6 × 10⁶ m³ between 1984 and 2011 (Fig. 6). The hydraulic conductivity values range from 10 to 20 m d^-1, and the aquifer specific yield is 0.16 to 0.23 (unconfined part), whereas storativity is 10^-3 in the semi-confined part of the system (El Sayed et al. 2004).

Fig. 6 Yearly extraction rate of groundwater (m³) from the Nubian aquifer between 1984 and 2011.

Development of the conceptual model

The development of a conceptual framework for the flow model requires the detail of the three main components: the hydrogeological framework, the nature of the flow system parameters and the water budget. Also the aquifer homogeneity, fracture system, boundary conditions and recharge are needed. The hydrogeological properties of the Nubian aquifer in southwest Egypt used in the flow model involve its lithology and boundary conditions. The aquifer is composed mainly of medium to very coarse grained sandstone with few shale intercalations. The base of the aquifer system was specified as a no-flow boundary as the Nubian sandstone overlies the Precambrian basement rocks. A thin Abu Ballas shale layer covers the aquifer in the north of the study area. The aquifer shows lateral and vertical facies changes and its thickness varies from south to north. The groundwater flows locally from southwest to northeast. The aquifer system has been in a transient-state condition since the current arid period began about 5000 years ago (Heinl and Thorweihe 1993, Churcher et al. 1999, Ebraheem et al. 2002, Gossel et al. 2004). Due to the prevalence of hyper-arid conditions, present day aquifer recharge is excluded; but the discharge rate has increased incrementally due to the more than 700 productive wells. The long-period simulation (10 000 years) of the impact of present and future groundwater extraction rates from the non-replenished Nubian sandstone in southwest Egypt indicates that the aquifer is still under the influence of the past humid period and has been in an unsteady depletion process (Ebraheem et al. 2002).

Appropriate space and time scales were chosen in relation to the heterogeneity of the system and the data available for calibration. The area hosts a thick, water bearing sandy horizon which is often intercalated with clay lenses occurring at variable depths. As evident from the previous sections, the groundwater flow is dominated by regional and sub-regional flow patterns (Fig. 4). The selected model area covers around 6500 km² (Fig. 7(a)). As the physical boundaries of the aquifer stretch in all directions and there is not enough information available about the boundaries, the model’s hydraulic boundaries were chosen according to the hydraulic information gathered from the drilled wells, such as total water head, water flow direction and magnitude and the hydraulic parameters of the aquifer. The hydrogeological units are based on the stratigraphic units, but the method of De Marsily et al. (1986) was adopted to represent the complex geological formation within the model area. The different hydrostratigraphic units of similar properties were identified by using the stratigraphic boundaries and the total hydraulic head. The thick sandstone section of the Nubian aquifer was subdivided into hydrogeological sub-units separated by silts and shale horizons; this hydrogeological configuration took into account the validity of the governing equations for multi-layered aquifers (Connorton 1985). The total range of integration of the aquifer is split over z (vertically), while it is assumed that qx and qy (horizontally) are sufficiently smooth over each sub-interval of integration.

Fig. 7 (a) Location map of the groundwater model area. (b) Schematic diagram showing horizontal and vertical flows in the model layer.

Based on the integration of hydraulic head, transmissivity values and lithology, a conceptual model of the flow system was constructed. Flow patterns of the selected model layers show different vertical and horizontal flow directions (Fig. 7(b)). The assumptions taken into consideration during component quantification of the flow model are as follows:

1. The flow occurs under confined and leaky conditions,

2. The aquifer was in steady-state conditions prior to 1984 that is suggested as a base year,

3. The aquifer is treated as in an unsteady state post-1984 due to overpumping activities,

4. The water level drops slightly to deal with it as in a horizontal attitude, to which its depth is averaged, and

5. The aquifer can be subdivided into horizons of continuous lithological homogeneity.

Estimation of the water budget involves all possible recharge and discharge components, including inflow and outflow within the model domain. The suitable water balance equation must be satisfied at the zero time of the simulation. The base year (1984) is considered for the balance calculations. The aquifer is mainly recharged by inflow of groundwater across the southern and western boundaries. The discharge components include pumping and outflow of groundwater to the northern and eastern parts of the study area.

Development of the numerical model

Hydrogeological numerical modelling simulates the average behaviour of an aquifer and can provide a range of potential future conditions that can satisfy the management objectives of any developmental phases. These models use finite differences to treat 3D saturated and unsaturated flow. The US Geological Survey’s finite-difference groundwater model, MODFLOW, performs such groundwater simulation. It has been assumed that the groundwater system in the modelled area occurred in steady-state conditions, with negligible head variations prior to 1984. So a steady-state model can be used to simulate the initial conditions. The governing equations for flow are described by:

(1)

where K_xx, K_yy and K_zz are values of hydraulic conductivity [L T^-l], along the x, y and z coordinate axes, respectively; h is the potentiometric head [L]; W is the volumetric flux per unit volume and represents sources and/or sinks of water [T^-l]; S_s is the specific storage of the porous material [L^-1]; and t is time [T].

Equation (1) implies a long-term balance between natural recharge and discharge in the area. Therefore, the term δh/δt is neglected. A tentative water balance was calculated with the help of software recommended by Stoertz and Bradbury (1989) and Anderson and Woessner (2002). In assigning the model grid, all finite different cells were given specified head values. The 2D and 3D meshes and the graphical interface of MODFLOW were used to calculate the total inflow and outflow through the boundaries, and to define the recharge and discharge. However, MODFLOW has easier tools to make a good conceptual model, and check and calibrate the model. The water budget is very dependent on the estimated average transmissivity of the different hydrogeological units and on the interpolated heads taken from the water-level contour map (Fig. 4(a)). The total inflow and outflow in the calibrated model is accepted at within 10% of the estimated value. The groundwater recharge values have resulted in a significant vertical flow component (return flow). The upward flux (-ve) represents the capture zone around the pumped well of interest where boundary inflow and outflow are calculated.

The study area was discretized into a non-uniform square grid, composed of 161 rows and 140 columns, with an average grid spacing of 1520 m outside the major pumping field and 760 m inside the major pumping fields. The grid spacing was determined in view of the available data, computational time and computer storage requirements. The ground surface elevation of the study area was calculated from the digital elevation model (Shuttle Radar Topography Mission, SRTM 90m DEM) and corrected according to detailed topographic maps. The outline of the grid is determined by the available geological and hydraulic information. The boundaries were placed sufficiently far from the effective pumping wells that the pumping effect does not reach the boundary within the time period used in transient simulation. Using the schematic stratigraphic cross section (Fig. 2), the model layers were discretized according to the hydro-stratigraphic units.

All boundaries of the model were assigned as specified head boundaries to allow inflow to and outflow from the system. A specified head was also justified because the general head boundary conditions in MODFLOW assume a constant head away from the model boundary. The constant head (Dirichlet head) and no-flow boundaries were assigned (Fig. 7(a)). The constant head boundary was used where the groundwater levels were known or could be interpolated (southern and northern boundaries). The no-flow boundary was set when the flow line is almost parallel to the border of the model (eastern and western boundaries). Additionally, singular point sources were assigned at borehole locations to simulate pumping rates. It is assumed that impermeable basement rocks represent the base of the model.

The initial conditions in MODFLOW define the initial status of the pressure head and the moisture content. The initial conditions were set to cold start with pressure head computed from the constant heads, where these initial conditions for the flow equation are given by: h = h_i (x, y, z) in R, where R is the region of interest and h_i is the prescribed initial condition, which can be obtained by field measurement. For the fluid properties, the density of water was set as 1000 kg m^-3, the viscosity was set to 4.68 kg m^-1 h^-1. The compressibility of water was set to 3.4 × 10^-17 m h^-2 kg^-1, and the acceleration due to gravity was set as 1.27 × 10⁸ m h^-2. The aquifer hydraulic properties for the beginning of the calibration were estimated from pumping test interpretation (Nour 1996, Ghoubachi 2004). The adopted mathematical model (governing equations, boundary conditions and initial conditions) was solved using the point-wise iterative matrix solver.

Regional models at best simulate the average head or the drift (Christensen and Cooley 1999), so the variance of the residuals between the simulated and the measured head lies within the range of the head variance. For the steady-state simulation, the horizontal and vertical hydraulic conductivity of the aquifer’s aquitard and aquifuge layers are needed. The MODFLOW model calculates the transmissivities and vertical leakage terms from the thickness of the model layers implicitly, given by the top and bottom elevation of each layer. Steady-state conditions were calibrated by means of the head measurements of 1984. The homogenous aquifer sub-units were assigned uniform flow and constant hydraulic conductivity. In areas where no data exist, the hydraulic conductivity was assigned within values gathered from the available hydrogeological information. Calibration of steady-state conditions began with adjustment of the boundary conditions and return flow recharge rates, while hydraulic conductivity values were kept constant, until a scatter-plot, roughly coinciding with a 45° line (1:1 fit) was obtained. Under steady-state flow conditions, the estimated head variance provides an estimate of the expected range of model residuals (h measured–h simulated). The head loss in the model area was estimated to be 17 m. Model calibration involves the adjustment of the vertical leakage rates, boundary flow rate and hydraulic conductivity distribution. Model fit to the calibration target was assessed using the root mean square (RMS) error and scatter plot of measured vs simulated heads close to a 45° sloping line. The calibration process involves calibrating both steady-state and transient conditions. Head calculations for the steady-state run showed a great difference between the calculated and the observed heads (eight observation wells), so the hydraulic conductivity and transmissivity were changed by the trial-and-error method (Doherty 1990) until the calculated heads matched the observed heads. Under transient conditions, the heads resulting from the steady-state simulation were used as initial hydraulic heads. The model was run under these new conditions and the heads were checked, against the actual measurements in 2003, 2006 and 2011. The hydraulic parameters and boundary conditions are changed until the observed and calibrated heads were comparable. The model calibrated over the period 1984 to 2006 is justified and so can be used to demonstrate impacts of different scenarios on the Nubian sandstone aquifer up to 2041.

RESULTS

Steady-state conditions

The calculated boundary inflow and outflow were estimated as 9.4 × 10⁴ m³ d^-1. The water budget from the numerical model lies within the tolerable range of the previously estimated total budget (8.9 × 10⁴ m³ d^-1). After calibration of the boundary conditions and the return flow rate, the reducing head variance was adjusted by modification of the hydraulic conductivity distribution (Fig. 8(a)). The hydraulic conductivity varied between 1 × 10^-2 and 5 × 10^-5 m s^-1, which is consistent with values calculated from pumping tests carried out on sandstone water bearing formations (Nour 1996, Ghoubachi 2004), and sandstone values in the literature. The total head varies between 242 and 258 m, while the RMS error is 1.05 m (Fig. 8(b)). It is evident that the total head of groundwater in the study area correlated well with the initial heads of 1984 (Fig. 4(a)).

Fig. 8 (a) Hydraulic conductivity distribution and (b) total head distribution (m) of the pre-development calibrated steady simulation flow model (1984), and (inset) calculated vs observed head of the calibrated steady flow model (1984).

The pressure head was computed by subtracting the node elevations from total heads and values varied between −58.5 and 660.6 m (standard deviation 205.5 m). Negative values of pressure head are recorded in the southwest of the model where cells are located in the aeration zone above the water table; the ground elevation exceeds 470 m and the water level is at about 260 m. Higher pressure heads characterize the northeastern part, where the cells are below the water level. The computed nodal velocity varies between 3.2 × 10^-6 m d^-1 and 0.01 m d^-1, with a standard deviation of 0.04 m. The flow velocity increases southwestward as a consequence of higher hydraulic gradients.

Transient conditions

Based on the pattern of the aquifer parameters obtained during steady-state calibration, the model was subjected to transient calibration for a period of 27 years between 1984 and 2011. The total simulation period for transient calibration is 9855 d. It was divided into three stress periods and 30 time steps. During each stress period all external stresses (well extractions and return flow) and boundary conditions were kept constant. The transient run was simulated from the calibrated steady flow model. It was calibrated against the measured heads obtained from field observations. At this stage of calibration, variation in storativity values (specific storage) and hydraulic conductivity are of primary importance for transient model calibration. The simulated piezometric levels at 2003 and 2006 are shown in Fig. 9(a) and (b), respectively.

Fig. 9 Total head distribution (m) of the calibrated transient flow model: (a) 2003, (b) 2006 and (c) 2011.

The calibrated transient model resulted in higher natural inflow and lower outflow rates, compared to the steady-state flow model. The calibrated model for 2003 shows that the inflow and outflow rates are higher than the steady-state flow model (1984). The inflow attains 7.8 × 10⁵ m³ d^-1, and the extraction rate increased to 2.7 × 10⁵ m³ d^-1. The calibrated transient model for 2006 implies that the inflow and outflow exceed that of the transient model for 2003 (Fig. 9(a) and (b)). In 2006, the inflow attains 1.2 × 10⁶ m³ d^-1, and the wells abstraction rate increased to 6.9 × 10⁵ m³ d^-1. The Egyptian Government plans to drill about 383 water wells in the East Owienat area (Robinson et al. 2007), which might cause severe pumping rates and a likelihood of increase in the size of the groundwater troughs, consequently affecting the hydraulic gradient of the model area.

A sensitivity analysis of the model indicated that it is less sensitive to changes in specific storage, while hydraulic conductivities and the boundary conditions were identified as the most sensitive parameters for the cases considered. The best scenario for model verification is one that uses the calibrated model to simulate the aquifer under stressed conditions. The processes of model verification required further calibration refinement of the model. When the model has successfully reproduced measured changes in field conditions for both the calibration and historical time periods, it is ready for predictive simulations. The model verification was carried out with the final run of the transient calibration stage which resulted in prediction of the potentiometric surface at the end of 2011 (Fig. 9(c)); a scatterplot compares the simulated and observed. The model for 2011 shows that the inflow and outflow rates are greater than in the transient models of 2003 and 2006; the inflow reaches 2.1 × 10⁶ m³ d^-1, and the wells abstraction rate is increased to 1.9 × 10⁶ m³ d^-1.

Model predictions and aquifer potential

Unfortunately, human activities, such as groundwater extraction for irrigation, have resulted in excessive groundwater level declines and aquifer overdraft. A successfully calibrated and verified simulation model can be used for formulating various different water development schemes. The schemes can be compared in terms of their feasibility for efficient utilization of the available groundwater resources. The selected development schemes must, when implemented, meet the future water demand at a minimum cost commensurate with legal, organizational, political and environmental considerations. Three scenarios were considered to predict the consequences of future development in the model area. A planning period of 30 years (2011–2041) was used for the three scenarios, related to development plans of the Ministry of Agriculture and Land Reclamation (MALR 2007). The planning period is longer than the period for which the model was calibrated and validated (27 years). The starting conditions in each scenario were those obtained from the transient simulation at the end of 2011. Estimation of abstraction rates for the planning period (2011–2041) was used in establishing reliable alternative schemes. Three scenarios, which might be realistic as an increase in pumping rate is directly proportional to the increase in population and other developments, are as follows:

Scenario I The first scenario assumes that the groundwater abstraction observed during 2011 continues without modification, i.e. 1.9 × 10⁶ m³ d^-1 (6.9 × 10⁸ m³ year^-1), and the irrigated area is nearly 110 000 acres. The resulting drawdown contour map (Fig. 10) shows a cone of depression with a maximum drawdown of about 15 m at the center of the modelled area where agricultural development activities prevail (Table 1).

Fig. 10 Calibrated transient flow model for Scenario I—total head distribution (m) for: (a) 2021, (b) 2031 and (c) 2041; and (d) drawdown distribution (m) for 2041.

Scenario II Abstractions from the Nubian sandstone aquifer are assumed to increase by 150% from those of Scenario I, as the first development plan option of the MALR (2007) assumes an increase in cultivated areas by 150% (147 000 acres). This implies an increase in groundwater pumping from 1.9 × 10⁶ in 2011 to 2.9 × 10⁶ m³ d^-1 (1.1 × 10⁹ m³ year^-1) by the end of 2041. Fig. 11 shows the groundwater drawdown contour maps up to the end of 2041; the maximum drawdown in the agricultural area is about 23.3 m (Fig. 11(d) and Table 1).

Fig. 11 Calibrated transient flow model for Scenario II—total head distribution (m) for: (a) 2021, (b) 2031 and (c) 2041; and (d) drawdown distribution (m) for 2041.

Scenario III Abstractions are assumed to increase by 200% from those of Scenario I; the second development plan option of the MALR (2007) assumes an increase in cultivated areas by 200% (220 000 acres), which implies an increase in groundwater pumping from 1.9 × 10⁶ m³ d^-1 in 2011 to 3.8 × 10⁶ m³ d^-1 (1.4 × 10⁹ m³ year^-1) by the end of 2041. The model prediction shows that the maximum drawdown exceeds 30 m at the exploitation sites (Fig. 12 and Table 1).

Fig. 12 Calibrated transient flow model for Scenario III—total head distribution (m) for: (a) 2021, (b) 2031 and (c) 2041; and (d) drawdown distribution (m) for 2041.

CONCLUSION

Numerical groundwater modelling is a powerful tool for assessment, development and management of groundwater resources. A numerical groundwater model for the East Owienat arid area in Egypt was developed using MODFLOW software to simulate regional groundwater changes in the Nubian Sandstone aquifer under transient conditions. Hydrogeological data from different sources and previous studies were integrated in the model. Initially a conceptual model was developed comprising identification of hydrogeological units, parameterization of these units, identification of the system boundaries and identification of the water budget components. These were transferred into the mathematical groundwater model. All available measurements of water level, well discharges and hydraulic parameters were used for model calibration. The model was optimized through improved definition of the boundary conditions and other relevant parameters, such as well discharges and return flow. When the model calibration and validation were complete, different groundwater management scenarios were evaluated to find an optimal management solution to satisfy future needs. Three different development schemes were assessed to test the Nubian aquifer under long-term water pumping stresses. Scenario II, with a gradual increase of pumping for agricultural purposes over the next 30 years, to 150% of current rates, should be adopted for the protection and better management of the aquifer. This is the base level for a safe yield that maintains a balance between groundwater inflow and outflow in the Nubian aquifer of the modelled area with advancing agricultural activities. Much irrigation water is lost due to the low efficiency of the managed irrigation programmes, thus it is recommended to encourage modern irrigation techniques, which could easily lead to a major reduction in agricultural water use.

Acknowledgements

The author thanks the Deanship of Scientific Research and the Research Center at the College of Science at King Saud University for encouragement and support for this project. Also the author would like to thank the referees and editors for their thorough review and very useful comments that helped improve the clarity and the relevance of this paper.