2.1. Problem description

We study the soluble matter dispersing in a turbulent open channel flow, as shown in figure 1. By assuming that the mass and momentum transport in turbulent flows can be treated analogously (Reynolds’ analogy) and that isotropic turbulent diffusion can be well modelled using an eddy diffusivity $K^\ast$ that varies with depth (Taylor Reference Taylor1954), a simple yet pedagogically significant and theoretically useful model for turbulent open channel flows is obtained (Elder Reference Fischer1959; Chatwin Reference Chatwin1970; Barton Reference Barton1983a ). In dimensional form, the Reynolds-averaged convection–diffusion equation is given by

(2.1) \begin{eqnarray} \frac {\partial C^\ast }{\partial t^\ast } + u^\ast (z^\ast ) \frac {\partial C^\ast }{\partial \xi ^{\ast }} = K^\ast (z^\ast ) \frac {\partial ^2 C^\ast }{\partial \xi ^{\ast 2}} + K^\ast (z^\ast ) \frac {\partial ^2 C^\ast }{\partial z^{\ast 2}}, \end{eqnarray}

where $C^*$ is the concentration of solute, $t^*$ is the time, $u^*$ is the flow velocity, $\xi ^*$ is the longitudinal coordinate, and $z^*$ is the vertical coordinate. Dimensionless parameters are introduced as

(2.2) \begin{align} C=\frac {C^{*} h^{*2}} { Q^{*} }, \quad t=\frac {t^* \bar {K}^*}{h^{*2}}, \quad \textit{Pe} =\frac {u_\tau ^* h^{*}/\kappa }{\bar {K}^*} , \quad u=\frac {\kappa u^*}{u_\tau ^*}, \quad K=\frac {K^*}{\bar {K}^*}, \quad \xi =\frac {\xi ^*}{ h^{*}}, \quad z=\frac {z^*}{h^{*}}, \end{align}

where $h^*$ is the height of the channel, $Q^*$ is the total mass per unit length, $\bar {K}^*$ is the average diffusivity, $\kappa$ is the von Kármán constant, and $u_\tau ^*$ is the friction velocity. The average turbulent diffusivity $\bar {K}^*$ is defined through the dimensional eddy diffusivity $K^*=\kappa h^* u_\tau ^* (1-z^*/h^*) z^*/h^*$ as

(2.3) \begin{eqnarray} \bar {K}^* \triangleq \frac {1}{h^*} \int _0^{h^*} {K^*(z^*)}\, \mathrm{d} z^*= \frac {\kappa h^* u_\tau ^*}{6}, \end{eqnarray}

and the dimensional average velocity $\bar {u}^*$ is defined by

(2.4) \begin{eqnarray} \bar {u}^* \triangleq \frac {1}{h^*} \int _0^{h^*} {u^*(z^*)}\, \mathrm{d} z^*. \end{eqnarray}

Note that the eddy diffusivity $K$ is generally a function of both space and time, as it reflects the underlying eddy structures, which can evolve temporally and vary spatially. Different assumptions regarding the form of $K$ lead to different modelling approaches, each involving trade-offs between empirical observations and theoretical formulation (Pope Reference Refugio Cabañas-Mendoza, Olguín, Sánchez-Galván, Melo and Barrientos2001). The dimensionless velocity and turbulent diffusivity are adopted as (Elder Reference Fischer1959; Fischer Reference Fischer1966, Reference Fischer1976, Reference Follett, Hays and Nepf1979; Chatwin Reference Chatwin1970; Sayre Reference Sayre1975; Barton Reference Barton1983a ; Guo et al. Reference Hanjalić and Launder2021)

(2.5) \begin{eqnarray} u(z) = \frac {\kappa \bar {u}^*}{u_\tau ^*} + 1+ \ln z, \quad K(z)=6z(1-z). \end{eqnarray}

Despite its simplicity, the log-law turbulent velocity profile combined with a basic eddy diffusivity formulation remains a classical and widely used model for turbulent dispersion. This approach has a long history of use in both theoretical and applied studies (Taylor Reference Taylor1954; Elder Reference Fischer1959; Fischer Reference Fischer1966, Reference Fischer1976, Reference Follett, Hays and Nepf1979; Chatwin Reference Chatwin1970; Sayre Reference Sayre1975; Barton Reference Barton1983a ; Guo et al. Reference Hanjalić and Launder2021; Barik & Tsai Reference Barik and Tsai2025). The moment theory and diffusive formulation developed hereafter are general and can readily accommodate alternative or more accurate eddy diffusivity profiles (Ghisalberti & Nepf Reference Gill2005; Nepf & Ghisalberti Reference Nepf2008; Liang & Wu Reference Liu and Wang2014; Follett, Hays & Nepf Reference Ghisalberti and Nepf2019). In this sense, the present choice serves primarily as a canonical case that highlights the dominant physical mechanisms while allowing straightforward extensions. These idealisations not only reflect dominant physical mechanisms, but also serve to make the problem analytically tractable and facilitate comparisons with canonical models. Specifically, in a two-dimensional coordinate system (where $\xi$ and $z$ denote the dimensionless longitudinal and vertical coordinates, respectively), the governing equation of concentration $C$ becomes

(2.6) \begin{eqnarray} \frac {\partial C }{\partial t }+ \textit{Pe} u(z)\, \frac {\partial C }{\partial \xi }= K(z) \frac {{{\partial }^{2}}C }{\partial {{\xi }^{2}}}+\frac {{\partial } }{\partial {{z }}} \left ( K(z) \frac {{{\partial }}C }{\partial {{z }}}\right )\!. \end{eqnarray}

Figure 1.

Sketch for streamwise dispersion of soluble matter with eddy diffusion in a fully developed turbulent open channel flow following von Kármán’s logarithmic law from an initial line source discharge.

The initial condition reads

(2.7) \begin{eqnarray} \left . C\right |_{t=0}=\mathcal{C}^0(\xi ,z), \end{eqnarray}

where $\mathcal{C}^0$ is the initial concentration distribution corresponding to

(2.8) \begin{eqnarray} \int _{-\infty }^{\infty }\mathrm{d}\xi \int _0^1 \,\mathrm{d}z\,\mathcal{C}^0=1. \end{eqnarray}

Boundary conditions are employed as

(2.9) \begin{eqnarray} &C \rightarrow 0\quad \text{as } \quad |\xi | \rightarrow \infty , \end{eqnarray}

(2.10) \begin{eqnarray} & {{\left . \left (K(z)\frac {\partial C }{\partial z }\right ) \right |}_{z =0}}={{\left . \left (K(z)\frac {\partial C }{\partial z }\right ) \right |}_{z =1}}=0. \end{eqnarray}

2.2. Solutions of local and average moments

Since our focus is on the streamwise dispersion of a solute, an efficient way to characterise its evolution is through the local and average moments of the concentration. The moment method offers an analytical framework to track the centre of mass, the variance and higher-order structures of the concentration distribution over time. This approach generalises Taylor’s classical theory by incorporating spatial variations and depth-dependent diffusivity profiles. Building on the mathematical equivalence between moments and cumulants, we will later formulate a streamwise dispersion model in terms of local statistics. The analytical procedures for deriving these moment solutions are presented in this subsection.

Local moments of concentration can be defined as

(2.11) \begin{eqnarray} C_n (z,\,t) =\int _{-\infty }^{\infty } \xi ^n\, C(\xi ,z,t) \,\mathrm{d} \xi , \quad n=0,1,2,\ldots . \end{eqnarray}

Assuming that the concentration and its derivatives decay sufficiently rapidly as $\left |\xi \right |\rightarrow \infty$ (e.g. exponentially), the boundary terms arising from integration by parts vanish, i.e.

(2.12) \begin{eqnarray} \xi ^n \frac {\partial ^n C}{\partial \xi ^n} \rightarrow 0\quad \text{as } \quad \left |\xi \right |\rightarrow \infty ,\quad \text{for all } n. \end{eqnarray}

Then $C_n$ possesses analytical solutions of the following problem (Aris Reference Aris1956; Barton Reference Barton1983a , Reference Bartonb ):

(2.13) \begin{equation} \left .\begin{array}{c} \dfrac {\partial C_n }{\partial t } - \dfrac {{\partial } }{\partial {{z }}} \left (K(z) \dfrac {{\partial }C_n }{\partial {{z }}} \right ) = n(n-1)\, K(z)\, C_{n-2}+n\, \textit{Pe} u(z)\, C_{n-1},\\[9pt] \left .C_n\right |_{t=0}=\mathcal{C}_n^0,\\[9pt] {{\left . \left (K(z)\dfrac {\partial C_n }{\partial z }\right ) \right |}_{z =0}} ={{\left . \left (K(z)\dfrac {\partial C_n }{\partial z }\right ) \right |}_{z =1}}=0, \!\! \end{array}\right \} \end{equation}

where $C_{-2}=C_{-1}=0$ , and $\mathcal{C}_n^0 (z)= \int _{-\infty }^{\infty } \xi ^n \mathcal{C}^0\, \mathrm{d} \xi$ is the initial vertical distribution of the moment of concentration. The eigenvalue problem is

(2.14) \begin{eqnarray} \frac {\mathrm{d}}{\mathrm{d}z} \left (K \frac { \mathrm{d} f}{ \mathrm{d}z} \right )=-\eta f, \end{eqnarray}

with $\eta$ and $f$ denoting the eigenvalue and eigenfunction, respectively. The boundary conditions are

(2.15) \begin{eqnarray} {{\left . \left (K\frac {\partial f }{\partial z }\right ) \right |}_{z =0}} ={{\left . \left (K\frac {\partial f }{\partial z }\right ) \right |}_{z =1}}=0. \end{eqnarray}

Note that the above conditions are automatically satisfied because the eddy diffusivity is zero at the boundaries; cf. (2.5). Hence this is a singular Sturm–Liouville problem, and the problem can be transformed equivalently into Legendre’s equation by a simple transformation $\zeta =2z-1$ , with corresponding boundary conditions given as finite conditions (Barton Reference Barton1983b ). The eigenfunctions can be obtained with the Legendre polynomial $P_n(\zeta )$ as

(2.16) \begin{eqnarray} f_n (z)= P_n(2z-1), \quad n=0,1,2,\ldots , \end{eqnarray}

and corresponding eigenvalues are

(2.17) \begin{eqnarray} \eta _n = 6n (n+1), \quad n=0,1,2,\ldots . \end{eqnarray}

With the discrete set of eigenvalues and eigenfunctions that are normalised and orthogonal, the solution for $C_n$ from instantaneous initial releases reads

(2.18) \begin{eqnarray} C_n = \sum \limits _{i=0}^{\infty } p_{ni}(t) \exp {(-\eta _i t)}\,f_i(z),\quad n=0,1,2,\ldots , \end{eqnarray}

where the coefficients $p_{ni}$ are determined by inner products of initial distributions and eigenfunctions. The definitions of average moments are given by

(2.19) \begin{eqnarray} M_n (t) =\int _0^1 C_n\,\mathrm{d}z=\int _{-\infty }^{\infty } \int _0^1 \xi ^n\, C(\xi ,z,t)\,\mathrm{d}z\,\mathrm{d} \xi . \end{eqnarray}

Specifically, $M_n$ is the solution of the problem

(2.20) \begin{eqnarray} &\dfrac {\mathrm{d}M_n}{\mathrm{d}t}=n(n-1) \int _0^1 \textit{KC}_{n-2} \,\mathrm{d}z +n\textit{Pe} \int _0^1 uC_{n-1} \,\mathrm{d}z, \end{eqnarray}

noting that $C_{-1}=C_{-2}=0$ , and initial conditions read

(2.21) \begin{eqnarray} &M_n(0)= \int _0^1 \mathcal{C}_n^0 \,\mathrm{d}z. \end{eqnarray}

The analytical expressions for both local and average moments are detailed in Appendix A.

3.1. Justification for diffusive description

Taylor’s dispersion model provides an accurate description of the mean concentration in the asymptotic long-time limit. However, the characterisation of the local concentration distribution has traditionally relied on diffusive approximations (cf. Gill Reference Guan and Chen1967; Chatwin Reference Chatwin1970; Smith Reference Smith1982b ), which, despite their practical utility, lack rigorous theoretical justification in the transient regime. In this work, we emphasise that the streamwise dispersion model, formulated directly from the convection–diffusion equation without invoking early averaging or diffusive assumptions, provides a comprehensive and self-consistent solution that remains valid across all time scales. We begin by clarifying the fundamental distinction between local and average moments, which underpins the statistical complexity of transient dispersion processes. Subsequently, we analyse the asymptotic behaviour of phenomenological dispersion coefficients derived from both local and average perspectives, thereby justifying the emergence of diffusive characteristics in Taylor dispersion from a more general streamwise perspective.

Specifically, the first-order average moment $M_1$ is given by

(3.1) \begin{eqnarray} M_1= \textit{Pe} t + \textit{Pe} \sum \limits _{i=1}^{\infty } \overline {\mathcal{C}_0^0 f_i}\,\overline {u f_i}\frac {1-\exp {(-\eta _i t)}}{\eta _i}. \end{eqnarray}

Here, the overline denotes the cross-sectional average. The leading term $\textit{Pe} t$ reflects Taylor’s classical result that the centre of mass ultimately propagates at the discharge velocity, whereas the remaining terms represent a small, transient displacement relative to the initial centroid.

The mean centroid position is defined as $\bar {\nu }_1= {M_1}/{M_0}$ , whereas the local centroid position, dependent on transverse distributions, is given as $\gamma _1={C_1}/{C_0}.$ Accordingly, the mean drift velocity $\bar {U}_d$ is defined as the temporal rate of change of the mean centroid position,

(3.2) \begin{eqnarray} \bar {U}_d = \frac {\mathrm{d}\bar {\nu }_1}{\mathrm{d} t} = \frac {1}{M_0} \frac {\mathrm{d}M_1}{\mathrm{d} t} - \frac {\bar {\nu }_1}{M_0}\frac {\mathrm{d}M_0}{\mathrm{d} t} =\frac {\mathrm{d}M_1}{\mathrm{d} t} =\textit{Pe}+ \textit{Pe} \sum \limits _{i=1}^{\infty } \overline {\mathcal{C}_0^0 f_i}\,\overline {u f_i}\exp {(-\eta _i t)}, \end{eqnarray}

while the local drift velocity is expressed as

(3.3) \begin{eqnarray} {U}_d^{\star } = \frac {\mathrm{d}\gamma _1}{\mathrm{d} t} = \frac {1}{C_0} \frac {\mathrm{d}C_1}{\mathrm{d} t} - \frac {\gamma _1}{C_0}\frac {\mathrm{d}C_0}{\mathrm{d} t} . \end{eqnarray}

For asymptotically long times, we have

(3.4) \begin{eqnarray} \lim _{t \to \infty }C_0= M_0 =1. \end{eqnarray}

Hence the transient difference between local and average centroid position is

(3.5) \begin{eqnarray} \lim _{t \to \infty } \left (\gamma _1-\bar {\nu }_1\right )= \textit{Pe}\sum \limits _{j=1}^{\infty } \frac {{{f}_{\!j}}}{{\eta }_{\!j}} \,\overline {u{{f}_{\!j}}}, \end{eqnarray}

which is independent of time and of the initial conditions. This implies that in the asymptotic limit, the mean drift velocity equals the asymptotic local drift velocity, both approaching the discharge velocity:

(3.6) \begin{eqnarray} \lim _{t \to \infty }\bar {U}_d = \lim _{t \to \infty } U_d^{\star } = \textit{Pe}. \end{eqnarray}

Note that differences among these three quantities exist for earlier periods.

Explicit expressions for the eigenfunctions and eigenvalues used in these calculations are listed in table 1. For the line source (i.e. $\mathcal{C}^0=\delta (\xi )$ and $\mathcal{C}_0^0=\mathcal{C}_1^0=0$ ), with the understanding that $M_0=C_0=1$ and $M_1=\textit{Pe}t$ , the first-order local moment $C_1$ is given by

Table 1.

Explicit expressions of normalised eigenfunctions $f_i(z)$ and eigenvalues $\eta _i$ ( $i=1,2,\ldots ,6$ ) for streamwise dispersion in a turbulent open channel flow $0 \leqslant z \leqslant 1$ .

(3.7) \begin{align} C_1 &=\textit{Pe}t+ \textit{Pe}\sum \limits _{j=1}^{\infty }{{{f}_{\!j}}}\, \overline {u{{f}_{\!j}}}\dfrac { 1-\exp ( -{{\eta }_{\!j}}t )}{{\eta }_{\!j}}\nonumber \\ &=\textit{Pe}t-\frac {1}{8} \left [1-\exp {(-12t)}\right ](2z-1) -\frac {5}{216}\left [1-\exp {(-36t)}\right ](1-6z+6z^2) +\cdots . \end{align}

Thus the difference between local and mean drift velocities is

(3.8) \begin{eqnarray} U_d^{\star }-\bar {U}_d= \textit{Pe}\sum \limits _{j=1}^{\infty }{{{f}_{\!j}}}\, \overline {u{{f}_{\!j}}}\exp {(-\eta _j t)}, \end{eqnarray}

which decays exponentially with time. Consequently, the local drift velocity is always greater than the average drift, while both tend towards the discharge velocity in the long-time limit. Consistent with Guan & Chen (Reference Guan, Jiang, Tao, Chen and Lee2024), this local velocity is referred to as the ‘streamwise drift velocity’ because this parameter precisely depicts the drift along each streamline, in contrast to the bulk average quantity defined purposely for asymptotic descriptions. We emphasise that only the cross-sectional average of the streamwise drift accurately represents the actual transport behaviour of the solute cloud, whereas its average value coincides with the mean drift only in the asymptotic limit.

Once again, with the initial conditions $\mathcal{C}_0^0=\mathcal{C}_1^0=\mathcal{C}_2^0=0$ for the line source, the mean variance $\bar {\sigma }_x^2=M_2/M_0-M_1^2/M_0^2$ reduces to

(3.9) \begin{eqnarray} \bar {\sigma }_x^2 = 2t+ 2t\textit{Pe}^2+\textit{Pe}^2 t^2 \sum \limits _{i=1}^{\infty }\frac {\overline {u f_i f_i}}{\eta _i} -2\, \textit{Pe}^2 \sum \limits _{i=1}^{\infty }\overline {u f_i f_i} \,\frac {1-\exp {(-\eta _i t)}}{\eta _i^2}. \end{eqnarray}

Note that this variance differs from the local variance, prompting an investigation into the precise nature of this difference. The difference $\varDelta$ between the average local variance $\int _0^1 {\sigma _x^{\star }}^2\,\mathrm{d}z = \int _0^1 (C_2/C_0-C_1^2/C_0^2 )\,\mathrm{d}z \equiv \overline {{\sigma _x^{\star }}^2}$ and the mean variance $\bar {\sigma }_x^2$ is given by

(3.10) \begin{eqnarray} \varDelta &=& \overline {{\sigma _x^{\star }}^2}- \bar {\sigma }_x^2 =\overline {\left (\frac {C_2}{C_0}- \frac {C_1^2}{C_0^2}\right )} - \big (M_2-M_1^2\big ) \nonumber \\ &=& \overline { C_1^2}- M_1^2 = \textit{Pe}^2 \sum \limits _{i=1}^{\infty } (\overline {u{{f}_{i}}} )^2 \left [\frac {1-\exp \left ( -{{\eta }_{i}}t \right )}{{\eta }_{i}}\right ]^2 \!, \end{eqnarray}

resulting from the order in which vertical averaging and streamwise integral are performed. Equation (3.10) is derived based on the normalisation of eigenfunctions and the known analytical forms of the moments under the line source condition. Precisely, $\overline {{\sigma _x^{\star }}^2}$ quantifies the average variance over the vertical direction, thereby indicating the extent of dispersion along each streamline, while the mean variance $\bar {\sigma }_x^2$ corresponds to the variance that one would obtain under the assumption of complete vertical mixing. As time tends to infinity, the discrepancy converges to

(3.11) \begin{eqnarray} \lim _{t \to \infty } \varDelta = \lim _{t \to \infty } \big (\overline { C_1^2}- M_1^2 \big ) = \textit{Pe}^2 \sum \limits _{i=1}^{\infty }\, \frac {(\overline {u{{f}_{i}}})^2}{\eta _i^2} - \textit{Pe}^2, \end{eqnarray}

which is constant. Note that $M_0=1$ for this instantaneous pulse release at long times. Further, $\varDelta \to 0$ as $t \to \infty$ in this moving coordinate at the speed of mean flow velocity. In other words, $\bar {\sigma }_x^2$ is only accurate for asymptotically long times, when the drift velocity has reached its uniform value and the solute is vertically well mixed. Yet for short-time dynamics, accurately capturing local drift and dispersivity requires resolving the detailed local structure of the concentration field.

Following the pioneering work of Taylor (Reference Taylor1954), we recall that the longitudinal coordinate is replaced by $x=\xi -\bar {u} t$ , and streamwise velocity by $u^\prime =u-\bar {u}$ , where $\bar {u}$ denotes the dimensionless mean flow velocity. The Taylor–Aris dispersion model is expressed as

(3.12) \begin{eqnarray} \frac {\partial \bar {C}}{\partial t} = \bar {D}_T \frac {\partial ^2 \bar {C}}{\partial x^2}, \end{eqnarray}

possessing a well-established Gaussian solution:

(3.13) \begin{eqnarray} \bar {C}= \frac {1}{\sqrt {4\pi \bar {D}_Tt}} \exp {\left (-\frac {x^2}{4\bar {D}_Tt}\right )}. \end{eqnarray}

In the long-time limit, $\bar {D}_T(\infty )$ becomes Taylor’s dispersivity as constant. By multiplying both sides of (3.12) by $x^2$ , integrating over $x\in (-\infty ,\infty )$ and applying integration by parts, we obtain the mean dispersivity as

(3.14) \begin{eqnarray} \bar {D}_T(t)= \frac {1}{2} \frac {\mathrm{d}\bar {\sigma }^2_x}{\mathrm{d}t} =\frac {1}{2}\frac {\mathrm{d}M_2}{\mathrm{d}t}-M_1\frac {\mathrm{d}M_1}{\mathrm{d}t}, \end{eqnarray}

which precisely recovers Taylor’s result in the long-time limit. In parallel with the local variance, a local dispersivity is defined by

(3.15) \begin{eqnarray} D_T^{\star }(z,t)= \frac {1}{2} \frac {\partial {{\sigma _x^{\star }}^2}}{\partial t} = \frac {1}{2} \frac {\partial }{\partial t}\left [\frac {C_2}{C_0}-\frac {C_1^2}{C_0^2}\right ]\!, \end{eqnarray}

and the average over the cross-section gives rise to the average local dispersivity:

(3.16) \begin{eqnarray} \overline {D_T^{\star }}(t)= \int _0^1 D_T^{\star }\,\mathrm{d}z=\int _0^1 \frac {1}{2} \frac {\partial {{\sigma _x^{\star }}^2}}{\partial t}\,\mathrm{d}z. \end{eqnarray}

Equation (3.15) is consistent with the definition of streamwise dispersivity used by Vilquin et al. (Reference Vilquin, Bertin, Soulard, Guyard, Raphaël, Restagno, Salez and McGraw2021).

By evaluating the local variance and dispersivity, we show that solute particles do not uniformly sample vertical positions in open channel flows. Moreover, we demonstrate that the average Taylor–Aris dispersivity is a conceptually asymptotic quantity, valid only as a long-time limit of the full solute cloud dynamics. We stress that solely the local dispersivity precisely captures the axial displacement along each streamline and distinctly showcases the notable enhancement of axial spreading. Nevertheless, both dispersivity definitions converge to the same constant in the long-time limit.

Having established the distinction between local and bulk average statistics, it is clear that local moments, rather than their averaged counterparts, are the appropriate quantities for constructing diffusive series expansions of the transient concentration distribution. We now provide further justification for the diffusion-type model, which, extensively employed in practice as an empirical model, lacks solid theoretical foundation in the pre-asymptotic regime.

Following the approach of Guan & Chen (Reference Guan, Jiang, Tao, Chen and Lee2024) in a cylindrical coordinate system, we introduce a convected coordinate system that moves with the local fluid velocity along each streamline, which is devised for the present configuration,

(3.17) \begin{equation} {x}^{\prime }=\xi -\textit{Pe} u(z)\,t, \end{equation}

as a convection transformation to decouple convection from diffusion, to give a purely diffusion equation

(3.18) \begin{equation} \frac {\partial C}{\partial t}=\boldsymbol{\nabla }\boldsymbol{\cdot } (\unicode{x1D63F}\, \boldsymbol{\nabla }C). \end{equation}

Using the chain rule, the equivalent diffusivity tensor is expressed as

(3.19) \begin{equation} {\unicode{x1D63F}}=K(z)\left ( \begin{matrix} 1+{{ \left ( {\textit{Pe}{t}} \right )}^{2}}{u_z^{\prime 2}} & -\textit{Pe}{{t}}\,{ u_z^\prime } \\[3pt] -\textit{Pe}{{t}}\,{ u_z^\prime } & 1 \end{matrix} \right )\!, \end{equation}

where $K(z)$ is exactly the diffusivity introduced in (2.5).

This transformation naturally segregates the diffusivity tensor $\unicode{x1D63F}$ into three components, each scaling with a different power of the product of the Péclet number and the dimensionless time:

(3.20) \begin{equation} \unicode{x1D63F}={\unicode{x1D63F}^{(\unicode{x1D64A})}}+{\textit{Pe}{t}}\,{\unicode{x1D63F}^{(\unicode{x1D63E})}}+{{\left ( {\textit{Pe}{t}} \right )}^{2}}\,{\unicode{x1D63F}^{(\unicode{x1D640})}}, \end{equation}

where

(3.21) \begin{equation} {\unicode{x1D63F}^{(\unicode{x1D64A})}}=\left ( \begin{array}{*{35}{c}} K & 0 \\[3pt] 0 & K \end{array} \right )\!,\quad {\unicode{x1D63F}^{(\unicode{x1D63E})}}=\left ( \begin{matrix} 0 & K u_z^\prime \\[3pt] K u_z^\prime & 0 \end{matrix} \right )\!,\quad {\unicode{x1D63F}^{(\unicode{x1D640})}}=\left ( \begin{matrix} K u_z^{\prime 2}& 0 \\[3pt] 0 & 0 \end{matrix} \right )\!. \end{equation}

The first tensor $\unicode{x1D63F}^{(\unicode{x1D64A})}$ represents the isotropic turbulent diffusion of soluble matter; the second tensor ${\unicode{x1D63F}^{(\unicode{x1D63E})}}$ denotes a pair of symmetric cross-diffusivities due to the velocity shear; and the third tensor $\unicode{x1D63F}^{(\unicode{x1D640})}$ quantifies a streamwise diffusion enhancement. With help of (3.21), the purely diffusion equation can be decomposed as

(3.22) \begin{equation} \frac {\partial {{C }}}{\partial {{t}}}={{\boldsymbol{\nabla }}} \boldsymbol{\cdot } \big ( {\unicode{x1D63F}^{(\unicode{x1D64A})}}\,{{\boldsymbol{\nabla }}}{{C}} \big )+\textit{Pe}{{t}}\,{{\boldsymbol{\nabla }}}\boldsymbol{\cdot } \big ( {\unicode{x1D63F}^{(\unicode{x1D63E})}}\,{{ \boldsymbol{\nabla }}}{{C }} \big )+{{\big ( \textit{Pe}{{t}} \big )}^{2}}{{\boldsymbol{\nabla }}}\boldsymbol{\cdot } \big ( {\unicode{x1D63F}^{(\unicode{x1D640})}}\,{{\boldsymbol{\nabla }}}{{C }} \big ). \end{equation}

For a given Péclet number, for a short time we obtain a conventional isotropic diffusion equation

(3.23) \begin{equation} \frac {\partial {{C }}}{\partial {{t}}}= K(z)\frac {\partial ^2 {{C }}}{\partial {{x}^{\prime 2 }}} +\frac {\partial {} }{\partial {{z}}} \left ( K(z) \frac {\partial {C} }{\partial {{z}}} \right )\!, \end{equation}

and at asymptotically long times, we retain a one-dimensional streamwise diffusion equation with a substitution of the log-law velocity profile as

(3.24) \begin{equation} \frac {\partial {{C }}}{\partial {{t}}}= K(z)\!\left (1+\frac {\textit{Pe}^2 t^2}{z^2}\right )\,\frac {\partial ^2 {{C }}}{\partial {{x}^{\prime 2 }}}, \end{equation}

parallel to the one-dimensional dispersion equation for mean concentration in Taylor’s theory. This coordinate transformation confirms that the dispersion process is fundamentally diffusive in nature when viewed through the lens of this space- and time-dependent effective diffusivity characterised by (3.19). Crucially, it shows that the enhanced spreading characteristic of Taylor dispersion is not a separate mechanism but emerges naturally from the $\unicode{x1D63F}^{(\unicode{x1D640})}$ tensor, whose influence grows quadratically with $\textit{Pe} t$ . This perspective not only justifies the use of diffusive models across all time scales, but also provides a pathway for constructing unified models that bridge the gap between the short-time and asymptotic regimes.

3.2. Derivation of concentration distribution in the moment space

Grounded in first principles of convection–diffusion dynamics, a streamwise dispersion model is formulated for turbulent open channel flows using local moments rather than bulk average statistics. This approach enables a rigorous derivation of the local concentration field via series expansions, capturing non-Gaussian features of the solute distribution. In probability theory, cumulants provide an alternative characterisation to moments; the two sets are equivalent in that they uniquely determine each other. Using this equivalence, the local concentration is analytically expressed in terms of local cumulants, with solutions obtained at various truncation levels.

We define the local probability density function $P$ of the longitudinal concentration distribution as

(3.25) \begin{eqnarray} P(\xi ,z,t)=\dfrac {C(\xi ,z,t)}{\int _{-\infty }^{\infty }C\,\mathrm{d}\xi }=\frac {C(\xi ,z,t)}{C_0(z,t)}, \end{eqnarray}

which gives the probability of finding a tracer particle at point ( $\xi ,z$ ) at time $t$ . The Fourier transform of $P$ yields the characteristic function with regard to the angular frequency $k$ (Kendall & Stuart Reference Kumar, Thomson, Powers and Harris1958; Chatwin Reference Chatwin1970; Jiang & Chen Reference Jiang and Chen2018; Guan & Chen Reference Guan, Jiang, Tao, Chen and Lee2024):

(3.26) \begin{eqnarray} \varphi (k,z,t) \equiv \int _{-\infty }^{\infty } P \exp (\mathrm{i} k \xi )\,\mathrm{d}\xi =\sum \limits _{n=0}^{\infty }\frac {\left (\mathrm{i}k\right )^n}{n!} \gamma _n, \end{eqnarray}

where $\gamma _n$ is the $n$ th-order local raw moment of $P$ , defined as

(3.27) \begin{eqnarray} \gamma _n(z,t) = \int _{-\infty }^{\infty } \xi ^n\, P(\xi ,z,t)\, \mathrm{d}\xi =\frac {C_n}{C_0}. \end{eqnarray}

Correspondingly, the local central moment of order $n$ is

(3.28) \begin{eqnarray} \mu _n(z,t) = \int _{-\infty }^{\infty } (\xi -\gamma _1)^n\, P(\xi ,z,t)\, \mathrm{d}\xi . \end{eqnarray}

The cumulant-generating function is the natural logarithm of the characteristic function, and the cumulants $\kappa _n$ are defined by the expansion

(3.29) \begin{eqnarray} \ln \varphi (k,z,t)=\sum \limits _{n=1}^{\infty }{\kappa _n}\frac {\left (\mathrm{i}k\right )^n}{n!}. \end{eqnarray}

Substituting (3.29) into the inverse Fourier transform of $\varphi$ gives

(3.30) \begin{equation} \frac {C(\xi ,z,t)}{\int _{-\infty }^{\infty } C(\xi ,z,t)\,\mathrm{d}\xi } =\frac {1}{2\pi }\int _{-\infty }^{\infty } \exp \left [{\sum \limits _{n=1}^{\infty }{\kappa _n}\frac {\left (\mathrm{i}k\right )^n}{n!}}\right ] \exp (-\mathrm{i} k \xi )\,\mathrm{d} k. \end{equation}

Using the identity

(3.31) \begin{equation} \int _{-\infty }^{\infty } \exp \left (\!-\mathrm{i} k\xi -\frac {1}{2} k^2\right )(\mathrm{i} k)^n \,\mathrm{d} k=\sqrt {2 \pi }\,\mathit{He}_{n}(\xi )\,\exp \left (\!-\frac {1}{2} \xi ^2\right )\!, \end{equation}

together with a Taylor expansion of the exponential involving higher-order cumulants in (3.30), we obtain the series representation

(3.32) \begin{equation} C(\xi ,z,t)=\frac {{C}_{0}}{\sqrt {2 \pi {{\kappa }_{2}}}}\exp\! \left [ -\frac {{{ (\xi - \kappa _1) }^{2}}}{{2! \, {\kappa }_{2}}}\right ] \left [1+ \frac {\kappa _3}{3!\,\left (\kappa _2\right )^{3/2}} \mathit{He}_{3}\! \left ( \frac {{{ \xi - \kappa _1}}}{{{\kappa }_{2}^{1/2}}} \right ) + \frac {\kappa _4}{4!\,\kappa _2^2} \mathit{He}_{4}\! \left ( \frac {{{ \xi - \kappa _1 }}}{{{\kappa }_{2}^{1/2}}} \right ) +\cdots \right ]\!, \end{equation}

where $\mathit{He}_{n}$ denotes the Chebyshev–Hermite polynomial of order $n$ . The cumulants are related to raw moments and central moments by

(3.33) \begin{eqnarray} &\kappa _1=\gamma _1, \end{eqnarray}

(3.34) \begin{eqnarray} &\kappa _2=\gamma _2-\gamma _1^2=\mu _2, \end{eqnarray}

(3.35) \begin{eqnarray} &\kappa _3=\gamma _3-3\gamma _2\gamma _1+2\gamma _1^3=\mu _3, \end{eqnarray}

(3.36) \begin{eqnarray} &\kappa _4=\gamma _4-4\gamma _3\gamma _1-3\gamma _2^2+12\gamma _2\gamma _1^2-6\gamma _1^4=\mu _4-3\mu _2^{ 2}. \end{eqnarray}

Note that cumulants coincide with raw moments only for the first order, and with central moments only up to third order; for higher orders, they differ. For ease of interpretation, we introduce phenomenological dispersion coefficients $\lambda _n$ :

(3.37) \begin{eqnarray} \lambda _1=\kappa _1,\quad \lambda _2=\kappa _2,\quad \lambda _3=\dfrac {\kappa _3}{ \left (\kappa _2\right )^{3/2}},\quad \lambda _4=\dfrac {\kappa _4}{\kappa _2^2}. \end{eqnarray}

Physically, $\lambda _1$ is the local mean, $\lambda _2$ is the local variance, $\lambda _3$ is the local skewness (quantifying asymmetry), and $\lambda _4$ is the local kurtosis (indicating peak sharpness). A normalised streamwise coordinate is defined as

(3.38) \begin{eqnarray} X \equiv \frac {\xi -\lambda _1}{\sqrt {\lambda _2} }, \end{eqnarray}

so that the probability distribution function in terms of $X$ has zero mean and unit variance. This coordinate naturally follows the local velocity at short times and the mean flow speed at long times, as discussed in the context of the two time limits (see e.g. (3.23) and (3.24)). In terms of $X$ , the series solution (3.32) becomes

(3.39) \begin{eqnarray} C(\xi ,z,t)=\dfrac {{C}_{0} \exp \left ( -{X^2}/{2} \right )}{\sqrt {2\pi {{\kappa }_{2}}}} \left [1+ \frac {1}{6} \lambda _3\, \mathit{He}_{3}(X) + \frac {1}{24} \lambda _4\, \mathit{He}_{4}(X) +\cdots \right ]\!. \end{eqnarray}

The truncation order of the streamwise model is set by the available phenomenological dispersion coefficients $\lambda _n$ . A second-order truncation (retaining only $\lambda _1$ and $\lambda _2$ ) exactly represents the streamwise drift and dispersivity, while a fourth-order model (including $\lambda _3$ and $\lambda _4$ ) captures skewness and kurtosis with full accuracy up to that order.

This approach is reminiscent of the method used by Aris (Reference Aris1958) to derive the Gram–Charlier series of type A, later refined by Chatwin (Reference Chatwin1970) through Edgeworth’s modification. Our results are distinguished by the explicit incorporation of local moments and cumulants. The phenomenological dispersion coefficients are time-dependent and locally defined. Consistent with Taylor (Reference Taylor1954), we employ a moving coordinate for visualisation, namely $x=\xi -\bar {u} t$ , corresponding to a reference frame travelling at the mean flow velocity $\bar {u}$ .