2.1. Governing equation

To illustrate the analytical solution procedure, we consider mass transport in a parallel channel with a time-dependent velocity field. For simplicity, we focus on the two-dimensional transport problem. A typical case is shown in figure 1.

Figure 1.

Sketch of solute transport problem in an oscillatory flow created by the oscillating channel wall. Several simulated particle trajectories are shown. The time variable in the oscillatory velocity field can be virtually viewed as an additional ‘transverse’ spatial variable.

The dimensionless form of the whole transport problem can be adopted as follows:

(2.1a) \begin{align} \frac {\partial C}{\partial t} + u (y, t) \frac {\partial C}{\partial x} & = \frac {1}{\textit{Pe}^2} \frac {\partial ^2 C}{\partial x^2} + \frac {\partial ^2 C}{\partial y^2}, \end{align}

(2.1b) \begin{align} \left . \frac {\partial C}{\partial y} \right |_{y = 0} &= \left . \frac {\partial C}{\partial y} \right |_{y = 1} = 0, \end{align}

(2.1c) \begin{align} C | _{t = 0} &= C_{{{ini}}} (y) \delta (x) . \end{align}

The range of the dimensionless variables is $- \infty \lt x \lt + \infty$ (the channel length is infinite), $0 \lt y \lt 1$ (the dimensionless channel width is $1$ ) and $t \gt 0$ (the time). Note that the following dimensionless parameters and variables have been chosen (the superscript asterisk denotes dimensional variables and parameters):

(2.2) \begin{align} x &= \frac {1}{\textit{Pe}} \frac {x^{\ast }}{L^{\ast }}, \quad y =\frac {y^{\ast }}{L^{\ast }}, \quad t = \frac {t^{\ast }}{{L^{\ast }}^2 /D^{\ast }}, \nonumber \\ u &= \frac {u^{\ast }}{U^{\ast }_{\textrm{c}}}, \quad C = \textit{Pe} \frac {C^{\ast }}{C^{\ast }_{{\text{c}}}}, \quad \textit{Pe} = \frac {U^{\ast }_{\textrm{c}} L^{\ast }}{D^{\ast }}, \quad \omega = \frac {\omega ^{\ast }}{D^{\ast } / {L^{\ast }}^2}, \end{align}

where $L^{\ast }$ is the channel width, $D^{\ast }$ is the solute diffusion coefficient, $U^{\ast }_{{\text{c}}}$ is a characteristic velocity and $\textit{Pe}$ is the corresponding Péclet number. Here, $C^{\ast }_{{\text{c}}}$ is a characteristic concentration constant such that the dimensionless concentration distribution satisfies

(2.3) \begin{equation} \int ^{\infty }_{- \infty } \int ^{1}_0 C (x, y, t) \, \textrm{d} y \,\textrm{d} x = 1. \end{equation}

Namely, $C$ is a well-defined probability density function.

For the advection term in (2.1a ), the dimensionless time-dependent flow velocity $u (y, t)$ is assumed to be oscillatory and the dimensionless angular frequency is $\omega$ . Namely, the dimensionless period of the flow velocity $T = {2 \pi }/{\omega }$ and

(2.4) \begin{equation} u (y, t + T) =u \left(y, t + \frac {2 \pi }{\omega }\right) = u (y, t). \end{equation}

Equation (2.1b ) is the non-penetration boundary condition. Equation (2.1c ) is the initial condition. The solute is assumed to be instantaneously released in a cross-section of the channel. Here, $\delta$ denotes the Dirac delta function and $C_{{{ini}}} (y)$ is initial transverse distribution.

2.2. Concentration moments

The goal of this study is to derive the complete time-dependent solution of concentration moments for the above-mentioned transport problem with unsteady flow velocity. First, we follow the work by Aris (Reference Aris1956) to define the $n$ th-order concentration moment for $C (x, y, t)$ as

(2.5) \begin{equation} C_n (y, t) \triangleq \int ^{\infty }_{- \infty } x^n C (x, y, t) \, \textrm{d} x \end{equation}

for $n = 0, 1, 2, \ldots$ . The moments of the cross-sectional mean concentration distribution are more practical in applications. Here, we use the bar to denote variables related to the cross-sectional mean concentration distribution and

(2.6) \begin{equation} \bar {C} (x, t) \triangleq \int ^{1}_{0} C (x, y, t) \; \textrm{d} y. \end{equation}

Similarly, the $n$ th order moment of $\bar {C}$ is defined as

(2.7) \begin{equation} \bar {C}_n (t) \triangleq \int ^{\infty }_{- \infty } x^n \bar {C} (x, t) \, \textrm{d} x \end{equation}

for $n = 0, 1, 2, \ldots$ .

The first four order moments are related to the most fundamental overall characteristics of the solute transport problem. Here, $\bar {C}_0$ represents the total amount of solute mass and $\bar {C}_0=1$ for the current transport problem. The averaged moving speed of the solute cloud, i.e. the drift, is

(2.8) \begin{equation} u_d \triangleq \frac {\textrm{d} \overline {\mu _x}}{\textrm{d} t} =\frac {\textrm{d}}{\textrm{d} t} \left ( \frac {\bar {C}_1}{\bar {C}_0} \right ) = \frac {\textrm{d} \bar {C}_1}{\textrm{d} t}, \end{equation}

where $\overline {\mu _x}$ is the mean position of the solute cloud (expected value).

The effective diffusion (apparent dispersivity) is related to the second-order moment $\bar {C}_2$ as

(2.9) \begin{equation} D_e \triangleq \frac {1}{2} \frac {\textrm{d} \bar {\sigma }^2}{\textrm{d} t} = \frac {1}{2} \frac {\textrm{d}}{\textrm{d} t} \left ( \frac {\bar {C}_2}{\bar {C}_0} - \frac {\bar {C}_1^2}{\bar {C}_0^2} \right ) = \frac {1}{2} \frac {\textrm{d}}{\textrm{d} t} (\bar {C}_2 - \bar {C}_1^2), \end{equation}

where $\bar {\sigma }^2$ is the variance of the mean concentration distribution.

The higher-order statistics including skewness and kurtosis are also of interest because they can basically reflect the deviation from a Gaussian distribution. Skewness is a measure of the asymmetry of the longitudinal concentration distribution, whose definition is related to the third-order moment $\bar {C}_3$ as

(2.10) \begin{equation} \textit{Sk} \triangleq \frac {\bar {\kappa }_3}{\bar {\sigma }^{3}} = \frac {\bar {C}_3- 3 \bar {C}_1 \bar {C}_2 + 2 \bar {C}_1^3}{\bar {\sigma }^{3}}, \end{equation}

where $\bar {\kappa }_3$ is the third-order cumulant.

Kurtosis (excess kurtosis) is a measure of the outliers of a distribution (the tail weight) compared with a Gaussian distribution. One can use the fourth-order cumulant $\bar {\kappa }_4$ to define the excess kurtosis as

(2.11) \begin{equation} \textit{Ku} \triangleq \frac {\bar {\kappa }_4}{\bar {\sigma }^{4}} = \frac {\bar {C}_4 - 4 \bar {C}_1 \bar {C}_3 - 3 \bar {C}_2^2 + 12 \bar {C}_1^2 \bar {C}_2 - 6 \bar {C}_1^4}{\bar {\sigma }^{4}}. \end{equation}

Therefore, we need the fourth-order moment $\bar {C}_4$ .

Typically, the next step of the method of concentration moments is to derive the governing equations of moments and then solve them order by order (see the paper by Ding et al. Reference Ding, Hunt, McLaughlin and Woodie2021 for details). The governing equation of moments $C_n$ for the current transport problem is given in Appendix A. However, due to the time-dependent velocity field, the tedious treatment of the time integral with the oscillation terms is a key issue (cf. Singh & Murthy Reference Singh and Murthy2023, with (D5)–(D8) numerically integrated). The corresponding solution procedure is much more complicated than that of the steady-flow case, making it difficult to obtain the complete explicit expressions for higher-order moments. Therefore, we will not solve the original moment equations. Instead, we introduce an associated two-time-variable transport problem to get rid of the time-dependent velocity field.

3.1. Two-time-variable system

First, we impose a two-time-variable system with

(3.1) \begin{equation} t_0 = t, \quad t_1 = \omega t. \end{equation}

While $t_0$ represents the original basic time, the new auxiliary time $t_1$ , referred to as oscillation time, is introduced directly associated with the frequency $\omega$ of the oscillating flow to characterise the inherent oscillation in the dispersion. Thus, the oscillatory velocity is only relevant to $t_1$ . For convenience, we define a new velocity function of $t_1$ as

(3.2) \begin{equation} u_1 (y, t_1) \triangleq u \left(y, \frac {t_1}{ \omega }\right)=u (y, t), \end{equation}

such that we can simply replace $u (y, t)$ with $u_1 (y, t_1)$ in the original transport equation (2.1a ). Note that $u_1$ is a periodic function of $t_1$ with a period of $2 \pi$ .

Next, the concentration distribution $C(x,y,t)$ needs to be extended to some virtual ‘concentration’ function in the two-time-variable system, denoted as $P(x, y, t_0, t_1)$ . The relationship between $C$ and $P$ is

(3.3) \begin{equation} C (x, y, t) = P (x, y, t_0=t, t_1=\omega t). \end{equation}

Note that there may exist many possible extensions from $C$ to $P$ in the higher-dimensional system. One of the most straightforward choices is to set $P$ periodic with respect to $t_1$ because $t_1$ is the oscillation time and, thus, it requires

(3.4) \begin{equation} P(x,y,t_0, t_1)= P(x,y,t_0, t_1 + 2\pi ). \end{equation}

For the time derivative, using the chain rule, we have

(3.5) \begin{equation} \frac {\partial C}{\partial t} = \frac {\partial P}{\partial t_0} + \frac {\partial P}{\partial t_1} \frac {\partial t_1}{\partial t} = \frac {\partial P}{\partial t_0} + \omega \frac {\partial P}{\partial t_1} . \end{equation}

Readers may find that the above-mentioned two-time-variable procedure is similar to the two-time-scale homogenisation method used in ‘previous studies (e.g. Ding et al. Reference Ding, Hunt, McLaughlin and Woodie2021). In fact, if one replaces $\omega$ with a so-called small perturbation parameter $\epsilon$ , then the formulations are the same. However, we would not dive into the perturbation analysis with the multi-scale expansion of concentration. More importantly, the angular frequency of the velocity in this study is not required to be confined as in the two-scale method. Neither a small parameter nor scale analysis is needed. Thus, we use the term ‘two-time variable’ instead of ‘two-time scale’. What we are looking for is the complete solution of concentration moments, not the long-time asymptotic concentration distribution (which is the focus of the homogenisation method). Our work finds that the two-time variable formulation could provide a new perspective for the analytical solution of concentration moments. We have compared the analytical result by the two-time variable method with numerical results for different magnitude values of parameters, see figures 3 and 9.

3.2. Governing equation for the two-time-variable system

Equation (3.3) offers the relationship between $C$ and $P$ . Substituting $P$ into the transport equation (2.1a ) and using the time derivatives in (3.5), we can obtain the governing equation for $P$ as

(3.6) \begin{equation} \frac {\partial P}{\partial t_0} + \omega \frac {\partial P}{\partial t_1} + u_1 (y, t_1) \frac {\partial P}{\partial x} = \frac {1}{\textit{Pe}^2} \frac {\partial ^2 P}{\partial x^2} + \frac {\partial ^2 P}{\partial y^2}. \end{equation}

Apparently, the second term on the left side reflects the auxiliary flux change due to periotic oscillation in the oscillation time dimension $t_1$ , actually a translation effect in the oscillation space, with the oscillation frequency standing for the translation speed in the oscillation time dimension. When $\omega =0$ , i.e. for the steady case, the form of (3.6) for $P$ is the same as the original (2.1a ) for $C$ .

Because $P$ is periodic with respect to $t_1$ as in (3.4), we can confine the whole problem within one periodic interval $[0, 2 \pi ]$ of $t_1$ . Namely, $t_1$ is with the modulus $2 \pi$ ,

(3.7) \begin{equation} t_1 = (\omega t) \mod 2 \pi , \end{equation}

then we can impose the following periodic conditions with respect to $t_1$ for $P$ :

(3.8) \begin{equation} P |_{t_1=0} = P |_{t_1 = 2 \pi }. \end{equation}

It is further required that $P$ represents the ‘probability density function’ in the two-time-variable system (within one period of $t_1$ ), namely,

(3.9) \begin{equation} \int ^{\infty }_{- \infty } \int ^{1}_0 \int ^{2 \pi }_0 P (x, y, t_0, t_1) \, \textrm{d} t_1\, \textrm{d} y\, \textrm{d} x = 1. \end{equation}

We can define the probability moment of $P$ in § 3.3.

The boundary condition of $P$ at the channel walls is the same as that of $C$ in (2.1b ),

(3.10) \begin{equation} \left . \frac {\partial P}{\partial y} \right |_{y = 0} = \left . \frac {\partial P}{\partial y} \right |_{y = 1} = 0. \end{equation}

To complete the transport problem of $P$ in the two-time-variable system, an ‘initial’ condition at $t_0=0$ is needed. According to the initial condition (2.1c ) of $C$ , the simplest choice for $P$ is

(3.11) \begin{equation} P | _{t_0 = 0} = C_{{{ini}}} (y) \delta (x). \end{equation}

Namely, the initial condition is independent of $t_1$ . In fact, we just set the initial condition for $P$ to satisfy

(3.12) \begin{equation} C (x, y, t=0) = P (x, y, t_0=0, t_1=0), \end{equation}

and $P | _{t_0 = 0}$ is periodic with respect to $t_1$ .

We have obtained the whole initial-boundary value problem (3.6), (3.8), (3.10), (3.11) of $P$ in the two-time-variable system. If a solution of $P$ is given, then through (3.3), the solution of $C$ can be immediately obtained because it satisfies all the equations in the transport problem (2.1) for $C$ . However, obtaining the exact solution of $P$ is challenging, as is the case for $C$ . Therefore, we consider the concentration moments instead, as in § 2.2.

3.3. Probability moments

Similar to the definition of concentration moment in (2.5), one can also define the moments of $P(x, y, t_0, t_1)$ and then obtain the corresponding moment equations. The $n$ th-order probability moment is

(3.13) \begin{equation} P_n (y, t_0, t_1) \triangleq \int ^{\infty }_{- \infty } x^n P(x, y, t_0, t_1) \, \textrm{d} x \end{equation}

for $n = 0, 1, 2, \ldots$ . According to the relationship (3.3) between $P$ and $C$ , we have the relationship between $P_n$ and $C_n$ ,

(3.14) \begin{equation} C_n (y, t) = P_n(y, t_0=t, t_1=\omega t). \end{equation}

Namely, once we solve $P_n$ , the solution for $C_n$ can be immediately obtained.

According to (3.6), (3.8), (3.10), (3.11), the whole initial-boundary value problem for $P_n$ is

(3.15a) \begin{align} \frac {\partial P_n}{\partial t_0} + \omega \frac {\partial P_n}{\partial t_1} - \frac {\partial ^2 P_n}{\partial y^2} & = n \, u_1 (y, t_1) P_{n - 1} + \frac {n (n-1)}{\textit{Pe}^2}P_{n - 2}, \end{align}

(3.15b) \begin{align} P_n |_{t_1 = 0} &= P_n |_{t_1 = 2 \pi }, \end{align}

(3.15c) \begin{align} \left . \frac {\partial P_n}{\partial y} \right |_{y = 0} &= \left . \frac {\partial P_n}{\partial y} \right |_{y = 1} = 0, \end{align}

(3.15d) \begin{align} P_n | _{t_0 = 0} & = \left \{\begin{array}{ll} C_{{{ini}}} (y), & n = 0,\\ 0 , & n = 1, 2, \ldots , \end{array} \right . \end{align}

with auxiliary terms $P_{-1}=P_{-2}=0$ .

3.4. Eigenvalue problem

In the two-time-variable system, the governing equation (3.15a ) of moments for the oscillatory flow is similar to the case of a steady flow because $t_1$ can be viewed as a spatial variable and the velocity is independent of $t_0$ . Therefore, now we are able to directly apply the generic solution expressions of moments given by Barton (Reference Barton1983, Reference Barton1984) for the case of steady flows. The key is to solve the corresponding eigenvalue problem to determine the eigenvalues and eigenfunctions in Barton’s expressions.

According to (3.15), the complete eigenvalue problem for probability moment $P_n$ is

(3.16a) \begin{align} - \omega \frac {\partial f}{\partial t_1} + \frac {\partial ^2 f}{\partial y^2} & = - \lambda f, \end{align}

(3.16b) \begin{align} f |_{t_1=0} &= f |_{t_1 = 2 \pi }, \end{align}

(3.16c) \begin{align} \left . \frac {\partial f}{\partial y} \right |_{y = 0} &= \left . \frac {\partial f}{\partial y} \right |_{y = 1} = 0, \end{align}

where $\lambda$ is an eigenvalue and $f=f(y, t_1)$ is the corresponding eigenfunction.

As mentioned in the introduction, we remark that Mukherjee & Mazumder (Reference Mukherjee and Mazumder1988, (3.1)) have introduced a similar ‘eigenvalue’ problem (time-dependent) to extend Barton’s method to oscillatory flows. Furthermore, Shapiro & Brenner (Reference Shapiro and Brenner1990, (16) and (17)) proposed the same form of the time-dependent problem. They pointed out that the eigenvalue problem is non-self-adjoint and gave the full solution form of the zeroth-order moment.

The fundamental difference we make here is the introduction of the oscillation time variable $t_1$ to make the oscillatory flow ‘steady’ with respect to time variable $t_0$ . The eigenvalue problem (3.16) is thus independent of time $t_0$ , intrinsically different from the time-dependent one discussed by Mukherjee & Mazumder (Reference Mukherjee and Mazumder1988) and Shapiro & Brenner (Reference Shapiro and Brenner1990). Thus, the introduction of $t_1$ not only enables one to use Barton’s result directly, but also can clarify the solution structure of moments (cf. § 6).

The eigenvalue problem (3.16) can be solved using the complex form of Fourier series. The solution of eigenfunction can be written as (see also Appendix B)

(3.17) \begin{equation} f_{m, k} (y, t_1) = A_m \cos (m \pi y) \frac {1}{\sqrt {2 \pi }} \textrm{e}^{\textrm{i} k t_1}, \end{equation}

where indices $k = 0, \pm 1, \pm 2, \ldots$ for non-zero $\omega$ , and $m = 0, 1, 2, \ldots$ . When $\omega =0$ for the steady case, we keep $k=0$ and then the whole solution framework reduces to Barton’s formulation. The imaginary unit is $\textrm{i}$ and $\textrm{i}^2=-1$ . The factor ${1}/{\sqrt {2 \pi }}$ for the Fourier series of $t_1$ and the undermined coefficient $A_m$ for the function of $y$ is used for the normalisation condition (3.20).

The corresponding eigenvalue is

(3.18) \begin{equation} \lambda _{m, k} = (m \pi )^2 + \textrm{i} \omega k. \end{equation}

Namely, the ‘eigenvalue’ can be a complex number instead of a real number, because the eigenvalue problem can be non-self-adjoint.

The inner product for the eigenvalue problem is defined as

(3.19) \begin{equation} \langle f_1, f_2 \rangle \triangleq \int ^{2 \pi }_0 \int ^1_0 f_1^{\star } (y, t_1) f_2 (y, t_1) \, \textrm{d} y \,\textrm{d} t_1, \end{equation}

where $f_1$ and $f_2$ are arbitrary functions. The superscript $\star$ denotes the complex conjugate. Note that due to the complex form of Fourier series, the orthonormal relationship (orthogonal and normalised) between eigenfunctions can be obtained:

(3.20) \begin{equation} \langle f_{m_1, k_1}, f_{m_2, k_2} \rangle = \delta _{m_1 m_2} \delta _{k_1 k_2} , \end{equation}

because

(3.21) \begin{equation} \int ^{2 \pi }_0 \left ( \frac {1}{\sqrt {2 \pi }} \textrm{e}^{\textrm{i} k_1 t_1} \right )^{\star } \frac {1}{\sqrt {2 \pi }} \textrm{e}^{\textrm{i} k_2 t_1} \; \textrm{d} t_1 = \frac {1}{2 \pi } \int ^{2 \pi }_0 \textrm{e}^{\textrm{i} (- k_1 + k_2) t_1} \; \textrm{d} t_1 = \delta _{0 (- k_1 + k_2)} = \delta _{k_1 k_2} . \end{equation}

Here, $\delta _{i j}$ is the Kronecker delta,

(3.22) \begin{align} \delta _{i j} = \left \{ \begin{array}{ll} 1, & i = j,\\ 0, & i \neq j. \end{array} \right. \end{align}

3.5. Solutions of moments

The final step of solving $P_n$ is to substitute the eigenvalues (3.18) and eigenfunctions (3.17) into the generic solution expressions of Barton (Reference Barton1983, Reference Barton1984). Some recent applications (Wang & Chen Reference Wang and Chen2017a ; Yang et al. Reference Yang, Jiang, Wu, Wang, Wu, Zhang and Zeng2021; Jiang et al. Reference Jiang, Zeng, Fu and Wu2022) have also given detailed expressions. Because the expressions are similar (only the notation is different), here, we just discuss $P_0$ and $P_1$ as illustrative examples, while the solutions of the higher-order moments are given in Appendix C. We mainly follow § 4.4 of Jiang et al. (Reference Jiang, Zeng, Fu and Wu2022) and reformulate the solution expressions using the current notation.

First, we can reindex the eigenvalues and eigenfunctions for convenience. Note that due to the introduction of $t_1$ , there are two indices $m$ and $k$ for the eigenfunctions ( $f_{m, k}(y, t_1)$ , with two dimensions). Hereafter, we use index $i$ for the list of eigenfunctions $\{f_i(y, t_1)\}_{i=0}^{\infty }$ and the list of eigenvalues $\{\lambda _i\}_{i=0}^{\infty }$ , sorted in the ascending order of the absolute value of $\lambda _{m,k}$ . The list will then be truncated at a finite number when evaluating the series summation in the following solution expressions. Notice that the index $i$ should not be confused with the imaginary unit $\textrm{i}$ .

For the zeroth-order moment, based on (3.1) of Barton (Reference Barton1983),

(3.23) \begin{equation} P_0(y, t_0, t_1) = \sum ^{\infty }_{i = 0} a_i \textrm{e}^{- \lambda _i t_0} f_i(y, t_1), \end{equation}

where $a_i$ is the coefficient related to the initial condition,

(3.24) \begin{equation} a_i \triangleq \langle f_i, P_0 | _{t_0 = 0} \rangle = \int ^{2 \pi }_0 \int ^1_0 f_i^{\star } (y, t_1) C_{{{ini}}} (y) \, \textrm{d} y\, \textrm{d} t_1, \quad i = 0, 1, \ldots . \end{equation}

Note that the initial condition (3.15d ) is independent of $t_1$ and thus the coefficient $a_i$ for the oscillatory term with $t_1$ is zero, namely, there is no oscillation in (3.23). Although we introduce a more complicated function base ( $f_i(y, t_1)$ ), $P_0$ is independent of $t_1$ . In addition, according to (3.14), the oscillatory term in $\textrm{e}^{- \lambda _i t_0}$ is cancelled out with the term of $t_1$ in $f_i(y, t_1)$ , namely, $\textrm{e}^{- \textrm{i} \omega k t_0} \textrm{e}^{\textrm{i} k t_1} =1$ . We obtain the same conclusion: no oscillation for $P_0$ . This is not surprising because the zeroth-order moment is the total amount of mass in each parallel streamline (cf. the concept of streamwise dispersion (Guan & Chen Reference Guan and Chen2024)), which is not affected by the oscillatory flow.

For the first-order moment $P_1$ , based on (2.9) by Barton (Reference Barton1984) (here, we recollect the terms with respect to $f_i$ and $t_0 f_i$ ),

(3.25) \begin{equation} P_1(y, t_0, t_1) = t_0 \sum ^{\infty }_{i=0} a_i B_{i, i} \textrm{e}^{- \lambda _i t_0} f_i(y, t_1) + \sum ^{\infty }_{\substack {i,j=0 \\ j \neq i}} \frac {a_j B_{i, j} (\textrm{e}^{- \lambda _j t_0} - \textrm{e}^{-\lambda _i t_0}) }{\lambda _i - \lambda _j} f_i(y, t_1), \end{equation}

where $B_{i, j}$ is the expansion coefficient related to the inner product with the velocity profile,

(3.26) \begin{equation} B_{i, j} \triangleq \langle f_i, u_1 f_j \rangle = \int ^{2 \pi }_0 \int ^1_0 f_i^{\star } (y, t_1) u_1 (y, t_1) f_j (y, t_1) \, \textrm{d} y\, \textrm{d} t_1 \end{equation}

for $i, j = 0, 1, \ldots$ . Similar to $P_0$ , note that the first term in (3.25) is independent of $t_1$ . The second term may be oscillatory because it contains $a_j B_{i, j} \textrm{e}^{- \lambda _j t_0} f_i(y, t_1)$ (when $i \neq j$ ).

Note that in practice, the oscillatory velocity profile can be generally expressed in the form of the Fourier series of $t_1$ (or $t$ in the original transport problem). Therefore, the following formula for the integration of the product of three eigenfunctions with respect to $t_1$ can be helpful,

(3.27) \begin{equation} \int ^{2 \pi }_0 \frac {1}{\sqrt {2 \pi }} \textrm{e}^{\textrm{i} k_1 t_1} \frac {1}{\sqrt {2 \pi }} \textrm{e}^{\textrm{i} k_2 t_1} \frac {1}{\sqrt {2 \pi }} \textrm{e}^{\textrm{i} k_3 t_1} \; \textrm{d} t_1 = \frac {1}{\sqrt {2 \pi }} \delta _{0 (k_1 + k_2 + k_3)} . \end{equation}

Namely, when $k_1 + k_2 + k_3 = 0$ , the integration result is $ {1}/{\sqrt {2 \pi }}$ .

The solutions for the second- and higher-order moments $P_n$ are given in Appendix C. Once we solve $P_n$ , the solution for $C_n$ can be immediately obtained through (3.14). Then, the dispersion characteristics can be calculated according to their definitions in § 2.2.

4.1. Velocity profile of oscillatory Couette flow

A sketch of the transport problem is shown in figure 1. See also the recent experiment by Ding et al. (Reference Ding, Hunt, McLaughlin and Woodie2021). The bottom wall ( $y^{\ast }=0$ ) is at rest, while the upper wall ( $y^{\ast }=L^{\ast }$ ) is executing an oscillatory motion with velocity $U_A^{\ast } \cos (\omega ^{\ast } t^{\ast })$ , where $U_A^{\ast }$ is the amplitude of the velocity oscillation and $\omega ^{\ast }$ is the angular frequency.

Based on the dimensionless parameters and variables in (2.2), we introduce two additional dimensionless parameters:

(4.1) \begin{equation} \textit{Sc} =\frac {\nu ^{\ast }}{D^{\ast }}, \quad \textit{Wo} = L^{\ast } \sqrt {\frac {\omega ^{\ast }}{\nu ^{\ast }}}. \end{equation}

Notice that $\omega = \textit{Sc} \textit{Wo}^2$ , where $\textit{Sc}$ is the Schmidt number and $\textit{Wo}$ is the Womersley number. Here, $\nu ^{\ast }$ is the kinematic viscosity of the fluid. For the choice of the characteristic velocity in (2.2), we use the velocity oscillation scale, namely $U^{\ast }_{c} = U_A^{\ast }$ .

The dimensionless Navier–Stokes equation and the boundary conditions for this case are

(4.2) \begin{equation} \frac {1}{\textit{Sc}} \frac {\partial u}{\partial t} = \frac {\partial ^2 u}{\partial y^2}, \quad u |_{y = 1} = \cos (\omega t), \quad u |_{y = 0} = 0. \end{equation}

The oscillatory velocity profile has a solution (Drazin & Riley Reference Drazin and Riley2006; Ding et al. Reference Ding, Hunt, McLaughlin and Woodie2021):

(4.3) \begin{equation} u (y, t) = {\textrm {Re}}\left [\textrm{e}^{\textrm{i} \omega t} \frac {\sinh \left ( \textrm{e}^{({ {\pi }/{4})} \textrm{i}} \textit{Wo} y \right )}{\sinh \left ( \textrm{e}^{({ {\pi }/{4})} \textrm{i}} \textit{Wo} \right )} \right ]. \end{equation}

Notice that $\textrm{e}^{({ {\pi }/{4})}\textrm{i}} = \sqrt {\textrm{i}} = ({\sqrt {2}}/{2}) (1+\textrm{i})$ . Here, $\textrm {Re}$ means taking the real part of the complex number.

For the two-time-variable system (3.6), we need to rewrite the velocity profile (4.3) in the form of the Fourier series of $t_1$ (see also the definition of $u_1$ in (3.2)) to calculate the expansion coefficient $B_{i, j}$ in (3.26):

(4.4) \begin{equation} u_1 (y, t_1) = \sum _{k = \pm 1} \frac {\sqrt {2 \pi }}{2} \left [ \frac {1}{\sqrt {2 \pi }} \textrm{e}^{\textrm{i} k t_1} \frac {\sinh \left ( \textrm{e}^{({ {\pi }/{4})} \textrm{i}} \sqrt {k} \textit{Wo} y \right )}{\sinh \left ( \textrm{e}^{({ {\pi }/{4})} \textrm{i}} \sqrt {k} \textit{Wo} \right )} \right ], \end{equation}

see also (5) of Ding et al. (Reference Ding, Hunt, McLaughlin and Woodie2021). Note that only one frequency is considered here. For the multiple-frequency case with more Fourier components (see also the paper by Vedel & Bruus Reference Vedel and Bruus2012), the solution procedure is similar but with more expansion coefficients for the ‘time’ variable $t_1$ .

Based on the experimental study by Ding et al. (Reference Ding, Hunt, McLaughlin and Woodie2021), we consider the transport problem of illuminated fluorescein dye in water. We perform the numerical simulation and analytical investigation based on the parameters listed in table 1. More cases of different parameter set-ups can be found in figure 9. In this case, the temporal evolution of the velocity profile is shown in figure 2. Due to the small Womersley number, the velocity profile is almost linear along the channel, namely, the boundary layer effect is weak.

Table 1.

Dimensional and dimensionless parameters in (2.2) and (4.1) for the transport of illuminated fluorescein dye in an oscillatory Couette flow of water. The parameters are partially taken from the experimental study by Ding et al. (Reference Ding, Hunt, McLaughlin and Woodie2021).

Figure 2.

Dimensionless velocity profile $u(y,t)$ of the oscillatory Couette flow with the parameters listed in table 1.

4.2. Comparison with numerical simulation

To verify our analytical solution of the concentration moments, the Brownian dynamics simulation method (also known as the Monte Carlo method or random walk simulation) is adopted to obtain the numerical solutions. This simulation approach is widely applicable in various fields (Elhanati et al. Reference Elhanati, Dror and Berkowitz2023; Wu et al. Reference Wu, Zhan, Zhu, Zhang, Chang and Fu2026), as discussed in the comprehensive review by Berkowitz et al. (Reference Berkowitz, Cortis, Dentz and Scher2006). The numerical algorithm can be found in the recent paper by Ding et al. (Reference Ding, Hunt, McLaughlin and Woodie2021, § 3) (see also Wang et al. Reference Wang, Jiang, Chen, Tao and Li2021, Appendix A; Jiang & Chen Reference Jiang and Chen2021, Appendix B).

First, we need to find the corresponding stochastic differential equations (SDEs) for the original transport problem (2.1):

(4.5a) \begin{align} \textrm{d} X_t &= u (Y_t, t) \; \textrm{d} t + \frac {\sqrt {2}}{\textit{Pe}} \; \textrm{d} W_1, \end{align}

(4.5b) \begin{align} \textrm{d} Y_t &= \sqrt {2}\; \textrm{d} W_2, \end{align}

where $X_t$ and $Y_t$ are the random horizontal and transverse coordinates of the position of a particle, respectively. Here, $W_1$ and $W_2$ are independent standard Brownian motions.

In our simulation, the Euler–Maruyama method is adopted to numerically solve the SDEs. The details are given in Appendix D. As shown in figure 8, after checking the independence of the time step ( $\Delta t$ ) and the number of particles ( $N$ ), we take $\Delta t = 10^{-3}$ and $N = 10^5$ . For the point-source release considered in figure 3, all particles are initially released at $x_0=0$ and $y_0 = 1$ , which corresponds to

(4.6) \begin{equation} C_{{{ini}}} = \delta (y-1) \end{equation}

in (2.1c ). The trajectories of particles are recorded. The positions of particles are used to statistically calculate concentration moments.

Figure 3.

Comparison between the analytical solution of concentration moments and the numerical result by Brownian dynamics simulation: (a) the first-order moment of mean concentration distribution; (b) the variance of $\bar {C}$ ; (c) the skewness of $\bar {C}$ ; and (d) the kurtosis of $\bar {C}$ . $N_y$ denotes the number of eigenfunctions of $y$ used in the truncated summation of the analytical solution. The initial condition is a point-source release at $y_0 =1$ . Other parameters are listed in table 1.

For the analytical solution, the infinite summations in the solution expressions (such as (3.23)) need to be truncated for calculation or numerical evaluation (Barton Reference Barton1983). Only the first few eigenvalues (3.18), with indices $k$ and $m$ , and eigenfunctions (3.17) are used for the summations. First, for the Fourier series of $t_1$ for the oscillation in (3.17), we take $| k | \leqslant 4$ because we only consider concentration moments up to the fourth order for the discussion of kurtosis and the velocity profile is of one frequency. Higher-order moments will require more Fourier components of $t_1$ .

Second, for the trigonometric functions of $y$ in (3.17), the first $N_y$ eigenfunctions of $y$ (i.e. $m \leqslant N_y-1$ , because $m$ starts from 0) are used. We check the accuracy of the truncated summations by comparing the results of $N_y=4, 6, 8, 10$ . Take $N_y=8$ as an example. The total number of eigenfunctions $f_{m,k}$ is $8 \times (2\times 4 + 1) = 72$ . Then, we recollect all the eigenfunctions and eigenvalues into single-indexed lists $\{f_i(y, t_1)\}_{i=0}^{71}$ and $\{\lambda _i\}_{i=0}^{71}$ , and substitute them into the solution expressions of concentration moments.

As shown in figure 3, the analytical solutions of the transport characteristics of the mean concentration distribution $\overline {C}$ agree well with the numerical results, showing the correctness of the proposed auxiliary-time extension procedure. More comparisons under different parameter set-ups can be found in figure 9.

The truncation error may be large when $t$ is small, because a point-source-release initial condition (4.6) is considered. As shown in figure 3(c, d) for skewness and kurtosis, the curves of the results with $N_y=4$ and $N_y=6$ deviate from the numerical results in the initial stage. Namely, more eigenfunctions may be required for the higher-order moments. Generally speaking, the Dirac delta function in the initial condition (4.6) requires an infinite number of eigenfunctions of $y$ for the series summations in § 3.5 when $t\rightarrow 0$ . Fortunately, the truncation error will exponentially decrease as time increases (Barton Reference Barton1984). Note that the real part of the eigenvalue in (3.18) is $(m \pi )^2$ . Thus, the series summations in (3.23), (3.25), (C1) converge exponentially (with $\exp [-(m \pi )^2 t]$ ), and the truncation error diminishes with a large enough $m$ . The results of $N_y=8$ and $N_y=10$ are close to the numerical results, as shown in figure 3. Therefore, we use $N_y=8$ in the following discussions. The linear–log plot (with a logarithmic scale of time) is used to highlight the alignment of curves in the initial short-time period, showing that the current configurations for the analytical and numerical solution procedures are sufficient.

We also compare the result of the Taylor dispersion coefficient ( $D_T = \lim _{t \rightarrow \infty } \bar {\sigma }^2 / (2 t)$ ) with the value obtained by Ding et al. (Reference Ding, Hunt, McLaughlin and Woodie2021, table 2, Trial 1). For $N_y=8$ , the dimensional $D_T$ is $5.38\times 10^{-9}\textrm {m}^2\,\textrm {s}^{-1}$ , which is consistent with the theoretical value of $5.39\times 10^{-9}\textrm {m}^2\,\textrm {s}^{-1}$ given by Ding et al. (Reference Ding, Hunt, McLaughlin and Woodie2021), although both values deviate from the experimental result of $4.26\times 10^{-9}\textrm {m}^2\,\textrm {s}^{-1}$ .

6.1. Solution procedure to address phase shift

First, we follow Ding et al. (Reference Ding, Hunt, McLaughlin and Woodie2021) and add a phase shift $s$ in the oscillatory velocity field as

(6.1) \begin{equation} u_s (y, t) \triangleq u (y, t + s). \end{equation}

Then, replace the original velocity field $u$ in (2.1a ) with the shifted velocity field $u_s$ .

Typically, one needs to re-solve the whole problem of the concentration moments to impose the phase shift. However, with the auxiliary-time extension method, we can directly obtain the solution with a simple transformation of the time variable $t_1$ . Notice that $t_1$ is imposed in (3.1) as $t_1 = \omega t$ . To address the phase shift, we impose the following time variables instead:

(6.2) \begin{equation} t_0 = t, \quad t_{1 s} = \omega (t + s). \end{equation}

Notice that $t_{1 s} = t_1 + \omega s$ . Then, similar to (3.2), we define $u_{s1}$ to replace the oscillatory velocity field for the shifted two-time-variable system and

(6.3) \begin{align} u_{s1} (y, t_{1 s}) & \triangleq u_s (y, t) = u (y, t + s)\notag \\ & = u \left(y, \frac {t_{1 s}}{\omega }-s + s\right) = u \left(y, \frac {t_{1 s}}{\omega }\right)\\ & = u_1 \left(y, t_{1 s}\right).\notag \end{align}

Namely, the function $u_{s1}$ is in the same form as $u_1$ . Thus, the governing equation for the shifted two-time-variable system is the same as that of the original system in (3.6).

To sum up, there is no need to re-solve the whole problem for the phase shift in flow velocity. The solution of the moments for the shifted two-time-variable system can be obtained by replacing $t_1$ with $t_{1 s}$ in the solution of the original two-time-variable system given in § 3.5. Namely,

(6.4) \begin{equation} P_{s n} (y, t_0, t_{1 s}) = P_n (y, t_0, t_{1 s}), \end{equation}

where $P_{s n}$ is the $n$ th-order moment of the shifted two-time-variable system. Based on (3.14), the $n$ th-order concentration moments $C_{s n}$ of the shifted problem is

(6.5) \begin{equation} C_{s n}(y, t) = P_{s n} (y, t_0, t_{1 s}) = P_n (y, t, \omega (t + s)). \end{equation}

Note that one could not obtain $C_{s n}(y,t)$ through $C_{n}(y,t)$ by simply replacing $t$ with $t + s$ . Fortunately, this could be done by $ P_n$ with $t_1$ replaced by $t_1 + \omega s$ , which can be considered as a benefit of the auxiliary-time extension method. In addition, as discussed by Ding (Reference Ding2025), the decay rate to the Taylor dispersion regime can be characterised by the non-zero eigenvalue with the smallest real part, $1/\textrm {Re}(\lambda )$ . Replacing $t_1$ would not influence the eigenvalue (3.18). Therefore, the decay rate to the Taylor dispersion regime is not affected by the phase shift.

We also perform the numerical verification for the phase shift effect, as shown in figure 11. The analytical result by (6.5) agrees well with the numerical result.

6.2. Results of transport characteristics

To illustrate the effect of the phase shift on the transport characteristics, we consider the same oscillatory velocity profile as that in previous sections, but with different phase shift $s=0, \pi /4, \pi , \pi /2, 3\pi /4$ . A point-source release is considered ( $y_0 = 1$ ). Note that the case with $s=\pi$ is the reverse case of the case with $s=0$ , due to the symmetry of the transport problem. Thus, only four cases are considered.

For the first-order moment of $\bar {C}$ and the effective drift, it can be seen in figure 6(a, b) that the effect of the phase shift is remarkable. Note that the effective drift is the solute-weighted average velocity, while the overall transverse distribution is not affected by the velocity field. Thus, the phase shift of the velocity will directly be reflected in the effective drift. It then leads to the accumulative change in the mean position of the solute cloud as time increases.

Figure 6.

Temporal evolution of (a) the mean longitudinal position, (b) the effective drift velocity, (c) the variance of mean concentration distribution and (d) the effective diffusion coefficient for different phase shift $s$ . The initial condition is a point-source release at $y_0 = 1$ . Other parameters are listed in table 1.

Figure 6(c, d) shows the phase shift effect on the variance of $\bar {C}$ and the effective diffusion coefficient. Note that $D_e$ is initially zero because of the instantaneous release. The phase shift will change the time when $D_e$ reaches its maximum value. Then, $D_e$ will decay to the oscillatory state and reach the Reference TaylorTaylor dispersion regime. For the variance, similar to the effect of the initial condition discussed in § 5, the accumulative influence of the phase shift is weak and all the curves are nearly close to each other as $t$ increases.

As discussed by Ding et al. (Reference Ding, Hunt, McLaughlin and Woodie2021), the phase shift can change the sign of skewness, which is clearly observed in figure 7(a). Unlike the line-source case of Ding et al. (Reference Ding, Hunt, McLaughlin and Woodie2021), the point-source case discussed here has a large non-zero skewness in the initial stage of the transport process. Thus, the oscillation of skewness is not apparent. However, the phase shift of the velocity field can still make a qualitative change in the asymmetry of the mean concentration distribution. If we compare the cases with $s=0$ and $s=\pi$ (not shown here), then the curves will be flipped with respect to the $t$ -axis due to the symmetric property of the transport problem.

Figure 7.

Temporal evolution of (a) skewness and (b) kurtosis of the mean concentration distribution for different phase shift $s$ . The initial condition is a point-source release at $y_0 = 1$ . Other parameters are listed in table 1.

For kurtosis, the phase shift will change the sign of $\textit{Ku}$ , as shown in figure 7(b). The phase shift of the velocity field can be transferred to kurtosis, as the oscillation of $\textit{Ku}$ is apparent during the initial stage of the transport process. For the four considered cases, the temporal evolution curves have many intersections with each other even in the first period, showing the phase shift effect is complex.