Influence of Secondary Currents on Solute Dispersion in Curved Open Channels

Abstract

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 𝐷_𝑛𝑛 [3–9]. 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 𝐃_𝑐𝑑𝐶𝑑𝑡=∇⋅𝐃^𝑐∇𝐶=∇^𝑇𝐃_𝑐∇𝐶,𝐃_𝑐=⎡⎢⎢⎣𝐷^𝑠𝑠0𝐷0_𝑛𝑛⎤⎥⎥⎦.(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 𝐗^𝑇=𝐗⁻¹. According to the theory of orthogonal eigenvectors of a specified symmetric matrix, a certain diagonal matrix 𝐃^ is obtained by the following relationship:𝐃=𝐗𝐃^𝐗⁻¹=𝐗𝐃^𝐗^𝑇.(2.5) It should be noted that 𝐃^ is a similar matrix of 𝐃 by orthogonal transformation: 𝐃^=𝐗⁻¹𝐃𝐗. 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:𝜆₁,𝜆₂=𝐷^𝑠𝑠+𝐷_𝑛𝑛2±𝐷^𝑠𝑠−𝐷_𝑛𝑛2²+𝐷_2𝑠𝑛.(2.8) The obtained eigenvalues in (2.8) are the diagonal components of 𝐃^:𝐃^=⎡⎢⎢⎣𝜆¹0𝜆0₂⎤⎥⎥⎦(2.9) by the calculation of 𝐃^=𝐗⁻¹𝐃𝐗. 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 𝜆₁ and 𝜆₂ 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𝐷_𝑠𝑛²+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𝜋ℎ𝑡√𝜆₁𝜆₂exp−𝑥^24𝜆₁𝑡−4𝑦𝜆^2₂𝑡,(2.12) where 𝑀 is the total mass of tracer. With this solution, the equiconcentration curves are ellipses with √𝜆₁ and √𝜆₂ as the lengths of semimajor and minor axes. Using (2.12), the peak concentration at time 𝑡 is inversely proportional to √𝜆₁𝜆₂, which is given by√𝜆₁𝜆₂=𝐷_𝑠𝑠𝐷_𝑛𝑛−𝐷_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 𝜆₁,𝜆₂ 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:max𝜆₁,𝜆₂=𝐷^𝑠𝑠+𝐷_𝑛𝑛2+𝐷^𝑠𝑠−𝐷_𝑛𝑛2²min𝜆₁+𝐷_2𝑠𝑛≥max𝐷_𝑠𝑠,𝐷_𝑛𝑛,,𝜆₂=𝐷^𝑠𝑠+𝐷_𝑛𝑛2−𝐷^𝑠𝑠−𝐷_𝑛𝑛2²+𝐷_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 [11–15]. 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×10⁴ kg/m with constant isotropic diffusion coefficient 𝑘= 5 m²/s:𝐶(𝑥,𝑦,0)=𝑀4𝜋𝑘𝜏exp−4𝑥𝑘²−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 m²/s. With these longitudinal and transverse coefficients, the off-diagonal components of the SSFD coefficient were determined as 𝐷_𝑠𝑛=𝐷_𝑛𝑠=3.125 m²/s, following the dispersion tensor for the application example in Fischer [2]:𝐃=ℎ²𝜀⎡⎢⎢⎢⎢⎣𝑈^𝑜21205𝑈^𝑜𝑉_𝑜5𝑈_𝑜192𝑉_𝑜192𝑉12^𝑜2⎤⎥⎥⎥⎥⎦,(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} m²/s by (2.8) and (2.9), and the coefficient tensor of the corresponding CSFD model is 𝐃_𝑐=¹⁰⁰⁰¹ m²/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} m²/s were computed (2.11) to be oriented toward 17.4° and −72.6° with respect to the 𝑠 axis.

(a) CSFD model, 𝑡=6 day (in 0 degree flow)

(b) SSFD model, 𝑡=6 day (in 0 degree flow)

(c) CSFD model, 𝑡=6 day (in 30 degree flow)

(d) SSFD model, 𝑡=6 day (in 30 degree flow)

(e) CSFD model, 𝑡=6 day (in 45 degree flow)

(f) SSFD model, 𝑡=6 day (in 45 degree flow)

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.01−0.002−0.0020.001} m²/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 m²/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).

(a) CSFD

(b) SSFD

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 [20–22] 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.005−0.002−0.0020.001} m²/s was arbitrarily selected for the SSFD model, and 𝐃_𝑐=^0.005000.001 m²/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.

(a) CSFD

(b) SSFD

Figure 9 

Concentration distribution in Rozovski’s channel at 𝑡=4 s.

(a) CSFD

(b) SSFD

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.

Nomenclature

𝐶:	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
𝑡:	Time
Δ𝑡:	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.