Large-Scale Water Storage in Aquifers: Enhancing Qatar’s Groundwater Resources

Abstract

Qatar’s water resource has been largely overexploited, leading to the severe depletion of its aquifers and degradation of water quality due to saline intrusions. Qatar envisions employing regional aquifers to store water via forced injection of desalinated water and thus increase available from a few days to two months. A strategy for the implementation of forced injections is proposed based on a spatially distributed model of groundwater flow at the scale of the whole country. The model is based on calibration under steady-state flow conditions and for a two-dimensional single regional aquifer due to the lack of data. Injection scenarios include various mean injection rates at the scale of the whole system and are interpreted under the assumption that the additional storage should feed 2.7 M inhabitants for two months at a rate of 100 L/person/day. When this water supply stock is reached, the model is run to define the infiltration rate, which allows the stock to remain constant over time as a result of an even balance between infiltrations, withdrawals and also leaks or inlets through the boundary conditions of the system.

1. Introduction

Qatar is one of the driest countries in the world, with an average annual gross precipitation of 80 mm per year. There is no accumulation of water at the surface (rivers, lakes, ponds, etc.), except during rare but intense rainy episodes, with water running off and then accumulating over short periods before infiltration straight at the topographic depressions. The average annual rainfall is around 80 mm, which is not enough to recharge its aquifers. Groundwater is the only natural source of water suitable for consumption and was the only freshwater resource until the early 1960s, when desalination was introduced. This resource has been largely overexploited, leading to the depletion of aquifers and water quality degradation because of saline intrusions facilitated by a very weak piezometric gradient along the coastal areas. This is evident in high-salinity areas, especially on the eastern coast of northern Qatar around Al-Khor town. Despite an extremely arid climatic context, Qatar is one of the most water-consuming countries, with an average household consumption of more than 500 L/capita/day [1]. This is because of the high living standards and other socioeconomic factors [2].

To meet the ever-increasing demand for water, Qatar currently operates several desalination plants that produced a gross volume of 600 million m³ of water in 2016, becoming almost the sole drinking water resource for the country. These plants were produced just in time, and Qatar has only a few days of freshwater reserves [3].

Faced with this challenging situation, Qatar envisions forcing the natural reservoirs (aquifers) with injections of desalinated water and increasing the stocks and availability of water from a few days to two months. Unlike other means of storage, such as on-land tanks, aquifer recharge can accommodate large quantities of water and is less prone to damage or contamination. Compared to the mega storage tanks that Qatar has built that are capable of storing a one-week supply, aquifers can store water quantities for several months according to population needs. The objective of the present work is to propose a strategy for the implementation of recharge wells to meet at best the storage objectives.

The first part of this manuscript is devoted to the geological and hydrogeological context of the country. The second part deals with the hydrogeological model developed to build the recharge scenarios. The last part is dedicated to strategies fitting the volumes of water to be stored and the hydrogeological context.

2. Geological and Hydrogeological Context

2.1. Qatar’s Main Geological Structures

Qatar is a country located in the Persian Gulf. This country is a peninsula, sharing with the south a 70 km long land border with Saudi Arabia. The country is 180 km long and 80 km wide and has an area of 11,570 km². The topography is very flat, with altitudes varying between 0 m and 107 m.

The sedimentary formations visible at the outcrop in Qatar (Figure 1) are mainly limestone layers dating from the Neogene (23–2.5 Ma, BP) but with few occurrences of earlier formations dating from the Paleocene (65–55 Ma, BP) and Eocene (55–35 Ma, BP). Some Quaternary formations, such as beach deposits and Sabkha formations (evaporites), are observed along the southeast coast. In general, three main formations comprise the sedimentary pile visible at the outcrop in Qatar: the Umm er Radhuma formation, the Russ formation and the Dam and Dammam formation, from bottom to top.

2.1.1. Umm er Radhuma Formation

This geological layer inherits its name from a well “Umm Radmah” located in the locality of the same name in Saudi Arabia [5]. The first studies reporting on this formation date back to 1952 [6]. Significant outcrops are observed in the Arabian shield from Iraq to Saudi Arabia, whereas in Qatar, the Umm er Radhuma formation is only visible in a few topographic depressions. Data from the various drillings in Qatar indicate that Umm er Radhuma has an almost uniform thickness of approximately 300 m and mainly consists of chalk and limestone of the Paleocene age.

2.1.2. Rus Formation

This geological layer inherits its name from a hill located in Saudi Arabia, “Umm er Ru’us.” Powers et al. [7] used this terminology, replacing the usual name, “chalky zone,” which was initially assigned to the formation. In Qatar, it is possible to observe some outcrops, notably to the north of Doha (districts of Umm Salal, Abuthaylah and Al-Khor) and along the Dukhan Anticline Fold (Fhaihil and Dukhan domes), as illustrated in Figure 1. The thickness of this layer ranges from 20 m (in the center of the country) to 110 m (in the southeast). This layer consists of Eocene chalk, shale and some gypsum beds located in the south of the country.

2.1.3. Dam and Dammam Formation

This geological formation was first mentioned in 1956 by Thralls and Hasson, [8]. The formation is present over a large surface area of Qatar. The thickness of this layer varies from 30 to 50 m. It consists of limestones and dolomites dating from the Miocene.

2.2. Qatar’s Hydrogeology

Qatar’s aquifers are part of a larger system named the Eastern Arabian Peninsula. It originates in Saudi Arabia and is distributed from east to west to Bahrain and Qatar. These aquifers are not well documented, which renders difficult the characterization of their geometry and hydraulic properties. In the literature, the Qatar underground is often separated into three hydrogeological zones [5,9], mostly distinguished according to their difference in water salinity (Figure 2.).

The northern basin is the most exploited zone for drinking water resources. Groundwater can be described as “floating lenses” of freshwater over salty water [10]. Freshwater is usually located in areas of topographic depressions. The northern basin is bounded to the north by the shore with the Persian Gulf and to the south by an evaporite formation (from the southern basin). It is also the best-documented area, as it is the main natural freshwater resource of the country. Nowadays, this aquifer is almost exclusively exploited to support agricultural activity.

The southern basin is less documented than the northern part and is also less exploited because of the high salt concentrations in water (local concentrations reaching 5 g/L—Schlumberger Water Service [11]). More than half of the southern part of Qatar’s upper geological formations is dominated by evaporitic deposits that form discontinuous thin aquitards of lesser importance compared to the northern aquifer. These evaporites consist of thick, impermeable, compact gypsum beds overlain by a thin layer of microporous dolomitic limestone of the upper aquifer. The gypsum beds are underlain by a thick layer of carbonates representing the lower aquifer. This zone resembles for its permeable formations to a connected multi-layered system, located in the Rus formation.

The southwestern basin is the least documented hydrogeological zone and is often aggregated with the south basin. The literature reports on an aquifer of approximately 30 m thickness, occurring in the form of limestone, predominantly dolomitic, and interspersed with marl. The water is also very salty, making it unfit for consumption (local concentrations of 2 to 4 g/L, Al-Hajari [5]).

2.3. Groundwater Recharge

Qatar’s climate is desert-like, with very hot summers and mild winters. In summer, temperatures can vary from 35 °C to 45 °C, while in winter, temperatures fluctuate between 15 °C and 25 °C. The amount of rainfall is small, with a gross annual average precipitation established by the Department of Agriculture and Water Research (DAWR) at 82 mm/year between 1990 and 2008. Rainfall events in Qatar are not very intense and generally occur in winter. Due to the karstification of the carbonate formations, the aquifers are covered with multiple areas of topographic depressions and sinkholes that actively participate in groundwater recharge. The contribution of precipitation is the main source of natural water in the country, justifying that several studies tried to quantify this effective recharge, whose main estimated values are gathered in

Table 1. Groundwater recharge estimates (* averaged over the whole Qatar surface area).

	Recharge (Million m3/year)	Recharge (mm/Year) *
Eccleston et al. [12]	50	4.32
Harhash and Yousif [13]	42	3.54
Kimrey [14]	27	2.33
Lloyd et al. [10]	10	0.86
Schlumberger Water Service [11]	56	4.84
Baalousha [15]	25	2.16

.

2.4. Available Piezometric Data

Several piezometric head measurements have been carried out over the last 50 years. However, they were mainly measured only once at a given location and not at the same time. Therefore, to feed the task of groundwater modeling with a consistent set of data, we only selected 166 piezometer heads measured in 2017, during the campaign held by the Qatar General Electricity and Water Corporation (Figure 3). Out of the 166 values retained, only 15 points (or 9% of the total) are in the southern part of the country. This lack of data in the southern part is a severe limiting factor for the groundwater model calibration in this region.

3. The Groundwater Model: Concepts and Calibration

3.1. Groundwater Flow Model

Given the lack of piezometric data to characterize the multiple aquifers, the model was simplified to represent the three basins in Qatar (north, south and southeast) via a single-layer system.

The single measurement campaign carried out in 2017, with no duplication of measurements over time at the same location, does not allow for transient model calibration. Therefore, only the transmissivity values of a single-layer model in a two-dimensional approach to flow will be calibrated based on the 2017 measured piezometric heads. It is assumed that the system is under steady-state flow conditions.

The two-dimensional saturated groundwater flow model is based on the combination of Darcy’s law and the mass conservation equation, considering Dupuit–Forchheimer’s assumption of constant hydraulic head over depth:

$$\begin{gathered} \left\{\begin{array}{cc}\nabla ·T\nabla h\left(x\right)=F\left(x\right)& \\ h\left(x\right)=h_D\left(x\right)& x\in \partial {\mathsf{\Omega }}_D\\ T\nabla h\left(x\right)·n=q_N\left(x\right)& \end{array}\right. \\ x\in \partial {\mathsf{\Omega }}_N \end{gathered}$$

where h is the groundwater head (L), and T is the transmissivity tensor (L²T⁻¹). It is given by $T=Ke$, where e is the aquifer saturated thickness (L), and K is the hydraulic conductivity tensor (LT⁻¹). F is the sink-source term (LT⁻¹), which represents recharge and injection/pumping wells.

$\partial {\mathsf{\Omega }}_D$ and $\partial {\mathsf{\Omega }}_N$ are partitions of the domain boundaries that correspond to Dirichlet and Neumann conditions, respectively, and n is the unit vector normal to the boundary, counted positive outward. $h_D\left(x\right)$ is the prescribed head value at the Dirichlet boundaries, $q_N\left(x\right)$ is the prescribed flux at the Neumann boundaries. The mathematical model is solved by a two-dimensional nonconforming finite element method [16].

The boundary conditions are shown in Figure 4.

The piezometric head at the Dirichlet boundary corresponds to the seaside and was set to 0.0 m neglecting the effects of density-driven flow and eventual seawater intrusion. A Neumann boundary was set up at the border with Saudi Arabia [9], with a total inlet flux estimated at 1 Mm³/year [12].

Based on the different estimated values (Table 1), a recharge of 42 million m³/year was distributed over 6 zones according to a spatial distribution described in Baalousha [15] and given in Figure 5.

Pumping wells were inventoried by Schlumberger Water Service between 10 May 2008, and 9 May 2009 [11]. The 4000 identified wells were gathered into various subareas of the domain (Figure 5), their exact locations being unknown. Without any further information on local extraction rates in each well, the total amount of pumped water (which is roughly known) was uniformly distributed over these wells (more precisely over the subareas and proportionally to the number of wells enclosed), which lead to a flow rate of 161 m³/day and per well.

3.2. Groundwater Model Calibration

Model calibration consisted of estimating the transmissivity values. It was performed using an inverse procedure by seeking iteratively a spatial transmissivity distribution minimizing an objective function (OF) as the sum of the quadratic differences between measured and computed heads. Minimization was performed by relying upon the adjoint state equation [16,17] and a multiscale parameterization [18,19]. The multiscale parameterization discretizes here the distribution of transmissivity values over the domain with a mesh independent of the mesh used for the groundwater flow modeling. Each cell of the flow grid is assigned a transmissivity value calculated by interpolation of the seed values at the nodes of the parameter mesh. The latter can be refined (mainly for seeking a better representation of parameter spatial heterogeneity) after each optimization iteration (or set of iterations) until the OF reaches an acceptable value of 10⁻¹ m². This approach is an alternative to pilot-point methods applied to the northern part of Qatar [20].

Because of the aquifer’s heterogeneity and the limited amount of information, the model calibration may not lead to a unique solution. Therefore, the calibration was performed several times with different initial parameter values for each parameter mesh. In our case, the calibration procedure was run 90 times, each calibration being iteratively conducted until both the criterion on the OF value was reached and also no significant changes were found in mean and standard deviation of the transmissivity value distribution.

Figure 6 shows that, out of the 90 solutions obtained, the north and northeast region of Qatar conceal the highest local mean values of transmissivity (over the 90 solutions). Moreover, the standard deviation associated with this mean (Figure 6) in these areas is low, which means that this pattern is substantially common to each solution. Some extreme values have few physical meaning for aquifers with a saturated thickness of approximately 20 m. The values obtained during the calibration indicate that the northern and north-eastern regions may show average transmissivities 10 to 100 times higher than the usual range of 100 to 1000 m²/day reported by the literature [5,9,21].

In the various approaches used to identify hydraulic conductivity at the Qatar scale, the northern area is often considered as the most conductive [10,20,21]. These high transmissivities can also be explained by the presence of karst formations mainly located in the north of Qatar [22].

The uncertainty on transmissivity evaluated by the standard deviation of the 90 inverse solutions is the highest in the southern part of the country, probably because it is very poorly documented, since only 15 of the 166 piezometers used for calibration are available in the area.

The spatial distribution of mean piezometric heads averaged over the 90 solutions is presented in Figure 7 and its standard deviation in Figure 7.

The simulated piezometric heads appear reliable and stable over the 90 inverse solutions, as shown by the low values of the standard deviation. This is due to computed heads close to measured values and weakly sensitive to the various calibrated transmissivity. Note also that the computed head values are to some extent conditioned by a long perimeter of Dirichlet boundary conditions prescribing the heads along the coasts.

The highest uncertainties can be associated with the recharge areas (Figure 5), where the piezometric heads are more sensitive to the transmissivity values.

Intensive withdrawals that are significantly higher than the natural recharge (

Table 2. Water balance for the Qatar aquifer based on the 90 calibrated models.

Water Fluxes	Inflow (Million m3/Year)	Outflow (Million m3/Year)
Recharge	42.0	-
Neumann boundary	1.00	-
Pumping wells	-	235
Dirichlet boundary	203 ± 4	11.0 ± 4
Total	246	246

) also lead to hydraulic heads below sea level and, therefore, possible saltwater intrusion along two broad coastal strips (east and west) in the northern part of the country. This is consistent with the work of Al-Hajari [5], which reported saline intrusion in the western part of the country.

It can also be mentioned that the 90 different hydrogeological models only slightly differ inon their water fluxes at the Dirichlet boundaries. This feature agrees with the small differences in head values close to this boundary (see Figure 7) due to the prescribed head value at the boundary and also the small differences in transmissivity values between the 90 simulations (see also Figure 6) along the coastal boundary between the 90 simulations.

3.3. Groundwater Storage Strategy

On the basis of the United NationsNation recommendations [23], 3 levels of water consumption werehave been retained: 70, 100 and 200 L/capita/day. The level of 70 L/capita/day corresponds to the minimum value to ensure survival over a medium-term crisis period, i.e., the stress on water distribution envisioned for this study. To compensate for eventual losses during water distribution, the level of 100 L/capita/day was also used to define the storage strategy. Finally, a much higher value of 200 L/capita/day was also studied. This value, well beyond the minimum recommendation, would account for the doubling of the population at the time of crisis, or simply because Qatar, as a large consumer of water, cannot enforce a brutal and severe reduction of its consumption with a decrease of more than a factor of 5 on the ongoing values.

The artificial recharge was defined based on two objectives: (i) estimate an infiltration rate to store the requested amount of water within less than 5 years, and (ii) estimate the infiltration rate to maintain the groundwater water level when the required storage is reached without additional withdraw associated with the redistribution of stored water.

Many different strategies were analyzed using a steady-state calibrated model representative (with a transmissivity distribution close to the average solution) of the 90 calibrated models. The most appropriate strategy relies upon 60 infiltration wells grouped in 3 zones (20 wells/zone) located in the central northern part of the aquifer (Figure 8).

The simulations were performed using the following parameters and assumptions:

-

The amount of water pumped out of the aquifer is the value used for model calibration (235 million m³/year), the same hypothesis applies for groundwater recharge set at 42 million m³/year.

-

The population is estimated at 2.7 M inhabitants.

-

The water storage must meet 60 days of daily demand. Considering the daily water consumption of 70, 100 and 200 L/capita/day, the total amount of water to store is 11, 16 and 32 million m³, respectively.

-

The porosity (storage capacity) of the (unconfined) aquifer is assumed to be uniform.

-

The simulations were performed in a transient flow regime to estimate the time needed to store the required amount of water over 5 years, considered a reasonable time to reach the objectives. The initial conditions correspond to the piezometric situation observed in 2017.

Three different infiltration rates were studied to estimate the time needed to reach the required volume of injected water: 10 million m³/year (456.6 m³/day/well), 20 million m³/year (913.2 m³/day/well) and 30 Mm³/year (1369.9 m³/day/well), and three values of porosity (2%, 5% and 8%) were used, following field data from Dietrich et al. [24].

The water is mainly stored at the central part of the aquifer, as it can be observed by comparing the piezometric levels estimated for 2017 (Figure 7) and the new piezometric levels obtained when the requested storage is reached (Figure 8).

Figure 9 shows an almost linear evolution of the amount of stored water over time, the slope being proportional to the total infiltration rate for a given porosity. In the range of the used values, the porosity does not significantly affect the evolution of the amount of stored water over time.

The times required to reach the 3 objectives within 5 years are summarized in

Table 3. Time (year) needed to reach the required amount of water storage (* objective not reached within 5 years).

Infiltration Rate (Million m3)	Objective (L/Capita/Day)	Porosity: 0.02	Porosity: 0.05	Porosity: 0.08
10	70	2.2	1.8	1.7
	100	3.7	2.7	2.5
	200 *	-	-	-
20	70	0.9	0.8	0.8
	100	1.4	1.2	1.1
	200	3.7	2.7	2.5
30	70	0.6	0.5	0.5
	100	0.9	0.8	0.7
	200	2.0	1.6	1.5

. Most of the objectives are reached in less than 5 years, except for the highest storage value (32 million m³) and the lowest infiltration rate (10 million m³/year).

The objective of storing 11 million m³ (70 L/capita/day) in the underground can be easily reached, in less than one year for a recharge plan of 20 or 30 million m³/year. The objective of 100 Ll/capita/day for 60 days of supply (16 million m³) can be reached, according to the model, within one year if the recharge plans are of 20 or 30 Mm³/year. The objective of 200 L/capita/day, which corresponds to slightly more than a third of the current daily consumption of the inhabitants, represents a water mass of 32 million m³ for 60 days of supply. This ambitious objective is only reached after a few years for recharge plans greater than or equal to 20 million m³/year.

It is worth noting that water storage in the aquifer also increases the amount of water leaving the aquifer by its coastal boundaries. This water loss depends on both the total infiltration rate and the porosity (

Table 4. Total amount of infiltrated water (Inf.) and of outflow at the coasts in Mm³.

Infiltration Rate (Million m3)	Objective (Million m3)	Porosity: 0.02 (Inf.–Outflow)	Porosity: 0.05	Porosity: 0.08
10	11	22.0–11.0	18.0–7.0	17.0–6.0
	16	37.0–21.0	27.0–11.0	25.0–9.0
	32	-	-	-
20	11	19.0–8.0	16.0–5.0	15.0–4.0
	16	29.0–13.0	24.0–8.0	23.0–7.0
	32	74.0–42.0	53.0–21.0	49.0–17.0
30	11	17.0–6.0	15.0–4.0	14.0–3.0
	16	26.0–10.0	23.0–7.0	21.0–5.0
	32	60.0–28.0	49.0–17.0	45.0–13.0

). The total amount of water going back to the sea is larger for small porosity values (up to more than 50% of the total artificial infiltration). For the same amount of stored water, a small porosity value increases the mean water level in the central part of the aquifer and, therefore, the head gradients between this central part and the seashore, which explains the higher losses. The losses also increase with the infiltration rate.

Once the system has reached the target of 16 million m³ of stored water, the infiltration rate must be adjusted to maintain the amount of stored water. To estimate this conservation rate, new simulations were performed, and it was concluded that an infiltration rate of 5.1 Mm³/year (respectively 1.6 million m³/year) was necessary to maintain the groundwater levels for a porosity of 2% (respectively 8%). For both anthropic injection rates, the shape of the piezometric maps are similar (Figure 10) the stored water quantity being the same.

4. Conclusions

This study was mainly motivated by increasing the stocks of freshwater in Qatar by injecting desalinated water into the regional aquifer. This study is highly prospective and needs various forecasts that cannot be discussed without a model of the regional aquifer. This model was built and conditioned onto available data before proceeding with various scenarios of forced reinjections in the aquifer.

Unfortunately, the lack of data compromised building a sophisticated model. For example, the absence of history (duplication of measurements over time) of groundwater levels at various locations over the Qatar territory became a shadow area of the study. It only allowed us to set up a steady-state flow model within a two-dimensional single aquifer taken as a single-porosity medium. This model was inverted for its transmissivity values based on available head levels measured once at various wellbore locations.

The model showed that the central and northern areas of the Qatar territory are assigned very high transmissivity values that also go with the fact that the regional limestone aquifer is highly karstified in these regions. It could have been well advised to locally rely upon more sophisticated settings of the hydrogeological model; for example, by treating the highly conductive areas as a dual-continuum porous matrix + fractures or by posing some explicit draining networks over the single-porosity system. Notwithstanding, these more sophisticated approaches are only relevant for simulating transient flow conditions, as the physics of flow are mainly conditioned by the exchanges of fluxes between compartments. If the problem is not treated in the transient mode (with transient data for conditioning), the models come back to a steady-state flow regime in a single-porosity medium with contrasted transmissivity values.

Fortunately, the geological and geometrical settings of the regional aquifer are such that the system is inertial with very flat head gradients. These should naturally render smooth and weak flushes of freshwater from the central part of the aquifer toward the seashore. This makes it so that weak, average or high values of the storage capacity (a parameter not seen by steady-state flow) in the system do not influence the capability of artificially storing water. Small water storage would increase the head gradients from the central area of the aquifer to the seashore. This increase is small enough to simply favor the loss of injected water toward the sea without canceling out the efforts for storing water. It thus takes more time to reach the storage objectives, as there are more leaks. In opposition, a high storage capacity allows the head gradients to remain very flat but improves the yield of total injected water versus leaks, which quickens the process of storing water.

To feed a population of more than 2.5 M inhabitants, at the rate of 100 L per day and per capita over two months, with no other resource, approximately 20 Mm³ of stock is needed.million m³. This storage value is reached in less than two years for forced recharge rates of approximately 20 million m³/year, irrespective of the conjectured storage capacity of the aquifer. That being said, after reaching the requested stock, keeping it as such by simply compensating for the leaks would necessitate a continuous injection of less than 5 million m³/year. This appears a small value compared with the total amounts of desalinated water (600 million m³/year) produced by plants in Qatar.

Finally, the idea of storing desalinated water in aquifers in Qatar appears a feasible option in order to mitigate the case of a crisis stemming from the total failure of a desalinated water supply. The question of how a limestone aquifer subjected to injection of partly demineralized water (which is the case for desalinated water) behaves still holds, especially regarding the reactivity of the host rock. Incidentally, the model proposed in this study appears stable, meaning that it does not render awkward inverse solutions. This is a good indicator of reliability, but it does not mean that the model behavior over transient flow is a good one. Stated differently, and to improve the efficiency of reinjection scenarios, it is best suited to rehandle the basis of the present study by carrying out model calibrations in the transient mode and by accounting for preferential flow paths in a contrasted karstified system. In this projection, the key point is to be able to rely upon available data rendering information on the transient behavior of the actual system.