ABSTRACT

Currently, the mathematical code Modflow is widely used to simulate groundwater flow in aquifers. Due to the ease which exists today to create mathematical models through Modflow visual Interfaces, it is possible to obtain contaminant transport results which may not have much support, especially when simulating the transport of contaminants with little groundwater flow information. Domenico’s equation is an analytical solution for transport of contaminants in groundwater that can be used when not much groundwater flow information exists. The objective of this study is to model, using Domenico’s equation, a groundwater contaminant plume that discharges into a tidally influenced river. The study area was a wood treatment facility located on the bank of a river which is influenced by tides. Previous studies have found the presence of creosote in the subsurface and the formation of a groundwater plume that apparently discharges into the river. Domenico’s equation was selected to model this site because of the limited piezometric data available at the site to properly simulate the daily hydraulic gradient inversion due to the river tides. Domenico’s equation was successfully used to model this plume and reproduce the field distribution of naphthalene, benzene and 1- methyl-naphthalene. Two sources 40 m inland had to be defined to properly simulate the plume behavior. It was determined through modeling that biodegradation plays an important role on the plume’s behavior. These were key issues in the conceptual model understanding of the plume at this site.

1. Introduction

Currently, the mathematical code Modflow is widely used to simulate groundwater flow in aquifers [1]. Modflow is a powerful tool with several modules that help in modeling increasingly complex scenarios [2]. Among its features, it is able to model the movement of contaminants in the aquifer (MT3D) using as framework the previously obtained groundwater flow solution [3]. However, the simulation of the transport of contaminants in aquifers is not easy in the Modflow-MT3D environment. It can get more complicated than modeling groundwater flow itself, and it can carry the errors made during the development of the groundwater flow modeling [4]. Since the results of a model will never be better than the quality of the data used for modeling, a lack of an adequate conceptual groundwater flow model will carry errors when modeling the transport of contaminants. Due to the ease which exists today to create mathematical models through Modflow visual Interfaces such as Visual Modflow [5], it is possible to obtain contaminant transport results which may not have much support.

Domenico’s equation is an analytical solution for transport of contaminants in groundwater that can be used with a minimum of groundwater flow information from an aquifer [6]. In cases where the conceptual model of an aquifer is not completely known, Domenico’s equation can be used to easily determine the transport of contaminants in an aquifer [7]. The objective of this study is to model, using Domenico’s equation, a groundwater contaminant plume discharging into a tidally influenced river.

2. Methods

2.1. Site Description

The study area was a wood treatment facility that began operations in the 1930s. The facility is located on the bank of a river which is influenced by tides (Figure 1). The facility is no longer active, but for several years different substances were used as wood preservatives such as creosote. During the last 20 years there have been

several studies carried out in order to characterize and remediate the existing pollution site. These studies have found the presence of creosote in the subsurface and the formation of a groundwater plume that apparently discharges into the river.

The groundwater plume at the site (Figure 2) has already been characterized using a Waterloo Profiler [8-13]. These studies determined that the subsurface below the river consists of three major hydrostratigraphic units: 1) a silty clay upper aquitard, with a thickness between 1 and 8 m, which in places is absent; 2) a sandy aquifer with a thickness between 25 and 30 m with a hydraulic conductivity of 5.2 × 10^–4 m/s, formed by fluvial deposits containing discontinuous lenses of finer or coarser material; and 3) a very compact silty clay aquitard. Most of the groundwater that flows in the aquifer discharges into the river in areas where the upper aquitard is absent. At times during the day, groundwater flow is reversed due to the influence of the tides (flowing from the river into the formation), however the net groundwater flow is towards the river (Figure 1).

The contaminant plume consists of polynuclear aromatic hydrocarbons (PAHs), BTEX compounds (benzene, toluene, ethylbenzene and xylenes) and other hydrocarbons. Compounds that have been detected widely in groundwater at levels above the standards are naphthalene, pyrene and acenaphthene. Phenanthrene and anthracene have also been detected in some samples above the standard.

All these compounds belong to the group of PAHs, and the remaining organic compounds are detected at concentrations below regulatory limits. Because of the concentrations detected, the characterization of the plume for the study area has mainly focused on the distribution of naphthalene in groundwater.

The data obtained during the characterization of the contaminant plume suggests that the contaminant plume in a stationary state (Figure 2). The naphthalene plume discharges to the river at very low concentrations in areas where the top aquitard is absent. Detected inland naphthalene concentrations are in excess of 10,000 µg/l, while the highest concentration detected near the discharge area underneath the river is 64 µg/l (Figure 3).

2.2. Domenico Equation

The following partial differential equation describes the transport of solutes affected by advection dispersion, adsorption and biodegradation, in a saturated porous media in three dimensions and at steady state [7].

(1)

where:

C = concentration of the solute or contaminant in terms of x, y, z;

D_x, D_y, D_z = dispersion coefficients in x, y, z;

V_x = linear groundwater flow velocity;

R = retardation coefficient;

λ = first order biodegradation rate.

Domenico [6] developed an analytical solution of this equation:

(2)

where:

C_o = initial concentration of the dissolved contaminant in the border (source area);

x = distance downstream of the border in the direction of flow;

y = lateral distance perpendicular to the direction of flow;

z = vertical distance perpendicular to the direction of flow;

α_x = longitudinal dispersivity;

α_y = transverse dispersivity;

α_z = vertical;

Y = dimension of the source area in the y direction;

Z = dimension of the source area in the z direction.

This analytical solution was selected because of the limited piezometric data available at the site to properly simulate groundwater flow. The information necessary to use a model like Modflow would require piezometric data obtained at several intervals during the day to properly simulate the hydraulic gradient inversion, and thus reversing the direction of groundwater flow which occurs on a daily basis. The Domenico equation provides an alternative to this problem by not having the need to describe the full aquifer operation, since only an average net flow is necessary. Besides the plume can be modeled in three dimensions using Domenico’s equation and the solution can be easily programmed into a spreadsheet.

2.3. Data and Model Calibration

The initial concentration of naphthalene in the source zone (C_o) was set to 24,378 µg/l. This was the maximum concentration of naphthalene in water that Anthony [8] obtained by dissolving a sample of creosote found at the site and keeping it in balance with water. Anthony [8] also estimated a retarding factor from 2.5 to 20 in adsorption experiments, depending on the concentration of naphthalene.

In Domenico’s equation, the source zone is modeled as a vertical plane perpendicular to groundwater flow. To reproduce the observed distribution of naphthalene in the plume (Figures 2 and 3) it was necessary to define two different sources close to the river, 40 m from the edge of the river. The first defined source is 20 m long, centered with respect to the cross section E-E’. The second source is 70 m long, located 35 m east from the edge of the first source.

The net ground water flow rate was estimated to be 39.42 m/yr [8]. The dispersivities that provided the best fit with the field data during calibration were α_x = 3.75 m, α_y = 1.25 m (α_y = 1/3 α_x) and α_z = 0.05 m (α_z = 1/80 α_x). Charbeneau [14] states that transverse dispersivities field can range from 1/3 to 1/160 of α_x. The retardation factor and the rate of biodegradation that provided a good fit to the model were R = 10 and λ = 0.00063 d^–1.

3. Results and Discussion

Naphthalene field concentrations were closely matched by those of the model after the calibration process (Figure 4). The source zone, where free phase hydrocarbon is detected in the subsurface, is located in the central processing area (Figure 1). However, dissolved concentrations of naphthalene remain constant in groundwater, until they start quickly dropping near the river [11,13]. This shows a spatially variable biodegradation, however Domenico’s Equation does not support these variations in biodegradation. Therefore, during the model calibration process, different configurations were studied to

define the source area. Modeling with a source zone located 130 m inland from the river, in the area of the facility processes where free phase hydrocarbons are present, produced low concentrations of naphthalene near the limits of the river.

Modeling a single source zone 40 m inland produces a much wider plume with naphthalene concentrations in the core much lower than those observed in the field data. In addition, modeled naphthalene concentrations on the boundaries of the plume were much higher than the observed field concentrations. In the case of modeling a single source zone 40 m inland it was possible to obtain a reasonable calibration using extremely high transverse dispersivities (α_y = 300 m = 60 times α_x). No values were found in the literature related to the aquifer dispersivities in cases similar to the studied in here where a contaminant plume discharges into a tidally influenced river. Therefore it was determined that the configuration with two contaminant sources 40 m inland from the river, was the configuration that best simulated the observed field naphthalene concentrations.

A retardation factor of 20, the maximum value obtained by Anthony [8] requires a biodegradation rate (λ) of 1.3 × 10^–4 d^–1, to calibrate and match the observed field distribution of naphthalene. If the retardation factor is increased to 200, the rate of biodegradation necessary for calibration is negligible; however, there is no basis to increase the retarding factor to a high value. Similarly, a decrease in retardation factor would require an increased rate of biodegradation. A biodegradation rate (λ) obtained on microcosm experiments for this site [11] is 2.7 × 10^–3 d^–1, which would require a retarding factor of 2.3, which is at the lower end of possible retardation factors determined by Anthony [8]. Therefore it can be concluded that biodegradation of naphthalene plays a role in the plume’s behavior and distribution.

The behavior of other organic compounds found in the site (benzene and 1-methyl-naphthalene) was simulated using the calibrated model. The initial concentrations and retarding factors were obtained from Anthony [8], the remaining model variables in the model were kept unchanged from the naphthalene simulation. The model results produced reasonable matches to field data for benzene (Figure 5) and 1-methyl-naphthalene (Figure 6).

Both compounds required retardation factors calculated in the estimated average range [8] (8 for benzene and 12 for 1-methyl-naphthalene). The biodegradation rates were respectively 5.2 × 10^–4 and 4.4 × 10^–4 d^–1. These rates are only slightly lower than those obtained for naphthalene. The results obtained with these compounds suggest that the model was correctly implemented and that the framework developed not only reproduces the concentrations of a single compound in the plume, but it can be applied to model other compounds from the same plume.

• 4. Conclusions

• Domenico’s equation was used to model a groundwater contaminant plume discharging into a tidally influenced river.

• No complex groundwater flow scenarios were necessary to perform the simulation.

• Two sources 40 m inland had to be defined to properly simulate the plume behavior.

• The model was able to reproduce the field distribution of naphthalene, benzene and 1-methyl-naphthalene.

• It was determined through the modeling that biodegradation plays an important role on the plume’s behavior.

5. Acknowledgements

The author would like to thank Dr. James F. Barker (University of Waterloo) for his guidance and support during the development of this work.

REFERENCES