**1. Introduction**

48 Hydrodynamics – Natural Water Bodies

Thoman, R.V. & Segna, J.S. (1980). Dynamic phytoplankton-phosphorus model of Lake

Tucci, C.E.M. (1998). *Modelos Hidrológicos.* Coleção ABRH de Recursos Hídricos, UFRGS,

Van den Berg, M.S.; Coops, H.; Meijer, M.L.; Scheffer, M. & Simons, J. (1998). Clear water

Woodward, G. & Hildrew, A.G. (2002). Food web structure in riverine landscapes.

Wu, F.C.; Shen, H.W. & Chou, Y.J. (1999). Variation of roughness coefficients for

Wu, J. (1982). Wind-Stress Coefficients over Sea-Surface from Breeze to Hurricane. *Journal of* 

White, F.M. (1974). *Viscous Fluid Flow*, McGraw-Hill, ISBN 0070697124, New York.

Vol.125, pp. 934-942, ISSN Print 0733-9429, ISSN Online 1943-7900.

Publishers, ISBN 0250403323, Michigan.

ISBN8570258232, Porto Alegre, Brazil.

Verlag, ISBN 0387982841, New York.

2427.

Ontario: ten-year verification and simulations. In: *Phosphorus management strategies for lakes*, C. Loehr; C.S. Martin & W. Rast (Eds.), 153-190, Amr Arbor Science

associated with a dense Char a vegetation in the shallow and turbid Lake Veluwemeer, the Netherlands. In: *Structuring Role of Submerged Macrophytes in Lakes*, E. Jeppesen; M. Søndergaard & K. Kristoffersen (Eds.), 339-352, Springer-

*Freshwater Biology*, Vol. 47, pp. 777-798, ISSN Print 0046-5070. ISSN Online 1365-

unsubmerged and submerged vegetation. *Journal of Hydraulic Engineering-Asce*,

*Geophysical Research-Oceans and Atmospheres*, Vol.87, pp. 9704-9706, ISSN 0196-2256.

The common basis of the modeling activities is the numerical solution of the momentum and mass conservation equations in a fluid. For hydrodynamic modeling, the Navier-Stokes equations are usually simplified according to the specific water body properties, obtaining, for example, the shallow water equations, so called because the horizontal scale is much larger than the vertical. Therefore, in cases where the river has a relation width-depth of 20 or more and for many common applications, variations in the vertical velocity are much less important than the transverse and longitudinal direction (Gordon et al., 2004). In this sense, the equations can be averaged to obtain the vertical approach in two dimensions in the horizontal plane, which adequately describes the flow field for most of the rivers with these characteristics.

At the same time, the contaminant transport models have evolved from simple analytical equations based on idealized reactors to sophisticated numerical codes to study complex multidimensional systems. Since the introduction of the classic Streeter-Phelps model in the 1920 to evaluate the Biochemical Oxygen Demand and dissolved oxygen in a steady state current, contaminant transport and water quality models have been developed to characterize and assist the analysis of a large number of water quality problems.

This chapter presents the numerical solution of the two-dimensional Saint-Venant and Advection-Diffusion-Reaction equations to calculate the free surface flow and contaminant transport, respectively. The solution of both equations is based on a second order Eulerian-Lagrangian method. The advective terms are solved using the Lagrangian scheme, while the Eulerian scheme is used for diffusive terms. The specific application to the Coatzacoalcos River, Mexico is discussed, having as a main building block the water quality assessment supported on mathematical modelling of hydrodynamics and contaminants transport. The solution method here proposed for the two-dimensional equations, yields appropriate results representing the river hydrodynamics and contaminant behaviour and distribution when comparing whit field measurements.

In this work is presented the structure of a numerical model giving an overview of the program scope, the conceptual design and the structure for each hydrodynamic, pollutants transport and water quality modules that includes ANAITE/2D model (Torres-Bejarano and Ramirez, 2007). The numerical solution scheme is detailed explained for both Saint-Venant and the Advection-Diffusion-Reaction equation. To validate the model, some comparisons were made between model results and different field measurements.

A Study Case of Hydrodynamics and Water Quality Modelling: Coatzacoalcos River, Mexico 51

**Biochemical**: Dissolved Oxygen (DO), Biochemical Oxygen Demand (BOD), Fecal

**Eutrophication**: Ammonia Nitrogen (NH3), Nitrates (NO3), Organic Nitrogen (N\_org.),

 **HAPs**: Acenaphthene, Phenanthrene, Fluoranthene, Benzo(a)anthracene, Naphthalene. The transport and transformation of the different environmental parameters was carried out

> *CCC C C u v Ex Ey <sup>Γ</sup> t x yx x y y*

The reaction mechanism, *Γc*, is used to represent the water quality parameters, and it is

*C*

(4)

**Physics**: Temperature, Salinity, Suspended Solids, Electric Conductivity.

Inorganic phosphorous (phosphate, PO4), organic Phosphorous (P\_org.).

**Metals**: Cadmium, Chromium, Nickel, Lead, Vanadium, Zinc.

by applying the two-dimensional approach of A-D-R (eq. 4):

*C =* Concentration of any water quality parameter, (mg/L)

*Γc =* Reaction mechanism (specified for each parameter) (m-1)

*Ex =* Coefficient of longitudinal dispersion, (m2/s) *Ey =* Coefficient of transversal dispersion (m2/s)

Fig. 1. Flow diagram of ANAITE/2D numerical model

solved individually for each of them.

Coliforms (FC).

where:

### **2. The numerical model**

The developed model is a scientific numerical hydrodynamic and water quality model written in FORTRAN; the model has been named ANAITE/2D. The current version of this model solves the Saint-Venant equations for hydrodynamics representation and the Advection-Diffusion-Reaction (A-D-R) equation using a two-dimensional approach to simulate the pollutants fate.

### **2.1 The hydrodynamic module**

In order to obtain a better representation of the hydrodynamics reproduced by the ANAITE model (Torres-Bejarano & Ramirez, 2007), this work presents the solution to change from one-dimensional steady state approach to unsteady two-dimensional flow approximation, solving the two-dimensional Saint-Venant equations (eqs. 1, 2 and 3); these equations describe two-dimensional unsteady flow vertically averaged, representing the principles of conservation of mass and momentum and are obtained from the Reynolds averaged Navier-Stokes equations under certain simplifications (Chaudhry, 1993). These equations have a wide applicability in the study of free surface flow. For example, the flow in open channels with steep slopes (Salaheldin et al., 2000), flows over rough infiltrating surfaces (Wang et al., 2002), propagation of flood waves Rivers (Ying et al., 2003), dam break flow (Mambretti et al., 2008), among others.

$$\frac{\partial h}{\partial t} + \frac{\partial (hu)}{\partial x} + \frac{\partial (hv)}{\partial y} = 0 \tag{1}$$

$$\frac{\partial \mathbf{u}}{\partial t} = -\boldsymbol{u} \frac{\partial \mathbf{u}}{\partial \mathbf{x}} - \boldsymbol{v} \frac{\partial \mathbf{u}}{\partial \mathbf{y}} - \boldsymbol{g} \frac{\partial \mathbf{h}}{\partial \mathbf{x}} + \nu\_t \left( \frac{\partial^2 \mathbf{u}}{\partial \mathbf{x}^2} + \frac{\partial^2 \mathbf{u}}{\partial \mathbf{y}^2} \right) + \mathbf{g} \left( \mathbf{S}\_{\alpha \mathbf{x}} - \mathbf{S}\_{f\mathbf{x}} \right) \tag{2}$$

$$\frac{\partial v}{\partial t} = -\mu \frac{\partial v}{\partial \mathbf{x}} - v \frac{\partial v}{\partial y} - g \frac{\partial \mathbf{h}}{\partial y} + \nu\_t \left( \frac{\partial^2 v}{\partial \mathbf{x}^2} + \frac{\partial^2 v}{\partial y^2} \right) + g \left( \mathbf{S}\_{oy} - \mathbf{S}\_{fy} \right) \tag{3}$$

where:

*Sf =* friction slope, (·)


In this equations system, it was assumed that the effect of the Coriolis force and tensions due to wind at the free surface are negligible given the nature of the problems that focus this work, although their inclusion in the numerical scheme can be done without difficulty.

#### **2.2 The water quality module**

The water quality model, adapted to the main stream of Coatzacoalcos river, simulates the behaviour and concentration distributions for different water quality parameters. The water quality module solves the following parameters, grouped according to the chemical properties:


The transport and transformation of the different environmental parameters was carried out by applying the two-dimensional approach of A-D-R (eq. 4):

$$
\frac{\partial \mathbf{C}}{\partial t} + \boldsymbol{\mu} \frac{\partial \mathbf{C}}{\partial \mathbf{x}} + \boldsymbol{\nu} \frac{\partial \mathbf{C}}{\partial \mathbf{y}} = \frac{\partial}{\partial \mathbf{x}} \left( \mathbf{E} \mathbf{x} \frac{\partial \mathbf{C}}{\partial \mathbf{x}} \right) + \frac{\partial}{\partial y} \left( \mathbf{E} y \frac{\partial \mathbf{C}}{\partial y} \right) \pm \boldsymbol{\Gamma}\_{\mathbf{C}} \tag{4}
$$

where:

50 Hydrodynamics – Natural Water Bodies

The developed model is a scientific numerical hydrodynamic and water quality model written in FORTRAN; the model has been named ANAITE/2D. The current version of this model solves the Saint-Venant equations for hydrodynamics representation and the Advection-Diffusion-Reaction (A-D-R) equation using a two-dimensional approach to

In order to obtain a better representation of the hydrodynamics reproduced by the ANAITE model (Torres-Bejarano & Ramirez, 2007), this work presents the solution to change from one-dimensional steady state approach to unsteady two-dimensional flow approximation, solving the two-dimensional Saint-Venant equations (eqs. 1, 2 and 3); these equations describe two-dimensional unsteady flow vertically averaged, representing the principles of conservation of mass and momentum and are obtained from the Reynolds averaged Navier-Stokes equations under certain simplifications (Chaudhry, 1993). These equations have a wide applicability in the study of free surface flow. For example, the flow in open channels with steep slopes (Salaheldin et al., 2000), flows over rough infiltrating surfaces (Wang et al., 2002), propagation of flood waves Rivers (Ying et al., 2003), dam break flow (Mambretti et

> *h (hu) (hv) <sup>0</sup> tx y*

*uvg ν g S S*

*uvg <sup>ν</sup> gS S t xyy x y*

In this equations system, it was assumed that the effect of the Coriolis force and tensions due to wind at the free surface are negligible given the nature of the problems that focus this work, although their inclusion in the numerical scheme can be done without difficulty.

The water quality model, adapted to the main stream of Coatzacoalcos river, simulates the behaviour and concentration distributions for different water quality parameters. The water quality module solves the following parameters, grouped according to the chemical

*u u u h uu*

*t xyx x y* 

*v v v h vv*

*2 2*

*2 2*

*t 2 2 ox fx*

*t o 2 2 y fy*

(1)

(2)

(3)

**2. The numerical model** 

simulate the pollutants fate.

al., 2008), among others.

where:

properties:

*Sf =* friction slope, (·) *h =* water depth, (m)

*u =* longitudinal velocity, x direction (m/s) *v =* transversal velocity, y direction (m/s)

*t =* turbulent viscosity, (m2/s) *g =* acceleration due to gravity, (m/s2)

**2.2 The water quality module** 

**2.1 The hydrodynamic module** 

*C =* Concentration of any water quality parameter, (mg/L)

*Ex =* Coefficient of longitudinal dispersion, (m2/s)

*Ey =* Coefficient of transversal dispersion (m2/s)

*Γc =* Reaction mechanism (specified for each parameter) (m-1)

The reaction mechanism, *Γc*, is used to represent the water quality parameters, and it is solved individually for each of them.

Fig. 1. Flow diagram of ANAITE/2D numerical model

A Study Case of Hydrodynamics and Water Quality Modelling: Coatzacoalcos River, Mexico 53

*x*(1 ) *-p p x* ( )

*x x x x*

O

*U*

P

Where *f(P)*, *f1*, *f0* y *f2* are the values at points *P*, *x1*, *x0*, and *x2* respectively, *p* is a weight

second-degree interpolation, the three terms of equation (6) are used. Substituting the

 22 2 1 1 1 05 05 05 05 .. .. *nn n n p ii i*

As the value at point P is required in a two dimensional grid, the solution is expressed as

*p pq q q qq*

*p pq q q qq*

*i j i j*

\* \* , , () ; ()

11 11 11 11 1

 11 11 11 11 1 <sup>4</sup> *i j ij i j i j i j i j v vvvvv* \* , , ,,,,

*u uuuuu i j ij i j i j i j i j*

, , ,,,,

*i j ut vt*

*x y*

*nn n*

.. ..

05 05 1 05 05 05 05

. . .. ..

2 22 2

2 22 2

where *p* and *q*, are the Courant numbers in *x* and *y* directions, respectively. The calculation

05 05 1 05 05 05 05

. . .. ..

 *p pp pp* 

*u v*

4

*Cp Cq*

*p q qq qq*

22 2 2

, ,, ,

*i aj b i j i j i j*

1 1 05 05 05 05

*p*

*n*

*i*-2 *i*-1 *i*

*n n n*


*x*

Fig. 3. Notation used, shown in one-dimensional mesh

coefficient which positions the point *P* with respect to *φ<sup>i</sup>*

of the Courant numbers for *u* and *v* is as follows:

\*

bottom line

*i+*1

*n*

*<sup>n</sup>* and *φi-1n*. Since the polynomial is a

  (7)

 

1 2

 

 

1 1

*nn n i j i j i j nn n i j i j i j*

1 1 1 1 1

,, ,

1 1 1 1 1

(8)

,, ,

(9)

(10)

*n*

known values for the three points:

*x*

shown in equation (7).

where:

*t*

*n*+1

*t*

#### **2.3 Numerical solutions**

The Saint-Venant and A-D-R equations are numerically solved using a second order Eulerian-Lagrangean method; a detailed explanation of the solution is presented in this work. The method separates the equations into its two main components: advection and diffusion, which are solved using a combination of Lagrangian and Eulerian techniques, respectively. In this way, the entire equations are solved. Fig. 1 shows the flow diagram for the model general solution.

#### **2.3.1 Numerical grid**

A numerical grid Staggered Cell type is used (Fig. 2). In this grid the scalars are evaluated in the center of the cell and vector magnitudes on the edges.

Fig. 2. Notation for Staggered cell grid

#### **2.3.2 Advection (Lagrangian method)**

The advection solution uses a Lagrangian method whose interpolation/extrapolation principle is base on the characteristics method. In the characteristics method is assigned to each node at tn+1 a particle that does not change its value as it moves along a characteristic line defined by the flow. It locates its position in the previous time tn by the interpolation of adjacent values of a characteristic value, in this case the Courant number, which is assigned to node tn+1. For simplicity, the method is exemplified in one dimension, but is similar for two or three dimensions (Fig. 3).

Assuming that the value of the variable at point P (*φ<sup>p</sup> <sup>n</sup>*), can be calculated by interpolating between the values *φni-1* y *φ<sup>i</sup> <sup>n</sup>* from the adjacent points *x0* and *x1* respectively (Rodriguez et al., 2005).

If a particle at the point *P* travels at a constant velocity *U* will move a distance *x + U Δt* at time *t + Δt*, so it is:

$$
\varphi(\mathbf{x},t) = \varphi\left(\mathbf{x} + \mathbf{u}\Delta t, t + \Delta t\right) = \varphi\left(\mathbf{x} + p\Delta \mathbf{x}, t + \Delta t\right) \tag{5}
$$

if we apply the modified Gregory-Newton interpolation:

$$f\left(P\right) = f\left(\mathbf{x} + p\Delta\mathbf{x}\right) = f\_1 - p\left[f\_1 - f\_0\right] + \frac{p\left(p-1\right)}{2!} \left[f\_2 - 2f\_1 + f\_0\right] - \Lambda\tag{6}$$

Fig. 3. Notation used, shown in one-dimensional mesh

Where *f(P)*, *f1*, *f0* y *f2* are the values at points *P*, *x1*, *x0*, and *x2* respectively, *p* is a weight coefficient which positions the point *P* with respect to *φ<sup>i</sup> <sup>n</sup>* and *φi-1n*. Since the polynomial is a second-degree interpolation, the three terms of equation (6) are used. Substituting the known values for the three points:

$$\log\_p^n = \left(1 - p^2\right)\rho\_i^n + \left(0.5p^2 + 0.5p\right)\rho\_{i-1}^n + \left(0.5p^2 - 0.5p\right)\rho\_{i+1}^n$$

As the value at point P is required in a two dimensional grid, the solution is expressed as shown in equation (7).

$$\begin{split} \boldsymbol{\rho}\_{i-a,j-b} &= \left(1-\boldsymbol{p}^{2}\right) \left[\left(1-\boldsymbol{q}^{2}\right)\boldsymbol{\rho}\_{i,j}^{n} + \left(0.5\boldsymbol{q}^{2}+0.5\boldsymbol{q}\right)\boldsymbol{\rho}\_{i,j-1}^{n} + \left(0.5\boldsymbol{q}^{2}-0.5\boldsymbol{q}\right)\boldsymbol{\rho}\_{i,j+1}^{n}\right] \\ &+ \left(0.5\boldsymbol{p}^{2}+0.5\boldsymbol{p}\right) \left[\left(1-\boldsymbol{q}^{2}\right)\boldsymbol{\rho}\_{i-1,j}^{n} + \left(0.5\boldsymbol{q}^{2}+0.5\boldsymbol{q}\right)\boldsymbol{\rho}\_{i-1,j-1}^{n} + \left(0.5\boldsymbol{q}^{2}-0.5\boldsymbol{q}\right)\boldsymbol{\rho}\_{i-1,j+1}^{n}\right] \\ &+ \left(0.5\boldsymbol{p}^{2}-0.5\boldsymbol{p}\right) \left[\left(1-\boldsymbol{q}^{2}\right)\boldsymbol{\rho}\_{i+1,j}^{n} + \left(0.5\boldsymbol{q}^{2}+0.5\boldsymbol{q}\right)\boldsymbol{\rho}\_{i+1,j-1}^{n} + \left(0.5\boldsymbol{q}^{2}-0.5\boldsymbol{q}\right)\boldsymbol{\rho}\_{i+1,j+1}^{n}\right] \end{split} \tag{7}$$

where *p* and *q*, are the Courant numbers in *x* and *y* directions, respectively. The calculation of the Courant numbers for *u* and *v* is as follows:

$$\mathcal{C}\_{u}(p) = \frac{\boldsymbol{\mu}\_{i,j}^{\*}\Delta t}{\Delta \mathbf{x}\_{i}} \; ; \; \mathcal{C}\_{v}(q) = \frac{\boldsymbol{\upsilon}\_{i,j}^{\*}\Delta t}{\Delta y\_{j}} \tag{8}$$

where:

52 Hydrodynamics – Natural Water Bodies

The Saint-Venant and A-D-R equations are numerically solved using a second order Eulerian-Lagrangean method; a detailed explanation of the solution is presented in this work. The method separates the equations into its two main components: advection and diffusion, which are solved using a combination of Lagrangian and Eulerian techniques, respectively. In this way, the entire equations are solved. Fig. 1 shows the flow diagram for

A numerical grid Staggered Cell type is used (Fig. 2). In this grid the scalars are evaluated in

*Vi, j+1/2*

*Ci,j*

*Ci,j+1*

*<sup>y</sup> Ci+1,j*

*Vi, j-1/2*

The advection solution uses a Lagrangian method whose interpolation/extrapolation principle is base on the characteristics method. In the characteristics method is assigned to each node at tn+1 a particle that does not change its value as it moves along a characteristic line defined by the flow. It locates its position in the previous time tn by the interpolation of adjacent values of a characteristic value, in this case the Courant number, which is assigned to node tn+1. For simplicity, the method is exemplified in one dimension, but is similar for

Assuming that the value of the variable at point P (*φpn*), can be calculated by interpolating

If a particle at the point *P* travels at a constant velocity *U* will move a distance *x + U Δt* at

 

1 10 2 10 1

(6)

2

*Ui+1/2, j*

*<sup>n</sup>* from the adjacent points *x0* and *x1* respectively (Rodriguez et

*xt x u tt t x* ,, , *p xt t* (5)

2

edge *V*

edge *U*

**2.3 Numerical solutions** 

the model general solution.

the center of the cell and vector magnitudes on the edges.

*x*

if we apply the modified Gregory-Newton interpolation:

*p p fP fx px f pf f f ff* !

Fig. 2. Notation for Staggered cell grid

**2.3.2 Advection (Lagrangian method)** 

two or three dimensions (Fig. 3).

between the values *φni-1* y *φ<sup>i</sup>*

al., 2005).

time *t + Δt*, so it is:

*Ui-1/2, j*

**2.3.1 Numerical grid** 

$$
\mu\_{i,j}^\* = (1 - \alpha)\mu\_{i,j} + \frac{\alpha}{4} \left( \mu\_{i+1,j+1} + \mu\_{i-1,j+1} + \mu\_{i+1,j-1} + \mu\_{i-1,j-1} \right) \tag{9}
$$

$$
\psi\_{i,j}^\* = (1 - \alpha)\upsilon\_{i,j} + \frac{\alpha}{4} \left( \upsilon\_{i+1,j+1} + \upsilon\_{i-1,j+1} + \upsilon\_{i+1,j-1} + \upsilon\_{i-1,j-1} \right) \tag{10}
$$

A Study Case of Hydrodynamics and Water Quality Modelling: Coatzacoalcos River, Mexico 55

The 0.023 value was specifically obtained for the studied river. The dispersion terms in

*i j i j*

, ,

 

*i j i j*

This is the term that takes into account the external forces in the Saint-Venant equations, in this case the gravitational forces. It is solved with a centered difference of depth values in

Expanding the derivative and rearranging terms in equation (1), the continuity equation is

1 11 11

11 11 11 11

*n ij ij ij ij*

<sup>4</sup> *ij i j i j i j i j hh h h h* \*

Grouping the obtained terms, the following general equations for velocities are reached:

*<sup>n</sup> <sup>n</sup> <sup>n</sup> u Advec\_u t Dift\_u - Pres\_u tg S S i j ox fx* , ,

The first element of the last term is the free surface slope, which multiplied by the gravity

The second element of the last term represents the bottom stress, which causes a nonlinear

*n n n i j oy fy v Advec\_v t Dift\_v - Pres\_v tg S S* , ,

<sup>1</sup>

<sup>1</sup>

represents the action of gravitational forces. This term can be expressed as:

effect of flow delay and is calculated by the Manning formula:

, ,

*<sup>t</sup> Ey Ey C C CC*

*i i i*

*t Ex Ex Disp\_C C C CC x x x*

*j j j*

*y y y*

1

1

*Pres\_u g Pres\_v g*

*i j i j i j*

, ,,

*h Advec\_h t h h*

1

 

1 1 1

,, , ,

*i j ij ij i j*

*n*

*n*

*n*

*i j*

*i j*

(16)

(17)

(18)

1 1 1

 

*ij ij ij ij*

, , ,,

1 1 1 1 2 2

*x y* , , , , ; 

2 2

(20)

(21)

*uu vv*

*i j*

, ,,,, (19)

*S z xS z ox* 0 00 *<sup>y</sup>* /; / *y* (22)

*x y* \* \* ,, ,,

*ij ij ij ij h h h h*

*n n*

equation (4) are numerically solved as follow:

**2.3.4 The pressure terms** 

the calculation grid (Eq. 17).

**2.3.5 The continuity equation** 

**2.3.6 General solution for velocities** 

as follows:

where:

α is coefficient of relaxation with typical values between 0 – 1. In this work 0.075 was used. *u\*i,j* and *v\*i,j* are spatial velocities in *x* and *y* direction respectively.

The solution method is applied in a similar way to the advective terms presents in the continuity (Eq. 1), momentum (Eqs. 2 and 3) and A-D-R (Eq. 4) equations; *φ* = *h*, *u*, *v* y *C*, are the advective variables, obtained by the lagrangian method in its second order approximation.

#### **2.3.3 The diffusion (Eulerian method)**

*Turbulent diffusion*. The turbulent viscosity coefficient present in the Saint-Venant equations is evaluated with a cero order model o mixing length model in two dimensions (vertically averaged):

$$\nu\_t = l\_m^{-2} \sqrt{2\left(\frac{\partial u}{\partial \mathbf{x}}\right)^2 + 2\left(\frac{\partial v}{\partial y}\right)^2 + \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial \mathbf{x}}\right)^2 + \left(2.34\frac{u\_f}{\kappa h}\right)^2} \quad \therefore \quad l\_m \approx 0.267\,\kappa h \tag{11}$$

where:

*ℓm* = mixing length

*uf* = friction velocity, *uf ghS*

*κ* = von Kármán constant

Thus, the diffusion terms in *x* and *y* are solved respectively by the following formulas:

$$Diff\\_{\boldsymbol{u}} = \nu\_{\boldsymbol{f}} \left\{ \left[ \left( \frac{\boldsymbol{u}\_{i+1,j} - \boldsymbol{u}\_{i,j}}{\Delta \mathbf{x}\_{i}} \right) - \left( \frac{\boldsymbol{u}\_{i,j} - \boldsymbol{u}\_{i-1,j}}{\Delta \mathbf{x}\_{i-1}} \right) \right]^{n} + \left[ \left( \frac{\boldsymbol{u}\_{i,j+1} - \boldsymbol{u}\_{i,j}}{\Delta \mathbf{y}\_{j}} \right) - \left( \frac{\boldsymbol{u}\_{i,j} - \boldsymbol{u}\_{i,j-1}}{\Delta \mathbf{y}\_{j-1}} \right) \right]^{n} \right\} \tag{12}$$

$$Diff\\_\upsilon = \nu\_1 \left\{ \left[ \left( \frac{\upsilon\_{i+1,j} - \upsilon\_{i,j}}{\Delta x\_i} \right) - \left( \frac{\upsilon\_{i,j} - \upsilon\_{i-1,j}}{\Delta x\_{i-1}} \right) \right]^n + \left[ \left( \frac{\upsilon\_{i,j+1} - \upsilon\_{i,j}}{\Delta y\_j} \right) - \left( \frac{\upsilon\_{i,j} - \upsilon\_{i,j-1}}{\Delta y\_{j-1}} \right) \right]^n \right\} \tag{13}$$

Detailed analysis of turbulence, their interpretation and mathematical treatment can be found at Rodi, (1980), Rodríguez et al., (2005).

*Longitudinal and transverse dispersion in rivers*. The A-D-R equation evaluates the dispersion process by *Ex* and *Ey* coefficients.

In this work the longitudinal dispersion coefficient proposed by Seo and Cheong (1998) has been implemented (Eq. 14):

$$\frac{Ex}{hu\_f} = 5.915 \left(\frac{W}{h}\right)^{0.620} \left(\frac{u}{u\_f}\right)^{1.428} \tag{14}$$

The expression for estimating the transverse dispersion coefficient in a river is given by (adapted from Martin and McCutcheon, 1999):

$$Ey = 0.023h\Omega\text{I}^\*\tag{15}$$

The 0.023 value was specifically obtained for the studied river. The dispersion terms in equation (4) are numerically solved as follow:

$$\begin{split} \text{Disp\\_C} &= \frac{\Delta t}{\Delta \mathbf{x}\_i} \Bigg[ \frac{\text{Ex}\_{i+1,j}}{\Delta \mathbf{x}\_i} \big( \mathbf{C}\_{i+1,j} - \mathbf{C}\_{i,j} \big) - \frac{\text{Ex}\_{i,j}}{\Delta \mathbf{x}\_{i-1}} \big( \mathbf{C}\_{i,j} - \mathbf{C}\_{i-1,j} \big) \Bigg]^n \\ &+ \frac{\Delta t}{\Delta \mathbf{y}\_j} \Bigg[ \frac{\text{E} \mathbf{y}\_{i,j+1}}{\Delta \mathbf{y}\_j} \big( \mathbf{C}\_{i,j+1} - \mathbf{C}\_{i,j} \big) - \frac{\text{E} \mathbf{y}\_{i,j}}{\Delta \mathbf{y}\_{j-1}} \big( \mathbf{C}\_{i,j} - \mathbf{C}\_{i,j-1} \big) \Bigg]^n \end{split} \tag{16}$$

#### **2.3.4 The pressure terms**

This is the term that takes into account the external forces in the Saint-Venant equations, in this case the gravitational forces. It is solved with a centered difference of depth values in the calculation grid (Eq. 17).

$$\text{Press\\_}u = -g \left[ \frac{h\_{i+1,j} - h\_{i-1,j}}{2 \Delta x} \right]^n; \quad \text{Pres\\_}v = -g \left[ \frac{h\_{i,j+1} - h\_{i,j-1}}{2 \Delta y} \right]^n \tag{17}$$

#### **2.3.5 The continuity equation**

Expanding the derivative and rearranging terms in equation (1), the continuity equation is as follows:

$$h\_{i,j}^{n+1} = \text{Advec\\_h} - \Delta t \left[ h\_{i,j}^\* \left( \frac{u\_{i+1,j} - u\_{i-1,j}}{2 \Delta x\_i} \right) + h\_{i,j}^\* \left( \frac{v\_{i,j+1} - v\_{i,j-1}}{2 \Delta y\_j} \right) \right] \tag{18}$$

where:

54 Hydrodynamics – Natural Water Bodies

α is coefficient of relaxation with typical values between 0 – 1. In this work 0.075 was used.

The solution method is applied in a similar way to the advective terms presents in the continuity (Eq. 1), momentum (Eqs. 2 and 3) and A-D-R (Eq. 4) equations; *φ* = *h*, *u*, *v* y *C*, are the advective variables, obtained by the lagrangian method in its second order

*Turbulent diffusion*. The turbulent viscosity coefficient present in the Saint-Venant equations is evaluated with a cero order model o mixing length model in two dimensions (vertically

> 2 2 2 2 <sup>2</sup> 2 2 2.34 0.267 *<sup>f</sup>*

*l l h*

1 1

1 1

*n n*

*n n*

\* *Ey* 0.023*hU* (15)

(11)

1 11 1

*ii jj*

1 11 1

*ii jj*

1 428 0 620 . .

*i j ij ij i j ij ij ij ij*

*xx yy* ,, , , , , ,,

*v v vv v v vv*

Detailed analysis of turbulence, their interpretation and mathematical treatment can be

*Longitudinal and transverse dispersion in rivers*. The A-D-R equation evaluates the dispersion

In this work the longitudinal dispersion coefficient proposed by Seo and Cheong (1998) has

*f f Ex W u hu h u* 

The expression for estimating the transverse dispersion coefficient in a river is given by

5.915

*i j ij ij i j ij ij ij ij*

*xx yy* ,, , , , , ,,

*u u uu u u uu*

(14)

(12)

(13)

*t m m u v uv u*

*x y yx h*

Thus, the diffusion terms in *x* and *y* are solved respectively by the following formulas:

*u\*i,j* and *v\*i,j* are spatial velocities in *x* and *y* direction respectively.

approximation.

averaged):

where:

*uf* = friction velocity, *uf ghS κ* = von Kármán constant

*t*

*t*

found at Rodi, (1980), Rodríguez et al., (2005).

(adapted from Martin and McCutcheon, 1999):

process by *Ex* and *Ey* coefficients.

been implemented (Eq. 14):

*ℓm* = mixing length

*Dift\_u*

*Dift\_v*

**2.3.3 The diffusion (Eulerian method)** 

$$\boldsymbol{h}\_{i,j}^{\*} = \frac{1}{4} \left( \boldsymbol{h}\_{i+1,j+1} + \boldsymbol{h}\_{i-1,j+1} + \boldsymbol{h}\_{i+1,j-1} + \boldsymbol{h}\_{i-1,j-1} \right) \tag{19}$$

#### **2.3.6 General solution for velocities**

Grouping the obtained terms, the following general equations for velocities are reached:

$$\left[\boldsymbol{\mu}\_{i,j}\right]^{n+1} = \text{Advec\\_\mu} + \Delta t \left[\text{Diff\\_\mu} - \text{Pres\\_\mu}\right]^n - \Delta t \text{g} \left[\mathbf{S}\_{\alpha x} + \mathbf{S}\_{f\mathbf{c}\_{i,j}}\right]^n \tag{20}$$

$$\left[\boldsymbol{v}\_{i,j}\right]^{n+1} = \text{Advec\\_v} + \Delta t \left[\text{Diff\\_v - Pres\\_v}\right]^n - \Delta t \text{g}\left[\boldsymbol{S}\_{oy} - \boldsymbol{S}\_{f\circ\_{i,j}}\right]^n \tag{21}$$

The first element of the last term is the free surface slope, which multiplied by the gravity represents the action of gravitational forces. This term can be expressed as:

$$S\_{\alpha x} = -\partial \mathbf{z}\_0 \nmid \partial \mathbf{x} \; ; \; S\_{0y} = -\partial \mathbf{z}\_0 \; / \; \partial y \; \tag{22}$$

The second element of the last term represents the bottom stress, which causes a nonlinear effect of flow delay and is calculated by the Manning formula:

A Study Case of Hydrodynamics and Water Quality Modelling: Coatzacoalcos River, Mexico 57

Fig. 4. Location of Coatzacoalcos River

Fig. 5. Minatitlan-Coatzacoalcos industrial park and study zone

We have developed a numerical model that solves the two-dimensional Saint-Venant equations and the Advection-Diffusion-Reaction equation to study the pollutants transport.

**4. Results and discussion** 

$$S\_{fx} = \frac{n^2 \mu \sqrt{\mu^2 + v^2}}{h^{4/3}} \; ; \; S\_{fy} = \frac{n^2 v \sqrt{\mu^2 + v^2}}{h^{4/3}} \tag{23}$$

where:

*Sf* = friction slope, (·)

*n* = Manning roughness coefficient

#### **2.3.7 General solution for A-D-R equation**

Finally, grouping term for A-D-R equation the solution is obtained with:

$$\left[\mathbb{C}\_{i,j}\right]^{n+1} = \text{Advec\\_c} + \Delta t \left[\text{Discp\\_C}\right]^n \pm \Delta t \left[\Gamma\_C\right]^n \tag{24}$$

As mentioned in section 2.2, Γc is solved individually for each parameter.

#### **2.3.8 Stability requirements**

Because a finite difference scheme is used, should be considered the linear stability criteria. The selection of the time step and space must satisfy the condition of Courant-Friedrichs-Lewy (CFL) for a stable solution. The CFL condition for two-dimensional Saint-Venant equations can be written as (Bhallamudi and Chaudhry, 1992):

$$C\_n = \frac{\left(R + \sqrt{gh}\right) \cdot \Delta t}{\Delta x \cdot \Delta y} \sqrt{\Delta x^2 + \Delta y^2} < 1\tag{25}$$

where: *R* = magnitude of the resultant velocity, (m/s)

### **3. Study case: Coatzacoalcos River**

The last stretch of the Coatzacoalcos River located in the Minatitlan-Coatzacoalcos Industrial Park (MCIP), with about 40 km length, is part of an area of vast natural diversity, where the high population concentration creates important environmental changes, due to pressures arising mainly from consumption and industrial activities. Currently, insufficient information exists for MCIP regarding hydrodynamics and water quality at the river stretch. This is one of the most polluted rivers in Mexico and is consequently a critical area in terms of industrial pollution.

Coatzacoalcos is a commercial and industrial port that offers the opportunity to operate a transportation corridor for international traffic, the site is the development basis for industrial, agricultural, forestry and commercial in this region; by the volume of cargo is considered the third largest port in the Gulf of Mexico (Fig. 4). The area of Coatzacoalcos river mouth has had a rapid urban and industrial growth in the last three decades. In this area, the largest and most concentrated industrial chemical complex, petrochemical and derivatives has been developed in Latin-America.

The importance of this industrial park, formed mainly by the Morelos petrochemical complex, Cangrejera, Cosoloacaque and Pajaritos, is such that 98% of the petrochemicals used throughout the country are produced there. Fig. 5 shows the industrial and petrochemical facilities located in this area.

Fig. 4. Location of Coatzacoalcos River

*fx* 4 3 *fy* 4 3 *nu u v nv u v S S h h* / / ;

*<sup>n</sup>* <sup>1</sup> *<sup>n</sup> <sup>n</sup> C Advec\_c t Disp\_C t i j*, *<sup>C</sup>*

Because a finite difference scheme is used, should be considered the linear stability criteria. The selection of the time step and space must satisfy the condition of Courant-Friedrichs-Lewy (CFL) for a stable solution. The CFL condition for two-dimensional Saint-Venant

> 2 2 <sup>1</sup> *<sup>n</sup> R gh t <sup>C</sup> x y x y*

The last stretch of the Coatzacoalcos River located in the Minatitlan-Coatzacoalcos Industrial Park (MCIP), with about 40 km length, is part of an area of vast natural diversity, where the high population concentration creates important environmental changes, due to pressures arising mainly from consumption and industrial activities. Currently, insufficient information exists for MCIP regarding hydrodynamics and water quality at the river stretch. This is one of the most polluted rivers in Mexico and is consequently a critical area in terms

Coatzacoalcos is a commercial and industrial port that offers the opportunity to operate a transportation corridor for international traffic, the site is the development basis for industrial, agricultural, forestry and commercial in this region; by the volume of cargo is considered the third largest port in the Gulf of Mexico (Fig. 4). The area of Coatzacoalcos river mouth has had a rapid urban and industrial growth in the last three decades. In this area, the largest and most concentrated industrial chemical complex, petrochemical and

The importance of this industrial park, formed mainly by the Morelos petrochemical complex, Cangrejera, Cosoloacaque and Pajaritos, is such that 98% of the petrochemicals used throughout the country are produced there. Fig. 5 shows the industrial and

Finally, grouping term for A-D-R equation the solution is obtained with:

As mentioned in section 2.2, Γc is solved individually for each parameter.

equations can be written as (Bhallamudi and Chaudhry, 1992):

where: *R* = magnitude of the resultant velocity, (m/s)

derivatives has been developed in Latin-America.

petrochemical facilities located in this area.

**3. Study case: Coatzacoalcos River** 

of industrial pollution.

where:

*Sf* = friction slope, (·)

*n* = Manning roughness coefficient

**2.3.8 Stability requirements** 

**2.3.7 General solution for A-D-R equation** 

2 22 2 22

(23)

(24)

(25)

Fig. 5. Minatitlan-Coatzacoalcos industrial park and study zone
