International Journal of Environmental Research University of Tehran ISSN: 1735-6865 EISSN: 2008-2304 Vol. 3, Num. 2, 2009, pp. 317-326

International Journal of Environmental Research, Vol. 3, No. 2, Spring 2009, pp.317-326

Inverse Method to Estimate the Mass of Contamination Source by Comparing Analytical with Numerical Results

Ardestani, M. and Sabahi, M.S.^*

Graduate Faculty of Environment, University of Tehran, P. O. Box 14155-6135, Tehran, Iran

*Corresponding author E-mail:ssabahi@ut.ac.ir

Received 19 March 2008; Revised 15 Dec. 2008; Accepted 25 Dec. 2008

Code Number: er09034

ABSTRACT

The development of a conceptual model that best accounts for the different parameters influencing LNAPL fate and transport in groundwater is ultimate the key to a successful simulation of LNAPL concentration in groundwater. Characterization of hydrocarbon sources and identification of areas with heavy LNAPL loadings from point and non-point sources is important for land use planners and environmental regulators.Bistoon petrochemical site has discharged the light nonaqueous phase liquid (LNAPL) contamination for 10 months. Since the amount of contamination in this source is not clear and the only available data is the amount of contamination in the observed well, an attempt was made to solve this problem by using inverse method. The amount of source contaminant was found through observed data in the well. In this method the analytical results with numerical ones were compared. Upon computation of the contamination concentration, the mass of contaminants can be calculated by multiplying concentration by volume. The numerical method used is finite difference which is an engine in Modflow program.

Key words: Contamination, NAPLs, Benzene, Modeling, Numerical method, Analytical method

INTRODUCTION

Crude oil, refined petroleum products, as well as polycyclic aromatic hydrocarbons are ubiquitous in various environmental compartments (Onwurah, et al., 2007). Groundwater contamination by hazardous substance is commonly the result of accidental spill that occur during production, storage and transportation activities. The groundwater contamination constitutes of organic and inorganic contaminants. The most common classes of organic groundwater include aromatic hydro-carbons, chlorinated solvents and pesticides (Figs. 1, 2, 3). Common inorganic pollutant include nitrate (NO₃ ⁻ ), arsenic (As ⁺ ), selenium (Se ⁺ ), and toxic heavy metals such as lead (Pb ⁺² ), cadmium (Cd ⁺⁴ ), and chromium (Cr ⁺⁶ ). Petroleum hydrocarbons comprise a diverse group of compounds. Gasoline, which is a very common groundwater pollutant, is also a complex mixture. Gasoline constitutes are generally isopentane, p-xylene, n-propylbenzen, 2,3-dimethylbutane, n-butane, n-pentane, and toluene, which together make up over 50% of the mixture. Shahidi Bonjar (2007) reported about adverse inhibitory effect of V-Guard and E-Guard gasoline additives against soil beneficial Streptomyces. He concluded that gasoline additives contaminate soil and groundwater by fuel leaks and spills.

The most important constitute is benzene, which is highly soluble and as a result highly mobile in groundwater aquifers. The drinking water standard for benzene (5 μg / L ) is much more stringent than that of other mono-aromatic hydrocarbon such as toluene (1000μg / L ) and xylenes (10,000). Another important group of pollutants are the poly-nuclear aromatic hydrocarbons (PAHs), which are commonly found near coal conversion facilities and petroleum refineries.

These hydrophobic pollutants are of major concern to both public and environmental health because of their tendency to concentrate in food chain (Mcelory et al., 1989) and acute toxicity (Heitkamp and Cerniglia, 1988). PAHs are the principle constituents are creosote, which is a complex mixture of about 200 compounds also containing phenolic and heterocyclic pollutants. Hydrocarbons that exist as a separate, immiscible phase comes in contact with water and/or air – the so-called Non-aqueous phase liquids (NAPLs). Non-aqueous phase liquids (NAPLs) are differences in the physical and chemical properties of water and NAPL result in the formation of a physical interface between the liquids which prevent the two fluids from mixing. Non-aqueous phase liquids are typically classified as either light non-aqueous phase liquids (LNAPLs) which have densities less than water, or dense non-aqueous phase liquids (DNAPLs) which have densities greater than water. (Newell et al., 1997).

The development of a conceptual model that best accounts for the different parameters influencing LNAPL fate and transport in groundwater is ultimately the key to a successful simulation of LNAPL concentration in groundwater. Characterization of hydrocarbon sources and identification of areas with heavy LNAPL loadings from point and non-point sources is important for land use planners, environmental regulators. This is essential for developing fate and transport models. Accurate quantification of LNAPL leaching to groundwater is difficult due to the complex interaction between land use practices, on-ground LNAPL loadeing, groundwater recharge, and spil characteristics. THe objective of this paper is to develop a conceptual modeling framework that integrates the on-ground LNAPL loadings , and LNAPL fate and transport in groundwater. Such modeling framework accounts for constant concentration sources of LNAPL. This integration is of great importance to realistically account fo the different processes that LNAPL undergoes and in order to arrive at rational estimates of LNAPL concentrations in groundwater.

MATERIALS & METHODS

Bistoon is an area located in Kermanshah in westerd of Iran. Linear Alkyl Benzene (LAB) is one of the petrochemical productions that are produced in Bistoon plant. \the site is trapezoid with 63600m² area (Technical Reprt, Bistoon Petrochemical, 2006). This plant was build in 2005. There are some evidences implying that due to some problems, this plant has been the main source of contaminates for groundwater. Observations show that the period of leaching LNAPL in groundwater from plant is 10 months. Laboratory Observations also suggest that there are some kinds of LNAPL exist in the groundwater. However it has been mentioned already since toxicity of the benzene is higher than other mixture, this study focused on benzene as an indicator. Based on laboratory report benzene constitutes 2 percent of the mixture therefore the value that will be mentioned for concentration is only for benzene and not for the mixture. Most of the soils in the study area are categorized as well drained. The ground water depth is around 80 meters in the study area. Average annual precipitation is 448 mm. Coefficient infiltration to the aquifer is 0.6. The groundwater flow is from north east to south west (EIA Reports, 2005).

Fig. 4. depicts a pictorial representation of the proposed framework for modeling the impact of land use on LNAPL contamination of groundwater. The framework is a simplification to the itemized conception. The framework also incorporates the identification of the spatial distribution of the on-ground LNAPL sources and corresponding loadings, and the modeling of the groundwater flow system and the LNAPL fate and transport processes. The modeling framework relies on MODFLOW (Harbaugh and McDonald, 1996) for the simulation of the groundwater flow model and on RT3D (Clement, 1998) in order to simulate the LNAPL fate and transport processes in groundwater. Both, MODFLOW and RT3D utilize the same finite-difference grid. There are different mechanisms that control the fate and transport of LNAPL in the soil and groundwater zones. In addition, the models utilized in the framework of different levels of complexities among them. However it is very important to keep in mind that these models and components to be executed in a sequential manner as depicted in (Fig.4). These linkages and procedures are summarized as follows:

A- On ground LNAPL loading is computed and all surface losses of LNAPL are accounted for. Subsequently the net loading is considered as net input of LNAPL to the soil zone.

B- This LNAPL input is subject to different transformations in the soil zone. The LNAPL that leaches to the groundwater is dealt with spatially where it is written into the sink and source mixing (SSM) input file (supported by RT3D). SSM file contains all the information about the concentrations (or mass distribution) of LNAPL that enter the aquifer corresponding to the point and non-point sources. Generally, the SSM file is immense in terms of its size. This is due to the fact that at each time step can be specified for each cell of the model domain the corresponding mass of LNAPL leaching to groundwater.To assure efficiency in framework implementation, it is appropriate to develop the SSM file automatically. This feature would enhance theimplementation of the modeling framework in highly distributed fine resolution situations.

C- Once the SSM file is developed automatically; RT3D can be executed after the development of the necessary input files.

Since this modeling framework is intended for regional analysis, fine time steps were avoided. Monthly time steps were considered for the RT3D models while MODFLOW was developed under steady-state conditions and consequent assumptions. A transient groundwater flow model would capture more adequately many of the dynamics of the drivers that ultimately dictate the transport of LNAPL in groundwater and hence the concentration distribution. Nevertheless, a transient model would imply that all input data to be variable with time which was unavailable for the study under consideration.The first step to develop the modeling framework is to characterize of the spatial variability of the on-ground LNAPL sources and to estimate the corresponding spatial and temporal distribution of the LNAPL loadings. Fig. 5.depicts Schematic view of part of the petrochemical site. The main source of the LNAPL contaminant related to the manholes and pipe lines that were shown in (Fig. 5).

The partial differential equation that governs the three-dimensional transport of a single chemical constituent in groundwater, considering advection, dispersion, fluid sinks/sources, equilibriumcontrolled sorption, and first-order irreversible rate reactions is described as follow; (Zheng and Bennett, 1995):

Where C is the dissolved concentration (ML ⁻³ );is the adsorbed concentration (MM ⁻¹ ); t is time (T); D_ij is the hydrodynamic dispersion coefficient tensor (L ² T ⁻¹ );V_i is thepore water velocity (LT ⁻¹ ); q_s is the volumetric flow rate per unit volume of aquifer and represents fluid sources and sinks (T ⁻¹ ); C_s is the concentration of the fluid source or sink flux (ML ⁻³ ); λ is the reaction rate constant (T ⁻¹ ); R is the retardation factor (L⁰); ρ_b is the bulk density of the porous medium (ML); and is the porosity (L⁰).As can be concluded from Eq. (1), the LNAPL fate and transport model requires the velocity of the groundwater flow. It is necessary to develop a groundwater flow model to obtain the velocity field. The following governing equation of the threedimensional groundwater flow has to be solved and the head distribution and subsequently the velocity are obtained and computed (Schwartz and Zhang, 2003):

Where K _xx, K_yy and K _zz are values of hydraulic conductivity (LT −1 ) along x, y, and z coordinate axes; h is the hydraulic head, W is a flux term that

accounts for pumping (LT), recharge, or other sources and sinks; Ss is the specific storage (L); and t is time (T). The solution to Eq. (2) provides a transient prediction of hydraulic head in a threedimensional domain for an (4)anisotropic hydraulicconductivity field (Schwartz and Zhang, 2003).Coefficients of hydrodynamic dispersion are given by the following equations:

D_L = α_Lv_i+ D^* (3)

D_T = α_Tv_i+ D^* (4)

Where D_L and D_T are the longitudinal and transverse hydrodynamic dispersion coefficients, respectively; D^* is the effective diffusion

coefficient (L² T⁻¹ ); and α_L and α_T are the longitudinal and transverse dispersive ties (L), respectively. In this model the effective diffusion coefficient is omitted due to high amount of Peclet number. The Peclet number is defined as follows (Perkins and Johnson, 1963):

Where v _x is the average linear velocity, d₀ is the mean grain diameter, and D^* is the diffusion coefficient. As the Peclet number increases, advective-dispersive processes (first term in equation 3 and 4) become increasingly more important than diffusion.There are two ways to solve equation 1: Numerical solution andAnalytical solution. The analytical solution was used in this paper is multidimensional transport from a finite, planar source of contaminant under transient conditions (Newell et al., 1996):

Where C= contaminant source, C₀= initial contaminant concentration at the source, x= distance down gradient of source, y= distance from the centerline of the source, z= vertical distance from the groundwater surface to the measurement point, Y= source withes, Z= source depth, α_x = longitudinal dipersivity, α y = horizontal transverse dispersivity, α_z = vertical transverse dispersivity λ = site specific first-order decay coefficient, t= time, v_c = contaminant velocity in groundwater, erf(x)= error function and erfc(x)= complementary error function= 1-erf(x) By solving equation 6, the amount of contaminant in the source zone can be estimated.The surficial geology of the study area shows that there is one major geologic layer. This layer is composed mainly of stratified sand and gravel. Typically, this layer is more than 80 m thick. Also the type of boundary conditions and the flow pattern are assumed as shown in (Table 1 & Fig. 6).

Table 1. Description of boundary conditions used in developing of the groundwater flow

The model domain was uniformly discretized into a finite-difference grid of 67×57 m² cells. Model calibration was conducted by altering the transmissivity values until the simulated potentiometric heads matched closely the observed value.The LNAPL leakage sources in the study are from the manholes as shown in Fig.5. As mentioned earlier, RT3D was used in developing the LNAPL fate and transport model using the same finite-difference grid as in MODFLOW. RT3D interfaces directly with MODFLOW. It retrieves the saturated thickness for each cell, fluxes across cell interfaces in three principle directions, and the locations of flow rates of the various sources and sinks. For temporal discretization, each year was divided into 12 stress periods where each period corresponds to one month during which all inputs are constant.In RT3D ‘Instantaneous Aerobic Decay of C₆H₆ Module’ is selected. In this module instantaneous reaction model is used to simulate aerobic degradation of C₆H₆. This module simulates the instantaneous biodegradation of fuel hydrocarbons under aerobic conditions. The overall aerobic reaction stoichiometry for a fuel hydrocarbon compound can be written as:

C₆ H₆ + 7.5O₂ → 6CO₂ + 3H₂O (7)

Equation (1) is solved both for benzene and oxygen simultaneously. At each time step, an instantaneous reaction algorithm is used to model the removal rates. According to this algorithm, either hydrocarbon or oxygen (which is limiting) will be reduced to zero within a grid cell, after a reaction time step. The reaction algorithm is expressed by (Rifai et al., 1990):

O (t+1) = O (t) -H (t).F,
H (t+1) = 0 when O (t) > H (t). F (9)

Where t refers to a particular time step and F is the stoichiometric ratio.

One of important factors in groundwater modeling is calibration. The purpose of calibrating the mathematical model of LNAPL fate and transport is to update the critical input parameters such that the simulated LNAPL concentrations are in close agreement with the field observed concentrations (Zheng and Bennett, 1995).

The most critical and uncertain parameters should be calibrated. Table 2 input parameters Since LNAPL contaminations occur after plant establishment initial concentration was assumed zero. Therefore only the longitudinal disprsivity is used for calibration. (Fig.7). depicts the algorithm that is used for solution in this study. The first step to solve equation 6 is to determine parameters constitute the equation. λ and v_c are the unknown(Table 3). Prior to determine v_c the value of retardation coefficient (R) should be known. Assuming the sorption can be described by linear partitioning between the water and the aquifer matrix, the retardation coefficient is given by equation 10 (Freeze and Cherry 1979):

Where ρ_b is the bulk density aquifer material, K_d is the distribution coefficient between groundwater and aquifer, and φ is the total porosityThe distribution coefficient (Kd) can be estimated on the soil adsorption coefficient for soil organic carbon (K_oc) and the fraction of soil organic carbon (ƒ_oc ) as in equation 11:

^Kd ^=ƒoc ^Koc (11)

The value of K_oc for benzene is 83.2 and the value is assigned for ƒ_oc is about 0.015 according to Bistoon EIA report (2003). It should be mentioned that this value is in the range of reports by Domenico and Schwartz (1998). By substituting the available data in equation 10 and 11, the results would be concluded: K_d is equal to 1.25 (Cm³/g) and. the value of R is about 7. Thus in this aquifer, benzene would travel 7 times slower than groundwater velocity which v_c is 0.3 (m/d). ë is both chemical-specific and site-specific, and is dependent on such site-specific conditions as electron-acceptor concentrations, soil chemistry, geology, water temperature, microbial populations, and concentrations of other constituents. Due to the great amount of uncertainty and difficulty in quantifying these conditions, the first-order decay rate is generally selected as a conservative literature value or simply neglected. (Nevin et al., 1997). The amount of set to 0.

RESULTS & DISCUSSION

The amount of the value of C₀ is used for numerical method obtained from the mean value of C₀ in well # 3 and C₀ in well # 100 out of analytical method by using equation 6. Therefore two values for well # 3 and well # 100 are available. By averaging these values a C₀for numerical method will be reached. Using this C₀ in Modflow and the result of the numerical will be compared with observation data. Because of assuming that some value is zero during the analytical method, the difference between numerical and observation will be occurring. This process is performed until the difference between numerical results and observation results is minus. (Fig. 8 & 9) depict the final results of the solution. Since in numerical results the decay coefficient is used by entering the amount of oxygen concentration the value of numerical results is estimated lower than analytical results. However in well # 100 the condition is different. In this well the value of C₀ for numerical results is larger than real value. This is because in this well the contamination takes a longer time than well # 3 to reach. Therefore all the concentration values are lower than well # 3. By averaging the two C₀ the value for well # 100 is larger than real one. However this difference is compensating with decay coefficient. Then in well # 100 the difference between observation data and numerical result is not as large as in well # 3.

1. Alvarez, P. J. J. and Illman, W.A. (2006). Bioremediation and Natural Attenuation. WILEY, INTERSCIENCE. 5-6.
2. Bistoon Petrochemical complex (2005). Technical reports. BPC, Kermanshah, Iran.
3. Charles J. Newell, Steven D. Acree, Randall R. Ross, and Scott G. Huling, (1997), Light Nonaqueous Phase Liquids, EPA Groundwater Issue, (EPA/540/S-95/500) Technology Innovation Office, Washington, DC.
4. Clement (1998). RT3D: A Modular Computer Code for Simulating Reactive Multispecies Transport in 3-Dimensional Groundwater Systems. PNNL-SA-11720
5. Domenico, P.A. and Schawartz, F.W. (1998). Physical and Chemical Hydrogeology. 2nd ed. John Wiley & Sons, Hoboken, NJ.
6. Freeze, R.A. and Cherry, J.A. (1979). Ground-water. Prentice Hall, Upper Saddle River, NJ.
7. Harbaugh, A.W. and McDonald, M.G. (1996). User’s documentation for MODFLOW-96, an update to the US Geological Survey modular finite-difference groundwater flow model: US Geological Survey Open-File Report 96-485, 56.
8. Heitkamp, M.A. and Cerniglia, C.E (1998). Mineralization of polycyclic aromatic hydrocarbon by a bacterium isolated from sediments below an oil field. Appl. Environ. Microbial. 54,1612-1614
9. Keith, L.H. and Telliard, W.A. (1979). Priority pollutants I- a perspective view. Environ. Sci. Technol. 13, 416-423.
10. Mashaldaran Petrochemical Company (2003). EIA report of Kermanshah LAB project. MPCo, Kermanshah, Iran.
11. McElroy,A. E., Farringtone, J.W. and Teal, J.M. (1989). Bioavailability of polycyclic aromatic hydrocarbons in the aquatic environment. In Metabolism of Polycyclic Aromatic Hydrocarbons in the Aquatic Environment, U. Varanasi (Ed.). CRC Press, Boca Raton, FL.
12. Nevin, J. P., Connor, J. A., Newell, C. J., Gustafson, J. B. and Lyons, K. A. (1997). FATE 5: A Natural Attenuation Calibration Tool for Groundwater Fate and Transport Modeling. NGWA Petroleum Hydrocarbons Conference, Houston, Texas.
13. Newell, C.J., Mcleod, R. K. and Gonzales, J. R. (1996). BIOSCREEN Natural Attenuation Decision Support System, User’ Manual Version 1.3, EPA/600/R-96/087, United States Environmental Protection Agency, Office of Research and Development Washington, DC.
14. Onwurah, I. N. E.; Ogugua, V. N.; Onyike, N. B.; Ochonogor, A. E. & Otitoju, O. F. (2007). Crude Oil Spills in the Environment, Effects and Some Innovative Cleanup Biotechnologies. Int. J. Environ. Res., 1(4),307-320.
15. Perkins, T.K. and Johnson. O.C. (1963). A review of diffusion and dispersion in porous media. J. Soc. Peteol. Eng., 3,70-83
16. Rifai, S. H. and Bedient, P. B. (1990). Comparison of biodegradation kinetics with an instantaneous reaction model for groundwater. Water Resources Res., 26(4), 637- 645.
17. Schwartz, F. and Zhang, H. (2003). Fundamentals of Ground Water. John Wiley and Sons, ISBN: 9780471137856.
18. Shahidi Bonjar, G. H. (2007). Potential Hazards of Gasoline Additives in Altering Soil Environment in Favor of Harmful Microorganisms. Int. J. Environ. Res., 1(1), 1- 4.
19. Zheng, C. and Bennett, G. D. (1995). Applied Contaminant Transport Modeling: Theory and Practice. Van Nostrand Reinhold, NewYork, 440.

Copyright 2009 - International Journal of Environmental Research

The following images related to this document are available:

Photo images