- Email: [email protected]
A 3D model for groundwater flow and seawater intrusion interface: Application to the Martil coastal aquifer system
539
A 3D m o d e l for g r o u n d w a t e r f l o w and s e a w a t e r intrusion interface: A p p l i c a t i o n to the Martil c o a s t a l a q u i f e r s y s t e m A. Aharmouch, A. Larabi and M. Hilali LIMEN, Ecole Mohammadia d'Ingenieurs, B.P. 765, Agdal, Rabat - Morocco
Coastal aquifers involve varying conditions in time and space, owing to the occurrence of free and moving boundaries, leading to upconing due to intensive pumping near the coast. A numerical procedure is developed for solving the location of these interfaces in steady and transient conditions with the Galerkin finite element technique. The conceptual approach developed to solve the finite element matrix equations is presented. Benchmark tests which involve confined and phreatic groundwater flow in coastal aquifers, are studied to demonstrate the capability of the present approach of predicting accurate results, including aquifers involving saltwater interface upconing. The model was also applied to the coastal aquifer of Martil in Morocco to simulate groundwater flow and to study effects of the seawater intusion interface under the ONEP (Office National de l'Eau Potable) pumping well rates. 1. INTRODUCTION The contact between freshwater and seawater in islands and coastal aquifers, requires special techniques of management, because of varying conditions in time and space, such as pumping near the coast. In response to pumping from a well in the fresh-water zone, the freshsaltwater interface moves vertically toward the well causing the upconing phenomenon. The problem of saltwater intrusion can be analysed by two methods. The first method considers both fluids miscible and takes into account the existence of transition zone [1]. The second method is based on an abrupt interface approximation [2] where it is assumed that (i) the freshwater and saltwater do not mix and (ii) are separated by a sharp interface. This method is considered as appropriate approximation in the case where the extent of the transition zone is relatively small compared to the aquifer thickness. The first attempts to solve the saltwater intrusion were based on the assumption of the existence of a sharp interface between freshwater and saltwater [3, 4, 5, 6, 7]. In this paper, a 3-D numerical model is presented which enables the prediction of the the fresh-saltwater interface displacement in response to varying conditions in time and space. It is assumed that an abrupt interface exists between freshwater and saltwater and its location, shape and extent must be determined. The Galerkin finite element approach is used and the numerical scheme is verified throughout benchmark problems and applied to the Martil coastal aquifer in Morocco to simulate regional groundwater flow and the transient interface displacement due to groundwater exploitation.
540 2. GOVERNING EQUATION The general form of the equation describing the flow in a partially saturated porous medium is:
Kxx
+
Kyy
+
Kzz
+ q : S(h)-~
(1)
where Kxx, K~, Kzz are hydraulic conductivity components, q is the sources or sinks per volume aquifer unit, S(h) = f(O, O, Ss) is the general storage term, 0 is the volumetric moisture content, ~bis the porosity, Ss is the specific storage, and h is the hydraulic head. 3. N U M E R I C A L FINITE E L E M E N T SOLUTION The model is based on a fixed domain discretization and includes the saturatedunsaturated groundwater zone and salt water region. So the free surface and the interface are not boundary conditions, but are parts of the problem. The position of these moving interfaces can be obtained by the pressure equilibrium that exists between groundwater and atmosphere or saltwater. The hydraulic head is expressed by:
h=z+
P = z +~ Pzg
(2)
= P is the pressure head, p is the interstitial pressure, pf is the freshwater density and g is
Pzg
the acceleration due to the gravity. The atmospheric pressure is often taken equal zero, so the free surface is the geometric location of points where h -- z. Because of the quasi stationary saltwater zone assumption, the hydrostatic pressure distribution is considered. Hence, in the saltwater region, the pressure is expressed by :
(3)
p : -p, gz
The interface must verify:
h=--~
,
6 = p ' -Pp' ~
(4)
In the present approach a fixed finite hexahedral element mesh is used to discretize the whole domain, including the saturated and unsaturated zones as well as the saltwater region. It is assumed that the groundwater flow occurs only in the freshwater region of the saturated zone and no flow is considered in the unsaturated and saltwater regions (the salt water is at rest). The water table and the interface are effectively hold impervious boundaries, but their location is unknown initially and must be obtained as part of the solution. For time discretization, a )~-weighted finite difference scheme is used ()~ = 0.5, Crank-Nicolson ; k = 1, backward Euler). The Galerkin finite element discretization yields the system of non linear algebraic equations: h ~ + I _
G(hk + ~)h ~+~ + C(h k+~)
z~;k + 1
h ~
=qt~( +~)
(5)
where h k + z = 2h k § 1 + (1 - 2)h k , h is the vector of nodal heads, superscript k denotes time step, G is the conductance matrix, C is the capacitance matrix and q is a vector containing the
541 boundary conditions. The coefficients of the conductance and capacitance matrices are given by: (6)
C~ = E ~SN, N/dVe
(7)
e
N~ and Nj are basis functions respectively for nodes i andj
3.1. Characteristic equations The flow occurs naturally in the saturated freshwater zone, however there is a small but negligible amount of water movement in the unsaturated and salt water zones. The flow in the unsaturated zone can be described [8] by :
o(~,) =
Kr(~/y,)
{0 r + (0 s - 0r)[1 + fll-m o,
_.
{I 1
+ /3)
.:mE< 2
1 if" /8)
m
-- ~ m
~,_<0
(8)
~u>0
j
2
~t'_< 0
(9)
~>0
Or is the residual moisture content, 0s is the saturated moisture content, fl =
~u
, ~s is the
capillary or air entry pressure head value, n can be interpreted as a pore size distribution index (mn = n -1, and (1.25 < n < 6)). A similar approach is also employed to describe flow in the saltwater region, which is transformed to an equivalent unsaturated zone where the pressure at any point of this sub-domain is defined by (~-p~). Hence, the following equations are used : Or + (0 s - Or)[1 +/3] -m o,
o(~,) =
K~(~) =
1 + fl) 2 (1 + fl)m _ tim
P'
r
(10)
P~
(11)
r
ps is the pressure head on the freshwater-saltwater interface, Or is the "residual water" in the equivalent unsaturated zone.
,8=
/ ~"- Ps / n ~s
The proposed relative conductivity in the saltwater zone (11) can be reduced in order to satisfy the hydrostatic assumption. The general storage term is expressed in the unsaturatedfresh saturated region as :
542 n-
(n - 1)(0,
dO
1
- or ~l~l
~<0
I
S =
dp'dO ,-~
with fl =
P'
Igsl"(l+~)m+l
(12) ~_>0
= S~
and in the freshwater-saltwater region by : dO
S=
( n - 1 ) ( 0 s - Or ) l ~ , -
~-
p
,,-1
~'
with fl =
I~i -
Iq/s[~(l§
Ps
n (13)
dO -d-~ = s , 9
~>p,
The storage term in the saltwater zone is also reduced, so that the hydrostatic assumption can still be valid. Note that the concept of approaching the flow in the 3 subregions as described in this section was also studied by [9], but using the so called the fast updating procedure. 3.2. Linearization technique Picard iteration scheme is used to solve (5) and gives : (m)
AG k+~ +
C
Ah
At~+l (m+l)
-
h~+~ -
+
qk+,~
(14)
A t ~+l (m)
where A h = h k +l _ h ~ +1 The resulting linear equations system is solved by the preconditioned conjugate gradients method with an incomplete factorization preconditioner based on M-matrix technique [ 10]. The updating process is repeated iteratively until there is no more significant change in the potential values over the fresh saturated flow domain. 4. VERIFICATION OF THE NUMERICAL MODEL
To verify the numerical model described above, several examples are considered, for which analytical solutions are available. Two examples are selected in this paper, which involve saltwater upconing and steady state flow with homogeneous isotropic hydraulic conductivity and using the Dupuits assumption. 4.1. Upconing in a confined aquifer overlain by an aquitard This problem is a field case [6], with the site-specific parameters shown in Table 1. Four simulations have been performed in steady state [11]. Both the numerical and the analytical position of the interface below the pumping well are plotted in Figure 1. The comparison between both solutions showed a good agreement except in the vicinity of the pumping well where the divergence starts to be relatively important for the last case. This is due to the horizontal flow approximation around the well assumed by [6].
543
Table 1 9Parameters used in example 4.1. Parameter Value -438.8 m Initial elevation of the saltwater-freshwater interface 10.9m Initial freshwater head 9.26x 10 "12 m 2 Horizontal intrinsic permeability k• kyy Leakance of the confining bed k'z/b' 1.2x 10 .3/day (in motz (1992)) Freshwater head in overlying constant head source bed 10.97 m Freshwater head at constant head boundary nodes 10.97 m Bottom elevation of well screen -260.0 m Length of well screen 221.0 m Distance from well bottom to the initial interface 178.8 m
4.2. Upconing beneath a pumping well in a coastal unconfined aquifer This example deals with seawater intrusion in an areal plane of an unconfined aquifer with a fully penetrating well near the coast line in steady state. The well withdraws freshwater only, causing upconing of underlying salt water [4]. The parameters used in this example are shown in Table 2.
Table 2 9Parameters for example 4.2. Parameter Value 1 m2/day (at x = 1000 m) Fresh water flux Qf 400 m3/day Well pumping rate Qw (6oo m, 0) Well location (xw,yw) 70 m/day Hydraulic conductivity, K 20m Depth from mean sea level to aquifer base The domain is composed of 1000 mx4000 m respectively along x and y directions. But, due to the symmetry upon the y axis, only the upper part of the domain was considered [ 11 ]. The simulated and analytical interfaces and water table positions are both plotted in Figure 2. The results show that there is a very good agreement between predicted water table positions and interfaces.
0
--" I
Dot
lines : Motz analytical
interface
Dashdot lines : Present
^
numerical
ution (1976) -s
Q = 112 Kg/s
'
D•
lillllllllll
f[ointhe well (,,,)
-20 Horizontal
Figure 1. Interface location for both analytical and numerical solutions for example 4.1.
diltance
(
m
)
Figure 2. Interface and water table positions for both analytical and numerical results for example 4.2.
544 4.3. Three-dimensional model for the Martil coastal aquifer The Martil aquifer is situated in the North-Western part of Morocco on the Mediterranean coast. It is the important local aquifer for potable water supply, irrigation and industry in the Martil region. In accordance with geology, the aquifer system consists of two aquifer units, separated by an aquitard. The upper unit is a Quaternary alluvial deposits of Martil river, the lower unit is composed of sandstone-limestone Pliocene formation, whereas the aquitard is generally marly and clayey. The aquifer thickness is very irregular in different directions. Figure 3a shows a cross-section from West to East. The conceptual model used in this study is a multi-layer aquifer system of variable thickness. The finite element mesh used for the numerical simulation of the Martil aquifer, shown in Figure 3b, is adjusted to fit the structure and the extension of the hydrogeological layers and the wells location. The 3-D mesh contains 56232 nodes, and 48288 hexahedral elements. The saturated hydraulic conductivity for different aquifers are Ks = 3.5 m/day in the alluvial material (upper), Ks = 4.5 m/day for the sandstone-limestone (lower), and Ks = 0.05 m/day for the marl/clay aquitard [9]. The other parameters are held constant (0s = 0.45, 0r = 0.1). West
Figure 3a. Martil West-East cross sectional view.
East
Figure 3b. Martil 3-D structured mesh.
The boundaries are treated as impervious, except at the eastern part, which is in direct contact with the sea (outflow face). Its extent is not known and constitutes a part of the problem. We also consider Dirichlet conditions at nodes belonging to the Martil and Alila rivers. In order to take into account the dry seasons these values are not set necessarily to be equal to the elevation. For the effective rainfall, we use a non-uniformly distribution taking into account the variably material distribution on the top layer. In a first step steady state simulations are performed under natural conditions such that natural groundwater flow can be reproduced efficiently. The purpose is to adjust conductivities and effective recharge. The obtained results are compared to observed piezometric map of 1966. The calibration is done using the trial and error technique. This procedure is repeated until a satisfactory pattern is found. Computed groundwater potential heads, shown in Figure 4, are compared to recorded values at the observation wells in 1966. The correspondence is relatively good. The 1966 interface location is also illustrated in Figure 5.
545
N
Figure 4. Computed (left) and observed (right) steady state groundwater potentials in the phreatic aquifer for 1966. A second set of simulations is performed starting from the 1966 head distribution, to predict the seawater intrusion in the aquifer. A transient simulation is made of over 30 years, considering variable pumping rates used for water public supply shown in Table 3. For irrigation, an amount of 7949 m3/day is pumped from several wells to irrigate 300 ha of agricultural area. The recharge is maintained constant along the simulation time. Table 3 : ONEP pumping rates Pumping rate (m3/day) Well number 992/2 780 994/2 1351.56 995/2 1181 996/2 418.75 1153/2 1260 1174/2 248.4 The results shown in Figure 6 show that the saltwater interface starts to move until 1982 and then remains stable along the time, although pumping wells withdraw upstream in the aquifer. This simulation scenario shows that there is no significant displacement of the interface due to the ONEP pumping wells upstream of the Martil river. West
East
T \
Figure 5. Plane view of the 1966 interface toe location.
Figure 6. Cross-sectional view of the moving interface positions from 1966 to 2014 (a new steady state regime is reached since 1982).
546 5. CONCLUSION A numerical model based on a 3-D finite elements is developed to predict the location, the shape and the extent of the moving interfaces in coastal aquifers involving unconfined or confined flow. Verification of this numerical model for upconing problems was made using the available analytical solutions to find the interface location and its dome extent. The numerical and analytical results are generally in good agreement. The model was also applied to the coastal aquifer of Martil in Morocco to simulate groundwater flow and to study the impact of the ONEP pumping well rates used for water public supply. The simulation results show that there is no significant displacement of the interface. ACKNOWLEDGEMENT
This work has been supported partially by the PARS SPI060 research project. REFERENCES
1. C. I. Voss and W. R. Souza, AQUIFE-SALT : A finite element model for aquifer containing a seawater interface. Water Resources Investigations Report 84-4369, US Geological Survey (1984) 37. 2. J. Bear and A. Veruijt, Modeling Groundwater Flow and Pollution, Reidel, Dordrecht, Holland, 1987. 3. M. A. Collins and L. W. Gelhar, Seawater intrusion in layered aquifers, Water Resour. Res., 7(4) (1971 ) 971-979. 4. O. D. L. Strack, A single-potential solution for regional interface problems in coastal aquifers, Water Resour. Res., 12(6) (1976) 1165-1174. 5. O. D. L. Strack, Groundwater Mechanics. Prentice Hall, Englewood Cliffs. NJ, 1987. 6. L. H. Motz, Saltwater upconing in an aquifer overlain by a leaky confining bed. J. Ground water 30(2) (1992) 192-198. 7. J. W. Bower, L. H. Motz and Durden D. W, Analytical solution for determining the critical condition of saltwater upconing in a leaky artesian aquifer. J. Hydrol. 221 (1999) 43-54. 8. M. T. Van Genuchten, A closed form Equation for predicting the Hydraulic Conductivity of Unsaturated Soils, Soil Sci. Soc. Am. J., Vol. 44 (1980) 892-898. 9. M. A. Sbai, F. De Smedt and A. Larabi, A Generalized Approch for Modeling 3D Transient Free and Moving Boundaries in Coastal Aquifers, First International Conference on Saltwater Intrusion and Coastal Aquifers-Monitoring, Modeling, and Management. Essaouira, Morocco, April 23-25 (2001). 10.A. Larabi and F., De Smedt, Solving 3-D hexahedral finite elements groundwater models by preconditioned conjugate gradients method, Water Resour. Res., 30(2) (1994) 509-521. 11.A. Aharmouch and A. Larabi, Numerical Modeling of Saltwater Interface Upconing in Coastal Aquifers, First International Conference on Saltwater Intrusion and Coastal Aquifers-Monitoring, Modeling, and Management. Essaouira, Morocco, April 23-25 (2001).
A 3D m o d e l for g r o u n d w a t e r f l o w and s e a w a t e r intrusion interface: A p p l i c a t i o n to the Martil c o a s t a l a q u i f e r s y s t e m A. Aharmouch, A. Larabi and M. Hilali LIMEN, Ecole Mohammadia d'Ingenieurs, B.P. 765, Agdal, Rabat - Morocco
Coastal aquifers involve varying conditions in time and space, owing to the occurrence of free and moving boundaries, leading to upconing due to intensive pumping near the coast. A numerical procedure is developed for solving the location of these interfaces in steady and transient conditions with the Galerkin finite element technique. The conceptual approach developed to solve the finite element matrix equations is presented. Benchmark tests which involve confined and phreatic groundwater flow in coastal aquifers, are studied to demonstrate the capability of the present approach of predicting accurate results, including aquifers involving saltwater interface upconing. The model was also applied to the coastal aquifer of Martil in Morocco to simulate groundwater flow and to study effects of the seawater intusion interface under the ONEP (Office National de l'Eau Potable) pumping well rates. 1. INTRODUCTION The contact between freshwater and seawater in islands and coastal aquifers, requires special techniques of management, because of varying conditions in time and space, such as pumping near the coast. In response to pumping from a well in the fresh-water zone, the freshsaltwater interface moves vertically toward the well causing the upconing phenomenon. The problem of saltwater intrusion can be analysed by two methods. The first method considers both fluids miscible and takes into account the existence of transition zone [1]. The second method is based on an abrupt interface approximation [2] where it is assumed that (i) the freshwater and saltwater do not mix and (ii) are separated by a sharp interface. This method is considered as appropriate approximation in the case where the extent of the transition zone is relatively small compared to the aquifer thickness. The first attempts to solve the saltwater intrusion were based on the assumption of the existence of a sharp interface between freshwater and saltwater [3, 4, 5, 6, 7]. In this paper, a 3-D numerical model is presented which enables the prediction of the the fresh-saltwater interface displacement in response to varying conditions in time and space. It is assumed that an abrupt interface exists between freshwater and saltwater and its location, shape and extent must be determined. The Galerkin finite element approach is used and the numerical scheme is verified throughout benchmark problems and applied to the Martil coastal aquifer in Morocco to simulate regional groundwater flow and the transient interface displacement due to groundwater exploitation.
540 2. GOVERNING EQUATION The general form of the equation describing the flow in a partially saturated porous medium is:
Kxx
+
Kyy
+
Kzz
+ q : S(h)-~
(1)
where Kxx, K~, Kzz are hydraulic conductivity components, q is the sources or sinks per volume aquifer unit, S(h) = f(O, O, Ss) is the general storage term, 0 is the volumetric moisture content, ~bis the porosity, Ss is the specific storage, and h is the hydraulic head. 3. N U M E R I C A L FINITE E L E M E N T SOLUTION The model is based on a fixed domain discretization and includes the saturatedunsaturated groundwater zone and salt water region. So the free surface and the interface are not boundary conditions, but are parts of the problem. The position of these moving interfaces can be obtained by the pressure equilibrium that exists between groundwater and atmosphere or saltwater. The hydraulic head is expressed by:
h=z+
P = z +~ Pzg
(2)
= P is the pressure head, p is the interstitial pressure, pf is the freshwater density and g is
Pzg
the acceleration due to the gravity. The atmospheric pressure is often taken equal zero, so the free surface is the geometric location of points where h -- z. Because of the quasi stationary saltwater zone assumption, the hydrostatic pressure distribution is considered. Hence, in the saltwater region, the pressure is expressed by :
(3)
p : -p, gz
The interface must verify:
h=--~
,
6 = p ' -Pp' ~
(4)
In the present approach a fixed finite hexahedral element mesh is used to discretize the whole domain, including the saturated and unsaturated zones as well as the saltwater region. It is assumed that the groundwater flow occurs only in the freshwater region of the saturated zone and no flow is considered in the unsaturated and saltwater regions (the salt water is at rest). The water table and the interface are effectively hold impervious boundaries, but their location is unknown initially and must be obtained as part of the solution. For time discretization, a )~-weighted finite difference scheme is used ()~ = 0.5, Crank-Nicolson ; k = 1, backward Euler). The Galerkin finite element discretization yields the system of non linear algebraic equations: h ~ + I _
G(hk + ~)h ~+~ + C(h k+~)
z~;k + 1
h ~
=qt~( +~)
(5)
where h k + z = 2h k § 1 + (1 - 2)h k , h is the vector of nodal heads, superscript k denotes time step, G is the conductance matrix, C is the capacitance matrix and q is a vector containing the
541 boundary conditions. The coefficients of the conductance and capacitance matrices are given by: (6)
C~ = E ~SN, N/dVe
(7)
e
N~ and Nj are basis functions respectively for nodes i andj
3.1. Characteristic equations The flow occurs naturally in the saturated freshwater zone, however there is a small but negligible amount of water movement in the unsaturated and salt water zones. The flow in the unsaturated zone can be described [8] by :
o(~,) =
Kr(~/y,)
{0 r + (0 s - 0r)[1 + fll-m o,
_.
{I 1
+ /3)
.:mE< 2
1 if" /8)
m
-- ~ m
~,_<0
(8)
~u>0
j
2
~t'_< 0
(9)
~>0
Or is the residual moisture content, 0s is the saturated moisture content, fl =
~u
, ~s is the
capillary or air entry pressure head value, n can be interpreted as a pore size distribution index (mn = n -1, and (1.25 < n < 6)). A similar approach is also employed to describe flow in the saltwater region, which is transformed to an equivalent unsaturated zone where the pressure at any point of this sub-domain is defined by (~-p~). Hence, the following equations are used : Or + (0 s - Or)[1 +/3] -m o,
o(~,) =
K~(~) =
1 + fl) 2 (1 + fl)m _ tim
P'
r
(10)
P~
(11)
r
ps is the pressure head on the freshwater-saltwater interface, Or is the "residual water" in the equivalent unsaturated zone.
,8=
/ ~"- Ps / n ~s
The proposed relative conductivity in the saltwater zone (11) can be reduced in order to satisfy the hydrostatic assumption. The general storage term is expressed in the unsaturatedfresh saturated region as :
542 n-
(n - 1)(0,
dO
1
- or ~l~l
~<0
I
S =
dp'dO ,-~
with fl =
P'
Igsl"(l+~)m+l
(12) ~_>0
= S~
and in the freshwater-saltwater region by : dO
S=
( n - 1 ) ( 0 s - Or ) l ~ , -
~-
p
,,-1
~'
with fl =
I~i -
Iq/s[~(l§
Ps
n (13)
dO -d-~ = s , 9
~>p,
The storage term in the saltwater zone is also reduced, so that the hydrostatic assumption can still be valid. Note that the concept of approaching the flow in the 3 subregions as described in this section was also studied by [9], but using the so called the fast updating procedure. 3.2. Linearization technique Picard iteration scheme is used to solve (5) and gives : (m)
AG k+~ +
C
Ah
At~+l (m+l)
-
h~+~ -
+
qk+,~
(14)
A t ~+l (m)
where A h = h k +l _ h ~ +1 The resulting linear equations system is solved by the preconditioned conjugate gradients method with an incomplete factorization preconditioner based on M-matrix technique [ 10]. The updating process is repeated iteratively until there is no more significant change in the potential values over the fresh saturated flow domain. 4. VERIFICATION OF THE NUMERICAL MODEL
To verify the numerical model described above, several examples are considered, for which analytical solutions are available. Two examples are selected in this paper, which involve saltwater upconing and steady state flow with homogeneous isotropic hydraulic conductivity and using the Dupuits assumption. 4.1. Upconing in a confined aquifer overlain by an aquitard This problem is a field case [6], with the site-specific parameters shown in Table 1. Four simulations have been performed in steady state [11]. Both the numerical and the analytical position of the interface below the pumping well are plotted in Figure 1. The comparison between both solutions showed a good agreement except in the vicinity of the pumping well where the divergence starts to be relatively important for the last case. This is due to the horizontal flow approximation around the well assumed by [6].
543
Table 1 9Parameters used in example 4.1. Parameter Value -438.8 m Initial elevation of the saltwater-freshwater interface 10.9m Initial freshwater head 9.26x 10 "12 m 2 Horizontal intrinsic permeability k• kyy Leakance of the confining bed k'z/b' 1.2x 10 .3/day (in motz (1992)) Freshwater head in overlying constant head source bed 10.97 m Freshwater head at constant head boundary nodes 10.97 m Bottom elevation of well screen -260.0 m Length of well screen 221.0 m Distance from well bottom to the initial interface 178.8 m
4.2. Upconing beneath a pumping well in a coastal unconfined aquifer This example deals with seawater intrusion in an areal plane of an unconfined aquifer with a fully penetrating well near the coast line in steady state. The well withdraws freshwater only, causing upconing of underlying salt water [4]. The parameters used in this example are shown in Table 2.
Table 2 9Parameters for example 4.2. Parameter Value 1 m2/day (at x = 1000 m) Fresh water flux Qf 400 m3/day Well pumping rate Qw (6oo m, 0) Well location (xw,yw) 70 m/day Hydraulic conductivity, K 20m Depth from mean sea level to aquifer base The domain is composed of 1000 mx4000 m respectively along x and y directions. But, due to the symmetry upon the y axis, only the upper part of the domain was considered [ 11 ]. The simulated and analytical interfaces and water table positions are both plotted in Figure 2. The results show that there is a very good agreement between predicted water table positions and interfaces.
0
--" I
Dot
lines : Motz analytical
interface
Dashdot lines : Present
^
numerical
ution (1976) -s
Q = 112 Kg/s
'
D•
lillllllllll
f[ointhe well (,,,)
-20 Horizontal
Figure 1. Interface location for both analytical and numerical solutions for example 4.1.
diltance
(
m
)
Figure 2. Interface and water table positions for both analytical and numerical results for example 4.2.
544 4.3. Three-dimensional model for the Martil coastal aquifer The Martil aquifer is situated in the North-Western part of Morocco on the Mediterranean coast. It is the important local aquifer for potable water supply, irrigation and industry in the Martil region. In accordance with geology, the aquifer system consists of two aquifer units, separated by an aquitard. The upper unit is a Quaternary alluvial deposits of Martil river, the lower unit is composed of sandstone-limestone Pliocene formation, whereas the aquitard is generally marly and clayey. The aquifer thickness is very irregular in different directions. Figure 3a shows a cross-section from West to East. The conceptual model used in this study is a multi-layer aquifer system of variable thickness. The finite element mesh used for the numerical simulation of the Martil aquifer, shown in Figure 3b, is adjusted to fit the structure and the extension of the hydrogeological layers and the wells location. The 3-D mesh contains 56232 nodes, and 48288 hexahedral elements. The saturated hydraulic conductivity for different aquifers are Ks = 3.5 m/day in the alluvial material (upper), Ks = 4.5 m/day for the sandstone-limestone (lower), and Ks = 0.05 m/day for the marl/clay aquitard [9]. The other parameters are held constant (0s = 0.45, 0r = 0.1). West
Figure 3a. Martil West-East cross sectional view.
East
Figure 3b. Martil 3-D structured mesh.
The boundaries are treated as impervious, except at the eastern part, which is in direct contact with the sea (outflow face). Its extent is not known and constitutes a part of the problem. We also consider Dirichlet conditions at nodes belonging to the Martil and Alila rivers. In order to take into account the dry seasons these values are not set necessarily to be equal to the elevation. For the effective rainfall, we use a non-uniformly distribution taking into account the variably material distribution on the top layer. In a first step steady state simulations are performed under natural conditions such that natural groundwater flow can be reproduced efficiently. The purpose is to adjust conductivities and effective recharge. The obtained results are compared to observed piezometric map of 1966. The calibration is done using the trial and error technique. This procedure is repeated until a satisfactory pattern is found. Computed groundwater potential heads, shown in Figure 4, are compared to recorded values at the observation wells in 1966. The correspondence is relatively good. The 1966 interface location is also illustrated in Figure 5.
545
N
Figure 4. Computed (left) and observed (right) steady state groundwater potentials in the phreatic aquifer for 1966. A second set of simulations is performed starting from the 1966 head distribution, to predict the seawater intrusion in the aquifer. A transient simulation is made of over 30 years, considering variable pumping rates used for water public supply shown in Table 3. For irrigation, an amount of 7949 m3/day is pumped from several wells to irrigate 300 ha of agricultural area. The recharge is maintained constant along the simulation time. Table 3 : ONEP pumping rates Pumping rate (m3/day) Well number 992/2 780 994/2 1351.56 995/2 1181 996/2 418.75 1153/2 1260 1174/2 248.4 The results shown in Figure 6 show that the saltwater interface starts to move until 1982 and then remains stable along the time, although pumping wells withdraw upstream in the aquifer. This simulation scenario shows that there is no significant displacement of the interface due to the ONEP pumping wells upstream of the Martil river. West
East
T \
Figure 5. Plane view of the 1966 interface toe location.
Figure 6. Cross-sectional view of the moving interface positions from 1966 to 2014 (a new steady state regime is reached since 1982).
546 5. CONCLUSION A numerical model based on a 3-D finite elements is developed to predict the location, the shape and the extent of the moving interfaces in coastal aquifers involving unconfined or confined flow. Verification of this numerical model for upconing problems was made using the available analytical solutions to find the interface location and its dome extent. The numerical and analytical results are generally in good agreement. The model was also applied to the coastal aquifer of Martil in Morocco to simulate groundwater flow and to study the impact of the ONEP pumping well rates used for water public supply. The simulation results show that there is no significant displacement of the interface. ACKNOWLEDGEMENT
This work has been supported partially by the PARS SPI060 research project. REFERENCES
1. C. I. Voss and W. R. Souza, AQUIFE-SALT : A finite element model for aquifer containing a seawater interface. Water Resources Investigations Report 84-4369, US Geological Survey (1984) 37. 2. J. Bear and A. Veruijt, Modeling Groundwater Flow and Pollution, Reidel, Dordrecht, Holland, 1987. 3. M. A. Collins and L. W. Gelhar, Seawater intrusion in layered aquifers, Water Resour. Res., 7(4) (1971 ) 971-979. 4. O. D. L. Strack, A single-potential solution for regional interface problems in coastal aquifers, Water Resour. Res., 12(6) (1976) 1165-1174. 5. O. D. L. Strack, Groundwater Mechanics. Prentice Hall, Englewood Cliffs. NJ, 1987. 6. L. H. Motz, Saltwater upconing in an aquifer overlain by a leaky confining bed. J. Ground water 30(2) (1992) 192-198. 7. J. W. Bower, L. H. Motz and Durden D. W, Analytical solution for determining the critical condition of saltwater upconing in a leaky artesian aquifer. J. Hydrol. 221 (1999) 43-54. 8. M. T. Van Genuchten, A closed form Equation for predicting the Hydraulic Conductivity of Unsaturated Soils, Soil Sci. Soc. Am. J., Vol. 44 (1980) 892-898. 9. M. A. Sbai, F. De Smedt and A. Larabi, A Generalized Approch for Modeling 3D Transient Free and Moving Boundaries in Coastal Aquifers, First International Conference on Saltwater Intrusion and Coastal Aquifers-Monitoring, Modeling, and Management. Essaouira, Morocco, April 23-25 (2001). 10.A. Larabi and F., De Smedt, Solving 3-D hexahedral finite elements groundwater models by preconditioned conjugate gradients method, Water Resour. Res., 30(2) (1994) 509-521. 11.A. Aharmouch and A. Larabi, Numerical Modeling of Saltwater Interface Upconing in Coastal Aquifers, First International Conference on Saltwater Intrusion and Coastal Aquifers-Monitoring, Modeling, and Management. Essaouira, Morocco, April 23-25 (2001).