Table of Contents Author Guidelines Submit a Manuscript
Journal of Applied Mathematics
Volume 2012, Article ID 781695, 18 pages
Research Article

Influence of Secondary Currents on Solute Dispersion in Curved Open Channels

1Department of Civil and Environmental Engineering, Seoul National University, Seoul 151-742, Republic of Korea
2Department of Ocean Civil & Plant Construction Engineering, Mokpo National Maritime University, Mokpo, Jeollanamdo 530-729, Republic of Korea

Received 25 December 2011; Accepted 17 April 2012

Academic Editor: Shuyu Sun

Copyright © 2012 Myung Eun Lee and Gunwoo Kim. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


The dispersion coefficient tensor including off-diagonal components was introduced in the flow with secondary currents, which is called skewed shear flow dispersion (SSFD) coefficient tensor, in this paper. To observe the detailed effect of cross-dispersion terms in SSFD model on solute dispersion, mathematical analysis of eigenvalue problem with respect to the equation with SSFD coefficient tensor was performed. The analysis results show the several differences of SSFD model compared to CSFD (conventional shear flow dispersion) model: the oblique direction of principal dispersion with respect to the streamline, the increase of peak concentration, and the change in the eccentricity of elliptical tracer cloud. SSFD coefficient tensor in a streamwise curvilinear coordinate system of curved channel was transformed to those components of fixed Cartesian coordinate system, and 2D numerical model with finite element method was established in the Eulerian-Cartesian coordinate. Through this process, the transformation equation using the depth-averaged velocity field was derived. Several numerical tests were performed to assure the results obtained in the mathematical analysis and to show the applicability of the derived transformation equation on the flow with continuously changing flow direction.

1. Introduction

The advection and dispersion of passive solutes in open channels—which includes pollutant transport in artificial canals, natural streams, and rivers—is an important topic in environmental hydraulics. In open channels, once vertical mixing is completed in the initial period of solute transport, the vertical shear velocity profile increases the longitudinal spreading in the streamline direction [1]. As a result, in flows where the longitudinal flow is dominant, such as straight open-channel flows, solute spreading is commonly described with longitudinal shear dispersion and transverse turbulent diffusion:𝜕𝐶+𝜕𝜕𝑡𝑣𝑠𝐶=𝜕𝜕𝑠𝐷𝜕𝑠𝑠𝑠𝜕𝐶+𝜕𝜕𝑠𝜀𝜕𝑛𝑛𝜕𝐶𝜕𝑛,(1.1) where 𝑠 is the coordinate axis coinciding with the streamline direction; 𝑛 is local coordinate axis that is normal to the streamline; 𝐶 the depth-averaged concentration; 𝑣𝑠 the longitudinal depth-averaged velocity; 𝐷𝑠𝑠 the longitudinal dispersion coefficient; 𝜀𝑛 the transverse turbulent diffusion coefficient. In (1.1), the axis of the longitudinal dispersion always coincides with the streamline of the principal flow. Thus, the distribution of concentration shows axisymmetry with respect to the (𝑠,𝑛) axes, as shown in Figure 1.

Figure 1: Schematic diagram for CSFD model: 𝑥,𝑦: coordinate axes of the Eulerian-Cartesian coordinate system; 𝑣𝑥,𝑣𝑦: components of depth-averaged velocities on the 𝑥,𝑦 axes.

However, the secondary current around pronounced curvatures in many open channels introduces a large magnitude of transverse circulation combined with the principal longitudinal flow. Hence, the solute dispersion by the secondary current cannot be described by only the dispersion in the longitudinal direction; there is a dispersion effect in the transverse direction that is much more effective than the transverse turbulent diffusion. A flow with secondary currents, as that in Figure 2, has a structure with skewed shear profiles having different velocity profiles in two orthogonal directions.

Figure 2: Schematic diagram for SSFD model: 𝑣𝑠,𝑣𝑛: horizontal velocities on the 𝑠,𝑛 axes; 𝑧: coordinate axis of the vertical direction.

Fischer [2] proposed that the cross-dispersion terms should be included in the 2D depth-averaged dispersion equation, to deal with the effect of skewed vertical profiles on the horizontal dispersion process:𝜕𝐶+𝜕𝜕𝑡𝑣𝑠𝐶=𝜕𝜕𝑠𝐷𝜕𝑠𝑠𝑠𝜕𝐶𝜕𝑠+𝐷𝑠𝑛𝜕𝐶+𝜕𝜕𝑛𝐷𝜕𝑛𝑛𝑠𝜕𝐶𝜕𝑠+𝐷𝑛𝑛𝜕𝐶𝜕𝑛,(1.2) where 𝐷𝑠𝑠,𝐷𝑠𝑛,𝐷𝑛𝑠, and 𝐷𝑛𝑛 are components of the full dispersion coefficient tensor and 𝐷𝑛𝑛 denotes the transverse dispersion coefficient. The additional cross-dispersion terms 𝐷𝑠𝑛(𝜕𝐶/𝜕𝑠) and 𝐷𝑛𝑠(𝜕𝐶/𝜕𝑛) indicate mass transport in the longitudinal direction caused by the concentration gradient in the transverse direction and vice versa. These terms rotate a dispersing tracer cloud away from symmetry on the 𝑠 and 𝑛 axes. Once the open-channel flow with secondary currents is regarded as the skewed shear flow, the off-diagonal components 𝐷𝑠𝑛 and 𝐷𝑛𝑠 should be clearly considered on the solute dispersion in a curved flow. However, most studies and pollutant transport models related to the role of secondary currents in the dispersion process only focused on activation of the transverse dispersion that only increases 𝐷𝑛𝑛 [39]. Neglecting cross-dispersion terms, solute dispersion on a curved flow with secondary currents is still described with only the longitudinal and transverse dispersion coefficients in the widely used 2D environmental mixing models:𝜕𝐶+𝜕𝜕𝑡𝑣𝑠𝐶=𝜕𝜕𝑠𝐷𝜕𝑠𝑠𝑠𝜕𝐶+𝜕𝜕𝑠𝐷𝜕𝑛𝑛𝑛𝜕𝐶𝜕𝑛.(1.3) Hereafter, we call this kind of model as a conventional shear flow dispersion (CSFD) model.

In this study, we introduced the full dispersion coefficient tensor to the solute dispersion model with respect to the stream wise curvilinear frame of reference as described in (1.2), and we call this the skewed shear flow dispersion (SSFD) model. This model describes the skewed shear flow dispersion process in a flow with a secondary current that induces strong interaction between the longitudinal and transverse dispersions. For detailed study with respect to the effect of the off-diagonal terms on the passive solute dispersion, the dispersion equation with the SSFD coefficient tensor was mathematically analyzed by solving the eigenvalue problem and comparing the results with the CSFD model. The finite element formulation in fixed Cartesian coordinates was selected for computational modeling to test the SSFD model in curved open channels. To deal with the governing equation on a fixed Cartesian coordinate, a transformation equation for SSFD tensor, originally defined in the stream wise curvilinear coordinate system, was derived using the depth-averaged velocity field.

2. Analysis of SSFD Model

Two-dimensional advection-dispersion equation is simply expressed below in (2.1) with a substantial derivative in the stream-wise curvilinear frame of reference𝑑𝐶𝑑𝑡=(𝐃𝐶)=𝑇𝜕𝐃𝐶,=𝜕𝜕𝑠𝜕𝑛𝑇.(2.1) By the above expression, the full dispersion coefficient tensor 𝐃 of the SSFD model is expressed in a 2×2 matrix to correspond to (1.2) as follows:𝐷𝐃=𝑠𝑠𝐷𝑠𝑛𝐷𝑛𝑠𝐷𝑛𝑛.(2.2) However in the CSFD model represented by (1.3), 𝐃 in (2.1) is replaced by the diagonal dispersion coefficient tensor 𝐃𝑐𝑑𝐶𝐃𝑑𝑡=𝑐𝐶=𝑇𝐃𝑐𝐶,𝐃𝑐=𝐷𝑠𝑠00𝐷𝑛𝑛.(2.3) In order to analyze the effect of off-diagonal components on the full dispersion coefficient tensor of the SSFD model, we performed a mathematical analysis to solve eigenvalue problems related to the dispersion tensor and equation.

Because 𝑇𝐃 in (2.1) is in quadratic form, we may replace 𝐃 with the symmetric matrix 𝐃 by taking off-diagonal components together in pairs and writing the result as a sum of two equal terms [10]:𝑑𝐶𝑑𝑡=𝑇𝐃𝐶,𝐷𝐃=𝑠𝑠𝐷𝑠𝑛𝐷𝑠𝑛𝐷𝑛𝑛,(2.4) where 𝐷𝑠𝑛=(𝐷𝑠𝑛+𝐷𝑛𝑠)/2. Symmetric dispersion coefficient matrices such as 𝐃 have an orthonormal basis of eigenvectors. Thus, if we take these as column vectors, we obtain a matrix 𝐗 that is orthogonal—so that 𝐗𝑇=𝐗1. According to the theory of orthogonal eigenvectors of a specified symmetric matrix, a certain diagonal matrix 𝐃 is obtained by the following relationship:𝐃=𝐗𝐃𝐗1=𝐗𝐃𝐗𝑇.(2.5) It should be noted that 𝐃 is a similar matrix of 𝐃 by orthogonal transformation: 𝐃=𝐗1𝐃𝐗. Substituting (2.5) into (2.4) transforms the quadratic form 𝑇𝐃 to the principal axes form:𝑑𝐶𝑑𝑡=𝑇𝐗𝐃𝐗𝑇||𝐶(𝑠,𝑛)𝐃=||𝐶(𝑥,𝑦),(2.6) where 𝐗𝑇|(𝑠,𝑛)=|(𝑥,𝑦)=[𝜕/𝜕𝑥𝜕/𝜕𝑦]𝑇 denotes the gradient in the coordinates of the principal axes (𝑥,𝑦). As a result, (2.1) is finally transformed to the equation with the principal dispersion coefficient tensor 𝐃 in the rotated (𝑥,𝑦)-coordinate system. Because 𝐃 is a diagonal matrix, the symmetric axes of the equiconcentration line should be parallel to the directions of the 𝑥 and 𝑦 axes. Because the directions of the 𝑥 and 𝑦 axes are rotated from the original 𝑠 and 𝑛 axes by the amounts of orthogonal transformations, the axisymmetric concentration should also be rotated from the 𝑠 and 𝑛 axes. We estimated the magnitude of the rotated angle by solving the following eigenvalue problem:𝐃𝐶=𝜆𝑖𝐶,(2.7) where 𝜆𝑖(𝑖=1,2) is eigenvalue of the corresponding eigenvector in 𝐗, and each value can be derived as follows:𝜆1,𝜆2=𝐷𝑠𝑠+𝐷𝑛𝑛2±𝐷𝑠𝑠𝐷𝑛𝑛22+𝐷2𝑠𝑛.(2.8) The obtained eigenvalues in (2.8) are the diagonal components of 𝐃:𝐃=𝜆100𝜆2(2.9) by the calculation of 𝐃=𝐗1𝐃𝐗. Because the principal dispersion coefficients always have to be positive, the following is the relationship between each component of 𝐃:𝐷𝑠𝑠𝐷𝑛𝑛>𝐷𝑠𝑛.(2.10) By (2.10), we can find the eigenvectors corresponding to 𝜆1 and 𝜆2 that are the column vectors of 𝐗, and the directions of eigenvectors whose directions coincide with the directions of the 𝑥 and 𝑦 axes can be derived. The angle of the counterclockwise rotation 𝜓 of the 𝑥 and 𝑦 axes from the 𝑠 and 𝑛 axes is obtained by𝐷𝜓=arctan𝑠𝑠𝐷𝑛𝑛2𝐷𝑠𝑛±𝐷𝑠𝑠𝐷𝑛𝑛2𝐷𝑠𝑛2+1.(2.11) Therefore, the angles of 𝜓 represent orthogonal rotation of the axes of symmetry: (𝑥,𝑦) from the axes of (𝑠,𝑛)-coordinates.

Another remarkable effect of the off-diagonal terms in the SSFD model is the difference in the resultant peak concentration. When we present the analytical solution of the instantaneously dumped point mass with respect to the (𝑥,𝑦)-coordinates, it is given by𝐶𝑥,𝑦=𝑀,𝑡4𝜋𝑡𝜆1𝜆2𝑥exp24𝜆1𝑡𝑦24𝜆2𝑡,(2.12) where 𝑀 is the total mass of tracer. With this solution, the equiconcentration curves are ellipses with 𝜆1 and 𝜆2 as the lengths of semimajor and minor axes. Using (2.12), the peak concentration at time 𝑡 is inversely proportional to 𝜆1𝜆2, which is given by𝜆1𝜆2=𝐷𝑠𝑠𝐷𝑛𝑛𝐷2𝑠𝑛𝐷𝑠𝑠𝐷𝑛𝑛.(2.13) According to (2.12) and (2.13), the peak concentration with off-diagonal terms 𝑀/(4𝜋𝑡𝐷𝑠𝑠𝐷𝑛𝑛𝐷2𝑠𝑛) is always larger than the peak concentration neglecting off-diagonal terms 𝑀/(4𝜋𝑡𝐷𝑠𝑠𝐷𝑛𝑛). Thus, if the off-diagonal components in the SSFD condition are ignored, the dilution degree of the pollutant is likely to be overestimated. Assuming that the corresponding CSFD and SSFD models have the same diagonal components for each dispersion coefficient tensor, the ratio of the overestimated dilution degree due to neglect of the off-diagonal components by applying the CSFD model is𝐶𝑝𝐃𝑐,𝑡𝐶𝑝=𝐃,𝑡𝑀/4𝜋𝑡𝐷𝑠𝑠𝐷𝑛𝑛𝑀/4𝜋𝑡𝐷𝑠𝑠𝐷𝑛𝑛𝐷2𝑠𝑛=𝐷𝑠𝑠𝐷𝑛𝑛𝐷2𝑠𝑛𝐷𝑠𝑠𝐷𝑛𝑛<1,(2.14) where 𝐶𝑝 is the peak concentration. Although (2.14) is valid only in the case of instantaneous release of mass at 𝑡=0 in a uniform flow field, the above analytical analysis provides the knowledge of possible overestimation in dilution degree when the CSFD model is wrongly applied to a skewed shear flow field of secondary currents.

The final characteristic of the SSFD model compared to the CSFD model is the change in the eccentricity of ellipses. The dispersive scale in the direction of the symmetric axes of concentration depends on the pair of principal dispersion coefficients, which are 𝜆1,𝜆2 for the SSFD model and 𝐷𝑠𝑠,𝐷𝑛𝑛 for CSFD. Therefore, the comparison between the magnitudes of the principal dispersion coefficients accounts for the difference in the shapes of the ellipses by 𝐃 and 𝐃𝑐. The ellipses of the concentration in the SSFD model have larger eccentricity than that of the corresponding CSFD model by the relation:𝜆max1,𝜆2=𝐷𝑠𝑠+𝐷𝑛𝑛2+𝐷𝑠𝑠𝐷𝑛𝑛22+𝐷2𝑠𝑛𝐷max𝑠𝑠,𝐷𝑛𝑛,𝜆min1,𝜆2=𝐷𝑠𝑠+𝐷𝑛𝑛2𝐷𝑠𝑠𝐷𝑛𝑛22+𝐷2𝑠𝑛𝐷min𝑠𝑠,𝐷𝑛𝑛.(2.15)

Including (2.15), three characteristics of the SSFD model that contrast with the CSFD model have been pointed out by the eigenvalue problem solved in this section; the rotation of the major dispersion axis about the streamline, the larger peak concentration, and the larger eccentricity of the elliptical concentration. These results show that the application of the CSFD model to flows with a secondary current is not accurate because skewed vertical profiles clearly exist in the secondary current combined with the principal flow along the curves. The oblique direction of the principal dispersion with respect to the streamline and other characteristics related to the peak concentration and shape of the equiconcentration curves were demonstrated in a numerical experiment that used the computational model established in the Eulerian coordinate system, as presented in the next section.

3. Coordinate Transformation

To deal with diverse flow directions in natural streams and rivers with irregular boundaries, conventional river hydrodynamics and mass transport models are usually established in a fixed Eulerian coordinate system, where implementing a horizontal unstructured grid is convenient. In computational models established for such curved channels with continuously changing flow directions, the principal direction of anisotropic dispersion is usually not parallel to the axes of the Cartesian coordinates. Therefore, in commonly used CSFD models, components of 𝐃𝑐 constantly defined in a stream-wise curvilinear frame of reference (𝑠,𝑛) are transformed into nodal dispersion parameters with respect to the Eulerian-Cartesian coordinates [1115]. The dispersion coefficient tensor 𝐃𝑐 of (𝑠,𝑛)-coordinates is related to nodal parameters in (𝑥,𝑦)-coordinates through the Jacobian matrix expressed with a nodal velocity vector.𝐉=𝜕𝑥𝜕𝑠𝜕𝑥𝜕𝑛𝜕𝑦𝜕𝑠𝜕𝑦=𝜕𝑛𝑣𝑥||𝑣𝑠||𝑣𝑦||𝑣𝑠||𝑣𝑦||𝑣𝑠||𝑣𝑥||𝑣𝑠||,(3.1) where 𝑣𝑥,𝑣𝑦 are depth-averaged velocity components in the 𝑥 and 𝑦 directions, respectively. Because (3.1) is orthogonal matrix, the inverse of 𝐉 is equal to 𝐉𝑇, which results in |(𝑠,𝑛)=𝐉𝑇|(𝑥,𝑦). Thus, both (2.1) and (2.3) can then be transformed into equations with respect to fixed Cartesian coordinate system; by applying |(𝑠,𝑛)=𝐉𝑇|(𝑥,𝑦) to (2.3),𝑑𝐶=𝐉𝑑𝑡𝑇𝑇𝐃𝑐𝐉𝑇𝐶=𝑇𝐉𝐃𝑐𝐉𝑇𝐶.(3.2)𝐃𝑐 is transformed into the nodal dispersion coefficients with respect to Eulerian-Cartesian coordinates: 𝐉𝐃𝑐𝐉𝑇, where each component of 𝐉𝐃𝑐𝐉𝑇 can be written as𝐷𝑥𝑥=𝐷𝑛𝑛+𝐷𝑠𝑠𝐷𝑛𝑛𝑣𝑥2𝑣𝑠2𝐷,(3.3a)𝑥𝑦=𝐷𝑦𝑥=𝐷𝑠𝑠𝐷𝑛𝑛𝑣𝑥𝑣𝑦𝑣𝑠2𝐷,(3.3b)𝑦𝑦=𝐷𝑛𝑛+𝐷𝑠𝑠𝐷𝑛𝑛𝑣𝑦2𝑣𝑠2,(3.3c)where 𝐷𝑥𝑥,𝐷𝑥𝑦,𝐷𝑦𝑥, and 𝐷𝑦𝑦 are Cartesian components of the nodal dispersion coefficient tensor. As we describe for the commonly used CSFD model, (3.3a), (3.3b), and (3.3c) are the conventional way to determine nodal dispersion parameters by the longitudinal and transverse dispersion coefficients 𝐷𝑠𝑠 and 𝐷𝑛𝑛 constantly specified in a global domain.

The full dispersion coefficient tensor introduced by Fischer [2] also include 𝐷𝑛𝑠 and 𝐷𝑠𝑛 as cross-dispersion coefficients that can be specified by physical considerations as well as 𝐷𝑠𝑠 and 𝐷𝑛𝑛:𝐷𝑠𝑠1=0𝑣𝑠𝑧01𝜀(𝑧)𝑧0𝑣𝑠𝐷𝑑𝑧𝑑𝑧𝑑𝑧,𝑠𝑛1=0𝑣𝑠𝑧01𝜀(𝑧)𝑧0𝑣𝑛𝐷𝑑𝑧𝑑𝑧𝑑𝑧,𝑛𝑠1=0𝑣𝑛𝑧01𝜀(𝑧)𝑧0𝑣𝑠𝐷𝑑𝑧𝑑𝑧𝑑𝑧,𝑛𝑛1=0𝑣𝑛𝑧01𝜀(𝑧)𝑧0𝑣𝑛𝑑𝑧𝑑𝑧𝑑𝑧,(3.4) where is flow depth; 𝜀 is the vertical turbulent diffusion coefficient; 𝑣𝑠,𝑣𝑛 the vertical deviations of the point velocities with respect to depth-averaged velocities in the 𝑠 and 𝑛 directions, respectively. Because the tensor dispersion coefficients proposed by Fischer [2] is expressed in fixed coordinates, the same transformation as did in (3.3a), (3.3b), and (3.3c) is needed for components of SSFD tensor to be applied in a flow where the flow direction continuously changes as in a curved stream: the components of 𝐉𝐃𝐉𝑇 in SSFD model, which is analogous to 𝐉𝐃𝑐𝐉𝑇  of CSFD model, are derived as follows:𝐷𝑥𝑥=𝐷𝑠𝑠𝑣𝑥2𝑣𝑠2𝐷𝑠𝑛+𝐷𝑛𝑠𝑣𝑥𝑣𝑦𝑣𝑠2+𝐷𝑛𝑛𝑣𝑦2𝑣𝑠2,𝐷𝑥𝑦=𝐷𝑠𝑠𝐷𝑛𝑛𝑣𝑥𝑣𝑦𝑣𝑠2+𝐷𝑠𝑛𝑣𝑥2𝑣𝑠2𝐷𝑛𝑠𝑣𝑦2𝑣𝑠2,𝐷𝑦𝑥=𝐷𝑠𝑠𝐷𝑛𝑛𝑣𝑥𝑣𝑦𝑣𝑠2𝐷𝑠𝑛𝑣𝑦2𝑣𝑠2+𝐷𝑛𝑠𝑣𝑥2𝑣𝑠2,𝐷𝑦𝑦=𝐷𝑠𝑠𝑣𝑦2𝑣𝑠2+𝐷𝑠𝑛+𝐷𝑛𝑠𝑣𝑥𝑣𝑦𝑣𝑠2+𝐷𝑛𝑛𝑣𝑥2𝑣𝑠2.(3.5) Obviously, (3.5) reduces to (3.3a), (3.3b), and (3.3c) if 𝐷𝑠𝑛,𝐷𝑛𝑠 are assumed to be zero. Through (3.5), the coefficients of the SSFD tensor in (𝑠,𝑛)-coordinates are transformed into nodal dispersion coefficients in Eulerian (𝑥,𝑦)-coordinates (Figure 2).

Using the transformed component of the nodal CSFD and SSFD coefficient tensor with respect to the global coordinate system, we present expanded the Cartesian forms of (2.1) and (2.3) as follows:𝜕𝐶+𝜕𝜕𝑡𝑣𝑥𝐶+𝜕𝜕𝑥𝑣𝑦𝐶=𝜕𝜕𝑦𝐷𝜕𝑥𝑥𝑥𝜕𝐶𝜕𝑥+𝐷𝑥𝑦𝜕𝐶+𝜕𝜕𝑦𝐷𝜕𝑦𝑦𝑥𝜕𝐶𝜕𝑥+𝐷𝑦𝑦𝜕𝐶.𝜕𝑦(3.6) For complex geometries of curved streams with irregular boundaries, an unstructured grid with the finite element or finite volume method is more useful than a numerical method with a structure grid. Therefore, we solve (3.6) by the finite element model established in the Eulerian coordinate system [16]. In this model, the Petrov-Galerkin approximation with a bilinear shape function is applied for spatial discretization, and the Crank-Nicolson type of discrete time marching is used for transient term. For each nodal point, the components of the dispersion tensor are determined by (3.3a), (3.3b), and (3.3c) in the CSFD model and by (3.5) in the SSFD model.

4. Test for the Direction of Principal Dispersion

In order to observe the oblique direction of the principal dispersion with respect to the longitudinal streamline, a solute mixing in uniform oscillatory flow was simulated by the established CSFD and SSFD models. The direction of the oscillatory flow was varied at 𝜃=0, 30°, and 45° counterclockwise with respect to the 𝑥-axis, to test the applicability of the derived (3.5) for various flow directions. Each result of the SSFD model was compared with that of the corresponding CSFD model with the same 𝐷𝑠𝑠 and 𝐷𝑛𝑛. The velocity of the uniform oscillatory flow was defined by a cosine curve with respect to time as𝑣𝑥,𝑣𝑦=𝑣𝑠cos𝜙𝑡cos𝜃,𝑣𝑠cos𝜙𝑡sin𝜃,(4.1) where 𝑣𝑠=0.25 m/s and 𝜙=2π/(12 h). In this flow field, the concentration distribution with a two-dimensional Gaussian profile was given as the initial concentration at the center of a 20 km × 20 km rectangular computational domain; this is the conventional test condition for the advection-dispersion problem [17]. The initial Gaussian distribution was obtained by the analytic solution after 1 day of pure diffusion of point mass 𝑀=5×104 kg/m with constant isotropic diffusion coefficient 𝑘= 5 m2/s:𝐶𝑀(𝑥,𝑦,0)=𝑥4𝜋𝑘𝜏exp2𝑦4𝑘2||||4𝑘𝜏=1day.(4.2) By (4.2), the initial peak concentration was 9.21 ppm. A rectangular element was chosen as a finite element grid of the test case due to its simplicity. The grid and time step sizes were determined as Δ𝑥=Δ𝑦=1km and Δ𝑡=900sec. The maximum Courant number was (𝑣𝑠Δ𝑡)/Δ𝑥=0.225, which guaranteed a stable solution. The initial concentration was spread under oscillatory flow conditions over 960 time steps, which covered a period of 6 days. Diagonal components of the SSFD and CSFD coefficient tensors of this case were arbitrarily determined as 𝐷𝑠𝑠=10 and 𝐷𝑛𝑛=1 m2/s. With these longitudinal and transverse coefficients, the off-diagonal components of the SSFD coefficient were determined as 𝐷𝑠𝑛=𝐷𝑛𝑠=3.125 m2/s, following the dispersion tensor for the application example in Fischer [2]:𝐃=2𝜀𝑈𝑜21205𝑈𝑜𝑉𝑜1925𝑈𝑜𝑉𝑜𝑉192𝑜212,(4.3) where 𝑈𝑜,𝑉𝑜 are magnitude of the longitudinal and transverse velocity deviation defined in Figure 3. The dispersion coefficient in (4.3) was derived from the approximated mean velocity profile of a continental shelf, as given in Figure 3. When we assume a skewed shear flow structure for secondary currents as shown in Figure 3, the intensity of the secondary current is maximized for the same transverse velocity deviation. The diagonal dispersion coefficient tensor with eigenvalues of 𝐃 is 𝐃=10.979000.021 m2/s by (2.8) and (2.9), and the coefficient tensor of the corresponding CSFD model is 𝐃𝑐=10001 m2/s.

Figure 3: Representation of the shear flow on the continental shelf of the middle Atlantic bight [2].

In Figure 3, the concentration distributions of the SSFD model 𝐶(𝐃,𝑡) at 𝑡=6 days are compared to those of the corresponding CSFD model 𝐶(𝐃𝑐,𝑡). In the results from the CSFD model, which are given in Figures 4(a), 4(c), and 4(e), the direction of the principal dispersion coincided with the oscillatory flow direction of each case. However, the differences in the symmetric axes of the concentration by SSFD model and the corresponding oscillatory flow direction were found to be 𝜓=17.4, as shown in Figures 4(b), 4(d), and 4(f). This difference confirmed the eigenvalue analysis results given in the previous section. The directions of the eigenvectors of 𝐃=103.1253.1251 m2/s were computed (2.11) to be oriented toward 17.4° and −72.6° with respect to the 𝑠 axis.

Figure 4: Concentration distribution in uniform oscillatory flow (unit: ppm).

The results in Figure 4 also show the applicability of (3.5); the coordinate transformation of dispersion tensor by (3.5) successfully introduced the full SSFD coefficient tensor in the numerical grid of Eulerian coordinates. When we consider solute transport in complex flow with continuously changing flow direction such as meandering streams, the whole problem is effectively solved in one domain based on the Eulerian frame of reference using (3.5). As presented in the next section, two flows with secondary currents case studies for application are considered. Through those examples, other aspects of the SSFD model—the increase in peak concentration and eccentricity of the tracer ellipse—were examined.

5. Application in Flows with Secondary Currents

In order to investigate the performance of the SSFD model in a flow field with secondary currents, an example case similar to the classic teacup experiment was considered. First, we assumed a solid-body rotation of water in a coaxial cylindrical container, as in Figure 5, to produce a so-called forced vortex. In this kind of fluid rotation, there is a pressure gradient from the perimeter toward the center. When we stop the rotation of the container abruptly, this pressure gradient coupled with the slower speed near the bottom boundary layer causes the secondary flow that makes the boundary layer spiral inward to the axis of circulation. Except for near the side wall and bottom, the fluid continues to rotate as before. Thus, a rotating flow field of vortex with a constant angular velocity can be assumed for a while in the container. We assume that a certain passive solute is dropped as an instantaneous point source as soon as the container is stopped. For this case, both dispersion models were applied to observe the difference in modeling results: SSFD model accounted for the effect of the secondary flow and CSFD neglected the secondary flow effect.

Figure 5: Rotating water in coaxial cylindrical container.

The outer and inner cylinders had radii of 10 and 3 m, respectively, and the container rotated at an angular speed of 𝜔=2𝜋/(3min). The plan view of the container and grid for the simulation are presented in Figure 6; the grid sizes in the radial and tangential directions were 0.5 m and 𝜋𝑟/40, respectively, where 𝑟 is the distance from the axis of rotation. The simulation was performed by the SSFD model with an arbitrarily determined 𝐃=0.010.0020.0020.001 m2/s of which the magnitude of the off-diagonal component was smaller than those of (4.3). The negative sign in the off-diagonal components of 𝐃 was determined by (3.3b) and (3.3c) and the circulating direction of the secondary current, which was negative 𝑛-direction in the upper part of the flow and positive 𝑛-direction in the lower part. The dispersion coefficient tensor of the corresponding CSFD model was 𝐃𝑐=0.01000.001 m2/s. An initial concentration of 100 was dropped at the nodal point (𝑥,𝑦)=(6.5m,0).

Figure 6: Grid system for coaxial cylindrical container.

The simulation results at 𝑡=36,72,,180sec for both CSFD and SSFD model are given in Figure 7. Figure 7(a) shows that the curvilinear axis of the principal dispersion direction in 𝐶(𝐃𝑐,𝑡) coincided with the circular streamline of the rotating flow. However, in Figure 7(b), the result of 𝐶(𝐃,𝑡) shows that the direction of the major dispersion axis spiraled outward to the perimeter of the outer cylinder due to the rotation of the symmetric axis of the tracer cloud. As pointed out in Section 2, the peak concentration of 𝐶(𝐃,𝑡) was computed to be larger than that of 𝐶(𝐃𝑐,𝑡); when the tracer cloud rotated for one complete round at 180 s, 𝐶𝑝(𝐃𝑐,𝑡) and 𝐶𝑝(𝐃,𝑡) were founded to be 3.305 and 3.851, respectively. Along with the larger peak concentration, the SSFD model increased the eccentricity of the ellipse of the tracer cloud. The tracer clouds in Figure 7(b) are longer and narrower than those in Figure 7(a).

Figure 7: Concentration distribution in coaxial cylindrical container.

For another example case of flow with secondary currents, the dispersion problem in a strongly curved channel with secondary flow was solved by the SSFD and CSFD modes. The famous experiments provided by Rozovski [18] were chosen as an example problem, in which the velocity field was measured along the bend of a U-shaped laboratory channel with a rectangular cross-section. The width of the channel was 0.8 m with inner and outer bend radii of 0.4 and 1.2 m, respectively, and the water flowed on the zero-slope bottom with a mean velocity of 0.25 m/s for his experimental case Run 8.

The velocity field in the whole domain was reproduced by the commonly used depth-averaged flow analysis model RMA-2 in TABS-MD [19]. Figure 8 compares the computed velocities against the measured data of Rozovski [18] at five cross-sections in the bend and shows that the computational result described the flow field along the channel bend well. When water flowed along the bend of the rectangular channel, the usual path of velocity maxima is located near the inner bend at the entrance to the bend and shift to near the outer bank at the exit as shown in the measurement of longitudinal velocity in curved channels [2022] as observed in Figure 8. This is due to the increase of the water elevation at the outer bank and the decrease at the inner bank by centrifugal forces in curved channels [23]. The principal flow following the velocity maxima along the bend is accompanied by the secondary flow by which the path line is partly downstream and partly across the channel from the outer bank toward the inner bank at the bottom. To observe the solute spreading around the bend, the initial concentration of 100 was defined at the centered point of bend entrance (𝑥,𝑦)=(6m,5.75m). A dispersion coefficient tensor with 𝐃=0.0050.0020.0020.001 m2/s was arbitrarily selected for the SSFD model, and 𝐃𝑐=0.005000.001 m2/s was used for the corresponding CSFD model.

Figure 8: Comparison between computed and measured flow fields in Rozovski’s channel.

Concentration distributions at 𝑡=4 and 8 s were presented in Figures 9 and 10, respectively. According to Figures 9 and 10, advection by a lateral nonuniform flow coupled with the skewed shear flow dispersion causes the remarkable change in the shape of tracer cloud compared to the CSFD model: the elliptical curves of the equiconcentration computed by the CSFD model maintained the elongated shape along the streamline until the total mass exited the bend, whereas those of the SSFD model were much shorter. These results seem to be opposite to the elongated shape of the tracer cloud by the SSFD model shown in the previous example of a force vortex. The mechanism for the generation of the clustered concentration by the SSFD model in this example is explained as follows. When the SSFD coefficient tensor is applied, the direction of the major dispersion axis becomes oblique clockwise with respect to the streamlines; this is the typical mechanism for skewed shear flow dispersion as well as the previous example. However, in contrast to the force vortex, the nonuniform longitudinal velocities of the channel bend shift each particle in the skewed tracer cloud with different speeds along each streamline; the particles near the inner bank move faster than the particles near the outer bank. As a result, the dispersed particles in the tracer cloud skewed about streamlines are centered into the mean displacement of the moving fluids. In particular, at 𝑡=4 s, the shape of the equiconcentration curve appeared to be rounded and concentrated.

Figure 9: Concentration distribution in Rozovski’s channel at 𝑡=4 s.
Figure 10: Concentration distribution in Rozovski’s channel at 𝑡=8 s.

The rotation of the major dispersion axis in the SSFD model and the nonuniform advection along the curved streamline had a combined effect on the increase in the peak concentration. The peak concentrations at 𝑡=8 s, as given in Figure 10, were found to be 0.60 and 0.78 for the CSFD and SSFD models, respectively. This shows that the application of the CSFD model in the curved channel with a strong secondary current—instead of the SSFD model—may underestimate the peak concentration of the transported pollutant from upstream. Although the secondary current is known to activate transverse dispersion and increase the dilution effect, the results in this study indicate that the skewed vertical shear profile of the secondary current may offset the enhancement of dilution caused by the large transverse dispersion.

6. Conclusion

In this study, it was proposed that the SSFD coefficient tensor should be applied for 2D passive solute transport modeling in the flow with secondary current because of its vertically skewed shear flow structure. Mathematical analysis of eigenvalue problem pointed out several significant effects of the off-diagonal terms of dispersion tensor: the rotation of principal direction of dispersion with respect to the streamline, the increase of peak concentration, and the change in eccentricity of elliptical concentration. To apply full dispersion coefficient tensor defined in a stream-wise curvilinear coordinate system to the numerical model on the Eulerian-Cartesian coordinates, transformation relationship was derived with given depth-averaged velocity field. With the derived transformation equation, 2D numerical model was established with finite element method on the Eulerian coordinate system. Numerical tests show that the coordinate transformation relationship derived in this study successfully introduced the SSFD coefficient tensor in the numerical grid of Eulerian coordinates. It was also shown that there is a possibility of overestimation in dilution of pollutant if CSFD model was applied instead of SSFD model in the dispersion process affected by secondary currents. The conventional 2D solute mixing modules equipped in the various hydrodynamic modeling packages are expected to predict more reliable mixing patterns of pollutants by including off-diagonal terms as in SSFD model, when it is applied to flow field with secondary currents.


𝐶:Depth-averaged concentration
𝐶𝑝:Peak concentration
𝐷𝑠𝑠,𝐷𝑠𝑛,𝐷𝑛𝑠, and 𝐷𝑛𝑛:Components of the full dispersion coefficient tensor defined in stream-wise curvilinear coordinates
𝐷𝑠𝑛:Mean value of 𝐷𝑠𝑛 and 𝐷𝑛𝑠
𝐷𝑥𝑥,𝐷𝑥𝑦,𝐷𝑦𝑥, and 𝐷𝑦𝑦:Components of the nodal dispersion coefficient tensor in the Eulerian-Cartesian coordinates
𝐃:Matrix notation of SSFD coefficient tensor in curvilinear coordinates
𝐃𝑐:Matrix notation of CSFD coefficient tensor in curvilinear coordinates
𝐃:Symmetric version of SSFD coefficient tensor, which takes 𝐷𝑠𝑛 as off-diagonal entries
𝐃:Principal dispersion coefficient tensor, which takes eigenvalues 𝜆1,𝜆2 as entries
:Flow depth
𝐉:Jacobian matrix for coordinate transformation
𝑘:Isotropic diffusion coefficient
𝑀:Total mass of tracer
𝑛:Axis normal to the streamline in the stream-wise curvilinear coordinate system
𝑠:Axis along the streamline in the stream-wise curvilinear coordinate system
Δ𝑡:Discretization size in time
𝑈𝑜,𝑉𝑜:Magnitude of longitudinal and transverse velocity deviations in Figure 3
𝑣𝑠,𝑣𝑛:Horizontal velocities on the 𝑠-,𝑛-axis
𝑣𝑠:Longitudinal depth-averaged velocity
𝑣𝑠,𝑣𝑛:Vertical deviations of the point velocities with respect to depth-averaged 𝑣𝑠,𝑣𝑛
𝑣𝑥,𝑣𝑦:Ddepth-averaged velocity components in the 𝑥 and 𝑦 directions
𝑥,𝑦:Axes of the Eulerian-Cartesian coordinate system
𝑥,𝑦:Axes with identical directions as the principal axes of 𝐃
Δ𝑥:Discretization size in the 𝑥 direction
Δ𝑦:Discretization size in the 𝑦 direction
𝐗:Matrix that takes the eigenvectors of 𝐃 as column vectors
𝑧:Axis of the vertical direction
𝜀𝑛:Transverse turbulent diffusion coefficient
𝜀:Vertical turbulent diffusion coefficient
𝜙:Angular frequency of oscillatory flow
𝜃:Angle of oscillatory direction with respect to the 𝑥-direction
𝜆1,𝜆2:Eigenvalues of 𝐃
𝜔:Angular speed of rotation of coaxial cylindrical container
𝜓:Angle of counterclockwise rotation from the 𝑠 axis of the principal dispersion axes.


  1. G. I. Taylor, “Dispersion of soluble matter in solvent flowing slowly through a tube,” Proceedings of the Royal Society A, vol. 219, pp. 186–203, 1953. View at Publisher · View at Google Scholar
  2. H. B. Fischer, “On the tensor form of the bulk dispersion coefficient in a bounded skewed shear flow,” Journal of Geophysical Research, vol. 83, no. C5, pp. 2373–2375, 1978. View at Publisher · View at Google Scholar
  3. C. W. Almquist and E. R. Holley, “Transverse mixing in meandering laboratory channels with rectangular and naturally varying cross sections,” Tech. Rep. CRWR-205, University of Texas, Austin, Texas, 1985. View at Google Scholar
  4. K. O. Baek, I. W. Seo, and S. J. Jeong, “Evaluation of dispersion coefficients in meandering channels from transient tracer tests,” Journal of Hydraulic Engineering, vol. 132, no. 10, pp. 1021–1032, 2006. View at Publisher · View at Google Scholar · View at Scopus
  5. J. B. Boxall and I. Guymer, “Analysis and prediction of transverse mixing coefficients in natural channels,” Journal of Hydraulic Engineering, vol. 129, no. 2, pp. 129–139, 2003. View at Publisher · View at Google Scholar · View at Scopus
  6. J. B. Boxall, I. Guymer, and A. Marion, “Transverse mixing in sinuous natural open channel flows,” Journal of Hydraulic Research, vol. 41, no. 2, pp. 153–165, 2003. View at Publisher · View at Google Scholar · View at Scopus
  7. Y. C. Chang, Lateral mixing in meandering channels, Ph.D. thesis, University of Iowa, Iowa City, Iowa, USA, 1971.
  8. B. G. Krishnappan and Y. L. Lau, “Transverse mixing in meandering channels with varying bottom topography,” Journal of Hydraulic Research, vol. 15, no. 4, pp. 351–371, 1977. View at Publisher · View at Google Scholar · View at Scopus
  9. A. Marion and M. Zaramella, “Effects of velocity gradients and secondary flow on the dispersion of solutes in a meandering channel,” Journal of Hydraulic Engineering, vol. 132, no. 12, pp. 1295–1302, 2006. View at Publisher · View at Google Scholar · View at Scopus
  10. E. Kreyszig, Advanced Engineering Mathematics, John Wiley & Sons, 9th edition, 2006.
  11. V. Alavian, “Dispersion tensor in rotating flows,” Journal of Hydraulic Engineering, vol. 112, no. 8, pp. 771–777, 1986. View at Publisher · View at Google Scholar · View at Scopus
  12. D. Liang, X. Wang, R. A. Falconer, and B. N. Bockelmann-Evans, “Solving the depth-integrated solute transport equation with a TVD-MacCormack scheme,” Environmental Modelling and Software, vol. 25, no. 12, pp. 1619–1629, 2010. View at Publisher · View at Google Scholar · View at Scopus
  13. M. Piasecki and N. D. Katopodes, “Identification of stream dispersion coefficients by adjoint sensitivity method,” Journal of Hydraulic Engineering, vol. 125, no. 7, pp. 714–724, 1999. View at Publisher · View at Google Scholar · View at Scopus
  14. R. W. Preston, “The representation of dispersion in two-dimensional shallow water flow,” Report TPRD/L/278333/N84, Central Electricity Research Laboratories, Leatherhead, England, 1985. View at Google Scholar
  15. A. Younes and P. Ackerer, “Solving the advection-diffusion equation with the Eulerian-Lagrangian localized adjoint method on unstructured meshes and non uniform time stepping,” Journal of Computational Physics, vol. 208, no. 1, pp. 384–402, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  16. I. W. Seo, M. E. Lee, and K. O. Baek, “2D modeling of heterogeneous dispersion in meandering channels,” Journal of Hydraulic Engineering, vol. 134, no. 2, pp. 196–204, 2008. View at Publisher · View at Google Scholar · View at Scopus
  17. R. T. Cheng, V. Casulli, and S. N. Milford, “Eulerian-Lagrangian solution of the convection-dispersion equation in natural co-ordinates,” Water Resources Research, vol. 20, no. 7, pp. 944–952, 1984. View at Publisher · View at Google Scholar · View at Scopus
  18. I. L. Rozovski, Flow of Water in Bends of Open Channels, Israel Program for Scientific Translation, Jerusalem, Israel, 1965.
  19. I King, B. P. Donnell, J. V. Letter, W. H. McAnally, W. A. Thomas, and C. LaHatte, “RMA2 users guide 4.5x,” US Army, Engineer Research and Development Center, Waterways Experiment Station, Coastal and Hydraulics Laboratory, 2006.
  20. Y. C. Chang, Lateral mixing in meandering channels, Ph.D. thesis, University of Iowa, Iowa City, Iowa, USA, 1971.
  21. H. J. DeVriend, “Flow measurements in a curved rectangular channel,” Internal Report number 9-79, Delft University of Technology, Department of Civil Engineering, Laboratory of Fluid Mechanics, 1979. View at Google Scholar
  22. S. Fukuoka, Longitudinal dispersion in sinuous channels, Ph.D. thesis, University of Iowa, Iowa City, Iowa, USA, 1971.
  23. A. O. Demuren and W. Rodi, “Calculation of flow and pollutant dispersion in meandering channels,” Journal of Fluid Mechanics, vol. 172, pp. 63–92, 1986. View at Publisher · View at Google Scholar · View at Scopus