Numerical Modelling of Heavy Metals Transport Processes in Riverine Basins

There has been a growing concern in the international community and an increased awareness of riverine pollution problems, particularly with regard to water pollution (Falconer & Lin, 2003). Human and aquatic life is often threatened by the transport of pollutants through riverine systems to coastal waters and it is therefore not surprising to find that, from a water quality point of view, rivers have been studied very extensively and for longer than any other water bodies (Thomann & Mueller, 1987). This is probably due to the fact that people live close to, or interact with, rivers and streams. Many rivers and estuaries have suffered environmental damage due to discharges from manufacturing processes and wastewater from centres of pollution over several decades. In recent years these environmental concerns have made the development of computer models that predict the dispersion of pollutants in natural water systems more urgent. The main attraction of such models, in contrast with physical models, is their low cost and the fact that they easily adapt to new situations. Thus the widespread popularity of mathematical modelling techniques for the hydrodynamic and pollutant transport in rivers justifies any attempt to develop new models based on novel and rigorous approaches (Nassehi & Bikangaga, 1993).


Introduction
There has been a growing concern in the international community and an increased awareness of riverine pollution problems, particularly with regard to water pollution (Falconer & Lin, 2003). Human and aquatic life is often threatened by the transport of pollutants through riverine systems to coastal waters and it is therefore not surprising to find that, from a water quality point of view, rivers have been studied very extensively and for longer than any other water bodies (Thomann & Mueller, 1987). This is probably due to the fact that people live close to, or interact with, rivers and streams. Many rivers and estuaries have suffered environmental damage due to discharges from manufacturing processes and wastewater from centres of pollution over several decades. In recent years these environmental concerns have made the development of computer models that predict the dispersion of pollutants in natural water systems more urgent. The main attraction of such models, in contrast with physical models, is their low cost and the fact that they easily adapt to new situations. Thus the widespread popularity of mathematical modelling techniques for the hydrodynamic and pollutant transport in rivers justifies any attempt to develop new models based on novel and rigorous approaches (Nassehi & Bikangaga, 1993).
This chapter describes numerical modelling of heavy metals in a riverine basin. It would be necessary to recognize and introduce the heavy metals behaviours and different processes during their transportation along the rivers; for example their sources, chemical and physical reactions, and also introducing the environmental conditions affecting the rate of concentration variability of these substances. The one-dimensional (1D) partial differential governing equations (PDE) of hydrodynamic and water quality will be fully described with the corresponding numerical solution methods. As a part of water quality PDE equations the 1D Advection-Dispersion Equation (ADE) will be described and it will be shown that how the dissolved heavy metals, i.e. lead and cadmium, may be numerically modelled through the source term of this equation. Details of the development of a modelling approach for predicting dissolved heavy metal fluxes and the application of the model to the Karoon River, located in the south west of Iran are also provided in this chapter. The Figure 1 illustrates the dissolved heavy metal transport process in a riverine basin. This process is very complicated, since presence and mobility of the heavy metal is highly depended on the environmental conditions, e.g. bed and suspended sediments. A quick review in the literature for the couple of recent years shows that the main attention was focused mostly on the measurements of heavy metals in alluvial rivers with or without sediments. For example, Rauf et al. (2009), Akan et al. (2010) and Kumar et al. (2011) investigated the effect of sediments on the transport of heavy metals; the seasonal variations of the heavy metals in rivers were investigated by Papafilippaki et al. (2008) and Sanayei et al. (2009). The effect of heavy metals on the river self purification process was studied by Mala & Maly (2009).
It should be noted that much attention should be given to the study of heavy metal transport dynamics. In fact, the main factors closely related to the heavy metal pollution transport-transformation in natural bodies are: water flow, sediment motion, pH value, water salinity, water temperature, sediment size distribution, sediment concentration, mineral composition for sediments, degree of mineralization of water and time. Thus theoretically the mathematical model of heavy metal transport dynamics should include equations describing all factors mentioned above (Haung et al., 2007;Haung, 2010). One way to numerically model heavy metals in riverine basins is to assume a reaction coefficient in the source term of the ADE (Advection-Dispersion Equation), which will be described later, showing the presence of the desired substance in the solution. In many studies (such as: Nassehi & Bikangaga, 1993;Shrestha & Orlob, 1996;Wu et al., 2001Wu et al., & 2005 etc.) the researchers assumed a constant reaction coefficient with time, 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 the river. Roshanfekr et al. (2008aRoshanfekr et al. ( & 2008b found that pH and EC play an essential role in adsorption and desorption of heavy metals by the particles in solution.

Theoretical background
Numerical models provide a valuable tool for predicting the fate and transport of dissolved heavy metals in river environments and are increasingly used for such hydro-environmental management studies of river waters. However, computer-based tools used for predicting such heavy metal concentrations are still used infrequently, even though they can support decision-making by the regulatory authorities, marine environment agencies and industry (Ng et al., 1996).
The use of computers has provided the opportunity to better understand and assess our water resources through comprehensive numerical model simulations and testing of various schemes or options. The numerical model allows the user to assess the hydraulic conditions in the river basin and thus, establish a better understanding of human impacts upon a natural or modified river system. Any numerical model used to predict the flow and dissolved heavy metal transport processes in rivers depends primarily on solving the governing hydro-environmental equations. In most riverine systems, the basin is regarded as a 1D system, with longitudinal flow dominating throughout the system. Any type of hydro-environmental model www.intechopen.com Numerical Modelling 94 commonly used by environmental engineers and water managers to predict the dissolved heavy metals concentrations in rivers generally involves solving the hydrodynamic and water quality equations as described below.

Equations for hydrodynamic modelling
In order to dynamically model the heavy metals in riverine basins, the governing hydrodynamic partial deferential equations must be numerically solved. The governing equations and their numerical solutions for modelling the river hydrodynamics (i.e. velocity and water elevation at any point and time) are therefore presented in this section.
For unsteady flow and hydrodynamic modelling the velocity and the water-level at any point of the basin and any time is of interest. The velocity component in rivers is usually assumed as a one-dimensional vector. The one-dimensional governing hydrodynamic equations describing flow and water elevations in rivers are based on the well-known St. Venant equations, applicable to 1D unsteady open-channel flows. Various forms of the St. Venant equations have been formulated in the field for unsteady open-channel flows since the 1950s, when numerical model simulations were first developed. The most widely used form in practice is generally written as (Cunge et al., 1980;Wu, 2008): The individual terms in the momentum equation can be defined as: (1) local acceleration, (2) advective acceleration, (3) pressure gradient and (4) bed resistance, where: T =top width of the channel;  =water elevation above datum; Q =discharge;  =momentum correction factor due to non uniform velocity over the cross-section; A =wetted cross-section area; RAP /  =hydraulic radius; P =wetted parameter of the cross-section; L Q =lateral inflow or outflow (positive for inflow and negative for outflow); Z C =Chezy coefficient; x  =longitudinal distance between two consecutive nodes; g =acceleration due to gravity; xt , =river flow direction and time respectively.
Equations (1) and (2) are solved numerically to provide the varying values of discharge and water elevations.

Equations for water quality modelling
Water quality modelling involves the prediction of water pollution using mathematical simulation techniques. The following sections describe the governing equations and their numerical solution for heavy metals modelling in riverine basins. The model was then applied to Karoon River for lead and cadmium modelling.
The transport of heavy metals in the dissolved phase can be described by the following onedimensional advection-dispersion equation (ADE) (Kashefipour, 2002): The individual terms in the advection-dispersion equation refer to: (1) local effects, (2) transport by advection, (3) longitudinal dispersion and turbulent diffusion, (4) sources or sinks of dissolved heavy metals and (5) transformation term defining absorbed and desorbed particulate fluxes to or from sediments (source term), where: S =cross-sectional averaged dissolved heavy metal concentration; x D =longitudinal dispersion coefficient; d S 0 =source or sink of dissolved heavy metal; d t S =transformation term defining absorbed and desorbed particulate fluxes to or from sediments (source term).
Sources or sinks of dissolved heavy metals can be defined as: where: L Q =lateral inflow or outflow discharge; L S =lateral inflow or outflow dissolved heavy metal concentration; x  =distance between two consecutive cross-sections which can be either constant or variable.
The longitudinal dispersion coefficient ( x D ) in natural rivers is dependent upon many hydrodynamic parameters including: depth, width, velocity and shear velocity (Fischer et al. 1979). There are many empirical and/or semi-empirical equations describing this very important dynamic coefficient.  presented two empirical equations based on applying the dimensional analysis procedure to more than 80 data sets in 30 natural rivers, to estimate longitudinal dispersion coefficient in natural channels and showed that these equations performed relatively better than the other existing equations. In this chapter the   where: H =averaged depth over the cross-section; U = cross-sectional average velocity; U * =local shear velocity.
In addition Tavakolizadeh (2006) used this dispersion coefficient for water quality modelling in Karoon River and achieved acceptable results for different water quality parameters.

Heavy metals modelling
Heavy metals can exist in both the dissolved and adsorbed particulate phases in rivers. The distribution between these two phases may be expressed by a partition coefficient. In recent years much effort has been focused on correlating the partitioning rate of heavy metals in particulate and dissolved phases to several environmental factors and water properties. It would be possible to numerically model each type of heavy metals in the water column, separately (Wu et al., 2001(Wu et al., & 2005. However, it is sometimes important to the environmental managers to have a good understanding of the ratio of these two types of heavy metal presence in water body. A few researchers assume a reaction coefficient in the transformation term, i.e. d t S , in the form of Equation (6) to model either dissolved or particulate heavy metal in the water column.
Since in 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 d t S in Equation (3) is equal zero. On the other hand, the transformation of heavy metal from dissolved phase to the particulate phase and vice versa is assumed to be equal. A review of the literature has shown that a number of researchers include this type of assumption in their models, such as Wu et al. (2005). However, another group of researches, for example Nassehi & Bikangaga (1993), assumed a decay term having a form of Equation (6) with a constant coefficient.
The fate and decay of toxic subtances can result from physical, chemical, and/or biological reaction. Transformation processes are those in which toxic subtances are essentialy irreversibly destroyed, changed, or removed from the water system. These transformation processes are often described by kinematic equations. Most decay processes are expressed as first-order reactions. Therefore, in this chapter the first-order chemical reaction was used as the transformation parameter in Equation (3) for dissolved heavy metal modelling and is written as follows (Zhen-Gang, 2008):  is a reaction coefficient rate, which may have a positive or negative value as the dissolved heavy metals disappears or accumulates in a given river section.
Since the exchange of the heavy metal substance between particulate and dissolved phases is a chemical process and is highly dependent on the environmental conditions, it seems that www.intechopen.com

97
assuming zero value or a constant value for the reaction coefficient  may not provide an accurate simulation. Therefore, the following sections describe the procedure for calculating varied reaction coefficients for dissolved lead and cadmium modelling using pH and EC changes in the water column. The key point is that the chemical characteristics of the flow, such as pH and EC, can affect the dissolved heavy metals from sorption and desorption, to or from the sediments, and these characteristics can have an important effect on the dissolved heavy metal concentrations. For more accurate heavy metal modelling, varied reaction coefficients has been suggested linking pH and EC to the kinetic processes.
Based on the different characteristics of each heavy metal (such as lead, cadmium and etc.) the varied reaction coefficient should be computed and the corresponding relation of the reaction coefficient should be used separately for each metal.
Reaction coefficient, also known as the decay coefficient, is the ratio for the number of atoms that decay in a given period of time compared with the total number of atoms of the same kind present at the beginning of that period (Zhen-Gang, 2008). There are different environmental parameters, such as: temperature, pH, salinity and etc., which generally affect the reaction coefficient in heavy metals. Therefore,  can be defined as: The  value may be related with to temperature as given by the following equation (see Orlob, 1983): where: 20  = reaction coefficient at 20 o C; TEMP = temperature of water; O = temperature coefficient which it can vary from 1.047 to 1.135 (Orlob, 1983).
Theoretically the best mathematical model for heavy metal reaction coefficient should consist of all factors affecting the heavy metals concentration. Roshanfekr et al. (2008aRoshanfekr et al. ( & 2008b found that pH and EC play a more essential role in adsorption and desorption of heavy metals by the particles in solution. Therefore assuming a variable reaction coefficient seems more reasonable. The first step to calculate a variable reaction coefficient is to select the parameters that are most likely affecting the dissolved heavy metal concentration. Then the efforts can be focused on finding suitable functions to represent the reaction coefficient rate for dissolved heavy metals (e.g. lead and cadmium) in rivers. In calibrating the model against measured dissolved lead and cadmium data, five approaches for each dissolved metal can be used: 1. No rate of reaction for dissolved heavy metal (used by some researchers and models for equilibrium conditions). 2. A constant reaction coefficient for the rate of reaction during the whole simulation time (the general practice in dissolved heavy metals modelling used by many researchers). 3. A time varying reaction coefficient for the rate of the reaction using pH as a variable. 4. A time varying reaction coefficient for the rate of the reaction using EC as a variable. 5. 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 and cadmium concentrations with the corresponding measured values at sites and for the times of measured values. Final values of the reaction coefficients for each indicator 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 (20) to (25)) were added to the model as a part of the numerical solution of the ADE (see Equation (3)). The model was then validated using the corresponding measured data for different time series at the survey site.

Numerical modelling
Following sections describe the numerical methods used to solve the hydrodynamic and water quality partial differential equations for heavy metals modelling.

Numerical methods for hydrodynamic equations solutions
There are many implicit and explicit numerical methods used for solving 1D hydrodynamic equations (i.e. Equations (1) and (2)), in which the stability, accuracy and consistency of the numerical solution are important. Almost all implicit methods are unconditionally stable, however the accuracy of model predictions is highly depended on the Courant number (i.e.   r CU t x /   ). Different methods for numerical solution of the above equations may be found in Abbott and Basco (1997).
In this study, the numerical model FASTER (Flow and Solute Transport in Estuaries and Rivers) (Kashefipour et al., 1999) was used. This model was first developed by Kashefipour (2002) and has since been extended and improved to predict the dissolved heavy metals concentrations for different reaction coefficients. The hydrodynamic module of FASTER model numerically solves the Saint Venant equations using Crank Nicolson with an implicit staggered scheme (Wu, 2008). This model uses the influenced line technique, enabling the model to remain implicit and thereby unconditionally stable and accurate over the whole domain, especially in river confluences. This model can be applied for complex channel networks with complex geometry and has been successfully applied to many research projects in Cardiff University, UK . In the numerical method used for this model, the hydrodynamic equations were formulated on a staggered grid to provide advantages in treating the typical hydrodynamic boundary conditions that are commonly used in such models.

99
A summary of numerical solution method is described here. The difference form of the continuity equation using the Crank-Nicholson central scheme around the node i ( Figure 2) can be written as: where: n and n+1=refer to time t and tt   , respectively;  =a weighting coefficient between 0 and 1 to split the spatial derivatives between the upper and lower time levels ( 0 1    ). The non-conservative form of the momentum equation may be discretised using the finite difference central scheme around the node (i+1/2) as shown in Figure 2, yields: By rearranging Equations (9) and (10) The staggered varying grid size with the numerical scheme is alternatively applied to the continuity and momentum equations to produce a set of linear algebraic equations (i.e. Equations (11) and (12)) for each three consecutive ξ and Q points. Applying Equations (11) and (12) simultaneously at all grid points in the discrete solution domain, from time nt to (n+1) t, yields a matrix system of linear algebraic equations based on n 1   and n Q 1  . A general equation, which contains all of the linear algebraic equations and may be solved by the Thomas algorithm or Gauss elimination procedure for the numerical solution of the governing partial differential equations, may be defined using the following equation: Due to the staggered method, the matrix   B is usually given as a tri-diagonal matrix and then Equation (13) can be generally solved using the well known Thomas algorithm. More information regarding the numerical solution and application of the influence line technique in FASTER model to keep the whole numerical solution in implicit form, specially in river confluences and junctions, can be found in Kashefipour (2002).

Numerical methods for ADE solution
In order to solve the ADE, an implicit algorithm has been developed and used in the FASTER model. This finite volume based solution procedure calculates the advection of a concentrate of solute, or suspended sediments at each face of any control volume, by means of a modified form of the highly accurate ULTIMATE QUICKEST 1 scheme (Lin & Falconer, 1997). As before, a space staggered grid system is used to solve the finite volume form of the ADE, in which the variable S is located at the center of the control volume .
Double integration of the one-dimensional ADE, i.e. Equation (3), with respect to time and volume over the control volume, as shown in Figure 2 gives: where: V= volume.
In Equation (14) the term (1) describes the local change of the solute concentration within the control volume from time (t) to time (t+t). Terms (2) to (5) refer to changes in the solute concentration due to: advection, diffusion, lateral inputs and transformation, respectively. The discrete forms of the terms in Equation (14) using finite volume method can be written as follows: More information regarding the FASTER model may be found in Kashefipour (2002) and Yang et al. (2002).
In the current chapter the dissolved heavy metal reaction coefficient comprises two parameters, including pH and EC. The coefficient was formulated using a linear regression relationship. These varied reaction coefficients were then added to the model for predicting dissolved lead and cadmium. The procedure of development and the equations added to the model can be followed in the modelling application and dissolved heavy metal results sections respectively.

Case study
The Karoon River is the largest and the only navigable river in south west of Iran (see Figure  3(a)). In this study the Mollasani-Farsiat reach of the Karoon River, a distance of 110Km was selected due to the high amount of heavy metal concentrations along this reach (see Figure  3(b)). 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 (see Figure 3(c)) (Diagomanolin et al., 2004). Hydrodynamic and water quality data were acquired via Khuzestan Water and Power Authority (KWPA). A set of six 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, dissolved lead and cadmium concentrations at the Mollasani and Shekare gauging stations (see Figure 3(c)). 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 conducted by Khuzestan Water and Power Authorities in 2000.
The time series 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 main model, respectively. Also concentrations of dissolved lead and cadmium were measured from more than fifteen outfalls and industrial locations along the Mollasani and Farsiat reach. Cross-sections No.1, 36, 49 and 113 corresponded to the crosssections at the gauging stations of Farsiat, Shekare, Ahwaz and Mollasani, respectively.

Application of hydrodynamic modelling
The hydrodynamic module of the FASTER model was calibrated against the data provided for the year 2004, starting from the month of March. The main hydrodynamic parameter used for calibration was the manning roughness coefficient. The river was separated into 4 parts, with the manning coefficient varying from 0.026 to 0.050. Good agreement was obtained between the predicted water levels and corresponding field data at the Ahwaz gauging station as the hydrometric survey site, with a difference in results being less than 3% (see Figure 4(a)) and also the model discharges agreed well with the field data obtained at the Ahwaz gauging station with the difference being less than 16% (see Figure 4(b)). The hydrodynamic module was then validated using another series of measured data (see Figures 5(a) and (b)). As can be seen from these figures the predicted data also gave relatively good correlation with the corresponding measured values. A summary of the statistical analysis of the model results is illustrated in Table 1.

Application of dissolved heavy metals modelling
As discussed above, the rate of reaction plays an important role in predicting the concentration distribution of the dissolved heavy metals for river, estuarine and coastal waters. In the current section effort was made to find suitable functions to represent the reaction coefficient rate for dissolved lead and cadmium metal modelling in rivers. These functions were established using a comparison of the predicted heavy metal concentrations with the corresponding measured values at the Shekare gauging station (see Figure 3). In calibrating the model against measured dissolved lead and cadmium levels, five approaches for each dissolved metal were used. The model with the adjusted rate of reaction was then validated using the corresponding measured data for different time series at the survey site (Kashefipour et al., 2006). In the following sections the equations and the results of the modelling for dissolved lead and cadmium, with the derived equations for the reaction coefficients are illustrated.

Results of dissolved lead modelling
For the first run a conservative dissolved lead was assumed, leading to a zero value for the rate of reaction coefficient. The fit between the predicted and measured data showed 25.2% and 33.3% errors for calibration and verification of the model respectively. As can be seen from Figures 6(a) and (b) the predicted dissolved lead in this case did not agree well with the corresponding measured data at the survey site (i.e. Shekare gauging station).
In the second run for predicting the dissolved lead concentration, the dissolved metal concentration was assumed to be non-conservative with the reaction coefficient in Equation (6) being constant. The best fit between the predicted and measured dissolved lead concentrations occurred for a reaction coefficient of 0.12 day -1 . This assumption led to a prediction error of 3.4% and 17.1% for calibration and verification of the model, respectively (see Figures 6(a) and (b)). 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 (i.e. as a function of salinity).
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 of the following form: where: pH =the mean pH of the river at the site for each time.
The predicted results, for which the reaction coefficient was calculated using Equation (20) in the model, were compared with the corresponding measured values for calibration and verification in Figures 6(a) and (b), respectively. The comparison showed that the error of simulation had reduced to 1.9% and 15% for calibration and verification of the model, respectively.
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 of the following form: where: EC = (Electrical Conductivity), the mean EC of the river at the site for each time (micro mhos/cm).
The predicted lead concentration, for which the reaction coefficients were calculated using Equation (21)  For the last run of the dissolved lead, 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 of the following form: The predicted results, for which the reaction coefficient were calculated using Equation (22) in the model, were then compared with the corresponding measured values for calibration and verification in Figures 6(a) and (b) with the errors of 0.4% and 8.3% respectively. These results showed another improvement in the predicted dissolved lead concentrations.
A summary of the statistical analysis for the different model results is shown in

Results of dissolved cadmium modelling
The same procedure was carried out for dissolved cadmium modelling. For the first run cadmium was assumed to be conservative for which the predicted data did not show reasonable agreement with the measured data and the error was estimated to be 71.1% and 76.4% for calibration and verification of the model respectively (see Figures 7(a) and (b)).
For the second run cadmium was assumed to be non-conservative, with a constant reaction coefficient and the best fit between the predicted and measured data occurred when a reaction coefficient of 0.38 day -1 was used in the model. This assumption significantly reduce the error to 7.7% and 8.5% for calibration and verification of the model, respectively (see Figures 7(a) and (b)).
In the third run the reaction coefficient was assumed to be varying with pH. The most suitable relationship between the reaction coefficient for dissolved cadmium and pH in the river was found to be of following form: The predicted dissolved cadmium for which the reaction coefficients were calculated using Equation (23) in the model, were compared with the corresponding measured values for calibration and verification in Figures 7(a) and (b), respectively. This comparison showed that the error of simulation had declined to 1.8% and 4.2% for calibration and verification of the model, respectively.
The fourth model run was carried out using the EC as a variable in computing the reaction coefficient. For the measured data of dissolved cadmium, the most suitable function for relating the reaction coefficient with EC in the river was found to be: This function was added to the model for predicting the results of dissolved cadmium with EC. The results showed that the error of simulation was 2.5% and 2.6% for calibration and verification, respectively (see Figures 7(a) and (b)).
For the last run for dissolved cadmium a number of simulations were carried out to find a formulation for time varying reaction coefficients for the rate of reaction using both the pH and EC as variables. With using the measured data of dissolved cadmium, the most suitable function for relating the reaction coefficient with pH and EC in river was found to be: The predicted dissolved cadmium for which the reaction coefficients were calculated using Equation (25)

Discussion
Salinity has been found by many investigators to be more influential on the reaction coefficient than any 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 a function of salinity. Nassehi & 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 modelling the partitioning coefficient of heavy metals in the Mersey estuary and concluded that the modelling results agreed well with the measured data.
It should be noted that the proposed method in this chapter is valid for rivers with large variations in salinity and pH. Therefore, this method could be used for rivers either close to the coastal waters, and thus affected by tides, or such rivers that have many agricultural inputs from saline soils draining into them. The chosen reach of the Karoon River in this research was an example of the second type of river. The average minimum and maximum EC for three years data collection (2002 -2004) at the Ahwaz hydrometric station (see Figure  3) were 707 and 2254 µΩ -1 /cm, respectively. The pH values also ranged from a minimum of 7.3 to a maximum of 8.5 at this station.
In deriving Equations (20) to (22) for lead and the similar ones for cadmium (Table 3), 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 and cadmium at the site this assumption was thought to be valid. However, there are some limitations in using these equations. Firstly, simultaneous measurements of dissolved lead and cadmium 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 to the reaction coefficient for dissolved lead and cadmium. Secondly, a one-dimensional model was used. Although one-dimensional models have been successfully used in riverine hydrodynamic and water quality studies, it seems that applying a two-or three-dimensional model may improve the derived equations. However, using two-or three-dimensional models needs extensive field-measured data. The importance of the models is to estimate the desirable variables as accurately as possible. Measuring some special environmental variables, such as heavy metals, in the field is sensitive and ideally needs extensive laboratory studies with sophisticated instruments and with large investments. Measuring pH and EC in riverine systems is relatively straightforward and can be done with even portable instruments. The main idea from this research work is therefore to introduce a procedure that relates the pH and EC values to reaction coefficients of heavy metal substances, such as lead and cadmium, for model predictions. Hence, for heavy metals modelling studies, measurements of pH and EC would be a suitable tool for relatively accurate estimation of these substances.
The results show an average improvement of 25% and 71.5% in error estimations of lead and cadmium, respectively, when using pH and EC as two variables affecting the dynamic processes of these heavy metals.

Summary and conclusions
Details are given of the hydro-environmental model to predict the dissolved heavy metals concentration along rivers using a varied reaction coefficient approach to the source term of the Advection-Dispersion Equation (ADE). The main purpose of this chapter was to describe the dissolved heavy metals modelling procedure and assess the impact of pH and EC on the reaction coefficient used in dissolved lead and cadmium modelling. The hydrodynamic module was first calibrated and validated using the field-measured data taken at a site located along the Karoon River, the largest river in the south west of Iran. In order to find the best equation between pH and EC with the reaction coefficient used in the ADE too, many model runs were carried out and the water quality module was subsequently calibrated by adjusting the reaction coefficient. For each measured lead or cadmium value at any time the most appropriate reaction coefficient was specified and from there for the considered heavy metals a few equation between pH and EC with the reaction coefficient were proposed and added to the water quality module of the model. The main findings from the model simulations can be summarized as follows: