The hydrogeological role of an aquitard in preventing drinkable water well contamination: a case study.

Groundwater pollution has become a worrisome phenomenon, mainly for aquifers underlying industrialized areas. In order to evaluate the risk of pollution, a model of the aquifer is needed. Herewith, we describe a quasi-tridimensional model, which we applied to a multilayered aquifer where a phreatic aquifer was coupled to a confined one by means of an aquitard. This hydrogeological scheme is often met in practice and, therefore, models a number of situations. Moreover, aquitards play and important role in the management of natural resources of this kind. The model we adopted contains some approximations: the flow within the aquifers is assumed to be horizontal, whereas leakage is assumed vertical. The effect of some wells drilled in these aquifers is also taken into account. In order to evaluate the leakage fluxes that correspond to different exploitation conditions, we numerically solve a system of quasilinear and time-dependent partial differential equations. This model has been calibrated by the hydrogeological data from a water supply station of the Milan Water Works, where water is polluted by some halocarbons. Our simulations account for several experimental facts, both from the hydrogeological and hydrogeochemical viewpoints. Maxima of computed downward leakage rates are found to correspond with measured pollutant concentration maxima. Other results show how the aquitard can help in minimizing the contamination of drinkable water.


Introduction
Groundwater pollution by nonbiodegradable compounds is a widespread phenomenon affecting most aquifers underlying industrialized areas where said compounds were or are being disposed of improperly. On the other hand, the deeper the aquifer layers, the less likely they are to be polluted by substances that percolate down from the surface. In order to roughly estimate or predict the risk of this kind of pollution, a model of the aquifer is needed. When this model is calibrated to a particular area, it becomes an invaluable simulation and resource management tool. In particular, it helps in determining the optimal water pumping policy and in evaluating the effectiveness of certain actions aimed at removing the pollutants.
The model we consider herewith relates to the coupling of a phreatic aquifer to a confined one by means of an aquitard. To the model equations a "quasi-tridimensional" scheme is applied. Time-dependent viscous flow in a nonhomogeneous porous medium is assumed to be two-dimensional (in the horizontal plane) in both aquifers and one-dimensional (along the vertical direction) across the aquitard. For drinking water wells that pierced into the aquifer, we evaluated how the leakage depends on aquitard conductivity, the rate at which water is drawn, and the locations of the wells. The resulting coupled differential equations were discretized by finite difference schemes, then numerically solved by means of a computer code, developed on a VAX-8600 computer.
This model is applied to a water supply station of the Water Works of Milan, which supplies drinking water to the city where hydrogeological data supports the model on a local scale; moreover, the phreatic aquifer there is known to be highly polluted by some chlorinated halocarbons that are also found in the confined aquifer at concentration levels high enough to make water non-PONZINI ET AL. potable. The values of all relevant hydrogeological parameters of the area appearing in the flow equations come from experiments, i.e., pumping tests, well specific capacities, etc.

Regional Geological Overview
The city of Milan is situated at the center of the large alluvial plain of the Po River. From the ground level downward the following geological formations are found: a) Continental and Transition Quaternary, which is 250 to 300 m thick. It is lithologically composed of a sequence of gravel, sand, and clay; the thicker terms prevail near the surface and the thinner terms at deeper levels. b) Marine Quaternary, about 400 m thick, mainly consists of fine-grained terms, e.g., clay, clay interspersed with turf, sands, and thin sands. c) Upper and Medium Pliocene, the top lays at about -700 m (hereafter, we shall denote the depth below ground surface level by a negative number) and is essentially clay.
Milan's water supply comes entirely from an aquifer pierced by 800 wells owned by the Water Works and by about 1000 private wells. The city territory extends over an approximate area 270 kin2, where 38 municipal water supply stations are scattered ( Fig. 1) Each station collects the water pumped by 4 to 25 wells, which are about 100 m deep and adequately spaced. Most wells draw their water from gravel and sand layers. The aquifer is very rich; the typical specific capacity of a well is 15 L/sec/m, whereas the average yearly consumption of the city is about 3 x 108 m3.

Local Hydrogeological and Hydrogeochemical Data and Their Interpretation
The Espinasse water supply station (Fig. 2) consists of 16 wells homogeneously distributed over an area of 0.20 km2. The water of each well has been contaminated for the last 10 years by 200 or more parts per billion (ppb) of chlorinated halocarbons. Therefore, all wells were disconnected from the drinking water supply network and have been used as purge wells, continuously discharging into the sewage drains for the last 10 years.
The data of interest to us mainly come from (1-7) and from an unpublished technical report (8). The data can be classified as follows: General well data * Location, depth and well top levels referred to the mean sea level * Field lithological boring logs * Construction and hydraulic data, e.g., screened interval, etc. The location, the lengths of screened intervals, and the field lithological boring log of a typical well of the Espinasse Water Supply Station are shown by Figure 3. Of great interest to us is the following: Fact 1: Each well contains more than one screened interval, i.e., it draws water simultaneously from more than one layer. This feature cannot be easily modeled, neither can hydrogeological and hydrogeochemical data from a well be uniquely ascribed to a single aquifer layer.
Hydrogeological data * More than 10 piezometric surveys relative to different aquifer exploitation conditions are available The piezometric surface obtained on May 7, 1980, when all wells had been shut down for 14 to 16 hr, is represented by Figure 4. Another piezometric surface, obtained on June 10, 1980, when all wells had been discharging 40 L/sec for 6 to 8 hr is represented by Figure  5. All plots have been produced by linear interpolation and subsequent smoothing. * Specific capacities Their mean value is 14.7 L/sec/m with a standard deviation value of 6.7 L/sec/m. * Two pumping tests During the first test, well 13 was discharging, whereas wells labb, 2a, 3, and 6 were used as piezometers. During the second test, well 6 was discharging, whereas wells 5, 9, and 10 acted as piezometers. Both tests were carried out for 48 hr with a constant pumping rate. The hydrogeological parameters computed by first interpolating the data then applying both the Hantush and Theis and Jacob methods (9)  5.0* 10-2 1.5. 10-' 1.6* 10-4 1.6. 10-6 Well 5 7.9 * 10-2 4.0 * 10-3 2.0 * 10-4 Well 9 8.5 * 10-2 4.3 * 10-3 6.7 * 10-5 Well 10 9.9. 10-2 5.1. 10-1.2* 10-* Phreatic surface data The parameters come from four sample points inside the Milan area, one of which is located at the Espinasse station. The latter is a 25-m deep well, drilled about 10 m from a deeper drinking water well. The free surface level and the head measured contemporaneously at the deep well are shown by Figure 6 as a function of time over an interval of some years. Both levels have been measured under steady-state conditions while no well was pumping.
Hydrogeochemical data Chemical analysis of water samples started in 1977. Hereafter, wells have been sampled at time intervals ranging from 10 to 30 days. In 1977, only 6 of the 16 sampled wells were significantly contaminated. The Milan Water Authority then decided to use these wells as purge wells, continuously discharging water at a rate of about 30 L/sec. Thereafter, all other wells became significantly contaminated and underwent the same treatment. The chlorinated halocarbons found are trichloroethylene, tetrachloroethylene, trichloromethane, carbon tetrachloride, and methyltrichloromethane. Concentration data of a given species measured at wells under different operating conditions have been interpolated by a logarithmic law then smoothed in order to yield contour maps. Figures 7, 8, and 9, respectively, show trichloroethylene, tetrachloroethylene, and methyltrichloromethane contour plots at a time when all wells were discharging. The following is also of interest to us: area where the well's density is higher. This same feature is found by similarly processing concentration data coming from other sites in the city (2,5).

Adopted Hydrogeological Scheme
The boring logs of the wells have been grouped in six lithological sections, three with a north-south, and three with a west-east direction; similar lithological terms have been correlated. Two of these lithological cross sections are shown in Figures 10 and 11; section WE includes wells 9, 13, 2, 14 and 1; and section NS includes wells 15, 2a, 3, and 4. Given the high well density and the relatively homogeneous lithological nature of the aquifer, we considered the correlation results to be very good. Furthermore, the site hydrogeology and hydro-geochemistry are based on the trend of phreatic and piezometric levels shown in Figure 6 and on all data and test interpretations listed in the previous section.
We have adopted an hydrogeological scheme where we identify four basic units: a) an upper leaky phreatic aquifer, b) an intermediate semipervious layer or aquitard, c) a lower leaky confined aquifer, and d) the bedrock.

Leaky Phreatic Aquifer
Its lithology varies from gravel to sand and gravel to clean sand; some rare, thin levels of silty and clayey sand are also found. The free surface is situated at -25 m, and the saturated zone is about 40 m thick. From the data described in the previous section and from Fact 1, we estimated the hydraulic gradient in steady-state conditions to have a modulus of 0.005 to 0.006 and a direction from NNW to SSE. Hydraulic conductivity ranged from 10-3 m/sec to 10-4 m/sec (Table 1). From the data shown in Figure 6, we realized that the phreatic surface lies 1 m or more above the piezometric head of the underlying aquifer to which the phreatic one is connected by the aquitard. Since an unsaturated zone lies above the free surface and toxic wastes from human activities are disposed of in the former, we expected the phreatic aquifer to be the most severely polluted. The findings by Airoldi and Plos (4) confirmed the expectation and showed that pollutant concentrations rapidly decrease towards the confined aquifer. Aquitard A continuous aquitard is present under the area at approximately -60 to -65 m. This semipervious layer is mainly composed of silt and silty clay. Its thickness, about 15 m, is fairly constant (Figs. 10,11). Its conductivity, evaluated from pumping tests and given in Table   1 ranges from 10-6 to 10'5 m/sec. The leakage flux across it is directed downward and is a function of the aquitard thickness and conductivity; moreover, the flux depends on the head difference across the two semipervious boundaries. The simulations described below give numerical values for the leakage.

Leaky Confined Aquifer
It consists of sand and gravels, sands, and fine sands. Levels of clay and silt are more abundant here than in the phreatic aquifer. The aquifer's thickness at well sites varies from 20 to 25 m. The aquifer material has been variously described by different authors; however, its transmissivity can be reasonably considered 10  In order to simulate the behavior of the above de-scribed aquifer we introduced a system of coupled partial differential equations, which were discretized and numerically solved. We assumed that Dupuit's hypotheses holds in both aquifers and that leakage flux across the aquitard is only vertical. In other words, we neglected the horizontal flux component across the latter. Solution values for the two aquifers are not significantly affected whenever the ratio between the conductivities of the aquifers and the aquitard conductivity is greater than 100 (10,11). confined aquifers respectively; Sos is the specific storativity of the aquitard and has the dimensions of inverse length; qp and qc represent the flows per unit area extracted from the two aquifers, whereas q. is the flow per unit volume extracted from the aquitard; Dp, Dp, I°and FV are given functions. Remark 1. The aquifers and the aquitard are assumed to be isotropic but nonhomogeneous. Remark 2. The thickness of the phreatic aquifer is given by lp = $pbp (4) where bp represents the height above the mean sea level of the aquifer bottom. As a consequence Eq. (la) is quasilinear. Remark 3. No term representing the accretion rate of the free surface has been included into Eq. (la) because the modeled area is entirely covered by buildings and roads. The only inflow term from above is represented by losses from the sewers and can be ignored. Remark 4. Some analytical solutions to the system of Eqs. (1) to (3) in this section, which correspond to simple geometries, constant coefficients and regular source terms, are described in the literature (11,12). In general, when the domain Q has no simple shape and the coefficients are not constant, no closed form solution to the system exists. A suitable approximation is then introduced that leads to a numerical solution procedure. With reference to the real aquifer we are modeling, the square domain fl is replaced by a lattice with square meshes and a total of 144 nodes; the side (Ax or Ay) of each mesh is 50 m long. The actual positions of wells are approximated by those of the nearest nodes.
The nonlinearity which appears in Eq. (la) is treated as follows. If we consider a one-dimensional homogeneous phreatic aquifer in the steady state condition, without accretion, we obtain: Kp(Op-p d P =q=cont -bP q. =cnt (5) By integrating (5) between x = 0 and x = 1 we obtain: This implies that the flow q through a section of the aquifer is given by: where the term in square brackets represents the mean aquifer thickness. As a consequence in our scheme, the discretized counterparts of lp shall be the arithmetic mean of two adjacent nodal values.
In order to discretize the aquitard along the z direction we assigned four nodes to the [0,l1] interval, which Eq. (2a) refers to. The first and last nodes lie respectively at the boundary with the phreatic aquifer and the confined aquifer. The distance between two nearest nodes, Az, is 5 m. This kind of discretization deals with transient processes through the aquitard in a satisfactory way.
Concerning the discrete scheme, we identified a lattice node by means of the triple (i,j,k), where the indexes correspond to the x, y and z coordinates respectively. Nodes having k = 1 (k = 4) belong to the phreatic (confined) aquifer. The value of the piezometric head at node (i,j,k) at the time t will be denoted by h(i,j,k,t). A similar notation holds for other function values.
Definition 3. The harmonic means K(i,j,k) of hydraulic conductivities between nodes (i,j,k) and (i + 1,j,k), at fixed j,k and along the i direction, reads: K(i, j, k) + K(i + l, j, k) (8) Discontinuous conductivity functions are more effectively dealt with by means of this definition. Definition 4. Let At denote the time step, A,h(i,j,k; t) := h(i + 1, j, k; t) -h(i,j, k; t) (9) denotes the forward increment of h(.,j,k;t) along the i direction, denotes the average of h(.,j,k;t) along the i direction. Remark 5. Definitions 3 and 4 are easily modified to hold for node pairs along the j and k directions.
The above defined nonlinear algebraic system is solved by a modified successive overrelaxation method. The terms in square brackets in Eq. (11) in this section at the m-th iteration are those computed one iteration earlier. A similar approach has been adopted by Gambolati et al. (13).

Model Calibration and Sensitivity Analysis Numerical Values of Relevant Quantities
We now need to assign some realistic numerical values to all functions and parameters mentioned. Let us begin with the boundary conditions. The available data show that the hydraulic gradient in the area we are considering has a modulus of about .005 and is directed toward S-SW. If the wells drew their water from one aquifer only, then measured head values could be easily extrapolated to obtain the above-mentioned Dirichlet boundary conditions for the phreatic and the confined aquifer separately. This is not the case, however (Fact 1). Therefore, we ascribe measured values to the confined aquifer, extrapolate them up to the boundary, and assume the phreatic head is 1.5 m higher than the confined one.
Hydraulic conductivities have been computed from the experimental results mentioned above in the third section (Table 1). Since these values vary from 1.5 x 10-3 to 5.1 x 10-3 m/sec, the aquifer is nonhomogeneous. In order to obtain the nodal values of hydraulic conductivities, we processed data as follows: a) The ratios between conductivity values obtained from pumping tests and specific capacities pertaining to the same wells are seen to be almost constant. b) This "constant ratio" rule is applied to those wells where no pumping test data are available. c) Once conductivities are known at each well, linear interpolation yields the values at the nodes of the 12 x 12 inner lattice. d) Again, because of Fact 1, all of these values are actually weighted averages of different hydrogeological units. According to the previous remarks, since the phreatic aquifer is about twice as thick as the confined one, we assume: K1 =5 (15) Conductivity values at peripheral nodes (Remark 7) are assumed to be constant and equal to the arithmetic means of the inner node values.
According to the results of pumping tests and to the previously mentioned hydrogeological considerations, the following constant values are assigned: SP = 2 i10-2; S, = 2.10-4; K. = 3. 10-6 ma-; So' = 10-6 8-1 (16) Values of specific discharge rates (per unit area) are obtained from well discharge rates after equally dividing them between the phreatic and confined aquifers. This is physically equivalent to having two collocated wells drawing from separated layers. The discharge rate of a given well is assigned to a single lattice node and is considered a concentrated sink term.

Results of Computer Simulations
Several runs have been carried out by means of this model in order to simulate the various exploitation conditions of the resource. A typical CPU time for a steadystate run is about 30 seconds. We have considered the comparison with three sets of experimental data: a) In the steady state with all wells turned off, the computed phreatic and confined heads differ by about 1 m (Figs. 12,13). This agrees with experimental data shown in Figure 6. b) The computed hydraulic gradients fit the measured ones, both by their moduli and directions (Fig. 4). c) Some pumping test data (drawdown vs. time) are shown in Figure 14 together with the numerical results obtained by the model both for the phreatic and confined aquifers. The quality of the fit between computed and measured quantities differs from well to well. However, we consider it satisfactory if we take into account all three previously described approximations.
We do not claim this model's calibration to be the best possible. However, we believe it is difficult, if not useless, to improve calibration at this stage because of the poor quality of the data.

Sensitivity Analysis
In order to evaluate the model calibration errors, we have carried out sensitivity analysis, as it is usually done in hydrogeological modeling (14,15). The sensitivity coefficient uh,G of the discretized piezometric head h at a node with respect to a known parameter (or boundary value G at the same or at another node is defined by: For a detailed account of our tests, we refer the reader to Giudici and Ponzini (6). We summarized herewith the most relevant feature. Dirichlet boundary conditions far away from the area of interest considerably reduce the influence of possible errors in the assumed head values. If a value on the northern boundary is changed by 1   For simplicity, we assume that contaminants are solutes carried by water flowing through the porous me-dium without either diffusing, decaying, chemically reacting, or being adsorbed by clay. In fact, most of these processes have not yet been studied in a satisfactory way, at least in the area we considered. By our model we can compute leakage fluxes under any condition of aquifer exploitation. Given these hypotheses and the quoted experimental results (4) showing that the leaky phreatic aquifer is heavily polluted over the city area, the solute concentration of water extracted by wells is expected to be related to the magnitude of leakage flux which reaches the confined aquifer.
The steady state leakage fluxes have been computed with all wells switched off and on and each of them discharging 30 L/sec when all wells are used as purge wells. A map ofthe ratios between these values is shown in Figure 15. The dimensionless numbers mean that FIGURE 15. Ratios between computed steady-state leakage fluxes corresponding to the conditions of all wells switched on vs. all wells switched off. 91 leakage increases strongly when all wells are turned on, particularly at the center of the domain where the well density is highest. As long as each well draws equal amounts of water from both aquifers, the leakage fluxes affect the resulting water contamination only for a low amount. However, the aquitard actually acts as a filter, adsorbing part of the contaminants coming from the phreatic aquifer. Pollutants must therefore have accumulated in the aquitard for the past 10 years, because all wells were used as purge wells. As a consequence, the semipervious layer is highly polluted, especially where the maximum leakage fluxes take place. This area of the aquitard acts as a contaminant source for the confined aquifer. These results are in complete agreement with Fact 2 and the trends of chlorinated halocarbons concentration contour maps of Figures 7, 8, and 9. Our model explains why maxima of trichloroethylene, tetrachloroethylene, and methyltrichloromethane concentrations occur at the center of the water supply station domain near wells 2, 3, and 13. Leakage flow rates there are three or more times as high as those near the boundary.
These remarks together with the results listed in the preceding section lead us to believe our model to be realistic. This is why we have carried out a number of simulation tests which stress the role of the aquitard in preventing contamination of the confined aquifer. More precisely, the purpose of the following examples is to show the effects of different aquitard conductivity values and pumping patterns on computed leakage fluxes which govern the contaminant flows towards the confined aquifer.
In order to represent our results in a better way, we selected a typical cross section of the domain, which corresponds to the sixth column of the lattice and includes wells 2a, 2, and 3. Example 1. Two different values of aquitard conductivity have been considered; the former ten times higher, and the second ten times lower than the one assumed in the model. Figure 16a represents the piezometric heads of the two aquifers when all wells are switched on and KS = 3 x 10-6 m/sec (the assumed value). shows the low sensitivity of the phreatic aquifer head and the high sensitivity of the confined one w.r. to aquitard conductivity variations. Figure 17 represents the computed leakage fluxes at each node of the said column, corresponding to the three above mentioned values of aquitard conductivity. The leakage fluxes are higher at those nodes where wells are located and reach values of almost 4 L/sec when Ks = 3 x 10-5 m/sec, whereas if Ks = 3 x 10-7 m/sec these values are considerably lower, i.e., about 0.2 L/ sec. We realize that the leakage rate contributes from 0.6 to 13 % ofthe total discharge rate coming from wells; thus, it does not significantly affect the resulting water contamination under the current exploitation scheme.
On the other hand, a more sensible exploitation policy that differentiates between users suggests that they draw water for industrial purposes mainly from wells drilled into the phreatic aquifer and draw drinking water from wells drilled into the confined aquifer, without altering the present structure of the distribution network. In this case, the leakage rate would play an important role in determining the contamination level of the confined aquifer; the need'of careful resource management is thus apparent.
In order to show how critical the value of K8 is in the model, let us look at Figure 18, which shows the total leakage flux, integrated over the area, vs. aquitard conductivity. At low values of K, (10-7 m/sec of less) leakage flux is highly sensitive, more precisely directly proportional to K., since the two aquifers are almost independent and hydraulic gradient across the aquitard does not change significantly. However, leakage flux values are negligible and are not expected to affect contaminant transport. At high values of K8 (10-5 m/sec and higher), sensitivity is lower, but leakage flux values are important especially as contaminant carriers. Namely, in the whole water supply station area, the contaminated flux eventually exceeds 100 L/sec. This is why precisely determining conductivity of the aquitard must be an important target of future experimental tests.
Example 2. Some numerical experiments aimed at showing the effects of different exploitation patterns on leakage fluxes have also been performed. The findings are summarized by Table 2, stressing only the most important results: a) Using all area wells as purge wells causes an extremely strong increase of leakage flux over the area ( Table 2, patterns 1 and 2), in agreement with Figure  8.
b) The attempt to redistribute the wells where hydraulic conductivity is higher ( Table 2, pattern 3) does not significantly decrease the leakage because the aquifer is almost homogeneous. c) Better results can be obtained by drawing water from the phreatic aquifer for industrial purposes and drawing drinking water from the confined one. The actual fraction of water supplied by Milan Water Works to industrial users is unknown. We therefore conservatively assumed that fraction to be two-thirds and distributed the well discharge rates between the aquifers accordingly ( Table 2, pattern 5). The total leakage flux towards the confined aquifer reduces to about 100 L/sec. As expected, better results are obtainable-not so much by increasing the rate drawn from the phreatic aquifer as by reducing the flow drawn from the confined one-because of the difference between the transmissivity values of the two aquifers.
We have also produced strong evidence against the widespread and often unjustified remedy of letting all wells in a contaminated area continuously discharge for a long time-even years-in order to expel pollutants from the aquifer. This method is shown to be either counterproductive or ineffective in the coupled aquifers system we have been dealing with.
Another advantage of our model appears from Example 2, whereby we can determine the optimal pumping pattern and schedule, which minimizes the transport of contaminants towards the confined aquifer, given the needed total discharge rate.