Abstract: In the present research, the dissolved lead is modeled using a linked 1D and 2D hydro-environmental model. Details are given for the governing equations and solution methods used in these numerical models, together with a discussion of a new development in dissolved heavy metals modeling using varied reaction coefficients in the model. It is found that pH and EC play an essential role in adsorption and desorption of heavy metals by the particles in solution. Therefore, in this study it has been tried to find the best relationship between pH and EC with the reaction coefficient. Relatively close agreement between predicted results and field measured dissolved lead concentrations were obtained for different varied reaction coefficients using the linked 1 and 2D model. Finally, the best relationships for the reaction coefficients for dissolved lead were introduced and the results were successfully compared with the corresponding measured values.
INTRODUCTION
Nowadays river, estuaries and seas around the world have received heavy loads of toxic, persistent and bio-accumulative organic and inorganic compounds over many years. This has caused a major impact on the aquatic environment and the ecosystems for the receiving waters. Fine sediments act as a sink (or source) for the organic chemicals and heavy metals entering aquatic systems. Re-suspension of heavy contaminated bed sediments during flood periods may release a great amount of heavy metals into the water. This desorption of contaminants from their particulate phase can impact significantly on the aquatic environment and its ecosystems. Therefore, numerical tools for accurately predicting dissolved heavy metals transport and the chemical processes of heavy metals between their dissolved and particulate phases in rivers and estuaries are of vital importance for aquatic environmental management.
The behavior of heavy metals in the aquatic environment is strongly influenced by adsorption to organic and inorganic particles. The dissolved fraction of heavy metals may be transported via the process of Advection-Dispersion (Wu et al., 2005). These pollutants are non-conservative in nature and their concentrations depend on salinity and pH, which may vary with time and along a river. As a result, the dissolved metal may come out of solution or even redissolve depending on conditions along the time or channel (Nassehi and Bikangaga, 1993). In many studies (Such as: Nassehi and Bikangaga, 1993; Shrestha and Orlob, 1996; Wu et al., 2001, 2005) the researchers assumed a constant reaction coefficient with time in the source term of the ADE, whereas in the field this coefficient may vary according to the rate of pH, salinity, temperature or even other chemical substances and other hydraulic characteristics of river. Based on to the effect of different substances such as pH and salinity on dissolved heavy metal concentrations in rivers, the necessity of heavy metal modeling with more accuracy in predicting the concentration distributions is inevitable.
In this research it is assumed that the movement of heavy metals from/to the solution can be defined as a first order decay function with a variable reaction coefficient. This assumption was already made by a few researchers with a constant value. Details are also given of the development of a linked 1- and 2- dimensional modeling approach for predicting dissolved heavy metal fluxes and the application of the model to the Karoon River, located in the south west of Iran. Therefore, the main focuses here are: (1) Developing a 1D and 2D linked model and (2). Developing a new approach in the source term of the ADE for 1D and 2D models.
MATHEMATICAL MODELS
In modeling estuarine and riverine processes, the modeling domain often covers areas of different physical characteristics, e.g., large water basins with a two-or three-dimensional flow structure and narrow channels with a predominately one-dimensional flow structure. When a two-dimensional numerical model is used for such cases the detailed bathymetric features of a narrow channel may not be well represented unless a very fine grid system is used. In this case the CPU computation time will significantly increase. Similarly, if a one-dimensional model is used for the wider part of an estuary or river with a two dimensional nature may not be well resolved. For many engineering problems these physical features are prevalent in many estuarine and riverine basins.
The hydrodynamic model used to predict the water elevations and velocity fields in coastal, estuarine and riverine waters initially involves the solution of the governing equations of fluid flow. The 2D hydrodynamic equations are generally based on the depth-integrated 3D Reynolds equations for incompressible and unsteady turbulent flows, with the effects of the earth`s rotation, bottom friction and wind shear being included to give (Falconer, 1993; Falconer et al., 2005):
(1) |
(2) |
(3) |
| H | = | ξ+h = Total water column depth |
| ξ | = | Water elevation above (or below) datum |
| h | = | Water depth below datum |
| U, V | = | Depth-averaged velocity components in x, y directions |
| qy | = | VH, |
| qx | = | Unit width discharge components (or depth-integrated velocities) in x, y directions |
| β | = | Momentum correction factor |
| f | = | Coriolis parameter |
| g | = | Gravitational acceleration |
| τxw, τyw | = | Surface wind shear stress components in x, y directions |
| τxb, τxb | = | Bed shear stress components in x, y directions |
| = | Depth average eddy viscosity |
The 1D governing hydrodynamic equations describing flow and water elevations in rivers are based on the Saint Venant equations, applicable to 1D unsteady open channel flows and are written as (Cunge et al., 1980):
(4) |
(5) |
| T | = | Top width of the channel |
| ξR | = | Water elevation above datum |
| QR | = | Discharge |
| A | = | Wetted cross-sectional area |
| R =A/P | = | Hydraulic radius |
| P | = | Wetted perimeter of the cross-section |
| Cz | = | Chezy coefficient |
Heavy metal is an important factor in the spatial dynamics of estuarine processes. This solute plays an important role in the level of water pollution. It may also affect on the sediment and air pollution as well. As a non-conservative tracer, dissolved heavy metal distributions in rivers can be modeled by the advection-dispersion equation. The 2D depth-integrated advection dispersion equation for dissolved heavy metals modeling is given as (Falconer et al., 2005):
(6) |
| C | = | Depth-averaged dissolved heavy metal concentration |
| = | Sources or sinks of dissolved heavy metals | |
| = | Transformation term defining absorbed and desorbed particulate fluxes to or from sediments (source term or rate of reaction) | |
| Dx, Dy | = | Depth-averaged dispersion coefficients in x, y directions, respectively |
The 1D cross-sectional averaged advection dispersion equation is generally written as (Kashefipour, 2002):
(7) |
| C | = | Cross-sectional averaged dissolved heavy metal concentration |
| Dl | = | Longitudinal dispersion coefficient |
MODEL DEVELOPMENT
To predict both predominantly 2D and 1D hydrodynamic, solute and particulate transport processes in estuarine and riverine waters, the two dimensional model DIVAST (Depth Integrated Velocity and Solute Transport) and the one dimensional model FASTER (Flow And Solute Transport in Estuaries and Rivers) were applied and statically linked together to produce an integrated model for feasible and practical problems and to use both models` advantages. The ULTIMATE QUICKEST scheme was used to solve the ADE in these models (Kashefipour, 2002; Yang et al., 2002).
Since outfalls a proportion of a pollutant that is added to the water column
generally decays and settles according to the chemical and hydraulic characteristics
of the flow, it can be concluded that the pollutant may also be added from
or to the sediments. Therefore, for water bodies close to outfalls the conditions
are not generally consistent with equilibrium conditions. For equilibrium
conditions it can be assumed that the parameter
(8) |
κ is the reaction coefficient rate, which may have a positive or negative value as the dissolved heavy metals disappears or accumulates in a given river section. The reaction coefficient can be affected by several environmental factors such as: temperature, pH and salinity (Kashefipour et al., 2006). Therefore, κ can be defined as:
(9) |
The κ value may be related to temperature as given by the following equation (Orlob, 1983):
(10) |
where, κ20 = Reaction coefficient in 20°C, O = The temperature coefficient which may vary from 1.047 to 1.135 and TEMP = temperature.
In the current study effort was focused on finding suitable functions to represent the reaction coefficient rate for dissolved lead in rivers. In calibrating the model against measured dissolved lead data, five approaches for each dissolved metal were used: (i) No rate of reaction for dissolved heavy metal (used by some researchers and models for equilibrium conditions), (ii) A constant reaction coefficient for the rate of reaction during the whole simulation time, (iii) A time varying reaction coefficient for the rate of the reaction using pH as a variable, (iv) A time varying reaction coefficient for the rate of the reaction using EC as a variable, (v) and A time varying reaction coefficient for the rate of the reaction using both pH and EC variables. For each one of these five cases a number of simulation calibration runs were carried out and the initial reaction coefficient was subsequently adjusted by comparing the predicted dissolved lead concentrations with the corresponding measured values at sites and for the times of measured values. Final values of the reaction coefficients were adopted when the best fit occurred between the series of data. The adjusted rate of reaction coefficients were then correlated with pH, EC and both to find the best relationships for κ as a function of pH and/or EC. These equations (i.e., equations in Table 2) were added to the model as a part of the numerical solution of the ADE (Eq. 8). The model was then validated using the corresponding measured data for different time series at the site. With adding these equations to the model and according to the changes of pH and/or EC with time, the reaction coefficient (κ) in the source term of the ADE will change with time.
FIELD DATA COLLECTION
Karoon River is the largest and only navigable river in south west of Iran (Fig. 1a). In this study the Mollasani-Farsiat reach of the Karoon River, a distance of 110 km was selected for the 1D model due to the high amount of heavy metal concentrations along this reach (Fig. 1b). The Karoon River basin has a network of gauging stations and there are several effluent inputs to the river between gauging stations at Mollasani and Farsiat, including industrial units such as: piping, steel, paint making, agriculture, paper mill, fish cultivation and power plant industries draining from wastewater works into the river (Fig. 1c) (Diagomanolin et al., 2004) and due to the municipal wastes of Ahwaz city and high levels of dissolved heavy metals in this area of river the Zergan-Omoltomeyr reach was selected for the 2D model run which the inputs of the 2D model were obtained from the April 14, 2008 outputs of 1D model (Fig. 1).
| Fig. 1: | (a) Site Location, (b) The Karoon River Network and
gauging Stations and (c) Outfalls, the gauging stations of Karoon
River and the cross-sections used in the model between Mollasani and
Farsiat reach |
Hydrodynamic and water quality data were acquired via Khuzestan Water and Power Authorities (KWPA). A set of five field measured data were available from March 2004, including discharge and water levels measurements at the Mollasani, Ahwaz and Farsiat gauging stations and pH, EC and dissolved lead concentrations at the Mollasani and Shekare gauging stations (Fig. 1c). Also concentrations of dissolved lead were measured from more than fifteen outfalls and industrial locations along the Mollasani and Farsiat reach. Furthermore 113 cross-sections were used as initial topology inputs of the 1D model for the river and a 438x475 grid with a 50 m distance between the meshes was built using the easting and northing of the cross-sections for the 2D model. Cross-sections No.1, 5, 36, 49, 70 and 113 corresponded to the cross-sections at the Farsiat, Omoltomeyr, Shekare, Ahwaz, Zergan and Mollasani, respectively.
MODEL SET-UP
The numerical 1D and 2D models were setup to simulate the flow field and dissolved lead concentrations in the Karoon River between the Mollasani and Farsiat stations using a linked 1- and 2-D model. The water elevations recorded at the Farsiat hydrometric station were chosen as the downstream boundary and the measured discharges and heavy metal concentrations at the Mollasani station were used as the upstream boundary conditions for flow and water quality modules of the 1D model. The water elevations at Omoltomeyr were chosen as the downstream boundary and the discharges and heavy metal concentrations at the Zergan were used as the upstream boundary conditions for flow and water quality modules of the 2D model which the inputs of the 2D model were gained from the 1D model outputs.
The 1D grid, covering the region from Mollasani to Farsiat, was represented using 113 segments, with extensive bathymetric data at each cross-section being collected during the most recent bathymetric survey conduced by Khuzestan Water and Power Authorities in 2000 and also the 2D grid bathmetry was built using the easting, northing and bed elevations of the cross-sections in the 1D grid.
CALIBRATION AND VERIFICATION OF THE HYDRODYNAMIC MODEL
The hydrodynamic module of the FASTER and DIVAST models were calibrated against the data provided for March 2003. The main hydrodynamic parameter used for calibration of the models was the manning roughness coefficient. For the 1D model the river was separated into 4 parts, with the manning coefficient being varied from 0.026 to 0.050. Good agreement was obtained between the predicted water levels and field data at the Ahwaz gauging station, with a difference in results being less than 3% and also the model discharges agreed well with the field data obtained at the Ahwaz gauging station with the difference being less than 17%. For the 2D model the river was selected as 1 part, with the manning coefficient of 0.028. Good agreement was obtained between the predicted water levels and field data at the Ahwaz gauging station, with a difference in results being less than 4% and also the model discharges agreed well with the field data obtained at the Ahwaz gauging station with the difference being less than 17% The hydrodynamic module was then validated using another series of measured data. The predicted data also gave relatively good correlation with the corresponding measured values. A summary of the statistical analysis of the model results is shown in Table 1.
DISSOLVED LEAD MODEL RESULTS
As it was mentioned, a first order decay function was assumed to model the dissolved heavy metals. The reaction coefficient in this function varies due to some environmental factors such as salinity and pH. The functions between this coefficient and those factors were established using a comparison of the predicted heavy metal concentrations with the corresponding measured values at the Shekare gauging station (Fig. 1). In calibrating the model against measured dissolved lead levels, five approaches were used.
| Table 1: | A summary of the hydrodynamic model results |
| Table 2: | The reaction coefficient that were found and used in the 1D and 2D models |
The model was then validated using the corresponding measured data for different time series at the site. In the following sections the equations and the results of the simulation for dissolved lead with the derived equations for the reaction coefficients are shown in Table 1.
For the first run a conservative dissolved lead was assumed, leading to a zero value for the rate of reaction coefficient. In the second run for predicting the dissolved lead concentrations, the dissolved metal concentration was assume to be non-conservative with the reaction coefficient in Equation 8 being constant. However, some research results suggest that the reaction coefficients for different pH and salinity conditions were not constant. A more detailed investigation is being planned to determine the rate of reaction coefficient for different pH and EC.
According to the above findings, it seems that using a variable reaction coefficient, which can be adjusted automatically within a numerical model, depending on the pH, EC or pH and EC values may give better calibration results. A number of simulations were carried out to find a formulation for describing the relationship between the reaction coefficient and the pH value. Using the measured dissolved lead concentrations, it was found that the most suitable relationship between the reaction coefficient for dissolved lead and pH of the river was found (Table 2). Based on the fact that the reaction coefficient relates to the EC value, a number of simulations were also carried out to find a suitable formulation for describing the reaction coefficient with the EC value. Using the measured dissolved lead concentration, it was found that the most suitable relationship between the reaction coefficient for dissolved lead and the EC of the river was found which is shown in Table 2. For the last run, a number of simulations were carried out to find a formulation for describing the relationship between the reaction coefficient and both the pH and EC variables. Using the measured dissolved lead concentrations, it was found that the most suitable relationship between the reaction coefficient for dissolved lead and pH and EC as variables for the river was found and it is shown in Table 2.
The best fit between the predicted and measured data showed 25.2, 33.3, 26.8 and 36.4% errors for calibration and verification of the 1D and 2D models, respectively, when zero reaction coefficient considered. In the second run the dissolved metal concentration was assumed non-conservative with a constant reaction coefficient. The best fit between the predicted and measured dissolved lead concentrations occurred for a reaction coefficient of 0.12 day-1 for 1D and 2D models. This assumption led to a prediction error of 3.4, 17.1, 7.9 and 19.2% for calibration and verification of the 1D and 2D models, respectively. The predicted lead concentrations, when the variable reaction coefficients as a function of pH considered agreed relatively well with the corresponding measured values with the error being (1.9 and 15%) and (9 and 18.7%) for 1D and 2D models and for calibration and verification stages, respectively. The predicted lead concentrations when a variable reaction coefficient as a function of EC used, were compared with the corresponding measured values and the calculated errors were less than the previous run for both the 1D and 2D models. The errors for 1D model were calculated 0.8 and 8.3% for model calibration and verification, respectively. The corresponding errors for 2D model were calculated 8.8 and 15.4%, respectively.
The predicted lead concentrations for the case with a variable reaction coefficient as a function of pH and EC agreed much better than the other cases with the corresponding measured values with the calculated error percentages being (0.4 and 8.3%) and (8.7 and 12.9%) for the 1D and 2D models and for calibration and verification stages, respectively.
As it is clear from Table 3 the predicted dissolved lead concentrations improved, giving lower errors when varying reaction coefficients were applied as a part of ADE.
| Table 3: | A summary of the dissolved lead model results |
Salinity has been found by many investigators to be more influential on the absorption and desorption of heavy metals than other environmental or water properties in riverine and estuarine waters. The results published by Turner et al. (2002) showed that the trace metal distribution coefficient in estuarine waters is primarily function of salinity. Nassehi and Bikangaga (1993) calculated the value of the reaction coefficient for dissolved zinc in different elements of a river. Wu et al. (2005) used salinity for modeling the partitioning coefficient of heavy metals in the Mersey estuary and concluded that the modeling results agreed well with the measured data.
In deriving the equations in the Table 2, it was assumed that the environmental factors and water properties remained constant during the whole simulation period. Since the model was calibrated using measured dissolved lead at the site this assumption was thought to be valid. However, there are some limitations in using these equations. As it was said simultaneous measurements of dissolved lead were only made at one site and for six months. More field measured data are needed to validate and improve the formulae, which relate the pH and EC values with the reaction coefficient for dissolved lead. The importance of the models is to estimate the desirable variables as accurately as possible. As it can be seen from the 1D and 2D model results the error of the 2D model is higher than the 1D model and that is because the inputs of the 2D model is obtained from the 1D model and the 1D errors can be transferred to the 2D model.
Table 3a and b show an average improvement of 25 and 19% in error estimations of the predicted lead concentrations for 1D and 2D models, respectively, when using pH and EC as two variables affecting the movement of heavy metal from/to the solution.
CONCLUSIONS
Details are given of hydro-environmental modeling study to predict the heavy metal concentrations using a linked 1 and 2D models. In this model, it was assumed that the movement of heavy metal from/or to the solution acts as one order decay function with some environmental parameters such as EC and pH being affecting on it. This link was statically and the necessary information was transferred from 1D model to the 2D model. The linked model was then applied to the limited existing dissolved lead data in Karoon River, located in south west of Iran. In calibrating and validating the water quality modules of the models, the predicted dissolved lead were compared against the corresponding measured field data at a specific site and the main findings from these simulations were:
| • | Five different procedures were used for estimating the rate of reaction coefficient for dissolved lead, including: a zero reaction coefficient, a constant reaction coefficient, a varied reaction coefficient with pH, a varied reaction coefficient with EC and a varied reaction coefficient with pH and EC. |
| • | Improvements were achieved in the predicted dissolved lead concentration distributions when varied reaction coefficients were used. |
| • | The best fit between the predicted and measured values for simulation with a constant reaction coefficient was obtained for 1D and 2D models when the coefficient was set to 0.12 day-1 for dissolved lead. |
| • | According to equations in Table 2 and the measured pH and EC values, the range of reaction coefficients for the 1D and 2D models were evaluated to be: (0.11-0.18, 0.10-0.29, 0.10-0.43) for three suggested procedures respectively. The error estimation was decreased from an average of 30 to 4% for the 1D model and 30 to 11% for the 2D model. |
ACKNOWLEDGMENTS
The authors are grateful to the Khuzestan Water and Power Authorities for provision of data.