JEQ Grow Your Career With ASA
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Related articles in JEQ
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via ISI Web of Science (1)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Gu, R. R.
Right arrow Articles by Chung, S.-W.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Gu, R. R.
Right arrow Articles by Chung, S.-W.
Agricola
Right arrow Articles by Gu, R. R.
Right arrow Articles by Chung, S.-W.
Related Collections
Right arrow Surface Water Quality
Right arrow Pesticides
Right arrow Coupled Flow/Transport Models
Right arrow Other Environmental Contamination
Journal of Environmental Quality 32:620-632 (2003)
© 2003 American Society of Agronomy, Crop Science Society of America, and Soil Science Society of America

TECHNICAL REPORTS
Surface Water Quality

A Two-Dimensional Model for Simulating the Transport and Fate of Toxic Chemicals in a Stratified Reservoir

Roy R. Gu*,a and Se-Woong Chungb

a Dep. of Civil Engineering, Iowa State Univ., Ames, IA 50011
b Water Management Center, Korea Water Resources Corporation, Seoul, Korea

* Corresponding author (roygu{at}iastate.edu)

Received for publication February 20, 2002.

   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 MODEL DEVELOPMENT
 MODEL APPLICATION
 RESULTS AND DISCUSSION
 CONCLUSIONS
 APPENDIX
 REFERENCES
 
A two-dimensional reservoir toxics model is essential to establishing effective water resources management and protection. In a reservoir, the fate of a toxic chemical is closely connected with flow regimes and circulation patterns. To better understand the kinetic processes and persistence and predict the dissipation of toxic contaminants in the reservoir during a spill or storm runoff event, a toxics submodel was developed and incorporated into an existing laterally integrated hydrodynamics and transport model. The toxics submodel describes the physical, chemical, and biological processes and predicts unsteady vertical and longitudinal distributions of a toxic chemical. The two-dimensional toxicant simulation model was applied to Shasta Reservoir in California to simulate the physico–chemical processes and fate of a volatile toxic compound, methyl isothiocyanate (MITC), during a chemical spill into the Sacramento River in 1991. The predicted MITC concentrations were compared with those observed. The effect of reservoir flow regimes on the transport and fate of the toxic substance was investigated. The results suggested that the persistence of MITC is significantly influenced by different flow regimes. Methyl isothiocyanate is more persistent in the reservoir under an interflow condition due to reduced volatilization from deep layers than under an overflow condition. In the overflow situation, the plume moved more slowly toward the dam and experienced greater dissipation. This analysis can assist in toxic spill control and reservoir management, including field sampling and closure of water intakes.

Abbreviations: MITC, methyl isothiocyanate


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
MODEL DEVELOPMENT
 MODEL APPLICATION
 RESULTS AND DISCUSSION
 CONCLUSIONS
 APPENDIX
 REFERENCES
 
THE USE OF agricultural chemicals such as pesticides (herbicides and insecticides) has produced various environmental problems. Accidental spills or storm runoff of toxic substances into rivers and reservoirs can severely affect fisheries and local water supply systems in a short time period ranging from hours to days (Capel et al., 1988; Chatterjee, 1991). The environmental and economic effects of a chemical spill into a reservoir can be fatal for a community, in particular if the community is heavily dependent on the reservoir for water supply. In midwestern U.S. reservoirs, frequent occurrence of high levels of pesticides during late spring and early summer is one of the most serious water quality issues (Goolsby et al., 1993). The levels of some pesticide concentrations occasionally exceeded their maximum contaminant level (MCL) for drinking water with potential adverse effects on human health (Eldridge et al., 1994; Munger et al., 1997). Therefore, a full understanding of the transport and fate of toxic substances in a reservoir during a spill or peak runoff loading period is essential to establishing an effective management strategy and protecting the scarce water resources.

Once a toxic chemical spill or runoff enters a reservoir, its fate and transport processes are governed by various factors including flow conditions, chemical and biological conditions of the waterbody, weather factors, and properties of the toxicant. The processes can be further complicated if the reservoir is stratified and a density current forms due to the temperature difference between river and ambient waters. The kinetic (decay) processes can be largely influenced by flow regimes (i.e., overflow, plunge flow, underflow, and interflow) (Fig. 1) . On entering the reservoir, river flow pushes water ahead to some distance until the buoyancy force, attributable to the density difference between the river and reservoir surface water, becomes sufficient to counteract the momentum of incoming water. An overflow occurs when the density of the river inflow is less than or identical to that of the receiving reservoir water, moving forward over the surface of the reservoir. The negative buoyancy force drives the inflow downward to sink below the surface when the river water is heavier than the ambient reservoir water, forming a plunge flow. After plunging, the river flow is submerged down the old river channel and moves along the riverbed, becoming an underflow. Both plunge flow and underflow are induced by initial momentum and negative buoyancy due to density difference. At a distance where the density of the underflow is reduced by mixing and dilution to that of the ambient reservoir water, the underflow separates from the bottom and penetrates horizontally into the reservoir, forming an intruding interflow. Underflow and interflow are controlled by reservoir stratification and mixing with ambient water. When the inflow becomes an overflow, some chemicals dissipate more readily through volatilization, photolysis, and oxidation because sufficient turbulence, sunlight, and oxygen are available near the water surface. Therefore, a water quality model should adequately link the kinetic processes of a toxic chemical with reservoir hydrodynamics for better prediction of the persistence and spatial distribution of the toxicant during a spill or runoff.



View larger version (23K):
[in this window]
[in a new window]
  		Fig. 1. Flow processes in a stratified reservoir.

 
A laterally integrated two-dimensional hydrodynamics and transport model has been used for modeling water temperature and conventional water quality constituents in reservoirs (Martin, 1988; Cole and Buchak, 1994). Chung and Gu (1998) employed the model to simulate and analyze the transport of a toxic pesticide that spilled into a reservoir by treating the toxicant as a conservative substance. Due to the lack of a toxic contaminants modeling component, the model's application has been limited to the transport and mixing processes of a tracer under certain flow regimes such as a deep interflow (Fig. 1), where the kinetic processes of the toxic substance may be insignificant. To understand the effectiveness of physico–chemical dissipation processes and behavior of the toxic contaminant under various flow regimes, the development of a toxics modeling component is necessary.

In this study, to help better understand the kinetic processes and persistence and predict the dissipation of toxic contaminants in a reservoir, a two-dimensional reservoir toxics model is developed and tested against field data. The aim of this modeling investigation is to provide information useful in the analysis of contaminated flows in a stratified reservoir. The model was constructed by developing a toxics submodel and incorporating it into an existing reservoir hydrodynamics and transport model. It simulates unsteady vertical and longitudinal distribution of a toxic chemical with respect to various transport and transformation processes. The model was validated against field data collected from the Shasta Reservoir, California, during the pesticide VAPAM (sodium methyl dithiocarbamate) spill after railroad tank cars derailed on a canyon bridge over the Sacramento River in 1991. Observed and simulated chemical concentrations were employed to analyze the fate and transport processes of a product of the spilled chemical, methyl isothiocyanate (MITC). Numerical experiments, using computational runs for selected scenarios, were conducted to examine the effect of two reservoir flow regimes, interflow and overflow, on the dissipation of the volatile toxic compound (MITC) in the reservoir during the spill event. The analysis and results of this study are intended to assist in contamination control and remediation after a toxic spill and guide field monitoring and water quality management.


   MODEL DEVELOPMENT
TOP
 ABSTRACT
 INTRODUCTION
 MODEL DEVELOPMENT
MODEL APPLICATION
 RESULTS AND DISCUSSION
 CONCLUSIONS
 APPENDIX
 REFERENCES
 
Governing Equations
A stationary sediment bed may be validly assumed in reservoirs and lakes where wind has little effect on bottom conditions. As toxic substances appear to decay very slowly (first-order rate constant of 0.001–0.1 d-1) in adsorbed particulate form or sediment environment (Schnoor, 1996; Thomann and Mueller, 1987; DiToro et al., 1981), sediment decay rates of toxic chemicals are assumed negligible compared with those in the dissolved phase. If sediment decay rates are significant (greater than 0.1 d-1), laboratory or field data are required to accurately determine the rates. Linear sorption–desorption kinetics were adopted because chemical concentrations in reservoir water are mostly environmentally relevant (i.e., less than one-half water solubility). An instantaneous local equilibrium was assumed for sorption processes because the time scale for sorption reactions is much smaller than that of other kinetic and macroscopic transport (advection and diffusion) processes in reservoirs.

With the aforementioned assumptions, total chemical concentrations in water column (Ct,w) and bed sediments (Ct,b) are formulated by the conservation of mass:

[1]

[2]
where subscripts t, d, and p denote the total, dissolved, and particulate phases of a chemical, respectively; subscripts a, w, and b denote air, water, and bed, respectively; fd and fp are the fractions of dissolved and particulate chemicals to the total chemical; t is time; x is longitudinal Cartesian coordinate (positive to the right); B is waterbody width; U and W are longitudinal and vertical flow velocities, respectively; Dx and Dz are longitudinal and vertical constituent dispersion coefficients, respectively; Kf is diffusive exchange rate between water column and pore water of the bed; {phi} is the porosity of the bed sediments; KP, KH, KO, and KB are rate constants for photolysis, hydrolysis, oxidation, and biotransformation, respectively; kl is air–water exchange rate due to volatilization; H is Henry's constant; Ca is vapor phase concentration; {Phi}NPS is the lateral mass flow rate of constituent per unit volume; vs is the net settling velocity of sorbed chemical; z is the depth of water from water surface; and y is the distance from reservoir bottom. The chemical kinetic reaction rates and model input parameters can be determined from field and laboratory experiments, estimation using chemical properties, and from the literature (Lyman et al., 1982; Schnoor et al., 1987).

Transport of sediments in the water column and reservoir bed was formulated using the two-dimensional laterally integrated advection–diffusion equations, assuming a stationary bed condition:

[3]

[4]
where Cs is sediments concentration including inorganic and organic (detritus) sediments in water column; ql is lateral mass flow rate of sediments per unit volume for water column; vss is net settling velocity of suspended solids; Css is suspended solids concentration; {Delta}z is grid cell thickness; Kdt is detritus decay rate; vdt is net detritus settling rate; Cdt is detritus concentration in water column; Cb is the total bed sediment concentration; {gamma}om is rate multiplier for temperature adjustment; Ks is organic bed sediment decay rate; and qb is lateral mass flow rate of sediments per unit volume for bed sediments.

Physical, Chemical, and Biological Processes
A partition coefficient (K) was used to determine the fractions of dissolved and particulate chemical forms to the total chemical based on the linear sorption–desorption kinetics. For the equilibrium phase partitioning, two phases and particle interaction models (Di Toro, 1985; Bierman, 1994) were adopted. In the initial stages of a chemical spill, dissolved chemical in the water column may transfer to interstitial water in the bottom sediment by diffusion processes, whereas it may be diffused from the bottom sediment to the water column during the recovery phase, depending on the concentration gradient. This mass transfer was estimated as a function of molecular weight (M) of a chemical and porosity ({phi}) of bed sediments (Di Toro et al., 1981) as follows:

[5]

The amount of chemical loss due to photolysis is a function of the quantity and wavelength distribution of incident light, the light absorption characteristics of the chemical, and the efficiency at which absorbed light yields a chemical reaction. A first-order reaction constant was used to estimate the rate of photolysis (Thomann and Mueller, 1987):

[6]
where Kdo is the direct near-surface photolysis rate; Io is the daily amount of incoming solar radiation at the water surface; Io' is the light intensity at which Kdo was measured; D is the radiance distribution function; Do is the radiance distribution function near the surface; and Ke({lambda}max) is the net light extinction coefficient at {lambda}max, the wavelength of the maximum light absorption. Extinction coefficients for water, inorganic, and organic sediments were used to calculate a net light extinction coefficient, Ke (Cole and Buchak, 1994).

Hydrolysis is a major dissipation process in which a chemical compound reacts with water molecules and results in a cleavage of chemical bond. Generally, hydrolysis is a second-order reaction because of dependence on the molar concentrations of [H+], [OH-], or water mediators. A pseudo–first order rate constant, KH = Kn + Ka[H+] + Kb[OH-], was used, where Kn is the neutral hydrolysis rate, Ka is the acid catalyzed hydrolysis rate, and Kb is the base catalyzed hydrolysis rate. The Arrhenius function was used to adjust the rate with temperature (Thomann and Mueller, 1987).

Toxic chemicals can be oxidized by a reaction with either dissolved oxygen or free radicals in natural waters (Schnoor, 1996). A pseudo–first order reaction rate constant was used to compute degradation of a chemical by oxidation assuming that the rate of free-radical formation (oxidant) is relatively constant (steady state). The rate was adjusted for water temperature with the Arrhenius equation.

Biotransformation is the microbially mediated decay process. The model was designed to treat the biotransformation process either by second-order or by pseudo–first order kinetic reactions based on the Monod equation:

[7]
where CB is the bacterial concentration; yB is the bacterial yield coefficient; µmax is the maximum specific growth rate; and Kx is the half-saturation constant. The rate was adjusted for water temperature.

A modified two-film theory (Mackay, 1982; O'Connor, 1983) was used to compute the gaseous transfer of chemical from air to water and water to air. The transfer rate was proportional to the concentration gradient between the chemicals in the water column and in the overlying atmosphere. The overall transfer coefficient, kl, was given as follows:

[8]
where Kl is the liquid film coefficient and Kg is the gas film coefficient. The kl value was computed as a function of the chemical characteristics (H, Kl, and Kg), water velocity, and wind speed. Since the transfer coefficient for the open bodies of water such as reservoir and lake are largely affected by wind, O'Connor (1983) and Mackay (1982) equations were incorporated into the model for the estimation of the liquid and gas film transfer coefficients.

Numerical Method and Computer Programming
The governing equations were solved using the finite-difference method as in the laterally integrated hydrodynamics and mass transport model, CE-QUAL-W2 (Cole and Buchak, 1994). The dependent variables are toxicant and sediment concentrations in the water column and bed sediments. The independent variables are longitudinal distance, flow depth, and time. The flow chart in Fig. 2 briefly describes the overall algorithms of the model. The water surface elevation and horizontal momentum are computed simultaneously in the model based on an implicit finite-difference scheme, which allows the use of a reasonable time scale for field application over entire stratification cycles. The equations for toxicant and sediment transport were solved using the explicit QUICKEST (Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms) finite-difference scheme (Leonard, 1979). Vertical turbulent transfer of toxic contaminants was determined from the vertical shear of horizontal velocity and a density gradient dependent Richardson number function (Chapra, 1997). The computer program for the toxics submodel was created and linked to the CE-QUAL-W2 code by using the Microsoft (Redmond, WA) Fortran Powerstation program, which can be found in Chung (1998). The physical, chemical, and biological properties of a toxic chemical and its kinetic reaction rates need to be provided through an input file. The model can compute kinetic decay processes either independently by providing individual kinetic reaction rates or collectively by providing a lumped value of transformation rate or half-life. The values of half-life and its variation with time can be specified in the model.



View larger version (33K):
[in this window]
[in a new window]
  		Fig. 2. Flow chart of the two-dimensional reservoir toxics model.

 
Boundary and Initial Conditions
The model generates vertical and longitudinal distributions of a toxic chemical in response to time-varying boundary conditions. Flows and contaminant loadings are specified at the upstream boundary and tributaries. The discharges of contaminant through outflow are used as downstream boundary conditions. At the surface of the waterbody, precipitation, evaporation, and meteorological changes over time are provided. The lateral boundary conditions include overland (surface) runoff, ground water gain or loss, water withdrawal, and concentrations of contaminants contained in the lateral flows. The time-varying flows and mass loadings can be defined as either a step function or linearly interpolated value between two data points.

The initial toxic substance and sediment concentrations in the water column can be specified as either a uniform value, a single vertical profile for all segments, or a vertical profile for each segment. In bed sediment, they can be specified as either a uniform value or longitudinal profile for each segment. If nonuniform initial conditions are preferred, grid-wide vertical or longitudinal profiles should be provided in an input file.


   MODEL APPLICATION
TOP
 ABSTRACT
 INTRODUCTION
 MODEL DEVELOPMENT
 MODEL APPLICATION
RESULTS AND DISCUSSION
 CONCLUSIONS
 APPENDIX
 REFERENCES
 
The model was applied to the 1991 Sacramento River spill (Gu et al., 1996) for investigating the physical processes and chemical reactions of MITC under various flow regimes in the Shasta Reservoir about 50 km downstream of the spill site. The reservoir is represented by a finite-difference grid system consisting of 32 segments with lengths of 200 to 1000 m in the longitudinal direction and 52 vertical layers with thickness of 0.5 to 4.0 m (Chung and Gu, 1998). The computational domain has finer grids in the plunging region and coarser grids in other regions. In the present application, dilution by ambient reservoir water and decay due to physico–chemical reactions were quantified using observed and simulated MITC concentrations. Simulation results were compared with those under the assumptions that MITC can be treated as a conservative tracer and that the level of MITC dissipation by kinetic processes (hydrolysis, photolysis, and volatilization) is insignificant in the deep and cold reservoir environments with low levels of turbulence, light, temperature, and oxygen.

Spill and Site Descriptions
On 14 July 1991, railroad tank cars derailed on a canyon bridge at the locally infamous Cantara Loop, California. Approximately 49 000 to 72 000 L of VAPAM liquid formulation were estimated to have spilled into the Sacramento River about 50 km upstream of the Shasta Reservoir (Fig. 3) (Rosario et al., 1994; Gu et al., 1996). VAPAM (sodium methyl dithiocarbamate) is a soil fumigant and decomposes into more stable products, primarily MITC (CH3N=C=S), which is more toxic in water. During transport along the Sacramento River, high levels of light and oxygen induced photolytic and oxidative transformation of sodium methyl dithiocarbamate into MITC (Wang et al., 1997). The LC50 (96 h) of MITC (i.e., the lethal concentration of 50% of test organisms during 96 h of exposure) for bluegill (Lepomis macrochirus) is 130 µg L-1 (Tomlin, 1994; Worthing, 1987). A great number of fish died during the spill. A large portion of the MITC was vaporized into the air and impaired human health in the vicinity of the 1991 spill site (Chatterjee, 1991). The residual MITC, eventually, entered into the Shasta Reservoir at 0000 h of 16 July 1991.



View larger version (41K):
[in this window]
[in a new window]
  		Fig. 3. Map of the spill site and sampling stations in the Shasta Reservoir.

 
Water sampling and field monitoring were conducted by the California Environmental Protection Agency and Central Valley Regional Water Quality Control Board to keep track of the spill after it entered the Shasta Reservoir (Rosario et al., 1994; Gu et al., 1996). The sampling stations (1–18) are located in the upper reach of the Shasta Reservoir (Fig. 3). The spill managers marked the contaminant plume with the red dye rhodamine WT to track its progress in the reservoir after the plume plunged into the reservoir and became an underflow. An artificial mixing device was installed at the distance of about 2.6 km downstream from the head of the reservoir (between Stations 8 and 9 as shown in Fig. 3) to mix the contaminant plume with the reservoir water in an attempt to lower the level of contamination.

The reservoir has a capacity of 1.9 x 109 m3. The average longitudinal river bed slope and reservoir side expansion angle are 0.3% and 1°, respectively. The reservoir is approximately 109 m deep at the dam wall. The inflow depth and width at the head of the reservoir (near Doney Creek) are 1.0 and 70 m, respectively. The reservoir was stratified during the spill period. Temperatures in the reservoir showed that vertical and longitudinal variations were significant (Gu et al., 1996), but lateral variations are generally small in the upper reach from the head of the reservoir to its confluence with the Squaw River arm (Chung, 1996). Field observations by the California Environmental Protection Agency (Rosario et al., 1994; Gu et al., 1996) and California Central Valley Regional Water Quality Control Board (1991) showed that water temperatures in the reservoir ranged from 19 to 27.5°C in the epilimnion and 7 to 16°C in the hypolimnion and that the river flow entering the reservoir averaged 7.5 m3 s-1 with a temperature of 18 to 24°C during the spill. Under these conditions, the contaminated river flow plunged after entering the reservoir and formed an underflow and interflow (Chung and Gu, 1998).

In this model application, the time-varying MITC boundary concentrations were specified at the upper end of the computational domain, Doney Creek. Data were obtained on an hourly basis from 0000 h of 16 July 1991, when the plume arrived at Doney Creek, to 1010 h the next day. Three additional measurements were taken until 1200 h, 19 July 1991. The MITC concentrations varied with time, from 2 mg L-1 at 0000 h on 16 July to 35 mg L-1 (peak) at 0500 h, and 5 mg L-1 at 1000 h on 17 July. They dropped below the detection limit (<0.001 mg L-1) on 19 July 1991. Initial and boundary conditions for flows, temperature, and weather conditions were set based on the observed data (Chung and Gu, 1998).

Physico–Chemical Processes of Methyl Isothiocyanate
The dissipation of MITC due to physico–chemical processes is dependent on the chemical properties and environmental conditions of the reservoir including the level of turbulence, temperature, dissolved oxygen, pH, and solar radiation. These parameters in turn are functions of water depth in a stratified reservoir. Methyl isothiocyanate is a volatile, reactive compound with high solubility and low adsorption potential (Table 1).


View this table:
[in this window]
[in a new window]
  		Table 1. Physical and chemical properties of methyl isothiocyanate (MITC) (Tomlin, 1994).

 
Previous studies concluded that the primary pathways of MITC dissipation in surface water are volatilization and hydrolysis (Rosario et al., 1994; Tomlin, 1994; Wang et al., 1997). The direct photolysis rate of MITC in surface water is a function of water depth (Zepp and Cline, 1977). Draper and Wakeham (1993) found that MITC was resistant to photodegradation in water after five hours of irradiation under laboratory conditions. The sorption of MITC onto suspended solids is negligible in surface water if suspended solids concentration (Css) is low (<1000 mg L-1) (Wang et al., 1997). The fraction of dissolved MITC, fd, was estimated as a function of the partition coefficient, Css, and the fraction of organic matter. Because Css in the Shasta Reservoir was far below 1000 mg L-1 after the spill, more than 97% of MITC was in the dissolved phase. Thus, the simulations assumed that volatilization and hydrolysis dominated the dissipation of MITC.

The rate of volatilization from water to air is influenced by physico–chemical properties, such as molecular weight, vapor pressure, water solubility, and phase transfer coefficient (Henry's constant). Volatilization is also affected by environmental conditions at the air–water interface (i.e., turbulence as a function of wind speed and current velocity) and water depth. Moreover, the artificial device installed to mix the river water and reservoir water could have increased turbulence and thus volatilization. The hydrolysis of MITC is a strong function of the pH level of water (Tomlin, 1994). The first-order hydrolysis rate constant (Kh) for MITC was obtained from its half-life (t1/2) value in neutral water (490 h at pH = 7) using the relation of Kh = 0.693/t1/2.


   RESULTS AND DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 MODEL DEVELOPMENT
 MODEL APPLICATION
 RESULTS AND DISCUSSION
CONCLUSIONS
 APPENDIX
 REFERENCES
 
Reservoir Flow Regimes
Flow velocities and water temperature were used to depict the reservoir circulation patterns during the spill event. Figure 4 shows modeled flow velocity vectors and contour of water temperature, respectively, at 0900 h on 23 July after the spill. It should be noted in the vector plot that only part of the reservoir (near inflow boundary) is plotted due to the large spatial variations in velocities between upstream and downstream. The inflow formed an underflow and interflow after plunging near the head of the reservoir (about 0.8 km downstream from the head) due to the temperature difference between inflow and ambient water. The vector plot obviously shows a considerable vertical motion of the plume within the underflow region driven by the density-induced negative buoyancy forces. When the plume detached from the reservoir bottom, an interflow was created at the depth of 8 to 10 m below the water surface. The temperature contour also shows that the inflow with lower temperature (18–24°C) than the ambient water (26–28°C) plunged into the reservoir and traveled along the reservoir bottom until it reached its equilibrium temperature (21–22°C) about 8 m below the water surface. The occurrence and locations of various flow regimes (plunge flow, underflow, and interflow) captured by the simulated flow and temperatures (Fig. 4) were verified by field measurements in the reservoir after the spill. Field observations (Rosario et al., 1994) indicated that a chemical plume plunged to 6 m below the water surface to form an underflow, which can be seen in the model results (Fig. 4) as well. The field data presented by California Central Valley Regional Water Quality Control Board (1991) also showed that the plume moved along the bottom and intruded into the reservoir, as an interflow, at the depth of approximately 8 to 12 m.



View larger version (39K):
[in this window]
[in a new window]
  		Fig. 4. The flow velocity vector field (a) and contour of reservoir water temperature (b) in the Shasta Reservoir during the spill at 0900 h on 23 July 1991.

 
Fate and Transport of Methyl Isothiocyanate
Simulated MITC concentrations (mg L-1) in the reservoir are presented in Fig. 5 in the form of iso-concentration contours. The MITC contours are 3, 30, and 60 h after the peak passed the upstream boundary. The levels of MITC decay shown in Fig. 5b were obtained by subtracting the concentrations of MITC resulting from both transport and kinetic processes from the MITC concentrations resulting from transport processes only (Table 2). The decay levels are used to determine the effectiveness of the kinetic processes in the total MITC reduction. The chemical decay was greater in the plunge flow and underflow compared with the interflow. This is mainly attributed to the strong turbulence near the surface of the reservoir induced by wind and flow, which accelerated volatilization of the contaminant. In the interflow, the strong reservoir stratification limited the gaseous transfer of the contaminant from water to air.



View larger version (49K):
[in this window]
[in a new window]
  		Fig. 5. Contours of simulated methyl isothiocyanate (MITC) concentrations in mg L-1 (a) and concentration difference between tracer and MITC that represents the level of decay (b).

 

View this table:
[in this window]
[in a new window]
  		Table 2. The effects of flow dilution and physico–chemical reactions on total methyl isothiocyanate (MITC) reduction versus time in various flow regimes.

 
The effectiveness of chemical reactions on the total concentration reduction, ER, was quantitatively analyzed using the predicted maximum chemical concentrations at different times and in different flow regimes (Table 2). In Table 2, Ctracer and Cmitc denote the predicted maximum tracer and MITC concentrations in the reservoir after time t (hours), respectively. The total concentration reduction after time t was obtained by subtracting Cmitc from the peak concentration at the inflow boundary, Co. The total reduction rate (RT), which was defined by 100 x (Co - Cmitc)/Co, increased exponentially with time. More than 80% of concentration reduction was accomplished within 30 h after the peak passed the upper boundary of the reservoir. A dilution index (ID) and a reaction index (IR) were used to separate the contributions of flow dilution by mixing processes and chemical decay by physical and chemical reaction processes. The ID and IR values varied from 0 to 99%, and their sum is always equal to 100%. An index value of 100% means complete dissipation of chemical by dilution (ID = 100%) or reactions (IR = 100%). A large ER value indicates a great effectiveness of chemical reactions and thus rapid dissipation rate at time t.

The IR value reached as high as 7% in the plunge flow region due to the high chemical concentrations and flow turbulence. Although the IR value constantly decreased to 2 to 3% in the underflow and to less than 1% after 81 h in the interflow, the ER value increased with time as the turbulent mixing decreased. The results suggest that the reduction of MITC was primarily attained by flow dilution (i.e., ID > 90%) in the early stage of the spill, but the influence of physical and chemical reaction processes increases with time and can be significant in the late stage of the spill. Thus, it is possible that the persistence of MITC is more likely dependent on the decay processes rather than the turbulent mixing processes in the late stage.

Observed and simulated MITC concentrations at Stations 1, 2, 3, 9, and 13 (Fig. 3) are presented in Fig. 6 . The large deviation between observed and simulated concentrations at Station 1 may be attributed to two causes: (i) inaccurate prediction of wind mixing effects in the epilimnion due to lack of wind direction data and (ii) a continuous breakdown of a small amount of residual sodium methyl dithiocarbamate near the reservoir surface. The sudden decline of measured concentrations at Stations 9 and 13 is due to the artificial aeration system located between Stations 8 and 9 (Fig. 3). The model is not capable of simulating the artificial mixing process and thus the concentration reduction. The comparison of observed and simulated concentrations suggested that the artificial aeration system efficiently contributed to the large drop of downstream MITC level. Using natural log scales for time, in Fig. 7 the simulated tracer and MITC concentrations are compared with the observed value versus time. The model results slightly improved after taking the chemical reactions into account. The slopes obtained from the linear regressions of the tracer and MITC concentrations (R2 = 0.991 and 0.993) represent the dilution rate (q) and the combined rate of dilution and kinetic reactions (q + k), respectively. Thus, the difference (0.025 h-1) between q and q + k is the overall kinetic reaction rate (k) of MITC. The ratio of the overall kinetic reaction rate to the dilution rate (k/q) was about 4% along the various flow regimes, which agrees well with the results obtained in Table 2 (IR = 0.9–6.9%).



View larger version (48K):
[in this window]
[in a new window]
  		Fig. 6. Observed and simulated methyl isothiocyanate (MITC) concentrations at selected sampling stations during the spill.

 


View larger version (21K):
[in this window]
[in a new window]
  		Fig. 7. The linear relationships of observed and simulated methyl isothiocyanate (MITC) concentrations versus time during the early stage of the spill.

 
Effect of Flow Regimes
Flow patterns in a reservoir may have a significant effect on the transport and fate of a toxic contaminant. In previous investigations, various types of models, ranging from zero-dimensional to three-dimensional, were applied to surface water quality problems. Ünlü and Demirekler (2000) employed a zero-dimensional reservoir (a complete-mix reactor) model to examine water quality effects of soil contamination due to oil spills in a reservoir catchment. A vertical two-dimensional model of estuarine hydrodynamics and water quality was used by Liu et al. (2001) to study dissolved oxygen, nitrogen, and phosphorus in a main channel and estuarine wetland. Kopmann and Markofsky (2000) conducted water quality modeling with TELEMAC-3D to investigate the influence of hydrodynamics on the water quality of a lake. A comparison between multidimensional calculations and directly calculated zero-dimensional model results by Kopmann and Markofsky (2000) demonstrated that flow patterns, including current distribution and circulations, have a significant effect on the water quality (dissolved oxygen, nitrogen, phosphorus, and phytoplankton).

The effect of reservoir flow regimes on the persistence of MITC during the late stage of the spill was quantified to examine the hypothesis that an interflow tends to delay the degradation of MITC and an overflow enhances it. An overflow regime was created by using an inflow river temperature equal to or greater than the reservoir surface temperature (27°C) while all other conditions were unchanged. On a summer day, water temperature in a shallow river is usually higher than in a deep reservoir and varies diurnally, peaking at about 1800 h (Gu et al., 1998). A river temperature of 28°C was recorded by the California Central Valley Regional Water Quality Control Board (1991) in the Sacramento River at Station 8 near the Sugarloaf Creek at 1500 h on 21 July 1991 (Gu et al., 1996). It is expected to be higher than 28°C in the late afternoon and could reach its daily maximum at 1800 h. Therefore, a river inflow temperature of 30°C was selected to generate the overflow, which could be a realistic flow scenario. The spatial distributions of a conservative tracer (top) and MITC (bottom) concentrations for the two regimes are displayed in Fig. 8 , showing concentration contours 10 d after the peak concentration passed the inflow boundary.



View larger version (51K):
[in this window]
[in a new window]
  		Fig. 8. The spatial distributions of tracer (top) and methyl isothiocyanate (MITC) (bottom) concentrations in the interflow (a) and overflow (b) after 10 d (26 July 1991).

 
In the interflow situation (Fig. 8a), the plume resided within the thermocline layer. The plume spread out in the longitudinal direction 6 to 10 m below the water surface while the head of the plume reached up to 15 km from the upstream boundary. The peak concentrations are diluted to 350 to 450 µg L-1 due to mixing effect of flow only, which is indicated in the tracer concentrations. The combined effects of dilution and kinetic reactions resulted in further reduction of MITC concentrations to 130 to 200 µg L-1. In the overflow case (Fig. 8b), the plume stayed near the surface of the reservoir and the plume head reached only up to 10 km from the upstream boundary. The peak concentration of tracer was about 350 to 450 µg L-1, indicating that the level of dilution in the overflow situation is similar to that in the interflow. However, in the overflow situation, the toxic plume moved more slowly toward the dam and experienced greater dissipation than the interflow did because more MITC was volatilized from the reservoir surface than from the deep water into the air. The concentrations were less than 60 µg L-1 in most parts of the reservoir. Because drinking water intake and recreation facilities are generally located near the dam, overflow near this location should help water quality management.

Figure 9 shows the maximum tracer and MITC concentrations versus longitudinal distance for the interflow and overflow situations. The concentration difference between tracer and MITC represents the amount of MITC dissipation. It is clearly shown that the effectiveness of physical and chemical reaction processes is greater in the overflow than in the interflow. During both flow regimes, flow dilution only reduces peak concentrations of MITC to 460 µg L-1. In the interflow case, however, the MITC peak concentration (195 µg L-1) was still greater than the LC50 (i.e., the lethal concentration of 50% of test organisms) for bluegill (130 µg L-1). The peak concentration in the overflow was 107 µg L-1, which is less than the LC50 for bluegill.



View larger version (57K):
[in this window]
[in a new window]
  		Fig. 9. The longitudinal profiles of maximum tracer and methyl isothiocyanate (MITC) concentrations in the interflow (a) and overflow regimes (b) after 10 d.

 
The effects of reservoir flow regime on the fate and transport of MITC are summarized in Table 3. Total reduction of MITC concentrations was mainly achieved by flow dilution (ID > 95%) in the early stage of the spill, but the effectiveness of chemical reactions (ER) increased with time as the turbulent mixing diminished. After 10 d, the ER values reached up to 57.6% and 76.7% for the interflow and overflow, respectively, indicating that the reduction of MITC in the late stage of the spill, after it has traveled a long distance along the river from the spill site to the reservoir (Fig. 3), is primarily attained through kinetic processes. Therefore, the persistence of the contaminant would be shorter if an overflow occurred during the spill. Under environmental conditions suitable for the formation of an overflow, it may be the flow regime desirable for management and remediation of a volatile toxic chemical in a reservoir after a spill.


View this table:
[in this window]
[in a new window]
  		Table 3. Effects of reservoir flow regime on the fate and transport of the toxic contaminant, methyl isothiocyanate (MITC).

 

   CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 MODEL DEVELOPMENT
 MODEL APPLICATION
 RESULTS AND DISCUSSION
 CONCLUSIONS
APPENDIX
 REFERENCES
 
A two-dimensional reservoir toxics submodel was developed and incorporated into a laterally integrated hydrodynamics and transport model. The model is capable of simulating the fate and transport processes of toxic contaminants, including sorption and desorption, photolysis, hydrolysis, oxidation, biotransformation, volatilization, diffusive exchanges between the bottom sediment and water column, and sediment transport and deposition in a reservoir. The development of the toxics model enhanced the versatility of the original reservoir model and provided an effective tool for describing and predicting the behavior of a toxic contaminant in a stratified reservoir during a spill or runoff.

The model was applied to Shasta Reservoir in California to investigate the effect of various flow regimes on the fate of MITC that was spilled into the reservoir. Predicted flow velocities, water temperature, and chemical concentrations clearly identified various reservoir flow regimes during the spill, including plunge flow, underflow, and interflow. Comparisons between the observed and predicted MITC concentrations showed that the model is reliable for simulating the fate of MITC in the reservoir. Modeled underflow and interflow regimes showed that slow kinetic dissipation processes resulted in long persistence of MITC during the spill. In the early stage of the spill, MITC concentrations were mainly reduced by flow dilution due to transport and mixing processes. The effectiveness of kinetic processes for the chemical dissipation increased with time as the turbulent mixing diminished in the late stage of the spill.

This study demonstrated that reservoir flow regime can substantially affect the persistence and transport of a volatile toxic contaminant in the late stages of a spill. The dilution levels in the interflow and overflow regimes were similar, while the overflow plume moved more slowly and experienced greater chemical loss. In comparison with the interflow regime, the overflow regime reduced MITC contamination, shortened plume length, and lengthened travel time. Because water intake and recreational activities are often located downstream near the dam, occurrence of overflow during river inflow should facilitate water quality management. The coupling of a toxics submodel to the laterally integrated hydrodynamics and transport model results in a tool that can assist field sampling plans, spill control, contaminant remediation, and reservoir management for water quality (e.g., optimization of timing for the closure of water intakes).


   APPENDIX
TOP
 ABSTRACT
 INTRODUCTION
 MODEL DEVELOPMENT
 MODEL APPLICATION
 RESULTS AND DISCUSSION
 CONCLUSIONS
 APPENDIX
REFERENCES
 
Symbols used in this paper: {phi}, bed sediment porosity; {lambda}max, the wavelength of the maximum light absorption (m); µmax, maximum specific growth rate (d-1); {Phi}NPS, lateral mass flow rate per unit volume (mg L-1 d-1); {gamma}om, organic matter rate multiplier for temperature adjustment; {Delta}z, finite difference grid cell thickness (m); B, waterbody width (m); Ca, vapor phase concentration (mg L-1); Cb, sediment concentration in the stationary bed sediment (mg L-1); CB, bacterial concentration (mg L-1); Cd, dissolved chemical concentration; Cdt, detritus concentration in the water body (mg L-1); Cmitc, maximum MITC concentration (mg L-1); Co, peak concentration at inflow boundary (mg L-1); Cp, particulate adsorbed chemical concentration (mg L-1); Cs, sediment concentration in water column (mg L-1); Css, suspended solids concentration (mg L-1); Ct, total chemical concentration (mg L-1); Ctracer, maximum tracer concentration (mg L-1); D, radiation distribution function; Do, radiation distribution function near water surface; Dx, longitudinal diffusion coefficient (cm2 s-1); Dz, vertical diffusion coefficient (cm2 s-1); ER, the effectiveness of chemical reactions for the total concentration reduction; fd, the fraction of the total chemical concentration that is dissolved; fp, the fraction of the total chemical concentration that is particulate; H, Henry's constant; ID, dilution index; Io', light intensity at which Kdo was measured (cal cm-2 d-1); Io, the average daily amount of incoming solar radiation at the water surface (cal cm-2 d-1); IR, reaction index; K, partition coefficient; Ka, acid catalyzed hydrolysis rate (d-1); Kb, base catalyzed hydrolysis rate (d-1); KB, biotransformation rate (d-1); Kdo, direct near surface photolysis rate (d-1); Kdt, detritus decay rate (d-1); Ke, extinction coefficient (m-1); Kf, diffusive exchange rate between water column and pore water (cm2 s-1); Kg, gas film coefficient; KH, hydrolysis rate (d-1); kl, volatilization rate (d-1); Kl, liquid film coefficient; Kn, neutral hydrolysis rate (d-1); KO, oxidation rate (d-1); KP, overall photolysis rate (d-1); Ks, organic bed sediment decay rate (d-1); Kx, half-saturation constant (mg L-1); M, molecular weight of a chemical (g mol-1); qb, the lateral mass flow rate of sediments in reservoir bed per unit volume (mg L-1 d-1); ql, the lateral mass flow rate of sediments in water column per unit volume (mg L-1 d-1); RT, total MITC reduction rate; t, time (s) U, longitudinal velocity (m s-1); vd, deposition velocity of sorbed chemical in the air (m s-1); vdt, net detritus settling velocity (m s-1); vs, the net settling velocity of sorbed chemical in water column (m s-1); vss, the net settling velocity of suspended solids (m s-1); W, vertical flow velocity (m s-1); x, longitudinal Cartesian coordinate (m); y, distance from reservoir bottom (m); yB, bacterial yield coefficient; z, depth of water from surface (m).


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 MODEL DEVELOPMENT
 MODEL APPLICATION
 RESULTS AND DISCUSSION
 CONCLUSIONS
 APPENDIX
 REFERENCES
 


Related articles in JEQ:

This Issue in Journal of Environmental Quality

JEQ 2003 32: 377-382. [Full Text]  




This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Related articles in JEQ
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via ISI Web of Science (1)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Gu, R. R.
Right arrow Articles by Chung, S.-W.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Gu, R. R.
Right arrow Articles by Chung, S.-W.
Agricola
Right arrow Articles by Gu, R. R.
Right arrow Articles by Chung, S.-W.
Related Collections
Right arrow Surface Water Quality
Right arrow Pesticides
Right arrow Coupled Flow/Transport Models
Right arrow Other Environmental Contamination

HOME HELP FEEDBACK