Impact of Demographic Growth on Seawater Intrusion: Case of the Tripoli Aquifer, Lebanon

Abstract

Water resources in Mediterranean coastal aquifers are subject to overexploitation leading to an increase in seawater intrusion. Based on the United Nations Environment Program, “UNEP” 75% of people in the world will live in coastal cities by 2020. This is having a major impact on the salinization process. This paper deals with the impact of demographic evolution on seawater intrusion and considers the case of the lower Tripoli aquifer in Lebanon. A numerical model based on the sharp interface approach is implemented using “Freefem++” to access the seawater intrusion. The model is verified against an analytic and a numerical solution. It is calibrated and validated against hydraulic head observations (RMSD = 0.819 m). Then several scenarios of pumping are applied based on available demographic growth rates to quantify the impact on seawater intrusion. The projection scenarios show that if the current pumping rates are maintained while maintaining the demographic evolution, the pumping wells will be salinized within 2 decades in the highly populated areas.

1. Introduction

Groundwater resources in coastal aquifers around the Mediterranean sea are over-exploited. Irrigation, domestic and economic activities are leading to groundwater degradation [1]. One of the most important impacts of over-pumping in a coastal aquifer is the intensification of seawater intrusion. Several studies have already shown the intensification of seawater intrusion across the world, in North America [2], Mexico [3], Australia [4], and Oman [5]. The intensification of seawater intrusion can be attributed to changes in supply and demand. The decrease in supply of freshwater into the aquifer can be due to reduction in rainfall or variability in the snow fall. The increased demand on groundwater resources is mainly due to increase in economic activities and demographic evolution. This paper deals with the impact of demography on the seawater intrusion. We consider the case of the aquifer of Tripoli-Lebanon. Since the 1990s, the uncontrolled urbanization and intensification of internal migration to the coast increased the pressure on freshwater use in this aquifer. The water quality in terms of mineral concentration in the aquifer has been greatly degraded. This decrease has been measured in 2008 and 2009 [6]. With the predicted impact of climate change on the south-Mediterranean region [7] and the aggravation of conflicts in the region, the situation is expected to worsen. Future evolution of seawater intrusion based on scenarios is needed to elaborate mitigation strategies. We choose numerical modeling to achieve this objective.

The physical approach to model seawater intrusion is to use variable density flow in porous media [8,9]. The approach relies on the fully coupled implicit equations of variable density flow in porous media (generalized Darcy 3D flow in saturated media [10]) and the solute transport in porous media (advection-dispersion in porous media). Even though this approach is well established [9] and used in several modeling tools ([11] (FEFLOW), [12] (SUTRA), [13] (OpenGeoSys)), it presents some drawbacks. The main drawback is that it is more requiring in terms of computational needs and numerical convergence [14]. This approach is inherently 3D (or 2D vertical profile). Finally it requires several parameters for the solute transport equation that are difficult to acquire like the transverse and longitudinal dispersion. A second approach to modeling saltwater intrusion is the sharp interface approach. In this case the saltwater and freshwater are considered as two immiscible fluids separated by an interface. The main approximation is provided from the Ghyben-Herzberg [15,16] works that were latter generalized in several models [17]. This approximation need to be carefully verified before application. Its validity has been discussed in [18,19]. This approach has been used by [20] in Batinah aquifer of Oman, where explicit expressions for the water table, sharp interface location and stored volume of fresh water are obtained. In this study, the needed conditions for the use of the sharp interface approach are verified over the lower Tripoli aquifer prior to running the model.

In the following sections, the finite element model is presented for the Tripoli aquifer based on previous hydro-geological studies. The model is verified against analytical and numerical solutions in basic configurations. Then it is applied to the Tripoli aquifer. The last part of the paper shows results for several scenarios of pumping and their impact on the freshwater/saltwater interface. Finally, a summary and conclusions are given.

The following

Table 1. List of used parameters.

Parameter	Symbol	Unit	Dimension
Ground level	a	m	-
Hydraulic head	h	m	-
Freshwater/saltwater interface position	g	m	-
Substratum level	b	m	0
Specific yield of the freshwater zone	${\varphi }_f$	$\frac{m^3}{m^3}$	$0.15$
Specific yield of the saltwater zone	${\varphi }_s$	$\frac{m^3}{m^3}$	$0.15$
Freshwater source term	S	$\frac{m^3}{day}$	33,600
Seawater source term	$S_s$	$\frac{m^3}{day}$	0
Sea level	$h_0$	m	200
Hydraulic permeability	$\mathcal{K}$	m·s${}^{-1}$	$1.388\times {10}^{-5}$
Well radius	r	m	$10;40$
Specific flux	ψ	m${}^3$·m·${}^{-2}$·s${}^{-1}$	-
Area of pump zone	A	m${}^2$	-
Coefficient for pumping rate	λ	m${}^3$·m${}^{-2}$·day${}^{-1}$	-
Root mean square deviation	RMSD	m	-
Aquifer length	$L_x$	km	2
Saltwater wedge penetration length	$L_0$	km	$1.1$
Density of freshwater	${\rho }_f$	kg·m${}^{-3}$	1000
Density of seawater	${\rho }_s$	kg·m${}^{-3}$	1025

shows the list of parameters used in this paper.

2. Materials and Methods

2.1. The Tripoli Aquifer

Tripoli is the second largest city in Lebanon and the regional center of North-Lebanon. It has an area of about 45 km${}^2$, and it is divided into lower (18 km${}^2$) and upper zone (27 km${}^2$) (Figure 1). The aquifer limits are inside the “Abou Ali” and “Abou Hala" watershed limits. In this study, only the lower Tripoli aquifer is considered. Its climate is hot and dry in summer, cool and wet in winter with annual rainfall of 700 mm. The median high and low temperatures vary between 1 ${}^{\circ }$C and 35 ${}^{\circ }$C. The groundwater is mainly recharged from rainfall and snow melt from “Bcharré ” region situated to the east of Tripoli aquifer (See Figure 2a). The rainfall season is short and dry season is long with little or no rainfall. Thus, during a large period of the year, the city relies on water coming from highly inconsistent snow melt [21]. The inhabitants rely highly on pumping the freshwater aquifer all year long due to the inefficient freshwater distribution network and the water shortages during the summer. In the current situation the water use is mainly domestic and industrial, a small amount is used for agricultural activities inside the lower aquifer.

Based on [22,23] the substratum of the aquifer in the study zone is the impermeable layer “C6” formed by marl and marl limestone. This layer is at a constant level of about 200 m below the sea level, see Figure 2b. At the Neogene and Quaternary period, many permeable layers were formed above to the “C6” layer. These layers are mainly formed by marly sand and have a thickness of about 200 m [24].

2.2. Observation Data Set

In the last few decades, Tripoli city has faced a large demographic change, which increased the freshwater demand. The official population growth rate given by Lebanese Ministry of Social Affairs, is $1.6\%$ per year [25]. There are about 82,000 inhabitants in the private wells region. Most of the areas originally covered by Lemon trees plantation were transformed into urban residential areas. These new buildings form a new neighborhood in Tripoli called “Dam and Farz”. Figure 1 shows part of the private wells that are operated in this region. In addition, there are six active public wells managed by Tripoli Water Authority in the study area. Four of them are in Bhsas neighborhood and the two other are in “El-Jisr” and “Malouleh”. Figure 3a shows the private wells region as a rectangle with 3,000,000 m${}^2$ area. Yellow marks represent the public wells.

Table 2. Pumping rate of public wells.

Well	Pumping Rate (m${}^3$/h)
Bhsas	195
Malouleh	136
El-Jisr	204

indicates the pumping rate of these wells [6]. The site is also characterized by the lack of complete observation data set due to unstable political situation. In the study all available data sources are considered the setup and validation of the model.

In 2008, water meters installed in 38 buildings showed that water consumption through well extraction is 250 L per capita per day [6]. This value exceeds the recommendation of the Lebanese Ministry of Energy and Water, which is 120 L per capita per day. In addition, when a private well starts salinizing, it remains operational for a certain amount of time for domestic uses (cleaning, toilets, etc.). A large amount of the drinkable water is acquired from external sources (bottled water, Tankers).

The total private wells pumping rate is about 20,828 m${}^3$ per day. This excessive well pumping has caused an increase in the seawater intrusion [6]. The electric conductivity has been measured in nine wells in “Dam and Farz” neighborhood in 2008 [6]. These measurements show that this region can be divided into four zones based on water ionic content, see Figure 3b. The electrical conductivity was also measured in the same wells in 2009 and showed a notable increase, from 800 μS/cm to 1000 μS/cm [6]. These limits exceed the limits for drinkable water which are below electrical conductivity of 1000 μS/cm [26].

2.3. Mathematical Formulation

In order to model the seawater intrusion in Tripoli aquifer, we consider the following general configuration:

• The Tripoli aquifer is an unconfined aquifer.
• The flow direction is mainly horizontal.
• The flow in the aquifer must be less than 100 m/day, enabling to consider the Darcy flow equation in porous media without the Ward-Forcheihmer quadratic terms ([27,28]).

Based on the four assumptions above, the 3D Darcy velocity equation can be integrated vertically resulting in the Dupuit-Boussinesq approximation [29]. In addition, we base our model on the Ghyben-Herzberg assumptions which means that :

• The interface between freshwater and saltwater is assumed to be a sharp interface (the two fluids are considered as immiscible). This is applicable when the thickness of the transition zone is small. The interface is considered as representative of the lower percentile (50%) of the advection-dispersion approach.
• The sharp interface approach is applicable considering the pumping rate. Based on [19] the sharp interface is more compatible to the realistic case when the pumping rate is higher than 0.07 m${}^3$ per day and up to 0.15 m${}^3$ per day also when the wells are relatively deep. This is the case in the lower Tripoli aquifer where the wells are deep (about 100 m) and the pumping rate is about 33,600 m${}^3$ per day justifying the application of the sharp interface approximation.
• The fluxes in the saltwater zone are considered as negligible and the saltwater is considered in hydro-static conditions.
• An average sea water level is considered. Thus tides effects are not considered.
• The pumping occurs in the freshwater zone. Saltwater adapts to the pressure in the freshwater, and pumping is stopped when the saltwater wedge is reached.

Figure 4 shows a schematic representation of an unconfined aquifer, in which $a\left(x\right)$ denotes the ground level, $b\left(x\right)$ the impermeable bedrock supposed constant in Tripoli case (see Figure 2b, $h(x,t)$ the hydraulic head, $g(x,t)$ the freshwater/saltwater interface, and the seawater height $h_0$. Hence we have the following inequalities:

$$[b\left(x\right)\le g(x,t)\le h_0\le h(x,t)\le a\left(x\right)\mathrm{for}\mathrm{all}(x,t)\in \Omega \times [0,+\infty )]$$

(1)

The aquifer Ω is divided into two parts, $\Omega ={\Omega }_s\left(t\right)\cup {\Omega }_f\left(t\right)$, where ${\Omega }_s\left(t\right)$ is the part of the aquifer which contains saltwater and ${\Omega }_f\left(t\right)$ that contains freshwater.

The general form of equations describing the flux of freshwater and saltwater in coastal aquifers under sharp interface approximation derive from combining the Darcy’s law with the mass conservation law for each fluid [17].

$$\left\{\begin{array}{c}{\varphi }_f\frac{\partial (h-g)}{\partial t}+S=\mathrm{div}\left(-(h-g)\mathcal{K}\nabla h\right)\hfill \\ \hfill \\ {\varphi }_s\frac{\partial (g-b)}{\partial t}+S_s=\mathrm{div}\left(-(g-b)\mathcal{K}\nabla (1-{ϵ}_0)h+{ϵ}_{0}g\right)\hfill \end{array}\right.$$

(2)

where ${\varphi }_f$, ${\varphi }_s$ are the specific yield for the freshwater and saltwater regions respectively, ${\displaystyle {ϵ}_0=\frac{{\rho }_s-{\rho }_f}{{\rho }_s}}$ is the ratio of mass densities, ${\rho }_f$ is the density of freshwater, ${\rho }_s$ that of saltwater, S is the freshwater source term, $S_s$ is the seawater source term and $\mathcal{K}$ is the hydraulic conductivity of porous media.

We assume a steady state due to long-term constant boundary conditions and no source terms in the saltwater. Then the Equation (2) can be written as follows:

$$\left\{\begin{array}{c}S=\mathrm{div}\left(-(h-g)\mathcal{K}\nabla h\right)\hfill \\ \hfill \\ 0=-\mathrm{div}\left((g-b)\mathcal{K}\nabla (1-{ϵ}_0)h+{ϵ}_{0}g\right)\hfill \end{array}\right.$$

(3)

The model incorporates the pumping zones that are 500 m far from the coast. Since the saltwater is considered to be hydro-static, then the “Ghyben-Herzberg” relationship ([15,16]) can then be introduced:

$$(1-{ϵ}_0)h+{ϵ}_{0}g=h_0$$

(4)

Therefore, the second equation of Equation (3) becomes trivial. The System (3) becomes:

$$S=\mathrm{div}\left(-(h-g)\mathcal{K}\nabla h\right)$$

(5)

Based on “Ghyben-Herzberg” relationship, g is given by

$$g=max\left(b,{\displaystyle \frac{h_0-(1-{ϵ}_0)h}{{ϵ}_0}}\right).$$

(6)

Accordingly, the Equation (5) can be written as follows:

$$S=\mathrm{div}\left(-\mathcal{K}min\left((h-b),\left({\displaystyle \frac{h-h_0}{{ϵ}_0}}\right)\right)\nabla h\right)\mathrm{in}\{h<a\}$$

(7)

Assuming that the substratum b is a horizontal layer, then the only variable in the Equation (7) is h. By changing the variable h to u, the Equation (7) can be written in the following form:

$$S=\mathrm{div}\left(-\mathcal{K}u\nabla u\right)\mathrm{in}\{h<a\}$$

(8)

with

$$u=\min \left((h-b),\left({\displaystyle \frac{h-h_0}{{ϵ}_0}}\right)\right)$$

(9)

Equation (8) is completed by boundary conditions in order to have a unique solution. In order to solve numerically Equation (8), we linearize it by the forward schema as follows:

$$S=\mathrm{div}\left(-\mathcal{K}{u}^n\nabla u^{n+1}\right)$$

(10)

where n is an iteration index.

$S=\sum S_i$, for each pumping zone, $S_i$, we consider ${\lambda }_i$, the extraction coefficient and $A_i$ its area. These two are related by the Equation (11)

$$S_i={\lambda }_i\times A_i$$

(11)

Equation (10) is solved numerically using finite element method in “Freefem++” software [31]. “Freefem++” is a free software to solve 2D and 3D non linear equations based on finite element method using C++ as programming language [32]. Computing u numerically is straightforward, using the “Freefem++” software. Then, by using Equation (9) and Equation (4), we can detect the freshwater/saltwater interface g. We call this implementation in “Freefem++”: “FreeSWIM”.

2.4. Model Verification in Schematic Domain

Prior to applying the model to the Tripoli aquifer, we conducted a validation exercise using a schematic rectangular domain. This exercise was performed for two configurations: with and without pumping wells. The model is verified against an analytical solution (only when no pumping is used) [17] and a numerical model : “BFSWIM” [33]. “BFSWIM” approach is a specific implementation of the “BIGFLOW” model [34] and was validated to the “Feflow” model using the Henry problem setup [35,36]. The domain geometry using “BFSWIM” is bound to be rectangular, so it is not adapted to model the lower Tripoli aquifer. Benchmarks experiment like the Elder problem [37] and the Goswami-Clement problem [38], were not used as their setup and configuration requires a variable density flow approach. The Elder problem may need a thermo-hydro variable density flow in porous media. The Goswami-Clement is a rectangular experimental setup of (53 cm × 26 cm × 27 cm) dimension thus the depth to thickness ratio is not commonly observed for coastal aquifers.

We launch the simulations in an homogeneous rectangle of 6 km × 2 km, and a pumping well in the middle of the domain (Figure 5). The substratum is considered as the reference level b = 0 m, $h_0$ = 200 m. The porous media has a hydraulic conductivity of $\mathcal{K}$ = 1.388 × 10${}^{-5}$ m·s${}^{-1}$. Concerning the boundary conditions, we consider Dirichlet condition, we have $h=h_0$ on ${\mathsf{\Gamma }}_1$ (the shore), on ${\mathsf{\Gamma }}_2$ we consider that the hydraulic head is 7 m above sea level.

Using Equation (9) and Equation (10) we have the following system:

$$\left\{\begin{array}{cc}\mathrm{div}\left(-\mathcal{K}{u}^n\nabla u^{n+1}\right)=\lambda {\mathbf{1}}_{B(S,r)}\hfill & \mathrm{in}\Omega \hfill \\ \\ u=0\hfill & \mathrm{on}{\mathsf{\Gamma }}_1\hfill \\ \\ u=207\hfill & \mathrm{on}{\mathsf{\Gamma }}_2.\hfill \end{array}\right.$$

(12)

where λ represents the extraction coefficient of the well placed at the position S, and $B(S,r)$ is a disk, used to model the pumping well, centered at S and of radius r. ${\mathbf{1}}_{B(S,r)}$ is the characteristic function, the pumping rate is equal to the extraction coefficient multiplied by the disk area. The validation has been done for a pump rate S = 2.5 m${}^3$·s${}^{-1}$, λ = 8 × 10${}^{-4}$ m·s${}^{-1}$ and r = 10 m.

In case of no pumping, $\lambda =0$, the freshwater/saltwater interface is also calculated and we considered the analytic solution of the sharp interface freshwater saltwater, which in this case is given by [17].

$$g=\left\{\begin{array}{cc}{h}_0-\sqrt{-2\times Q\times x/{ϵ}_0/\mathcal{K}/({ϵ}_0+1)}\hfill & \mathrm{if}x\le L_0\hfill \\ & \\ b\hfill & \mathrm{if}x\ge L_0.\hfill \end{array}\right.$$

(13)

where $Q=-2\times \mathcal{K}(h_1^2-({ϵ}_0+1)\times h_0^2)/L_x$, $L_0=-2\times \mathcal{K}\times Q\times {ϵ}_0\times ({ϵ}_0+1)\times h_0^2$, $L_x$ is the aquifer width, and $L_0$ is the point where the interface g reaches the substratum b. We detect also the interface g.

The interface g, calculated by the three methods has been compared by using the cross section in the mid of the domain. The “FreeSWIM” interface g, is very close to the analytic and “BFSWIM” interface as shown in Figure 6 which represents the freshwater/saltwater depth to the width of the schematic domain. The root mean square error for the exact analytic solution and for the numerical solution using “BFSWIM” code was calculated and shown in

Table 3. The root mean square error of “FreeSWIM” model in each case in meters.

Interface	Without Pumping		With Pumpinp
	RMSE (Analytic)	RMSE (BFSWIM)	RMSE (BFSWIM)
h	$0.028$	$0.015$	$0.055$
g	$1.13$	$0.601$	$2.19$

).

The difference between errors in pumping case is higher than the case of no pumping case, this is because of the assumptions inherent to the model. Indeed, in the model we assume that the flow is quasi-horizontal, but near the active wells, there is an impact of up-coning, then the flow is not horizontal.

2.5. Model Setup and Calibration

In order to build a model describing the current situation of seawater intrusion in Tripoli aquifer we should solve Equation (8) over a fixed finite element mesh with the suitable boundary conditions and with the actual pumping rate.

By crossing geological information indicated in Section 2.1 and Section 2.2, the numerical terrain model, and the remote sensing optical image of the Tripoli city (See Figure 3a), the boundary of the aquifer has been identified and digitized. The domain starts from the sea side on the shore and includes newly gained surfaces from the sea. The observed outcrops in the optical image and during field visits indicate the fracture between the lower and the upper zone of Tripoli city, see Figure 2. This fracture defines the hydro-geological limits. The border is digitized using Geographic Information Systems tools that combines all the available sources (topography, geology, remote sensing images). From the existing boundary, a triangular finite element mesh is generated. The mesh grid density is increased at the location of pumping wells and at the borders of the domain. This results in a well-defined mesh of 99,412 nodes and 196,822 elements, see Figure 7.

Recharge over the urban area is considered as negligible in the majority of the aquifer as a sewage network exists and the bare soil and vegetated areas are minimal. Also the rainfall during the dry season is very low. Boundary conditions in the study zone are divided into two parts: ${\mathsf{\Gamma }}_1$ which represents the shore, and ${\mathsf{\Gamma }}_2$, which represents the eastern inflow/recharge condition. On ${\mathsf{\Gamma }}_1$, we have a Dirichlet boundary condition, $h=h_0$ representing the height of the sea. By using Equation (9) we get $u=0$. ${\mathsf{\Gamma }}_2$ is the flux of freshwater entering the domain from the eastern border. It is a Neumann boundary condition. We have then $\mathcal{K}{\displaystyle \frac{\partial u}{\partial n}}=\psi$ on ${\mathsf{\Gamma }}_2$, where ψ is the specific freshwater entering into domain from ${\mathsf{\Gamma }}_2$ side and n is normal vector on ${\mathsf{\Gamma }}_2$ border.

Consequently, Equation (8) can be written as follows:

$$\left\{\begin{array}{cc}\mathrm{div}\left(-\mathcal{K}{u}^n\nabla u^{n+1}\right)=\sum _{i=1}^n{\lambda }_i{\mathbf{1}}_{A_i}\hfill & \mathrm{in}\Omega \hfill \\ \\ u=0\hfill & \mathrm{on}{\mathsf{\Gamma }}_1\hfill \\ \\ \mathcal{K}{\displaystyle \frac{\partial u}{\partial n}}=\psi \hfill & \mathrm{on}{\mathsf{\Gamma }}_2.\hfill \end{array}\right.$$

(14)

Table 4. Configuration parameters and pumping rate for the operated wells.

Pump Zone	${\lambda }_i$ (m·Day${}^{-1}$)	$A_i$  (m${}^2$)	$S_i$ (m${}^3$·Day${}^{-1}$)
Private wells	$7\times {10}^{-3}$	3,000,000	21,000
Bhsas	$0.9315$	5024	4680
Malouleh	$10.394$	314	3264
El-Jisr	$15.592$	314	4896

gives the ${\lambda }_i$, $A_i$ and $S_i$ values for the private wells, Bhsas, Malouleh and El-Jisr.

Each public well was modeled as a disk with a radius of $r=40$ m for Bhsas wells and $r=10$ m for the two other wells. For the region where private wells are densely distributed, we consider a mean pumping rate in the whole region for several reasons. First, the well positions and rates are not well identified, as they are not official. Second, the population and the density is well known from official sources. Thus, it is easier to infer the pumping rate for the whole area.

Based on [22,23], the hydraulic conductivity $\mathcal{K}$ has been estimated by a constant value of 1.388 × 10${}^{-5}$ m·s${}^{-1}$. The impermeable layer is at $b=0$ m and the sea level is at $h_0=200$ m.

At the border ${\mathsf{\Gamma }}_2$, there is a flux of freshwater entering into the aquifer. The calibration of this flux, ψ, is based on a subset of hydraulic head measurements in several wells performed in 2008. At first, we used the measured hydraulic head in the wells near the ${\mathsf{\Gamma }}_2$ boundary [39] to launch the model with conditions of Dirichlet on ${\mathsf{\Gamma }}_2$, then the associated flux, ψ, is computed and forced into the model. The flux ψ is then optimized based on the minimization of the quadratic difference between the measured and modeled hydraulic head inside the domain using the “SIMPLEX” algorithm. $\mathcal{K}$ is not modified in the calibration. After calibration we obtain an RMSD of 0.891 m.

2.6. Scenario Build-up

The demographic evolution is projected up to 2028 based on a 1.6% per year increment rate as shown in

Table 5. Projected number of inhabitants in Tripoli aquifer area.

Year	Number of Inhabitants
2008	82,000
2018	95,000
2028	110,000

. This evolution does not take into consideration the extreme events like conflicts in neighboring countries which increased the number of inhabitants in recent years. We have used the projections in Table 5 and modified the pumping rate in the unauthorized and official wells. The nominal rate of 250 L and the recommended rate of 120 L per capita were considered.

Four scenarios are considered :

• ”NOP”: refers to the no-pumping scenario.
• ”REC”: refers to scenario with a 120 L per capita as recommended by the water authorities.
• ”10 Y”: refers to a 10 years projection with the current pumping rate of 250 L per capita.
• ”20 Y”: refers to a 20 years projection with the current pumping rate of 250 L per capita.

The ”NOP” no-pumping simulation is the natural equilibrium conditions for seawater intrusion in the aquifer. None of the scenarios include special mitigation strategies like extraction of water further into the aquifer, injection of treated water into the aquifer, crystallization of the soil, or expansion of the land to the sea. The yearly projection is based on the number of inhabitants and doesn’t take into consideration the change in climate conditions and the change in strategies in water supply. The scenario results are compared to the results of the conditions for year 2008, for which the model was calibrated based on observations. This condition is labeled “CUR” meaning current condition.

3. Results and Discussion

3.1. Status of Saltwater Intrusion in Tripoli

The “FreeSWIM” model is run based on the calibration of the boundary conditions, hydraulic conductivity and wells position for year 2008. The flux on the land boundary conditions is calibrated using the hydraulic head at border ${\mathsf{\Gamma }}_2$ of the domain measured in year 2008. The mean hydraulic head is about 5.77 m meters above the sea level. After calibration, using a root mean square minimization in combination to a simplex algorithm the obtained flux ψ is estimated at 26.0235 m${}^3$·m${}^{-2}$·s${}^{-1}$. “FreeSWIM” model computes u numerically. Then, by using Equation (9) and Equation (4), we calculate the hydraulic head h and the interface g in Tripoli aquifer, as shown in Figure 8. These results are labeled as current saltwater intrusion conditions “CUR”.

The “CUR” simulation results are compared to the measured hydraulic head in observation wells in the aquifer. The comparison is done using the wells inside the domain and don’t include the wells used for calibration that are at the outskirts of the domain. The validation data set is constituted of 15 wells localized in the main area of interest of “Dam and Farz”. Figure 9a shows the scatter plot of the observed head vs. the modeled head. We notice that there is little variability in the head values. This can be explained by the high depth of the aquifer. The hydraulic head gradient is consistent with the expected value for a coastal aquifer with sandy marl soil [8]. The figures also show that the modeled head fits well with the measured one with a root mean square deviation of $RMSD$ = 0.819 m. Figure 9b shows the localization of the observation wells. The color represents the difference between mean observed head and modeled head. The position of the observation wells in the area of interest are well shown. Figure 9 shows that the model gives satisfactory results in terms of hydraulic head and that little place for improvement can be achieved through a calibration using existing hydraulic heads. The difference map does not show a high spatial correlation between bias values. No specific pattern can be depicted.

3.2. Impact of Demographic Growth on Saltwater Intrusion

The suggested scenarios in Section 2.6 are run using the same boundary conditions of “CUR” simulation. Figure 10 shows the position of the saltwater/freshwater interface for the different scenarios. The difference in terms of interface position between the current condition “CUR” and the other scenarios is shown in Figure 11.

The “NOP” scenario Figure 11a which considers no-pumping shows a significant reduction of the saltwater advancement into the interior of the aquifer. As expected the enhancement concerns primary the pumping zones. A reduction of 30 m of the position of the interface can be expected. We notice that the recommendation of 120 L per capita if respected gives half the improvement than the natural state “NOP”. On the other hand, the projections for 10 years and 20 years show that the level of the interface will raise to 14 m in the areas of the pumping. Figure 12 shows the inland advancement of the 100 m contour level of the interface g based on the “10 Y” and “20 Y” pumping rates (250 L/capita/day) with respect to the number of inhabitants and the years. The 100 m contour was selected here because the majority of the private wells are installed at this depth. It shows an advancement of 100 m within the next 2 decades.

The modeling results show if the current pumping rates are maintained that the saltwater intrusion in Tripoli aquifer will be aggravated, extending from the coast to the entire basin in 2 decades. Given that the level of freshwater will decrease by 14 m a large part of the private pumping wells will become unusable. One mitigation to this would be to enhance the water distribution system and to pump in the upper part of the aquifer or the “Naher Abou Ali” watershed sources. Part of the supply pumping wells are close to the upper zone of Tripoli, so increasing their pumping rate is recommended as they are far from the sea. Note that other elements need to be investigated before applying this recommendation like pumping yield, impact of karstic system and local pollution in the upper zone. Also it seems to be a reasonable solution to reduce the impact of saltwater intrusion in short term as the 120 L per capita per day seems to be an achievable objective using public awareness programs, educational programs and restriction policies accompanied by monitoring. This result is consistent with previous studies that showed the impact of pumping rates on saltwater intrusion in Tunisia [40]. The results of the scenarios need to be considered with care for several reasons. The local impact of pumping is not well modeled here and may be important. Also, a major contributor to the scenarios is the climate variability which is expected to play an increasing role in future decades [7]. The climate change has a well known impact on seawater levels confirmed since several decades at a mean rate of 4 to 6 mm/year [7]. Precipitation rates (rainfall and snow) are also expected to decrease . This should be the subject for future studies where the impact of each effect is modeled separately or jointly to determine the most probable scenarios. In this study we focused on demographic evolution as it is the most noticeable and directly implementable impact.

4. Summary and Conclusions

In this paper we present a study on the impact of demographic evolution on saltwater intrusion in the lower Tripoli aquifer in Lebanon. We present results based on numerical modeling using the sharp interface approach. We create a numerical implementation called “FreeSWIM” using the “Freefem++” tool to compute the interface position in steady state. “FreeSWIM” is verified against analytic solution and the “BFSWIM” software. We apply the model to the Tripoli aquifer. The model is calibrated and validated using piezometric measurements in the aquifer. Different pumping scenarios were presented in order to make projections of the seawater intrusion. The results show that the excessive pumping will cause the salinization of the aquifer’s freshwater in the region of private pumping. The results show clearly the expected aggravation of seawater intrusion in the next two decades if pumping rates are maintained. The Freshwater/Saltwater interface will raise to about 14 m in the region of pumping. The current paper thus provides a methodology for the assessment of saltwater intrusion evolution in the context of demographic change in weakly instrumented sites. The current study examines only the impact of demographic change on saltwater intrusion in the lower zone of Tripoli aquifer. Future works should simulate the impact of pump in the upper zone of the aquifer and the different hydroclimatic parameters (snow, rainfall, seawater level) and their variability in combination with demographic evolution. The salinization is expected to increase even more considering the fact that this area is considered as a hot-spot of climate change.