INTRODUCTION
Nowadays river, estuaries and seas around the world have received
heavy loads of toxic, persistent and bioaccumulative 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. Resuspension 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 AdvectionDispersion
(Wu et al., 2005). These pollutants are nonconservative 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 twoor threedimensional flow structure and narrow
channels with a predominately onedimensional flow structure. When a twodimensional
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 onedimensional 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 depthintegrated 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):
Where:
H 
= 
ξ+h = Total water column depth 
ξ 
= 
Water elevation above (or below) datum 
h 
= 
Water depth below datum 
U, V 
= 
Depthaveraged velocity components in x, y directions 
q_{y } 
= 
VH, 
q_{x} 
= 
Unit width discharge components (or depthintegrated
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):
Where:
T 
= 
Top width of the channel 
ξ_{R} 
= 
Water elevation above datum 
Q_{R} 
= 
Discharge 
A 
= 
Wetted crosssectional area 
R =A/P 
= 
Hydraulic radius 
P 
= 
Wetted perimeter of the crosssection 
C_{z} 
= 
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 nonconservative
tracer, dissolved heavy metal distributions in rivers can be modeled by
the advectiondispersion equation. The 2D depthintegrated advection dispersion
equation for dissolved heavy metals modeling is given as (Falconer et
al., 2005):
Where:
C 
= 
Depthaveraged 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 
= 
Depthaveraged dispersion coefficients in x, y directions, respectively 
The 1D crosssectional averaged advection dispersion equation is generally
written as (Kashefipour, 2002):
Where:
C 
= 
Crosssectional averaged dissolved heavy metal concentration 
D_{l} 
= 
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
in Eq. 6 and 7 is equal zero. On the
other hand, no source or sink term is assumed. A review of the literature
has shown that many researchers include this type of assumption in their
models, such as Wu et al. (2005). However, another group of researches,
for example Nassehi and Bikangaga (1993), assumed a decay term having a
form of Eq. 8 with a constant coefficient. In this case it is assumed the
value to be written as follows:
κ 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:
The κ value may be related to temperature as given by the following
equation (Orlob, 1983):
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 MollasaniFarsiat
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 ZerganOmoltomeyr 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 crosssections 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 crosssections 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 crosssections
for the 2D model. Crosssections No.1, 5, 36, 49, 70 and 113 corresponded
to the crosssections at the Farsiat, Omoltomeyr, Shekare, Ahwaz, Zergan
and Mollasani, respectively.
MODEL SETUP
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 2D 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 crosssection
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 crosssections
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 nonconservative 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 nonconservative
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 hydroenvironmental 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.110.18, 0.100.29, 0.100.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.