Longitudinal Dispersion in Straight Open Channels: Anomalous Breakthrough Curves and First-Order Analytical Solution for the Depth-Averaged Concentration

: A ﬁrst-order analytical solution is proposed for the actual depth-averaged concentration of tracers in shallow river ﬂows in the presence of large Peclet numbers (deﬁned as the ratio of section-averaged velocity times channel width to turbulent diffusion coefﬁcient). The solution shows how complete transverse mixing is never achieved due to the typical shape of the velocity and diffusion coefﬁcient proﬁle, which alternatively tend—depending on the downstream location of the cross-section—to concentrate the mass at the centre or at the boundaries of the cross-section itself. The ﬁrst-order analytical solution proves to be consistent with the results of Lagrangian numerical simulations based on real-ﬁeld input data, which show how the solute mass breakthrough curves always exhibit anomalous behaviour and a considerable and persistent delay when compared with those that are analytically obtained by assuming a truly one-dimensional process.


Introduction
Scalar mixing and spreading processes are fundamental in many geophysical systems and engineering applications [1,2]. The increasing need for scientific tools aimed at predicting, monitoring and controlling chemical contamination and sediment transport in natural channels has recently inspired extensive studies about different aspects of dispersion in stream flows, from the analysis of the effects of heterogeneous channel morphology and asymmetry to that of release conditions and backwater flows, and even to suspended load dynamics. Transient analytical solutions of the advection-dispersion equation (ADE) for dissolved chemicals and fine sediments in open channels were given by the authors: (1) in the case of deterministic depth-averaged velocity distributions and different initial conditions based on Aris' method of moments [3]; (2) in the case of randomly uniform flows in the presence of morphologic heterogeneity represented by short-range correlations based on a stochastic Lagrangian approach [4][5][6]; (3) and in terms of second-order concentration statistics based on a stochastic Eulerian approach [7]. An analytical-numerical model was elaborated to solve ADE in a stochastic Lagrangian framework in the case of backwater flows in large rivers, highlighting the occurrence of a global re-densification of the cloud and a persistent transverse non-uniformity [8].
A stochastic Lagrangian numerical model was developed to handle ADE in the case of depth-dependent turbulent diffusion coefficient, highlighting the occurrence of concentration boundary side-pockets [9].
In the last few years, several studies (e.g., [10]) have described solute mixing in terms of Shannon's information entropy [11][12][13][14]. A steady-state uniform-flow solution was recently given by the authors for the concentration of small-size suspended sediments [15] by modeling the cross-sectional velocity distribution based on the principle of maximum information entropy and by investigating the effects of velocity-dip position (e.g., [16,17]). The great interest in the investigation of pollutant/sediment dispersion and the correct quantification of mixing in river flows is justified by its crucial role in implementing effective natural or engineered remediation strategies and in performing risk assessment analysis (e.g., [18]).
As a matter of fact, the classic one-dimensional advection-dispersion equation with constant coefficients and its instantaneous point source solution were always found to be inadequate for the interpretation of field measurements of chemical dispersion. The typical explanation given to justify this discrepancy is that the tracer's travel time was not long enough, implicitly meaning that section-averaged concentration would have been well-described by Taylor's [19] model once the theoretical asymptotic regime of dispersion had finally been achieved.
Nevertheless, all the data available in the literature highlight a persistent skewness of the distribution and an increasing peripheral drift of the solute clouds that seem to be incompatible with a pre-asymptotic phase of the process [20][21][22][23][24], the reason for the failure of Taylor's asymptotic theory is likely related to the permanent cross-sectional non-uniformity of solute mass. Several attempts were made to relate the persistent concentration skewness to the presence of dead zones induced by morphological singularities as edges and corners or subsurface retention phenomena (e.g., [25][26][27]). Ahmad analysed the effect of variable transverse mixing on stream-wise dispersion in terms of longitudinal mixing length for different initial conditions [28] and a classic large-time Fickian behaviour. Later, based on 2-D Lagrangian numerical simulations and the investigation of the time behaviour of cloud moments, Pannone found that tracer dispersion in large rivers is singularly affected by (depth-dependent) transverse diffusion [9]. Specifically, the interplay between zero-bank velocity and the zero-bank transverse mixing coefficient leads to the occurrence of persistent solute side pockets and an anomalous transient behaviour of the spatial moments.
In the present paper, the peculiarities of longitudinal dispersion in large rivers in the presence of suitably variable transverse mixing is discussed through the numerical analysis of the breakthrough curves. Additionally, the permanent peripheral drift of solute mass and the persistence of cross-sectional incomplete mixing, even in the absence of morphological singularities and adsorption phenomena, are justified by a first order analytical solution of the depth-averaged advection-diffusion equation, which is pursued for large channel aspect ratios and large Peclet numbers. The scope of the study is: (a) showing that the incomplete transverse mixing is a ubiquitous feature of open channel flows; (b) providing a closed-form analytical solution that could serve as a fast predicting tool for the deviatory depth-averaged concentrations of chemicals subject to advection and dispersion along almost uniform reaches of steady wide river flows.

Formulation and First Order Analytical Solution
Let us consider the time-mean advection-diffusion equation that governs steady uniform turbulent flows in straight open channels, Equation (1): where x, y, and z are the orthogonal axes of a reference frame originating on the free surface at the right bank, c = c(x, y, z, t) is the time-mean concentration of tracer, u = u(y, z) is the time-mean longitudinal velocity, and ε x (y, z), ε y (y, z), and ε z (y, z) indicate longitudinal, transverse, and vertical turbulent diffusion coefficient, respectively. The longitudinal macro-dispersion coefficient, widely investigated in the literature to analyse the asymptotic macroscopic characteristics of the process (e.g., [3][4][5][6]9,[29][30][31]), appears in the one-dimensional transport equation derived from Equation (1) after section averaging, which is here performed according to Equation (2) [9]: where h(y) is the local depth, A the flow area, B the free-surface width, and H eq indicates the depth of the equivalent rectangular cross-section. When the variation of the flow depth along the transverse coordinate is very slow and the width is much larger than the equivalent depth, the depth-averaged diffusive deviatory fluxes are negligible if compared to those that are associated with the joint variability of u and c. Considering the no-flux condition at the banks, Equation (3) is obtained: where U is the section-averaged velocity, D x represents the longitudinal macro-dispersion coefficient expressed by Equation (4) combined with Equations (5) and (6): the overbar indicates depth-average and the prime the corresponding deviation. Equation (3) holds in the reasonable and well-documented hypothesis that the dispersive effect induced by the transverse heterogeneity of u and c dominates the corresponding effect induced by vertical heterogeneity and turbulent diffusion (e.g., [20,30,32]). Dispersion and dilution phenomena in large rivers can be analysed either by solving advection-diffusion equations like Equation (1) or Equation (3) (e.g., [7]), or by estimating the spatial moments of the equivalent discrete distributions of mass in terms of solute particle trajectories after [33]. The solution of Equation (3) for instantaneous point solute source at x = 0 reads: where t 0 indicates the initial longitudinal relaxation time, M total mass and H the cross-sectional average flow depth.
In the present study, in the hypothesis of shallow waters, and assuming that vertical mixing is much faster than its transverse and longitudinal counterparts [30], the problem is solved based on a 2-D numerical Lagrangian approach (e.g., [9]) in terms of downstream breakthrough curves, and based on a 2-D analytical Eulerian approach in terms of actual depth-averaged concentration. To our knowledge, this is the first attempt at finding an explicit analytical solution for the deviatory concentration (whose contribution is routinely analysed in an average sense by the corresponding variance).
The depth-averaged advection-diffusion equation is formulated as, Equation (8): where only the transverse diffusive term, whose order of magnitude is comparable with that of the stream-wise component [32], is maintained because of its crucial role in triggering longitudinal dispersion by the cross-sectional displacement of solute particles that enables the sampling of the whole velocity distribution. Subtracting Equation (3) from Equation (8) gives the governing equation for the deviatory concentration, Equation (9): An approximate solution of Equation (9) can be pursued in the case of straight reaches of large regular rivers in steady uniform-flow conditions. At the first order in the velocity/concentration deviations, Equation (10): where, Equation (11) [30]: u * indicates shear velocity, g is the acceleration due to gravity and i f channel slope. According to the above-mentioned hypotheses, Equation (10) will be solved here for large cross-sections and large related Peclet numbers, Equation (12): where, Equation (13): The adoption of suitable scales for length, velocity, time and concentration, respectively, Equation (14): leads to the following (Equation (15)) dimensionless version of Equation (10): where, Equation (16): and function Ψ has to be defined according to the assumed flow-section shape. When Pe is very large, its inverse δ can be adopted as the small parameter of a perturbation expansion of the solution, Equation (17): Substituting Equation (17) into Equation (15), one finds that, at order δ 0 , the equation to be solved is, Equation (18): while at order δ 1 one has, Equation (19): with initial condition (uniform cross-sectional strip injection at x = 0 and t = 0), Equation (20): Let us consider a parabolic flow-section shape, Equation (21): and a velocity profile expressed by the generalized Manning's formula [34] transformed into the equivalent depth-averaged velocity, Equation (22) [9]: where θ = θ(B/H) is a calibration coefficient accounting for the difference between average flow depth and hydraulic radius. In the case of large cross-sections and large aspect ratios (B>>H), θ~1 and H eq~H . Thus, Equation (23): and, Equation (24): In a longitudinal reference frame moving with the section-averaged velocity, Equation (25): Equations (23) and (24) respectively become, Equation (26): and, Equation (27): whose general solutions are, Equation (28): and, Equation (29): In dimensionless terms, the gradient of the section-averaged concentration (7) is given by Equation (30): In Equation (30) and in what follows, Pe M represents the macro-Peclet number Pe M = UB/D x . Equation (28) has the following closed-form solution, Equation (31): with, Equation (32): and, Equation (33): where "Erf" indicates Error Function (e.g., [35]). Equation (29) then yields Equation (34): with, Equation (35): and, Equation (36): where "Ei" indicates Exponential-Integral Function (e.g., [35]). See Figures 1-6 for a graphical illustration of the solution with reference to two different macro-Peclet numbers (1 and 0.01, respectively) assuming t 0 = 0.1. Specifically, Figure 1a,b display Y 0 ( y) and Y 1 ( y), respectively; Figure 2a,b display X 0 x * , t and X 1 x * , t , respectively, as a function of x * ; finally, Figure 3a,b display X 0 x * , t and X 1 x * , t , respectively, as a function of t. Note that Equations (31) and (34) consistently give zero at the section-averaged concentration maximum ( x * = 0) and at the section-averaged concentration tails ( x * = ±∞). Additionally, they would reproduce the typical two-peaked longitudinal structure detected in previous studies on concentration variance in heterogeneous velocity fields (e.g., [7,[36][37][38]). Indeed, since concentration spatial variance in constant-width stream flows by definition is, Equation (37): for example, at order zero, we have, Equation (38): or, Equation (39): Therefore, the longitudinal variability of the leading-order variance coincides with the longitudinal variability of function X 0 2 . As one can see from Figure 2a,b, function X 0 (like function X 1 ) is characterized by a maximum (positive) and a minimum (negative), with the same absolute value and the same distance from the origin. Its square is therefore characterized by two equal and symmetrical positive peaks. Their distance is, Equation (40): Both c 0 and c 1 are obtained as the product of a function of x * and t ( X 0 or X 1 ) and a function of y ( Y 0 or Y 1 ). Where X 0 or X 1 is positive, the sign of c 0 or c 1 as a function of y coincides with the sign of Y 0 or Y 1 . Where X 0 or X 1 is negative, the sign of c 0 or c 1 as a function of y is opposite of Y 0 or Y 1 's. Thus, for example, where X 0 is negative (and that happens for x * = (x − Ut)/B < 0, that is, in the raising branch of the section-average Gaussian bell), we have c 0 < 0 in the central part of the cross-section and c 0 > 0 near the side boundaries. That means that, at order δ 0 , while in the central part of the cross section total concentration is lower than the average, at the side boundaries it is higher. Vice versa for positive X 0 . That means that, while at the back of the section-average Gaussian wave (7) solute mass is concentrated at the boundaries (the slower part of flow field), at the front of the section-average Gaussian wave (7) solute mass is concentrated at the center (the faster part of flow field). Specifically, along the transverse direction, c 0 exhibits two (positive) maxima at the banks and a (negative) minimum at the channel axis for x * < 0, where X 0 < 0; it is > 0 for y > 0.5 + 1/2 √ 3 and < 0 for y < 0.5 − 1/2 √ 3. Conversely, c 0 exhibits two (negative) minima at the banks and a (positive) maximum at the channel axis for x * > 0, where X 0 > 0; it is > 0 for 0.5 − 1/2 √ 3 < y < 0.5 + 1/2 √ 3 and < 0 for y < 0.5 − 1/2 √ 3 and y > 0.5 + 1/2 √ 3 (see Figure 1a combined with Figure 2a). Unlike the leading order component, the first-order diffusive contribution is characterized by four zeros ( y = 0, y = 1, y ∼ = 0.28, y ∼ = 0.72). Thus, depending on the local sign of function X 1 (> 0 for x * > 0 and < 0 for x * < 0), it has its maximum at the centre of the cross-section or within more peripheral layers, for 0 < y < 0.28 and 0.72 < y < 1 (see Figure 1b combined with Figure 2b). This asymmetry would justify an anomalous delay in the breakthrough curves between first particle arrivals and last particle arrivals. It should also be noted that, whereas c 0 , viewed as a function of time, tends to a constant proportional to Pe M , c 1 tends to increase linearly. In term of the deviatory first order concentration solution (according to Equation (17)), this continuous increase is partially balanced by the small value of δ = 1/Pe (see Figure 3).

Numerical Lagrangian Simulations
A numerical investigation of dispersion within straight reaches of open channel flows was already discussed in previous authors' work [9]. The study focused on the analysis of solute cloud spatial moments reconstructed by 2-D particle tracking based on input experimental data collected along six rivers of southern Italy [39]. The six study cases were equally distributed among compact rectangular (CR), bi-rectangular (BR), and triangular cross-sections (T), characterized by width-to-depth ratio B/H ranging between~15 and~56. At the time of the field survey, the six channels were characterized by relatively large bed grain size, with no suspended load and no relevant bed-forms. In order to work on a more dense and regular grid, a larger number of depth values was derived from field data by 3 to 5 level linear interpolation. Each river was discretized into a variable number of transverse elements (from 80 to 320). Each element was characterized by the suitable value of u and ε y according to Equations (22) and (11), respectively. A high number of particles (NP = 200,000) was instantaneously released at x = 0 and y = B/2 and then tracked based on Equations (41) and (42): with ε y = ε y [Y(t)] (depth-dependent transverse mixing) or, for the sake of comparison, ε y = ε y = const (uniform mixing). In Equations (41) and (42), G ζ (0, 1) represents the generic element of a normal distribution (i.e., a zero-mean and unit-variance Gaussian), and the length of the time step ∆t was selected in such a way as to minimize the number of diffusive multiple-element jumps while keeping the computational time within reasonable limits. Indeed, to avoid tracking the particle with wrong velocity values, ∆Y must not exceed the transverse dimension of the computational grid ∆n, Equation (43): Thus, Equation (44): Since, Equation (45): we had, Equation (46): With ζ = 5, G 5 (0, 1) ∼ = 10 −6 . The total number of tracked particles was NP = 200,000. Setting, Equation (47): we were certain that only 0.2 of 200,000 displacements had the probability to be performed with the wrong velocity at each time step.
In the present study, we resorted to the same methodology for the numerical reconstruction of the breakthrough curves referred to in three of the six study cases (Table 1) in the presence of uniform and non-uniform turbulent mixing at four different downstream cross-sections (x 1 = 2Ut d , x 2 = 4Ut d , x 3 = 6Ut d , x 4 = 9Ut d , with t d = B 2 /ε y indicating the average diffusive time). It is well known that the breakthrough curves are a measure of the percentage solute mass that has crossed a given section at a given time and allow us to establish whether, and to what extent, concentration distribution is affected by tails and decelerations. As a matter of fact, their use by-passes the evaluation of the actual c(x, y, t), which is analytically accessible only in approximate forms (an example was shown above) and is numerically affected by unavoidable errors due to discretization issues and virtual diffusion. Following the mentioned Lagrangian approach, we evaluated the cumulate of concentration by counting the percentage number of particles that, at the generic time t, were characterized by X(t) > x i . Figures 4a, 5a and 6a show the breakthrough curves for the three cases under investigation in conditions of uniform turbulent mixing; Figures 4b, 5b and 6b do the same in the presence of non-uniform mixing (for R1, R2 and R3, respectively). The time on the abscissa corresponds in every case to 2 average diffusive times. Finally, for the sake of comparison, Figure 7a-c display the theoretical cumulative distributions corresponding to the one-dimensional Gaussian in the presence of uniform mixing, Equation (48):    The longitudinal macro-dispersion coefficients D x used in Equation (48) (110, 220, and 105 m 2 /s, respectively, for the three cases under investigation), were consistently obtained as one-half of the time-derivative of the cloud central inertia moment S x , Equation (49) [9]: where, Equation (50): indicates longitudinal average trajectory. Note that the reconstruction of the theoretical breakthrough curves based on Equation (48) implicitly assumes a transversally uniform distribution, with c = 0. The comparison between the curves numerically obtained for uniform and non-uniform mixing reveals that the uniform approximation (which is in qualitative agreement with the cumulate of the section-averaged concentration (48)) yields totally misleading results. Indeed, in the case of non-uniform mixing, and regardless of the morphological typology of the cross-section, even the shape of the curves is anomalous and is characterized by a considerable delay with respect to the theoretical one. This considerable and persistent delay can find a very reasonable explanation on the basis of the first-order analytical solution presented in the previous section. As a matter of fact, as above mentioned, both the leading order component and first order diffusion-related component are characterized by (positive) maxima near the banks and (negative) minima at the centre of the cross-section for x * < 0, that is, at the back of the C-wave. Conversely, they are characterized by a (positive) maximum at the centre of the cross section and (negative) minima in the more peripheral layers for x * > 0, that is, at the front of the C-wave. That means that solute mass is concentrated in the faster areas of the cross-section at the front and in the slower areas of the cross-section at the back, thus producing the detected delay in the breakthrough curves. Additionally, and unlike the leading-order component of the solution c 0 , which tends to a constant different from zero for any x * = 0, ±∞, the first-order diffusion-related component c 1 linearly increases in time for any x * = 0, meaning that its effect gets stronger and stronger as time elapses. Thus, the first-order analytical model proposed by the present study and based on the analysis of depth-averaged concentration behaviour in the presence of a symmetric parabolic cross-section characterized by a large aspect ratio and large Peclet number, seems to be able to grasp the basic features of tracer dispersion also under the more general conditions simulated by the numerical analysis. Additionally, since in the presence of shallow waters (for which the problem of monitoring and predicting the fate of pollutant substances is even more serious) a parabolic shape of the cross-section can be an acceptable approximation of reality, the solution may be easily utilized to assess the maximum cross-sectional value of concentration and its location.

Conclusions
The assessment of the value and position of chemical concentration peaks in river flows is a fundamental requirement for the monitoring and the control of contamination events, which can be very frequent along river reaches that cross urban or industrial areas. We have shown that, for large aspect ratios, at the first order in the deviatory velocity and at least for large Peclet numbers, complete transverse mixing is never achieved. Therefore, the commonly used one-dimensional Gaussian is not a suitable predictor. The proposed analytical model allows evaluating the time-dependent deviatory depth-averaged concentration at any location along the channel axis and within the cross section. Assuming a large parabolic shape of the flow area, we found that, at the front of the section-averaged wave, the mass is concentrated at the centre. Conversely, at the back of the section-averaged wave, it is relegated near the banks. Additionally, while the large-time leading component of the solution tends to a constant, the first-order diffusion-related component linearly increases in time. This peculiar behaviour fully justifies what emerges from Lagrangian numerical simulations applied to three real-field cases represented by more or less asymmetric cross-sections and finite Peclet numbers. Indeed, all the reconstructed downstream breakthrough curves exhibited a considerable and highly persistent delay as compared to those obtained from the one-dimensional Gaussian, thus validating the theoretical conclusion even in more general and less restrictive cases.
Author Contributions: M.P. conceived and designed the work, provided the analytical closed-form solution and performed the numerical simulations; D.M. graphically elaborated the output and analysed the state-of-art (numerical and experimental studies); A.D.V. performed the analysis of field data and transformed them in computational format; B.M. analysed the state-of-art (theoretical studies) and critically revised the paper.

Conflicts of Interest:
The authors declare no conflict of interest.
x dimensionless x y dimensionless y x * dimensionless x in a reference frame moving with the section-averaged velocity X(t) longitudinal particle position at time t <X(t)> longitudinal average trajectory Y(t) transverse particle position at time t