Locate the pollution source 355 Locate the pollution source en quan Yang Zhenyu He xiaofei Zhejiang University Hangzhou, China Advisor: Zhang Chong Summary We develop a model for a strategy to detect new pollution. Three processes govern the movements of pollutants in groundwater: advection, dispersion, and retardation. Information from the wells is used to determine the rate and direction of groundwater movement, determine the horizontal and vertical extent of the pollutants, and analyze the underground structure and characteristics Regarding the diversity and complexity of the given data, we employ a two-step data selection to determine the pollutants most likely to cause new pollution <uring this period of time. We refine the data to choose those chemicals that est represent the variation during this period of time. Then, by using a grid search algorithm, we write a computer program to simulate the movement process and identify the location and start time of the pollution source. The program is written in C and runs on a PC. Four kinds of new pollution sources are located. The graph resulting from our model is in a good agreement with the given data. Finally, we test parameter sensitiv Assumptions All soil and aquifer properties are homogeneous and isotropic throughout both the saturated zone and the unsaturated zone The aquifer consists of sand and gravel The UMAP Journal 20(3)(1999)355-368. @Copyright 1999 by COMAP, Inc. All rights reserved. Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice. Abstracting with credit is permitted, but copyrights for components of this work owned by others than COMAP must be honored. To copy otherwise to republish, to post on servers, or to redistribute to lists requires prior permission from COMAP
Locate the Pollution Source 355 Locate the Pollution Source Shen Quan Yang Zhenyu He Xiaofei Zhejiang University Hangzhou, China Advisor: Zhang Chong Summary We develop a model for a strategy to detect new pollution. Three processes govern the movements of pollutants in groundwater: advection, dispersion, and retardation. Information from the wells is used to • determine the rate and direction of groundwater movement, • determine the horizontal and vertical extent of the pollutants, and • analyze the underground structure and characteristics. Regarding the diversity and complexity of the given data, we employ a two-step data selection to determine the pollutants most likely to cause new pollution during this period of time. We refine the data to choose those chemicals that best represent the variation during this period of time. Then, by using a gridsearch algorithm, we write a computer program to simulate the movement process and identify the location and start time of the pollution source. The program is written in C and runs on a PC. Four kinds of new pollution sources are located. The graph resulting from our model is in a good agreement with the given data. Finally, we test parameter sensitivity. Assumptions • All soil and aquifer properties are homogeneous and isotropic throughout both the saturated zone and the unsaturated zone. • The aquifer consists of sand and gravel. The UMAP Journal 20 (3) (1999) 355–368. c Copyright 1999 by COMAP, Inc. All rights reserved. Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice. Abstracting with credit is permitted, but copyrights for components of this work owned by others than COMAP must be honored. To copy otherwise, to republish, to post on servers, or to redistribute to lists requires prior permission from COMAP
356 The uMAP Journal 20.3(1999) Steady, uniform water flow occurs only in the vertical direction throughout the unsaturated zone, and only in the horizontal (longitudinal) plane in the saturated zone in the direction of groundwater velocity Physical processes play the greatest role, while the chemical processes are negligible All the parameters describing the characteristic of both zones are constant hroughout the monitoring period All the sources of the pollutants are point sources roblem one This problem is to estimate the location and start time of the source, so we consider the movement process of the pollution and the structure of the underground Data Analysis and Processing We assume that there is no interaction between pollutants so that we can process each pollutant separately. With the given data of the coordinate and water level of each well, we plot the water level map by using linear interpo- lation on the elevations of the monitoring wells, as in Figure 1. For simplicity of computation, we assume that all the underground water flows in the same direction o M4 MuD Mwz Direction of under ground water Figure 1. Water-level map
356 The UMAP Journal 20.3 (1999) • Steady, uniform water flow occurs only in the vertical direction throughout the unsaturated zone, and only in the horizontal (longitudinal) plane in the saturated zone in the direction of groundwater velocity. • Physical processes play the greatest role, while the chemical processes are negligible. • All the parameters describing the characteristic of both zones are constant throughout the monitoring period. • All the sources of the pollutants are point sources. Problem One This problem is to estimate the location and start time of the source, so we consider the movement process of the pollution and the structure of the underground. Data Analysis and Processing We assume that there is no interaction between pollutants so that we can process each pollutant separately. With the given data of the coordinate and water level of each well, we plot the water level map by using linear interpolation on the elevations of the monitoring wells, as in Figure 1. For simplicity of computation, we assume that all the underground water flows in the same direction. Figure 1. Water-level map
Locate the pollution source 357 Data selection Because we have thousands of data points about concentration of various pollutants, we must select data carefully. We do this in the following steps Because pollutants are strongly influenced by layers of different permeabil- ity, measurements of critical parameters and pollutant concentrations need to be done at intervals over the depth of the aquifer. We need a method for sampling at different depths in an aquifer. By analyzing the data set, we find that almost every pollutant affects only one part of a well( top, middle, or bottom). Thus, for each pollutant we need to consider only the effect on one layer of the well. Furthermore, the data from the bottom of each well (if any)remain constant or nearly so, hence we can neglect such data. We delete the data for some pollutants, such as tetrachloroethane, acrolein benzene bromomethane chlorobenzene cobalt, and so on because there are hardly any changes in concentrations of these pollutants in each well. We think that the pulse fluctuation about the pollutant concentration dur- ing a relatively stable period, such as for manganese, is caused by random factors. Thus, we eliminate these pollutants from the data set There is a particular constituent, the CarbonTotalOrganic, whose concen tration value decreases significantly, from more than 1000 to less than 1.5 hus, we eliminate it. Now only four pollutants remain: calcium, chloride, magnesium, and TDS Reselection F maining pollutant, to accurately reflect the tendency for the centration to change we reselect its data as follows For each well, we choose two concentration values for each year, one from the first half of the year and the other from the second half Because we do not know the locations of mw-27 and MW-33 and. moreover the concentration changes in these two wells are small, we do not consider their data According to the groundwater flow direction, the average concentration value of MW-9 should not be higher than that of MW-3 and MW-12, which contradicts the given data for calcium, chloride, and so on. This is also true for barium.(In 1997, concentrations in MW3M and MW12M vary from 50 to 85, whereas they vary from 80 to 95 in MW9M )Therefore, we think that MW-9 is a pumping well(Figure 2). Thus, we do not use the data from MW-9 in our analysis Finally, we list in Table 1 the data for calcium that we use to calculate the source location
Locate the Pollution Source 357 Data Selection Because we have thousands of data points about concentration of various pollutants, we must select data carefully. We do this in the following steps: • Because pollutants are strongly influenced by layers of different permeability, measurements of critical parameters and pollutant concentrations need to be done at intervals over the depth of the aquifer. We need a method for sampling at different depths in an aquifer. By analyzing the data set, we find that almost every pollutant affects only one part of a well (top, middle, or bottom). Thus, for each pollutant we need to consider only the effect on one layer of the well. Furthermore, the data from the bottom of each well (if any) remain constant or nearly so, hence we can neglect such data. • We delete the data for some pollutants, such as tetrachloroethane, acrolein, benzene, bromomethane, chlorobenzene, cobalt, and so on, because there are hardly any changes in concentrations of these pollutants in each well. • We think that the pulse fluctuation about the pollutant concentration during a relatively stable period, such as for manganese, is caused by random factors. Thus, we eliminate these pollutants from the data set. • There is a particular constituent, the CarbonTotalOrganic, whose concentration value decreases significantly, from more than 1000 to less than 1.5. Thus, we eliminate it. • Now only four pollutants remain: calcium, chloride, magnesium, and TDS. Reselection For each remaining pollutant, to accurately reflect the tendency for the concentration to change, we reselect its data as follows: • For each well, we choose two concentration values for each year, one from the first half of the year and the other from the second half. • Because we do not know the locations of MW-27 and MW-33 and, moreover, the concentration changes in these two wells are small, we do not consider their data. • According to the groundwater flow direction, the average concentration value of MW-9 should not be higher than that of MW-3 and MW-12, which contradicts the given data for calcium, chloride, and so on. This is also true for barium. (In 1997, concentrations in MW3M and MW12M vary from 50 to 85, whereas they vary from 80 to 95 in MW9M.) Therefore, we think that MW-9 is a pumping well (Figure 2). Thus, we do not use the data from MW-9 in our analysis. Finally, we list in Table 1 the data for calcium that we use to calculate the source location
358 The uMAP Journal 20.3(1999) the direction of the water direction of undergroyndw ater water velocity directs toward ping well the pumping well when near it ure 2. Groundwater movement near a pumping well Table 1 Data for calcium used in the model Date MW-3M MW-7M MW-1IT MW-12M 12/7/93 41 3/7/94 9/19/94 7/10/95 36.5 54.3 59.5 10/10/9519 43.2 3/6/96 507 824 10/9/ 619 3/18/97 125 12/15/97614 115 63.8 884 According to the data, there is some pollutant detected in an early year as 1990; we name it the background concentration(Cb). We think that the later pollutants concentrations consist of background concentration plus new injected concentration. According to Figure 1, MW-1 must be at the headwater level. Moreover, the data from its bottom hardly change during this period according to the data set. Thus, we estimate Cb using data from MW-IB as follows Cb= arithmetic mean of the concentration value from MW-1B during this period for a certain pollutant n Table 2 we collect the symbols used in this paper and their definitions
358 The UMAP Journal 20.3 (1999) Figure 2. Groundwater movement near a pumping well. Table 1. Data for calcium used in the model. Date MW-3M MW-7M MW-11T MW-12M 12/7/93 41 50 39 42 3/7/94 42 50 43 47 9/19/94 42 45 41 41 7/10/95 36.5 54.3 44.7 59.5 10/10/95 19.2 53 43.2 54.7 3/6/96 62.4 65.1 50.7 82.4 10/9/96 60.2 61.9 53.3 87.6 3/18/97 63.8 125 53.2 87.6 12/15/97 61.4 115 63.8 88.4 • According to the data, there is some pollutant detected in an early year such as 1990; we name it the background concentration (Cb). We think that the later pollutants’ concentrations consist of background concentration plus new injected concentration. According to Figure 1, MW-1 must be at the headwater level. Moreover, the data from its bottom hardly change during this period according to the data set. Thus, we estimate Cb using data from MW-1B as follows: Cb = arithmetic mean of the concentration value from MW-1B during this period for a certain pollutant In Table 2 we collect the symbols used in this paper and their definitions
Locate the pollution source 359 Table 2 Symbols used horizontal dispersion coefficient(m) cal dispersion coefficient(m) pollutant concentration (mg/liter) background concentration(described above) Co concentration in the pollutant source(mg/liter pervasion coefficient(m2/s) H water level (ft) hydraulic conductivity(gal/day/ft2) L horizontal distance in the direction of water flow(ft) discharge rate of the pollutant(mg/day) effective porosity discharge rate of the pollutant(liter/day) retardation factor so comp und parameter start time of the pollution (yr) angle between the direction of underground water and the x-axis Va groundwater velocity(ft/day) hantush function pollution source coordinate Model design Model formulation The movement of pollutants consists of advection, dispersion, and retarda- tion. Furthermore, regarding the large scale of the area, the vertical movement is negligible. Thus, movement of pollutant in the soil(saturated and unsatu rated) can be described by the following two-dimensional equation: ac Vaal 0.z2+ ao a2C 8C a2C Rd Model explanation The model equation applies to steady uniform flow. An analytical solu- tion to the equation can be developed for both continuous(step-function) and pulsed inputs of pollutants as boundary conditions. A step function implie the input of a constant concentration pollutant for an infinite amount of time, while a pulse load is a constant concentration input for a finite amount of time The terms"infinite"and"finite"are relative to the time frame of the analysis Ne assume that the pollution source is applied as a step function(continu- ously) with the following boundary conditions (x,y,0)=0,(x,y)≠(0,0) C(0,0,t)=C t)=C(x,土∞,t)=0
Locate the Pollution Source 359 Table 2. Symbols used. αL horizontal dispersion coefficient (m) αT vertical dispersion coefficient (m) C pollutant concentration (mg/liter) Cb background concentration (described above) C0 concentration in the pollutant source (mg/liter) D pervasion coefficient (m2/s) H water level (ft) I hydraulic gradient K hydraulic conductivity (gal/day/ft2) L horizontal distance in the direction of water flow (ft) m discharge rate of the pollutant (mg/day) n effective porosity q discharge rate of the pollutant (liter/day) Rd retardation factor S compound parameter t0 start time of the pollution (yr) θ angle between the direction of underground water and the x-axis Vd groundwater velocity (ft/day) W hantush function (x0, y0) pollution source coordinate Model Design Model Formulation The movement of pollutants consists of advection, dispersion, and retardation. Furthermore, regarding the large scale of the area, the vertical movement is negligible. Thus, movement of pollutant in the soil (saturated and unsaturated) can be described by the following two-dimensional equation: Rd ∂C ∂t = VdαL ∂2C ∂x2 + Vdαt ∂2C ∂y2 − Vd ∂C ∂x . (1) Model Explanation The model equation applies to steady uniform flow. An analytical solution to the equation can be developed for both continuous (step-function) and pulsed inputs of pollutants as boundary conditions. A step function implies the input of a constant concentration pollutant for an infinite amount of time, while a pulse load is a constant concentration input for a finite amount of time. The terms “infinite” and “finite” are relative to the time frame of the analysis. We assume that the pollution source is applied as a step function (continuously) with the following boundary conditions: C(x, y, 0) = 0, (x, y) = (0, 0); C(0, 0, t) = C0; C(±∞, y, t) = C(x, ±∞, t)=0, t ≥ 0
360 The uMAP Journal 20.3(1999) Model solution The function is a second-order partial differential equation. Equations of this form apply to a wide variety of problems, including mass transport, fluid dynamics, and heat transfer For an instantaneous point source at time t=0, there is an analytical solt tion of the form C(z, 3, t)=Sexp(2ar )w(o, b)-W(t, b)] (2) q 4TVa(aLat) and w(u, b)is the hantush function xp-y dy with b y 402 alOT sin Before computing, we classified the parameter used according to our as- ove are unknown, and so is the value of m. Thus, ro, yo, to, and S are variable o During the data processing, the coordination and time of the pollution soure The parameters aL, ar, 0, and Va are constants The main task is to find the location and the start time of the pollutants. Hence we develop a grid-search optimization routine to get an optimized solution: We estimate the location of the pollutant source and transform coordinates as follows: Set the point of the pollutant source to be the new origin Set the new c-axis to be parallel to the direction of the underground water flow Set the new y-axis to be perpendicular to the new -axis We construct an equation to calculate the movement of the pollutant un der the ground. We calculate the concentration changes in each well and compare with the changes according to the data set. We repeatedly adjust the location of the pollution source, the value S, and the value to(detailed in the following)until there is a satisfactory agreement. The criterion for convergence is the sum of the squares of the residuals between the data and the model predictions. The objective function to be minimized ∑[C-Cb)-C where Ci is the pollutant concentration data value for well i, Ci is the model prediction, and C is the background pollution level
360 The UMAP Journal 20.3 (1999) Model Solution The function is a second–order partial differential equation. Equations of this form apply to a wide variety of problems, including mass transport, fluid dynamics, and heat transfer. For an instantaneous point source at time t = 0, there is an analytical solution of the form C(x, y, t) = S exp x 2αL [W(0, b) − W(t, b)] , (2) where m = C0q, S = m 4πVd (αLαT ) 1/2 , and W(u, b) is the hantush function W(u, b) = ∞ u exp −y − b2 2y y dy with b = x2 4α2 L + y2 4αLαT . Before computing, we classified the parameter used according to our assumptions above: • During the data processing, the coordination and time of the pollution source are unknown, and so is the value of m. Thus, x0, y0, t0, and S are variable. • The parameters αL, αT , θ, and Vd are constants. The main task is to find the location and the start time of the pollutants. Hence, we develop a grid-search optimization routine to get an optimized solution: • We estimate the location of the pollutant source and transform coordinates as follows: – Set the point of the pollutant source to be the new origin. – Set the new x-axis to be parallel to the direction of the underground water flow. – Set the new y-axis to be perpendicular to the new x-axis. • We construct an equation to calculate the movement of the pollutant under the ground. We calculate the concentration changes in each well and compare with the changes according to the data set. We repeatedly adjust the location of the pollution source, the value S, and the value t0 (detailed in the following) until there is a satisfactory agreement. The criterion for convergence is the sum of the squares of the residuals between the data and the model predictions. The objective function to be minimized is i [(Ci − Cb) − C i] 2 , where Ci is the pollutant concentration data value for well i, C i is the model prediction, and Cb is the background pollution level.
Locate the pollution source 361 Parameter estimation We estimate the parameters for the saturated zone as following: Hydraulic Conductivity K: We consider hydraulic conductivity, measured in gallons per day per square foot, only in thehorizontal direction. According to the literature, K=265 gPd/ft(1 gpd/ft2=4.72 x 10-5 cm/sec) Hydraulic Gradient: According to Figure 1, made by interpolation, we assume that the direction of the underground water is one-dimensional Ground-Water (Interstitial Pore Water) Velocity Vd: According to Darcys Law, va is defined as Vd=-ki/n where I is the hydraulic gradient, K is hydraulic conductivity, and n is effective porosity. We assume that the soil type of the saturated zone is sand with porosity 20%0, so we estimate Va=1.5 ft/da Dispersion Coefficient o: This coefficient incorporates two forms of disper- sive process: dynamic dispersion and molecular diffusion. According to the literature the horizontal dispersion coefficient and the vertical dispersion coefficient are approximately equal both have the estimated value 25 ft Retardation Factor: Retardation is based on pollutant characteristics and aquifer composition. Since its effect is not very significant, we estimate Ra=1. Concentration in Pollution Source: According to the literature when the water table is usually sufficiently high so that the pollutant directly enters ground water, the Co value is the estimate of the source concentration. Results There are four new pollutants: calcium, chloride, magnesium, and TDS The location and the start time for the pollution sources as predicted by our model, are in table 3 Table 3 Source and start time for pollutants Pollutant a-coordinate y-coordinate Start time (m/d/y) TDS 7077 6538 8/12/91 6423 1/1/ floride 6931 5/18/91 Calcium 6040 9/1/93 Finally, we mimic the movement process of the pollutant in reverse and compare with the given data set( Figure 3). From the graphs, we conclude
Locate the Pollution Source 361 Parameter Estimation We estimate the parameters for the saturated zone as following: • Hydraulic Conductivity K: We consider hydraulic conductivity, measured in gallons per day per square foot, onlyin the horizontal direction. According to the literature, K = 265 gpd/ft2 (1 gpd/ft2 = 4.72 × 10−5 cm/sec). • Hydraulic Gradient: According to Figure 1, made by interpolation, we assume that the direction of the underground water is one-dimensional. • Ground-Water (Interstitial Pore Water) Velocity Vd: According to Darcy’s Law, Vd is defined as Vd = −KI/n, where I is the hydraulic gradient, K is hydraulic conductivity, and n is effective porosity. We assume that the soil type of the saturated zone is sand with porosity 20%, so we estimate Vd = 1.5 ft/day. • Dispersion Coefficient α: This coefficient incorporates two forms of dispersive process: dynamic dispersion and molecular diffusion. According to the literature, the horizontal dispersion coefficient and the vertical dispersion coefficient are approximately equal. Both have the estimated value 25 ft. • Retardation Factor: Retardation is based on pollutant characteristics and aquifer composition. Since its effect is not very significant, we estimate Rd = 1. • Concentration in Pollution Source: According to the literature, when the water table is usually sufficiently high so that the pollutant directly enters ground water, the C0 value is the estimate of the source concentration. Results There are four new pollutants: calcium, chloride, magnesium, and TDS. The location and the start time for the pollution sources, as predicted by our model, are in Table 3. Table 3. Source and start time for pollutants. Pollutant x-coordinate y-coordinate Start time (ft) (ft) (m/d/y) TDS 7077 6538 8/12/91 Magnesium 6423 7461 1/1/94 Chloride 6931 5823 5/18/91 Calcium 7750 6040 9/1/93 Finally, we mimic the movement process of the pollutant in reverse and compare with the given data set (Figure 3). From the graphs, we conclude:
362 The uMAP Journal 20.3(1999) Well3 at Middle Well 2 at Middle t: year Well ll at Top Well l2 at Middle ba Figure 3. Calcium concentrations at four wells. The thick curves are data, the thin curves are model predictic Fornear-ideal conditions, the modelis suitable; forregular use, a more robust model is desired Even though the two curves do not fit very well, they show a similar change tendency Sensitivity analysis We conduct a rudimentary sensitivity analysis to explain the stability of our model. We separately vary the values of the constants aL, aT, 0, and Va by 10% and compute the corresponding changes in the values of the location and time of the pollution source table 4 The model demonstrates good stability, but 0 has a relatively significant influence on the result of the model. Thus, it is reasonable to consider the parameter 0 as a variable and repeat our grid-searching algorithm in a five- dimensional space of 0, o, yo, to, and S. For calcium, we get the comparative results shown in table 5 In the expanded model, the value for 0 is 7% larger. We think that there ome deflection of the direction of the groundwater flow, as shown in Figure 4
362 The UMAP Journal 20.3 (1999) Figure 3. Calcium concentrations at four wells. The thick curves are data, the thin curves are model predictions. • For near-ideal conditions, themodelis suitable; for regular use, amore robust model is desired. • Even though the two curves do not fit very well, they show a similar change tendency. Sensitivity Analysis We conduct a rudimentary sensitivity analysis to explain the stability of our model. We separately vary the values of the constants αL, αT , θ, and Vd by 10% and compute the corresponding changes in the values of the location and time of the pollution source (Table 4) The model demonstrates good stability, but θ has a relatively significant influence on the result of the model. Thus, it is reasonable to consider the parameter θ as a variable and repeat our grid-searching algorithm in a fivedimensional space of θ, x0, y0, t0, and S. For calcium, we get the comparative results shown in Table 5. In the expanded model, the value for θ is 7% larger. We think that there is some deflection of the direction of the groundwater flow, as shown in Figure 4
Locate the pollution source 363 Table 4 Effects of perturbations of the parameter values Parameter Change in location Change in time CL <10 Table 5 Comparison of 4- and 5-dimensional models ension to 07857750604093.752.1 0847750610093.602.2 6000Mw1 Real direction aMw13 Average direction Brw Mow 3000 Mw14 50006000700080009000 Figure 4. Suspected deflection of groundwater flor Problem two Local assumptions The storage tanks are located underground in the saturated zone The direction of the groundwater flow remains the same The saturated zone is semi-infinite The leak process is continuous, since the primary cause of leaks in steel underground storage systems is corrosion
Locate the Pollution Source 363 Table 4. Effects of perturbations of the parameter values. Parameter Change in location Change in time (ft) (yr) θ 70 0.2 αL <10 <0.1 αT 10 <0.1 Vd <10 <0.1 Table 5. Comparison of 4- and 5-dimensional models. Dimension θ x0 y0 t0 S (ft) (ft) (yr) ×106 4 0.785 7750 6040 93.75 2.1 5 0.84 7750 6100 93.60 2.2 Figure 4. Suspected deflection of groundwater flow. Problem Two Local Assumptions • The storage tanks are located underground in the saturated zone. • The direction of the groundwater flow remains the same. • The saturated zone is semi-infinite. • The leak process is continuous, since the primary cause of leaks in steel underground storage systems is corrosion
364 The uMAP Journal 20.3(1999) Model design To detect the pollutant rapidly and accurately, we develop a three-step method 1. according to the shape and size of the storage and the direction of the groundwater flow, we determine the number and location of the first group of wells. Provided that the storage is a square S m on a side, the number of the first group of wells is N=$/20. That is, we drill a well every m in a line perpendicular to the direction of the groundwater, as shown in Figure 5 We monitor the data from the wells Source Point Q Groundwater flow Direction *x represent storage tank well Figure 5. Locations of monitoring wells. Empty circles represent the initial monitoring wells; filled circles are the wells drilled after pollution is detected and found to affect well 3 most of all 2. Once there is some evidence of pollution, we determine which well is most affected by the pollutant. Near this well, we drill a series of wells(perhaps five or more)along the direction of the groundwater flow. Thus, we can construct a three-dimensional formulation to calculate the fluctuation of the pollutant concentration. Here the area occupied by the storage facility may not be very large(with side less than 1000 ft), so we cannot use(1).We employ the three-dimensional equation D Because the leaking is a continuous process, we assume that the pollution
364 The UMAP Journal 20.3 (1999) Model Design To detect the pollutant rapidly and accurately, we develop a three-step method. 1. According to the shape and size of the storage and the direction of the groundwater flow, we determine the number and location of the first group of wells. Provided that the storage is a square S m on a side, the number of the first group of wells is N = S/20. That is, we drill a well every 20 m in a line perpendicular to the direction of the groundwater, as shown in Figure 5. We monitor the data from the wells. Figure 5. Locations of monitoring wells. Empty circles represent the initial monitoring wells; filled circles are the wells drilled after pollution is detected and found to affect well 3 most of all. 2. Once there is some evidence of pollution, we determine which well is most affected by the pollutant. Near this well, we drill a series of wells (perhaps five or more) along the direction of the groundwater flow. Thus, we can construct a three-dimensional formulation to calculate the fluctuation of the pollutant concentration. Here the area occupied by the storage facility may not be very large (with side less than 1000 ft), so we cannot use (1). We employ the three-dimensional equation Rd ∂C ∂t + Vd ∂C ∂x = D ∂2C ∂x2 + ∂2C ∂y2 + ∂2C ∂z2 + m n . Because the leaking is a continuous process, we assume that the pollution