The location of deep salinity sources in the Israeli Coastal aquifer

Journal of Hydrology 250 (2001) 63±77

www.elsevier.com/locate/jhydrol

The location of deep salinity sources in the Israeli Coastal aquifer Uri Shavit*, Alex Furman 1 Faculty of Agricultural Engineering, Technion, Israel Institute of Technology, 32000 Haifa, Israel Received 3 May 2000; revised 12 March 2001; accepted 17 April 2001

Abstract The salinization process of the Israeli Coastal aquifer has led to an average concentration of about 200 mgCl/l with a signi®cant number of discrete salinity plumes in the middle and southern regions. The salinity of these plumes is high (500±1000 mgCl/l) and is increasing rapidly. Geochemical evidence has suggested that the salinity source in the Be'er Tuvia plume (in the south part of the aquifer) is at the bottom of the aquifer. This paper describes a solution of the source inverse problem and its application in the Be'er Tuvia plume. A transient two-dimensional ®nite element model was solved and the source terms were computed at each node in a 14 £ 14 km 2 area. An error analysis has shown that when no errors are introduced in the input data the reconstruction is perfect. The results of a sensitivity analysis are presented and the actual reconstruction errors are estimated. Applying the model in the Be'er Tuvia region indicates that a salinity source exists about 1 km to the west and 1.5 km to the north of the center of the salinity plume. This source is believed to be the plume source. q 2001 Elsevier Science B.V. All rights reserved. Keywords: Groundwater; Inverse problem; Finite elements; Salinity sources; Israeli Coastal aquifer

1. Introduction Water resources in Israel suffer from severe salinization (Hydrological Service Annual Report, 1998). The increasing salinity is of special concern in the Coastal Aquifer that provides about 20% of Israel's total use and serves as the largest storage volume. The Coastal Aquifer has an area of about 1900 km 2 underlying the coastal plain, where much of the population and economic activity are concentrated (see Fig. 1). Pumping water from the aquifer (up to 450 MCM), at times in excess of its natural recharge (about 340 MCM), has caused intrusion of polluted water into it * Corresponding author. Tel.: 1972-4-829-3568; fax: 1972-4822-1529. E-mail address: [email protected] (U. Shavit). 1 Present address: Department of Hydrology and Water Resources, University of Arizona, Tucson, Arizona 85721, USA.

(e.g. sea water). Agricultural, industrial, and other human activities have caused serious deterioration in its water quality. There are today whole zones of the aquifer whose water is saline and cannot be used without treatment. In recent years, the Coastal Aquifer has reached a chloride average concentration of about 200 mg/l (Hydrological Service Annual Report, 1998). Its salinity continues to rise in a mean annual rate of 2±3 mgCl/l. Fig. 1 shows that apart from the continuous increase of the background salinity, a signi®cant number of salinity plumes were discovered in the middle and southern parts of the aquifer. Pumping in these regions has been abolished and data indicates that the area covered by some plumes is spreading with time. The rate of salinization in the center of these plumes is about ®ve times higher than that of their surroundings and the chloride concentration is as high as 1000 mgCl/l. As sewage water is expected to

0022-1694/01/$ - see front matter q 2001 Elsevier Science B.V. All rights reserved. PII: S 0022-169 4(01)00406-1

64

U. Shavit, A. Furman / Journal of Hydrology 250 (2001) 63±77

Fig. 1. The Israeli Coastal aquifer. Map and data are based on a report prepared by the Israeli Hydrological Service, February 1998. The Be'er Tuvia plume is shown at the south part of the aquifer.

U. Shavit, A. Furman / Journal of Hydrology 250 (2001) 63±77

Fig. 2. Chloride concentration and water level history in wells around the Be'er Tuvia plume. Wells inside the plume and around it are marked in Fig. 3 (see also Fig. 8). Water level plot shows the decline prior to the mid 1960s followed by a recovery.

provide most of the water for agriculture and since groundwater has been exploited above its natural replenishment, the problem of salinization is becoming crucial and deserves more attention. Until recently the common thought was that salinization of the aquifer is mainly due to direct human activities. A recent geochemical study (Vengosh and Ben-Zvi, 1994; Vengosh et al., 1999) shows, however, that the rapid salinity increase in some of the plume areas is not anthropogenic. Using the chemical signature of a variety of potential sources and applying isotope techniques, it was shown that the groundwater

65

in the Be'er Tuvia plume is different from any known possible source (Vengosh et al., 1996, 1999). The water pollution sources that were eliminated include agricultural return ¯ows, industrial wastewater, arti®cial recharge from the Sea of Galilee, seawater intrusion, and lateral ¯ux from the eastern aquifer. It was suggested that the salinity increase is caused by intrusion of high salinity sources from subaquifers that lie under the active aquifer. Historical deep drills in nearby regions indicate that such a subaquifer exists and its chemical characteristics may be similar to that of the Be'er Tuvia plume (Vengosh and Ben-Zvi, 1994). The objective of the current study was to identify and quantify the salinization sources, reported by the geochemical studies, using an inverse modeling approach. In the paper we present a horizontal twodimensional numerical solution to the source identi®cation problem. The source terms that presumably introduce high salinity water from the subaquifer were computed. Groundwater level and chloride concentration data from 66 observation wells (for the years 1970±1990) were used as the input data. A Galerkin ®nite element scheme was used with triangular elements whose corner locations are at the observation wells. Following an error estimation and sensitivity analysis, both water and chloride source terms were calculated. The numerical results of our study are in an agreement with the conclusions of the previous geochemical results. The paper begins with a description of the Israeli coastal aquifer followed by the theory used for solving the inverse problem. Then, the procedure for the collection of the data needed for the model is described. The model was tested against an arti®cial case study and the level of accuracy was estimated. Finally, the results of the sources reconstruction are presented and discussed. 1.1. The hydrology of the Israeli Coastal aquifer (the Be'er Tuvia plume) Fig. 1 shows the layout of the Coastal Aquifer. It is composed of Pliocene-Pleistocene calcareous sandstone and layers of clay (Gavish and Friedman, 1969). Its thickness varies from 200 m near the coastline to a few meters in the east. The aquifer in its central, southern, and eastern parts is phreatic and

66

U. Shavit, A. Furman / Journal of Hydrology 250 (2001) 63±77

Fig. 3. A map of the Be'er Tuvia area showing the 66 grid nodes. 113 ®nite triangle elements were generated to connect these nodes. The 91 wells from which pumping and arti®cial recharge data was used are not marked. Fig. 2 shows water level and salinity at the marked nodes.

considered uniform. In general, the aquifer overlies thick (up to 2000 m) and impermeable units of the Saqiye Group. Vinokrurov et al. (1991) and Rosenthal et al. (1992) showed that in some parts in the eastern coastal plain the Yaffo Formation aquiclude is totally missing where in others it is deeply incised by trenchlike troughs ®lled with sand and sandstone. Along the eastern margin, the aquifer overlies the Cenomanian (Judea Group) aquifer in the central region and the Eocene (Shefela Group) aquitard in the northern and southern regions. The history of water level and chloride concentration in and around the Be'er Tuvia plume is shown in Fig. 2. It shows that over-exploitation caused a continuous drop in the piezometric levels during the 1950s and the 1960s. In 1964 pumping was reduced when the Israeli National Carrier started to provide water from the sea of Galilee to other regions of the country including the south of the coastal plain. As a result, water levels started to rise. The hydrological depressions did not disappear and groundwater

continued to ¯ow towards their center. In the professional community there were expectations that the recovery of the hydrological depressions would be followed by a decrease of the plumes salinity. These expectations were not ful®lled (Hydrological Service Annual Report, 1998). Fig. 2 shows the chloride concentration of four observation wells. The locations of wells 3 and 6 are within the salinity plume where the location of the others is outside the plume (Figs. 3 and 8). Fig. 2 shows that the rate of salinization inside the plume is larger (,10.9 mg/l/year) than that of its surroundings (,3.6 mg/l/year). Curve ®tting suggests that the plume salinization rate is increasing with time. Although some possible mechanisms were suggested to explain the continuous salinity increase (Artzy, 1999), there is no comprehensive explanation. Vengosh and Ben-Zvi (1994) utilized ion ratios such as Br/Cl, Na/Cl, B/Cl, Ca/Mg, K/Cl, Ca/ (SO4 1 HCO3) and showed that the chemical composition of the plume water is different from that of the water of the previously mentioned possible sources.

U. Shavit, A. Furman / Journal of Hydrology 250 (2001) 63±77

Recent geochemical and isotopic results indicate that saline plumes in other parts of the central and southern Coastal Aquifer are also derived from the underlying saline sources (Vengosh et al., 1999; Artzy, 1999). Chloride-rich water bodies which lie at a depth of 100±300 m below the aquifer bottom were discovered by oil drills in the mid and late 1980s. Some of these water bodies have high salinity (occasionally over marine salinity) and high pressure. In the Be'er Tuvia region, the closest drill hole in which a similar water body was discovered is about 8 km to the northwest. Vengosh and Ben-Zvi (1994) suggest that the chloride source is water from the underlying con®ned water bodies with chloride concentrations of 570± 2700 or 800±3350 mg/l. Calculation of the water/solutes mass balance of the Be'er Tuvia plume is complex. In the east, water level gradients are steep but aquifer thickness is thin. In the west, water level gradients are very small but aquifer thickness is thick. The region surrounding the high salinity Be'er Tuvia plume is characterized by sudden changes in the head gradients and salinity that are relatively small all around the plume. The complex local hydrology means that ®rst-order approximations may lead to an inaccurate mass balance calculation and a more careful and detailed analysis is needed.

2.1. The inverse problem scheme Solutions of the inverse problem for source reconstruction using indirect methods were published and reviewed in the past (Yeh, 1986; Sun, 1994). Here, we utilize a direct method using the Galerkin formulation and a ®nite element inverse model. The model computes the mass ¯ow rate of potential sources at every node as a function of time. The general formulation is two-dimensional and can be solved for either phreatic or con®ned aquifers. For the case of the Israeli Coastal Aquifer, a twodimensional horizontal phreatic formulation was used to compute the unknown water source term, Q 2h 2 7…T7h† 2 G P 2 G R 2t

and the unknown solute source term, QS

2…bC† ~ 2 7…ubD7C† 2 G p Cp 1 7…ubVC† 2t

2 G R CR

…2†

where h is the aquifer head, S the storage coef®cient, T ˆ Kb the aquifer transmissivity, b the aquifer thickness (head±bottom), K the hydraulic conductivity, C the chloride concentration, V~ the groundwater velocity vector, D the dispersion coef®cient, and u is the porosity. G P and G R are the known ¯ow rates of pumping and recharge (arti®cial and natural) with the corresponding concentrations CP and CR. A solution for Q and QS will provide the location, ¯ow rate, and solute concentration, QS/Q, of the salinity sources in the two-dimensional aquifer. The components of the dispersion coef®cient are de®ned as D x ˆ aL Vx 1 aT Vy and Dy ˆ aT Vx 1 aL Vy with Vx and Vy ^ and a L as the velocity components …V~ ˆ Vx õ^ 1 Vy j† and a T as the dispersion constants in the longitudinal and transverse directions. The inverse Galerkin procedure (Frind and Pinder, 1973) was implemented to solve both Eqs. (1) and (2) and is demonstrated here for Eq. (1). We ®rst de®ne a function L as LˆS

2h 2 7…T7h† 2 Q 2 G P 2 G R 2t

…3†

The head, h, and both known and unknown source terms Q, G P, and G R are approximated by the sum of products of their value at the observation wells (nodes), i, and the basis functions, f i as follows:

2. Theory and application

QˆS

QS ˆ u

67

…1†

h<

n X iˆ1

Q<

h i fi ; G P <

n X iˆ1

Qi fi

m X jˆ1

G Pj fj dj ; G R <

n X iˆ1

G Ri fi ; (4)

where n is the number of nodes, m the number of pumping and recharge wells, and d j is delta function. The Galerkin formulation is derived by substituting the approximations (Eq. (4)) into the function L (Eq. (3)) and writing the approximated solution over the domain R as ZZ Lfi dx dy ˆ 0 …5† R

The weight function f i equals 1 at node i (the

68

U. Shavit, A. Furman / Journal of Hydrology 250 (2001) 63±77

Table 1 Model coef®cients Coef®cient

Value

Range for sensitivity analysis

Units

K S u aL aT Natural recharge factor

10 0.25 0.4 10 1 0.174

^5 ^ 0.05 ^ 0.05 ^5 ^ 0.5 ^ 0.075

m/day ± ± m m ±

observation well) and zero at the other nodes. The ®nite element formulation uses bilinear weight functions …fi ˆ ai0 1 ai1 x 1 ai2 y† and triangle elements. As only 66 wells were available, the element layout was distributed manually. Time derivative was formulated using a ®nite difference backward scheme, with a one-year time step. Pumping and recharge were accumulated during the past full year. An integration by parts of Eq. (5), a collection of all elements, and incorporation of the appropriate boundary conditions result in a set of equations from which both Q and QS can be calculated at each node AQ Q ˆ RhsQ ;

AQS QS ˆ RhsQS ;

…6†

where A is the coef®cient matrix and Rhs is the righthand-side vector. The reconstruction of source terms results in an algebraic equation rather than differential equations. Hence, boundary conditions are not required. However, when a numerical formulation is used, and the source term on the boundary is unknown, external ¯ux information is necessary. An internal (forward or backward) derivative and a Dirichlet boundary condition are not suf®cient to complete a mass balance on the boundary. In the current formulation, data from wells around the boundary were used to de®ne a ¯ux boundary condition. 2.2. Data collection The hydrological, geological, and meteorological data were collected to enable an analysis in a square region of nearly 14 £ 14 km 2 shown in Fig. 3 (118± 132, 118±132, old Israel net). Meteorological data from six stations were provided by the Israeli meteorological service. The stations cover a larger region to allow correct interpolation. Interpolation, both linear and quadratic, has shown that spatial variability is

always below 10% and that a uniform rainfall can be used for each year. Lithological characterization was obtained using four geological cross-sections (Tolments, 1977), providing both the depth of the Saqiye Group (aquifer bottom) and an estimate of the aquifer transmissivity, T. The hydraulic conductivity was taken as constant following the observation that the aquifer formation is relatively uniform and more detailed information is not available (see Section 2.3). Information such as well locations, water level, chloride concentration, and pumping/recharge rate was provided by the Israeli hydrological service. This includes observation wells, pumping wells, recharge wells, and recharge basins. Water level and chloride concentration were monitored in more then 125 wells within the above-mentioned square region. However, within the years 1970±1990 only 66 wells were found to be adequate. Wells were screened out when too much data was missing. In a few cases, single data points were spatially and temporally interpolated to create a complete and continuous set of data during the period 1970±1990 in all 66 wells. According to the Israel Hydrological service database, 91 wells were used for pumping and arti®cial recharge. The provided data was incorporated into the ®nite element Galerkin scheme. The pumping and recharge values were presented at the 66 grid nodes by the sum of products of node values and basis functions as in Eq. (4). Each one of the 91 wells was divided between the three nodes of the element according to its relative location. Table 1 summarizes the hydrological coef®cients used in the model. The coef®cients were estimated using both general literature and local information. Mechanical dispersion coef®cients, a L and a T, were matched with the Coastal aquifer lithology using the review papers by Gelhar et al. (1992) and Beven et al.

U. Shavit, A. Furman / Journal of Hydrology 250 (2001) 63±77

69

Fig. 4. Reconstruction accuracy as a function of input errors. Errors in head and chloride concentration were arti®cially introduced. The reconstruction results of Q (left axis), and QS (right axis) are shown.

(1993). The natural recharge coef®cient was chosen based on a number of estimates published in the local literature. These estimates were developed using an aquifer calibration in the Be'er Tuvia region around a 15 km region to the north (Mercado, 1976), and a total mass balance for the whole Coastal aquifer (Hydrological Service Annual Report, 1998). Porosity, storability, and hydraulic conductivity were estimated based on the Coastal aquifer lithology and a large number of references (Averbach 1958; Lee et al., 1980; Adams and Gelhar 1992; among others. For a complete list see Furman, 1998). As explained later, the choice of hydraulic conductivity was further examined using a calibration process. 2.3. Level of accuracy In most inverse problems the number of unknowns is larger than the number of equations and the problem is ill de®ned. Direct solution methods are known to generate large errors when the transmissivity is to be calculated and the head data is noisy. In the current study, the reconstruction of the source terms was obtained by using an equal number of given observations and unknowns. In that sense the problem is well posed. However, input data such as water level, chloride concentration, and the geohydrology of the aquifer contain errors. Therefore, the reconstruction accuracy as a function of incorrect input data must be evaluated. Since sensitivity to errors depends on

problem conditions, the best one can hope for is to undertake realistic case studies and investigate the reconstruction error as a function of noise in the measurements (Yakowitz and Duckstein, 1980). Hence, the current uncertainty analysis is based on two different tests. The ®rst test examined the ability of the inverse methodology to reconstruct the source terms of a synthetic case. The test was implemented before the actual Be'er Tuvia problem was solved. However, the chosen geometry and numerical values of the hydrological parameters were similar to those of the Be'er Tuvia case. In the second test, the differentials dQ and dQS were used to represent the errors in computing Q and QS. A linear error analysis was applied using an order of magnitude representation of Eqs. (1) and (2) and the actual data used for the reconstruction of the sources in the Be'er Tuvia region. 2.3.1. Error evaluation using a synthetic aquifer To examine the accuracy of the proposed method, sets of groundwater levels were generated using a ®nite element forward numerical scheme. Model geometry and parameters were chosen to maintain a similarity with the hydrology of the Be'er Tuvia region. First, a source term in a limited area was arti®cially introduced. The generated groundwater levels and chloride concentration were used as an input for the inverse problem model, and the source terms were reconstructed with negligible errors

70

U. Shavit, A. Furman / Journal of Hydrology 250 (2001) 63±77

Fig. 5. Reconstruction accuracy as a function of input errors. Chloride concentration was calculated by applying an arti®cial source and using a hydraulic conductivity, KForward in a forward numerical model. Incorrect hydraulic conductivity, KInverse, was then introduced and the inverse model was solved to compute the source term, QS. The reconstruction result of the chloride source, QS, as a function of the ratio between these two hydraulic conductivites is shown.

(0.003% in Q and 0.0039% in QS). Similar results were found for non-steady sources and for a random distribution of source value and location (Furman, 1998). Then, a uniformly distributed random noise was arti®cially introduced to the input data. Water level and chloride concentration were varied randomly within a range ^0.5 cm and ^1 mgCl/l, respectively. The reconstruction error generated as a result of random noise is demonstrated in Fig. 4. Three arti®cial source levels of both water and chloride were reconstructed. It shows that the error for zero source is small. However, deviation increased when higher source value was introduced. The intermediate source level results in errors less then 5% for the water source and less then 15% for the chloride source. As will be shown, the value of the source terms that were computed for the Be'er Tuvia region are similar to that of the intermediate source level. As expected, it was found that for such random error, time averaging of a constant source results in a nearly zero deviation. The estimation of the hydraulic conductivity is considered to be the most important and dif®cult. The hydraulic conductivity was varied, therefore, within a range of two orders of magnitude and its in¯uence on the reconstruction accuracy was tested. Fig. 5 shows that under the test conditions, underestimation of the hydraulic conductivity generates a relatively small error (10% for the intermediate source

intensity). Results show that the error in the water source reconstruction is even smaller (Furman, 1998). As shown, when the model hydraulic conductivity is overestimated, reconstruction errors are very large. The hydraulic conductivity of the Be'er Tuvia region was carefully chosen to avoid an overestimated erroneous value. Note, that erroneous conductivity by itself will affect the reconstructed source intensity but will not affect the accuracy of the source location. Inaccurate storage coef®cient, S, might introduce reconstruction errors. However, its estimation is typically more accurate than that of the hydraulic conductivity coef®cient. It was found that reconstruction errors increase linearly with an increase in the discrepancy in S. When the storage coef®cient was over or underestimated by 20%, the computation of the intermediate source level was about 20% higher or smaller, respectively. As with the in¯uence of inaccuracy in the hydraulic coef®cient, estimation errors of the storage coef®cient did not change the computed location of the sources. Since reconstruction errors cannot be avoided, the location of the source terms is considered to be more reliable than their speci®c value. In addition to the sensitivity analysis of K, S, h, and C, incorrect location of wells, incorrect topography of the bottom of the aquifer, errors in the aquifer boundary and initial conditions were analyzed. Incorrect topography of the bottom of the aquifer, for example,

U. Shavit, A. Furman / Journal of Hydrology 250 (2001) 63±77

71

Table 2 A list of variables that are potential error contributors to Q and QS. The table contains the derivatives of the source term as a function of these variables, the given characteristic values and the maximum estimated errors. This data was used with Eq. (7) to generate the results presented in Fig. 6 Variable, fi

Source derivative, u2Q=2fi u or u2QS =2fi u

Units of fi

Characteristic value of fi

dfi

Dh K S b G P, G R DC u a CP, CR V~ Dt L ˆ Dx, Dy

2Q=2…Dh† ù S=Dt 2 Kb=L2 ; 2QS =2…Dh† ù …1 2 a=L .†bKDC=L2 2Q=2K ù bDh=L2 ; 2QS =2K ù …1 2 a=L†bDCDh=L2 2Q=2S ù Dh=Dt; 2QS =2S ù bDC=Dt 2Q=2b ù KDh=L2 ; 2QS =2b ù uDC=Dt 2 u…aV†DC=L2 1 uVDC=L 2Q=2G P ˆ 2Q=2G R ˆ 1; 2QS =2G P ˆ CP ; 2QS =2G R ˆ CR 2QS =2…DC† ù ub=Dt 2 u…aV†b=L2 1 uVb=L 2QS =2u ˆ 0; 2QS =2…u† ù bDC=Dt 2 DC…aV†b=L2 1 DCVb=L 2QS =2a ù ubVDC=L2 2QS =2CP ˆ G P ; 2QS =2CR ˆ G R V~ ù KDh=uL

m m/day ± m m/day g/m 3 ± M g/m 3 m/day day m

5 10 0.25 60 0.00025 100 0.4 10 1000, 300 0.125 365 1000

0.5 5 0.05 6 0.00025 5 0.05 10 500, 10

was tested by modifying its depth both locally and globally. In both cases, a modi®cation of 10% results in a small error (1.5% in the water source and 3% in the chloride source). Note, however, that such incorrect input might generate errors that increase with time. 2.3.2. Error estimation using the source differentials dQ and dQS Following the reconstruction results produced by the model (see Section 3), an error analysis as a function of the model parameters and input data was

obtained. The errors were estimated using linear differential representation of each of the model parameters. The differentials dQ and dQS represent the errors in computing Q and QS dQ ˆ

X 2Q df ; 2fi i i

dQS ˆ

X 2QS dfi 2fi i

…7†

where dfi represents the maximum error of a given variable fi. As shown in Table 2, the derivatives 2Q=2fi and 2QS =2fi were calculated by using an order of magnitude representation of Eqs. (1) and

Fig. 6. Errors computed by the products …2Q=2fi † dfi and …2QS =2fi † dfi :

72

U. Shavit, A. Furman / Journal of Hydrology 250 (2001) 63±77

Fig. 7. Model calibration using MODFLOW. The sum of squares of the difference between the measured water head, hm, and computed water head, hc, at 40 wells in the Be'er Tuvia region.

(2). The variables fi that contribute errors to the reconstruction of Q and QS, can be sorted into two groups. Measured variables such as Dh, DC, CP, CR, G P, and G R (arti®cial recharge only), and estimated (and calibrated) parameters such as K, S, b, u , and a . Table 2 shows a summary of these variables, the derivatives 2Q=2f i and 2QS =2fi ; their characteristic values and their maximum errors, dfi. Some variables contribute errors to both Q and QS, whereas others contribute errors to QS only. The last three rows of Table 2 provide characteristic values only. The individual errors computed by the products …2Q=2fi † dfi and …2QS =2fi † dfi are shown in Fig. 6. By accumulating all these individual contributions from Fig. 6 (Eq. (7)), total errors of dQ ˆ 0:003 m=day and dQS ˆ 2:52 g=m2 =day were obtained. The errors dQ and dQS are about 10 and 30% of their corresponding maximum values, respectively. 2.3.3. Model calibration Model calibration is used to reduce errors when hydraulic parameters are unknown (Anderson and Woessner, 1992). Calibration of the current inverse model is inherently different from the calibration process in forward modeling where a comparison between measured water head, hm, and computed water head, hc, is obtained. In the inverse model presented here, water head and chloride concentration are used as input data and could not be applied for direct comparison. To overcome this dif®culty, forward modeling was obtained by assigning a given value to the source term. Using MODFLOW (Groundwater Vistas V1.0, Environmental Simulations), two-

dimensional forward solutions were obtained twice, once assuming that the unknown source, Q, is zero, and once with a source value equal to that found by the ®nite element inverse model (see Fig. 10). Assigning S ˆ 0:25; u ˆ 0:4; and natural recharge coef®cient of 0.174, water head (hc) in the Be'er Tuvia region was computed for the years 1970±1990. Results were generated every month. In cases where data was not available on a monthly basis, linear interpolation was used and a comparison with the MODFLOW computed results was obtained. When water level data (hm) was interpolated, it was assumed that the end of summer head represents the lowest level and the end of winter head represents the highest level. Forty wells were chosen for a least square analysis. Fig. 7 shows that a hydraulic conductivity of K ˆ 10 m=day is the best choice. Furthermore, it was found that within the range and distribution of the current source term values, the choice of a hydraulic coef®cient is not sensitive to the intensity of the source value. 3. Results and discussion Figs. 8 and 9 show the spatial distribution of chloride concentration and water level in 1990. The high salinity plume is shown in Fig. 8 with a concentration

Fig. 8. The Be'er Tuvia plume in 1990. In 1998 the chloride concentration at the plume center has reached 1000 mgCl/l.

U. Shavit, A. Furman / Journal of Hydrology 250 (2001) 63±77

73

Fig. 9. Water levels in the Be'er Tuvia region, 1990.

Fig. 11. Average chloride source, QS, 1970±1990.

of above 700 mgCl/l. The plume is isolated with less than 400 mgCl/l to the east and less than 300 mgCl/l to the west. Records show that the plume location does not change much with time but its concentration has increased and reached in 1998 a value of 1000 mgCl/l (Vengosh, 1998; Vengosh et al., 1999).

Fig. 9 shows that the Be'er Tuvia plume is located at the foot of a high gradient head slope. South and east from this plume the head gradient is high while a small slope exists to the north and west directions. These gradients are affected by the regional topography, the local geohydrology, and by the pumping and recharge practices. Despite the small gradients, ¯ow rate through the north and west boundaries is not zero at all times. Flow rate calculations consider both the data presented in Figs. 8 and 9 and the significant change in the aquifer thickness (200 m in the west and just a few meters in the east). Fig. 3 shows the spatial layout of the 66 wells. A grid of 113 triangle elements was generated by connecting these wells. Although automated algorithms for grid generation are available (Thompson et al., 1999), the small number of grid points allows a manual grid formation (using common considerations such as shortest available connections, center of gravity, and symmetry). A signi®cant advantage of the ®nite element method is that spatial interpolation is not required before mass balance calculations are conducted. Both the interpolation and equation solution are performed simultaneously. In that sense, the ®nite element method is based on linear interpolation and handles successfully the non-Cartesian grids and irregular boundaries.

Fig. 10. Average water source, Q, 1970±1990.

74

U. Shavit, A. Furman / Journal of Hydrology 250 (2001) 63±77

Fig. 12. Temporal variations of the water source, Q, and the chloride source, QS, at the `Azrikam source'. The insert is a reproduction of Fig. 11 with an arrow pointing at the `Azrikam source'.

The reconstruction of the source terms is shown in Fig. 10 (water source) and Fig. 11 (chloride source). These two ®gures represent an average source value for each grid point during 1970±1990, excluding negative values. The water source term of Eq. (1) is signi®cant in Fig. 10 at three major locations (around 119.8/128.6, 122.9/128.6, and 127.8/128.1) with some lower intensity sources around the center region between Be'er Tuvia and Ein Zurim. The source around 122.9/128.6 is considered to be the source of the Be'er Tuvia plume (`plume source'). It contributes more than the other sources to the salinization of Be'er Tuvia plume. It is postulated that the other two major sources are not necessarily sub-aquifer sources. As will be shown, the `Azrikam source' (119.8/128.6) is more likely to be identi®ed as a surface source. The nature of the `eastern source' (127.8/128.1) is less understood. Although no clear evidence was found around the surface of its vicinity, it is assumed that this source is also a surface source (or an artifact which was generated due to a local low conductivity). As shown in Fig. 11, the chloride sources were found in similar locations (119.8/128.6, 122.9/128.6, and 127.8/128.1). Such a ®nding strengthens the reliability of the reconstruction results. Furthermore, it suggests that the source term is a result of a signi®cant

convective mass ¯ow rather than a rock water interaction. To further examine the correlation between the water source term, Q, and the chloride source term, QS, a comparison was made between the temporal behavior of these two terms. It was found that in locations where a signi®cant source was identi®ed, the temporal behavior of the water and chloride sources has a similar shape. This similarity is nearly perfect in the Azrikam source. Assuming that the Azrikam source is due to arti®cial recharge, such a similarity is expected. The temporal similarity of Q and QS, which were found in the other two signi®cant sources, show excellent correlation in the eastern source and good correlation (with some variations) in the plume source. The high temporal correlation in the eastern source supports the hypothesis that it is a surface source. The plume source has an `eight' shape and is more signi®cant in its south half (123.1/127.9) than at its north half. Furthermore, this part of the source coincides with the location of the north part of the salinity plume (Fig. 8). Such a match supports the conclusion that the source found at 123.1/127.9 is the actual plume source, however, the salinity plume is stretched to the southeast with no matching salinity sources. Possible explanations are related to the threedimensional ¯ow problem that was analyzed here by a

U. Shavit, A. Furman / Journal of Hydrology 250 (2001) 63±77

two-dimensional model. The necessary vertical integration may create deviations when a deep source is to be reconstructed using well data from the upper layer of the aquifer. It may also be claimed that the hydrological depression of this area (see Fig. 2) could affect the plume location by convection from the northwest to the southeast. An examination of Figs. 8, 10, and 11 reveals that although signi®cant sources were computed in the Azrikam area, the salinity of this region is equal to the aquifer background salinity. A year-by-year presentation of both Q, and QS around Azrikam shows that the source terms did not exists prior to 1978±1979. As shown in Fig. 12, the water and salinity sources started to appear in 1979±1980. Following these computation results, the Israeli Hydrological service provided us with new data regarding an arti®cial recharge basin near Azrikam. Data shows that the Azrikam reservoir was used for arti®cial recharge of water from the Sea of Galilee. The archive includes the chloride information since 1981 and volume data since 1985. We were unable to sort out the 4±5 year difference between our ®nding, the chloride, and the volume data. Despite the ®ve years lack of information, this sequence of events served as an excellent example of the model capability to reconstruct signi®cant ¯ow and transport source terms. The solution of the inverse model as presented in Figs. 10 and 11 agrees with conclusions made by some of the previous studies, in particular, in the geochemical studies by Vengosh and Ben-Zvi (1994) and later by Vengosh et al. (1999). As described earlier, the authors suggest that the chloride source is from underlying con®ned water bodies with a chloride concentration of 570±2700 mg/l (or 800± 3350 mg/l). In order to compare the geochemical studies with the current study, the concentration at the sources was calculated. In general, the reconstructed chloride source, QS, was divided by the reconstructed water source, Q. This is not always possible since the water source, Q, is occasionally small, resulting in a very high and unrealistic concentration. Therefore, two methods were used. First, source concentration was calculated node by node, locally. Then, the water source and the chloride source were integrated over the area of the plume source and the area of the Azrikam source. Integration over the area covered by the sources provides the total average

75

annual contribution of these sources. The computed contribution of the sub-surface plume source is about 2 £ 10 6 m 3/year water and 1200 t/year chloride. The Azrikam surface source is about 5 £ 10 6 m 3/year water and 1100 t/year chloride. The concentration at the source, based on the integral values, shows that the plume source has a concentration of only 600 mgCl/l and the Azrikam source has a concentration of 220 mgCl/l. These results are considered more reliable than the nodeby-node calculation. However, the low 600 mgCl/l found in the plume source deserves an explanation as the plume concentration has reached in 1998 a value of 1000 mgCl/l. As most of the well data are from the upper one-third of the aquifer cross-section, it is postulated that both lateral and vertical dilution has led to this result. The node-by-node calculation was obtained when only both Q and QS were signi®cant. It was found that in the plume region, the source concentration varies between 300 and 5000 mgCl/l (with only a few nodes having more than 1000 mgCl/l). The average local concentration of the source at Azrikam was found to be between 160 and 330 mgCl/l). The chloride concentration found in Azrikam, using both calculations, is similar to the concentration of the Sea of Galilee (C < 220 mgCl/l) and similar to the salinity data provided by the Hydrological Service for the Azrikam recharge reservoir (Caverage ˆ 207 mgCl=l; Cmin ˆ 181 mgCl=l and Cmax ˆ 224 mgCl=l). The concentration at the source was expected to be higher due to the contribution of the unsaturated zone. Vengosh and Ben-Zvi (1994) have also presented an integral calculation and claimed that the salinity source contribution is about 1 g/m 2/day distributed uniformly under the Be'er Tuvia plume. Given a smaller area of the plume source (see Figs. 10 and 11), a somehow higher value of the source contribution is expected (3±4 g/m 2/day). Vengosh and BenZvi (1994) did not calculate the source water ¯ow rate. However, a combination of salt contribution rate and concentration could lead to similar results. The plume of Be'er Tuvia appeared for the ®rst time during the 1950s around 124.0/127.0. Based on this information, Vengosh and Ben-Zvi (1994) suggested that the plume location is also the location of the salinity source. Our inverse solution found that the source location is at 122.9/128.6, about 1 km to the west and 1.5 km

76

U. Shavit, A. Furman / Journal of Hydrology 250 (2001) 63±77

to the north. As noted, both the two-dimensional analysis of a three-dimensional ¯ow problem and the streamline direction during the hydrological depression can serve as an explanation to this calculated shift in the source location. 4. Mechanisms A major question that was not answered is the salinization mechanisms. A theory is needed to explain the increase in chloride concentration following the recovery of the hydrological depression. Theories have suggested that a lower artesian aquifer is the source of high salinity water ¯ow and the reason for the continuous salinity increases (Artzy, 1999). Other theories have postulated that underground changes have occurred during the depression and caused higher permeability at the interface between the active aquifer and the sub-aquifer. Neither the above theories nor the current analysis were able to provide, at this stage, a comprehensive explanation. 5. Conclusions Geochemical analysis has shown that the source of salinity in the Be'er Tuvia plume is not due to anthropogenic activity or lateral ¯ow. The geochemical studies suggest that high salinity water penetrates through the generally thick and impermeable units of the Saqiye Group at the bottom of the aquifer. Such a penetration is likely to occur through parts where the aquiclude is totally missing or through regions where the Saqiye Group is replaced by sand and sandstone. To accept or reject the above conclusion, an inverse problem model was developed and applied. A time dependent, two-dimensional, inverse solution of both the ¯ow and transport differential equations provides the volume ¯ow rate and the mass ¯ow rate of water and chloride source terms at each model node in the Be'er Tuvia plume area. Hydraulic head and chloride concentration were given and the hydrogeological model parameters were estimated. Error and sensitivity analysis has shown that when input data are error free, the source reconstruction is exact. However, when errors were arti®cially introduced in either the input data (head and concentration) or the hydrogeological parameters,

reconstruction errors are unavoidable. Using an order of magnitude representation of the governing equations, a linear error analysis was developed. The estimated uncertainty is about 10% for the water source and about 30% for the salt source. Most of these errors affect the source intensity more than its spatial location. The solution of the inverse problem provides an estimate of the location and intensity of the salinity sources. In general, the location of the water sources were found at the same location of the reconstructed chloride sources. Three major sources were found at the north half of the study area. A surface source at the northwest was reconstructed and later con®rmed to be a result of an arti®cial recharge by data from the Israeli Hydrological Service. The second major source was identi®ed as the sub-surface source. It was found at coordinates 122.9/128.6, about 1 km to the west and 1.5 km to the north of the center of the salinity plume. Results show that the ¯ux in the center of the source is about 4 g/m 2/day, with a water ¯ux of about 1 cm/day. Integration over the area covered by the source provides the total average annual contribution of these sources. The computed contribution of the sub-surface plume source is about 2 £ 10 6 m 3/year water and 1200 t/year chloride. The computed contribution of the Azrikam surface source is about 5 £ 10 6 m 3/year water and 1100 t/year chloride. Although the current study and the geochemical studies used a different approach, results show that both studies agree reasonably well. It is concluded that when adequate data is available, the direct inverse model can be used for source reconstruction analysis. Our inverse modeling approach is potentially superior since quantitative data, including the location of the source, the area covered by the source, and water and solute source ¯uxes, can be computed as a function of time. Neither geochemical analysis, nor model calibration can directly provide such quantitative analysis. It was shown that a sub-surface salinity source exists in the Be'er Tuvia area and causes the existence and growth of the salinity plume. However, the study results and conclusions must be treated with some caution. Since head and concentration vertical pro®les are unavailable, a two-dimensional analysis was obtained. The assumption that the aquifer is uniform and a constant conductivity can be used was adopted due to lack of detailed spatial information. It is clear

U. Shavit, A. Furman / Journal of Hydrology 250 (2001) 63±77

that this assumption holds in general but probably fails in some locations. Finally, the Be'er Tuvia plume is probably a rare case where direct inverse solution can be applied. When this method is to be applied in other regions, a great attention has to be paid to the reliability of the data and the speci®c hydrological circumstances. References Adams, E.E., Gelhar, L.W., 1992. Field-study of dispersion in a heterogeneous aquifer. 2. Spatial moments analysis. Water Resour. Res. 28 (12), 3293±3307. Anderson, M.P., Woessner, W.W., 1992. Applied Groundwater Modeling. Academic Press, San Diego, California. Artzy, Y., 1999. Salinization mechanisms of ground water in the coastal aquifer of Israel: the Gan-yavne Gedera region. Ms thesis, Department of Geological and environmental sciences, Ben Gurion University (in Hebrew with an abstract in English). Averbach, S., 1958. Aquifer parameters characterization using the strip calculation method. Tahal Report, Appendix 7b (in Hebrew). Beven, K.J., Henderson, D.E., Reeves, A., 1993. Dispersion parameters for undistributed partially saturated soil. J. Hydrol. 143, 19±43. Frind, E.O., Pinder, G.F., 1973. Galerkin solution of the inverse problem for aquifer transmissivity. Water Resour. Res. 9 (5), 1397±1410. Furman, A., 1998. Identi®cation of formation processes of salinity plumes in the coastal plain aquifer of Israel and possible solutions. Ms thesis, Faculty of Agricultural Engineering, Technion (in Hebrew with an abstract in English). Gavish, E., Friedman, G.M., 1969. Progressive diagnosis in quaternary to late tertiary carbonate sediments: sequence and time scale. J. Sediment. Petrol. 39, 980±1006. Gelhar, L.W., Welty, C., Rehfeldt, K.R., 1992. A critical review of data on ®eld-scale dispersion in aquifers. Water Resour. Res. 28 (7), 1955±1974. Hydrological Service Annual Report, 1998. Ground Water

77

Resources in Israel, Fall 1997. Hydrological Service, ISSN0793-1093, Jerusalem, 282pp. (in Hebrew). Lee, D.R., Cherry, J.A., Pickens, J.F., 1980. Groundwater transport of salt tracer through a sandy laebed. Limnologic Oceanogr. 25 (1), 46±61. Mercado, A., 1976. Nitrate and chloride pollution of aquifers: a regional study with the aid of a single-cell model. Water Resour. Res. 12 (4), 731±747. Rosenthal, E., Vinokurov, A., Ronen, D., Magaritz, M., Moshkovitz, S., 1992. Anthropogenically induced salinization of groundwater: a case study from the Coastal Plain aquifer of Israel. J. Contaminant Hydrol. 11, 149±171. Sun, N.Z., 1994. Inverse Problems in Groundwater Modeling. Kluwer Academic Publishers, Dordrecht. Thompson, J.F., Soni, B.K., Weatherill, N.P., 1999. Handbook of Grid Generation. CRC Press, Boca Raton, FL. Tolments, Y., 1977. The Israeli hydrogeological atlas, geological cross section along the coast. Hydrological Service. Vengosh, A., 1998. Personal communication. Vengosh, A., Ben-Zvi, A., 1994. Formation of a salt plume in the coastal plain aquifer of Israel: the Be'er Toviyya region. J. Hydrol. 160, 21±52. Vengosh, A., Artzi, Y., Zirlin, I., Pankratov, I., Harpaz, I.H., Rosenthal, E., Ayalon, A., 1996. Monitoring water quality in the Coastal Plain aquifer of Israel (Givat Brener, Gedera, Yavne). Hydrological Service, Hydro. Report 4/96, ISSN0334-3367, Jerusalem, 34 pp. (in Hebrew). Vengosh, A., Spivack, A.J., Artzi, Y., Ayalon, A., 1999. Geochemical and boron, strontium, and oxygen isotopic constraints on the origin of the salinity in groundwater from the Mediterranean coast of Israel. Water Resour. Res. 35 (6), 1877±1894. Vinokrurov, A., Rosenthal, E., Ronen, D., Magaritz, M., 1991. The geological structure of the Petah-Tiqva Kfar Sava area and its possible in¯uence on groundwater salinization processes. Hydrological Service of Israel, Jerusalem, Rep#4/91. Yakowitz, S., Duckstein, L., 1980. Instability in aquifer identi®cation: theory and case studies. Water Resour. Res. 16 (6), 1045± 1064. Yeh, W.W-G., 1986. Review of parameter identi®cation procedure in groundwater hydrology: the inverse problem. Water Resour. Res. 22 (2), 95±108.