3.1 Moment equation and solution in the Laplace space

Firstly, following the method of moments developed by Aris (Reference Aris1956), multiply (2.7a ) by $x^{n}$ and integrate in the $x$ direction to obtain the equation for $c_{n}^{\ast }(y,t)$ ,

(3.3) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}c_{n}^{\ast }}{\unicode[STIX]{x2202}t}+\mathit{Pe}\,u\int _{-\infty }^{\infty }x^{n}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}x}\,\text{d}x=\int _{-\infty }^{\infty }x^{n}\frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}x^{2}}\,\text{d}x+\frac{\unicode[STIX]{x2202}^{2}c_{n}^{\ast }}{\unicode[STIX]{x2202}y^{2}},\end{eqnarray}$$

where $c_{n}^{\ast }=\int _{-\infty }^{\infty }x^{n}c\,\text{d}x$ is the $n\text{th}$ longitudinal moment of concentration in the filament through $y$ , which is not yet transversely averaged. The moments $m_{n}$ introduced above are the transverse averages of $c_{n}^{\ast }$ . After integration by parts and noting that the concentration and all of its derivatives vanish at infinity, we have

(3.4a ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}c_{n}^{\ast }}{\unicode[STIX]{x2202}t}-n\mathit{Pe}\,uc_{n-1}^{\ast }=n(n-1)c_{n-2}^{\ast }+\frac{\unicode[STIX]{x2202}^{2}c_{n}^{\ast }}{\unicode[STIX]{x2202}y^{2}},\end{eqnarray}$$

where $c_{-1}^{\ast }=c_{-2}^{\ast }=0$ . Similarly, the boundary conditions (2.7b ) and (2.7c ) give

(3.4b ) $$\begin{eqnarray}\displaystyle & \displaystyle -\frac{\unicode[STIX]{x2202}c_{n}^{\ast }}{\unicode[STIX]{x2202}y}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FE}_{n}^{\ast }}{\unicode[STIX]{x2202}t}=k_{a}c_{n}^{\ast }-k_{d}\unicode[STIX]{x1D6FE}_{n}^{\ast }\quad \text{at}~y=1, & \displaystyle\end{eqnarray}$$

(3.4c ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}c_{n}^{\ast }}{\unicode[STIX]{x2202}y}=0\quad \text{at}~y=0, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}_{n}^{\ast }$ is defined as the $n\text{th}$ longitudinal moment of the surface concentration.

Laplace transformation in time reduces (3.4) to a system of ordinary differential equations (ODEs) involving only the transformed variable,

(3.5) $$\begin{eqnarray}{\hat{c}}_{n}^{\ast }(y,s)=\mathscr{L}\{c_{n}^{\ast }\}(s)=\int _{0}^{\infty }c_{n}^{\ast }\text{e}^{-st}\,\text{d}t,\end{eqnarray}$$

because the transformed longitudinal moments of surface concentration $\hat{\unicode[STIX]{x1D6FE}}_{n}^{\ast }=\mathscr{L}\{\unicode[STIX]{x1D6FE}_{n}^{\ast }\}$ in the boundary condition can be eliminated. In Laplace space, (3.4) are given by

(3.6a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}^{2}{\hat{c}}_{n}^{\ast }}{\unicode[STIX]{x2202}y^{2}}=s{\hat{c}}_{n}^{\ast }-c_{n}^{\ast }(t=0)-n\mathit{Pe}\,u{\hat{c}}_{n-1}^{\ast }-n(n-1){\hat{c}}_{n-2}^{\ast }, & \displaystyle\end{eqnarray}$$

(3.6b ) $$\begin{eqnarray}\displaystyle & \displaystyle -\frac{\unicode[STIX]{x2202}{\hat{c}}_{n}^{\ast }}{\unicode[STIX]{x2202}y}=s\hat{\unicode[STIX]{x1D6FE}}_{n}^{\ast }-\unicode[STIX]{x1D6FE}_{n}^{\ast }(t=0)=k_{a}{\hat{c}}_{n}^{\ast }-k_{d}\hat{\unicode[STIX]{x1D6FE}}_{n}^{\ast }\quad \text{at}~y=1, & \displaystyle\end{eqnarray}$$

(3.6c ) $$\begin{eqnarray}\displaystyle & \displaystyle -\frac{\unicode[STIX]{x2202}{\hat{c}}_{n}^{\ast }}{\unicode[STIX]{x2202}y}=0\quad \text{at}~y=0. & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D6FE}_{n}^{\ast }(t=0)=0$

$\hat{\unicode[STIX]{x1D6FE}}_{n}^{\ast }$

(3.7) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D6FE}}_{n}^{\ast }=\frac{k_{a}}{k_{d}+s}{\hat{c}}_{n}^{\ast },\end{eqnarray}$$

so that (3.6b ) turns into a Robin-type boundary condition,

(3.8) $$\begin{eqnarray}-\frac{\unicode[STIX]{x2202}{\hat{c}}_{n}^{\ast }}{\unicode[STIX]{x2202}y}=\frac{k_{a}s}{k_{d}+s}{\hat{c}}_{n}^{\ast }\quad \text{at}~y=1.\end{eqnarray}$$

The $\unicode[STIX]{x1D6FF}$ -function initial distribution of solute leads to the following initial conditions:

(3.9a,b ) $$\begin{eqnarray}c_{0}^{\ast }=1,\quad c_{1}^{\ast }=c_{2}^{\ast }=0\quad \text{at}~t=0.\end{eqnarray}$$

Therefore, (3.6a ), together with boundary conditions (3.6c ) and (3.8), gives the following system of ODEs for ${\hat{c}}_{0}^{\ast }$ , ${\hat{c}}_{1}^{\ast }$ and ${\hat{c}}_{2}^{\ast }$ :

(3.10a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}^{2}{\hat{c}}_{0}^{\ast }}{\text{d}y^{2}}=s{\hat{c}}_{0}^{\ast }-1, & \displaystyle\end{eqnarray}$$

(3.10b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}^{2}{\hat{c}}_{1}^{\ast }}{\text{d}y^{2}}=s{\hat{c}}_{1}^{\ast }-\mathit{Pe}\,u{\hat{c}}_{0}^{\ast }, & \displaystyle\end{eqnarray}$$

(3.10c ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}^{2}{\hat{c}}_{2}^{\ast }}{\text{d}y^{2}}=s{\hat{c}}_{2}^{\ast }-2\mathit{Pe}\,u{\hat{c}}_{1}^{\ast }-2{\hat{c}}_{0}^{\ast }, & \displaystyle\end{eqnarray}$$

(3.11a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}{\hat{c}}_{n}^{\ast }}{\text{d}y}=-\frac{k_{a}s}{k_{d}+s}{\hat{c}}_{n}^{\ast }\quad \text{at}~y=1, & \displaystyle\end{eqnarray}$$

(3.11b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}{\hat{c}}_{n}^{\ast }}{\text{d}y}=0\quad \text{at}~y=0, & \displaystyle\end{eqnarray}$$

$n=0,1,2$

However, it is not the analytical solutions of ${\hat{c}}_{n}^{\ast }(y,s)$ but the transverse-averaged moments $\hat{m}_{n}(s)=\int _{0}^{1}\hat{c_{n}^{\ast }}\,\text{d}y$ that are of interest here. For instance, $\hat{m}_{0}$ has the form

(3.12) $$\begin{eqnarray}\hat{m}_{0}=\frac{1}{s}-\frac{k_{a}\sinh (\sqrt{s})}{\sqrt{s}(k_{a}s\cosh (\sqrt{s})+\sqrt{s}\sinh (\sqrt{s})(k_{d}+s))}.\end{eqnarray}$$

The analytical forms of $\hat{m}_{1}$ and $\hat{m}_{2}$ are complex (given in the supplementary materials available at https://doi.org/10.1017/jfm.2017.546), but both of them and $\hat{m}_{0}$ can be written in a general form as

(3.13) $$\begin{eqnarray}\hat{m}_{n}(s)=\frac{N_{n}(s)}{E(s)^{(n+1)}}\quad \text{for}~n=0,1,2,\end{eqnarray}$$

where the denominator $E(s)$ is given by

(3.14) $$\begin{eqnarray}E(s)=k_{a}s\cosh (\sqrt{s})+(k_{d}+s)\sqrt{s}\sinh (\sqrt{s}),\end{eqnarray}$$

which is a transcendental function of $s$ and includes all the singularities of the moments. The numerators $N_{n}(s)$ are complex functions of $s$ obtained by a computer algebra system (MathWorks 2012). The transcendental function $E(s)$ has two important properties:

1. (i) There are only first-order singularities in $E(s)$ , and thus $\hat{m}_{0}$ , $\hat{m}_{1}$ and $\hat{m}_{2}$ have only first-, second- and third-order singularities, respectively. This helps to employ the residue theorem for the inverse Laplace transform.

2. (ii) All the singularities of $E(s)=0$ fall along the negative axis, and thus substituting $s=-p^{2}$ , where $p$ is a real positive number, leads to a transcendental equation of $p$ in real space (where $\tanh (\text{i}p)=\text{i}\tan (p)$ is used),

   (3.15) $$\begin{eqnarray}\tan (p)(p^{2}-k_{d})-k_{a}p=0.\end{eqnarray}$$
   Equation (3.15) has an infinite number of roots $p_{k},~k=0,1,\ldots ,\infty$ . These roots $p_{k}$ correspond to characteristic decay rates of the moments, and the lowest-order term with the smallest root, i.e. $p_{0}=0$ , dominates the behaviour at late times.

3.2 Inverse Laplace transform by the residue theorem

The inverse Laplace transform of the moments can be written as the Bromwich integral,

(3.16) $$\begin{eqnarray}m_{n}(t)=\frac{1}{2\unicode[STIX]{x03C0}\text{i}}\int _{{\mathcal{C}}}\hat{m}_{n}(s)\text{e}^{st}\,\text{d}s,\end{eqnarray}$$

where $\text{i}=\sqrt{-1}$ and ${\mathcal{C}}$ is a contour chosen so that all the singularities of $\hat{m}_{n}(s)$ are to the left of it. Further, if we apply the residue theorem to the above integral, we have

(3.17) $$\begin{eqnarray}m_{n}(t)=\mathop{\sum }_{k=0}^{\infty }R_{k},\end{eqnarray}$$

where $R_{k}$ are the residues of $\hat{m}_{n}\text{e}^{st}$ and can be calculated as

(3.18) $$\begin{eqnarray}R_{k}=\frac{1}{(l-1)!}\lim _{s\rightarrow s_{k}}\frac{\text{d}^{l-1}}{\text{d}s^{l-1}}(\hat{m}_{n}\text{e}^{st}(s-s_{k})^{l}),\end{eqnarray}$$

where $l$ is the order of the $k\text{th}$ singularity or pole  $s_{k}$ .

Since $\hat{m}_{0}$ , $\hat{m}_{1}$ and $\hat{m}_{2}$ have only first-, second- and third-order singularities, respectively, we have

(3.19a ) $$\begin{eqnarray}\displaystyle & \displaystyle m_{0}(t)=\mathop{\sum }_{k=0}^{\infty }\lim _{s\rightarrow s_{k}}(s-s_{k})\hat{m}_{0}\exp (st)=\mathop{\sum }_{k=0}^{\infty }a_{k}\exp (-p_{k}^{2}t), & \displaystyle\end{eqnarray}$$

(3.19b ) $$\begin{eqnarray}\displaystyle & \displaystyle m_{1}(t)=\mathop{\sum }_{k=0}^{\infty }\lim _{s\rightarrow s_{k}}\frac{\text{d}}{\text{d}s}[(s-s_{k})^{2}\hat{m}_{1}\exp (st)]=\mathop{\sum }_{k=0}^{\infty }b_{k}^{(1)}\exp (-p_{k}^{2}t)+b_{k}^{(2)}t\exp (-p_{k}^{2}t), & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(3.19c ) $$\begin{eqnarray}\displaystyle \displaystyle m_{2}(t) & = & \displaystyle \mathop{\sum }_{k=0}^{\infty }\frac{1}{2}\lim _{s\rightarrow s_{k}}\frac{\text{d}^{2}}{\text{d}s^{2}}[(s-s_{k})^{3}\hat{m}_{2}\exp (st)]\nonumber\\ \displaystyle & = & \displaystyle \mathop{\sum }_{k=0}^{\infty }c_{k}^{(1)}\exp (-p_{k}^{2}t)+c_{k}^{(2)}t\exp (-p_{k}^{2}t)+c_{k}^{(3)}t^{2}\exp (-p_{k}^{2}t),\end{eqnarray}$$

(3.20a ) $$\begin{eqnarray}\displaystyle & \displaystyle a_{k}=\lim _{s\rightarrow s_{k}}(s-s_{k})\hat{m}_{0}, & \displaystyle\end{eqnarray}$$

(3.20b ) $$\begin{eqnarray}\displaystyle & \displaystyle b_{k}^{(1)}=\lim _{s\rightarrow s_{k}}\frac{\text{d}}{\text{d}s}[(s-s_{k})^{2}\hat{m}_{1}], & \displaystyle\end{eqnarray}$$

(3.20c ) $$\begin{eqnarray}\displaystyle & \displaystyle b_{k}^{(2)}=\lim _{s\rightarrow s_{k}}(s-s_{k})^{2}\hat{m}_{1}, & \displaystyle\end{eqnarray}$$

(3.20d ) $$\begin{eqnarray}\displaystyle & \displaystyle c_{k}^{(1)}=\frac{1}{2}\lim _{s\rightarrow s_{k}}\frac{\text{d}^{2}}{\text{d}s^{2}}[(s-s_{k})^{3}\hat{m}_{2}], & \displaystyle\end{eqnarray}$$

(3.20e ) $$\begin{eqnarray}\displaystyle & \displaystyle c_{k}^{(2)}=\lim _{s\rightarrow s_{k}}\frac{\text{d}}{\text{d}s}[(s-s_{k})^{3}\hat{m}_{2}], & \displaystyle\end{eqnarray}$$

(3.20f ) $$\begin{eqnarray}\displaystyle & \displaystyle c_{k}^{(3)}=\frac{1}{2}\lim _{s\rightarrow s_{k}}(s-s_{k})^{3}\hat{m}_{2}. & \displaystyle\end{eqnarray}$$

$\hat{m}_{0}$

$\hat{m}_{1}$

$\hat{m}_{2}$

(3.21a ) $$\begin{eqnarray}\displaystyle & \displaystyle a_{k}=\frac{N_{0}}{T_{1}}, & \displaystyle\end{eqnarray}$$

(3.21b ) $$\begin{eqnarray}\displaystyle & \displaystyle b_{k}^{(1)}=\frac{T_{1}N_{1}^{\prime }-2T_{2}N_{1}}{T_{1}^{3}}, & \displaystyle\end{eqnarray}$$

(3.21c ) $$\begin{eqnarray}\displaystyle & \displaystyle b_{k}^{(2)}=\frac{N_{1}}{T_{1}^{2}}, & \displaystyle\end{eqnarray}$$

(3.21d ) $$\begin{eqnarray}\displaystyle & \displaystyle c_{k}^{(1)}=\frac{(12T_{2}^{2}-6T_{1}T_{3})N_{2}-6T_{1}T_{2}N_{2}^{\prime }+T_{1}^{2}N_{2}^{\prime \prime }}{2T_{1}^{5}}, & \displaystyle\end{eqnarray}$$

(3.21e ) $$\begin{eqnarray}\displaystyle & \displaystyle c_{k}^{(2)}=\frac{T_{1}N_{2}^{\prime }-3T_{2}N_{2}}{T_{1}^{4}}, & \displaystyle\end{eqnarray}$$

(3.21f ) $$\begin{eqnarray}\displaystyle & \displaystyle c_{k}^{(3)}=\frac{N_{2}}{2T_{1}^{3}}, & \displaystyle\end{eqnarray}$$

$N_{n}^{\prime }=\text{d}N_{n}/\text{d}s$

$N_{n}^{\prime \prime }=\text{d}^{2}N_{n}/\text{d}s^{2}$

$s=s_{k}$

$T_{n}=E^{(n)}/n!$

$n\text{th}$

$E(s)$

$s=s_{k}$

$p_{k}$

$s_{k}=-p_{k}^{2}$

$a_{k}$

$b_{k}^{(1)}$

$b_{k}^{(2)}$

$c_{k}^{(1)}$

$c_{k}^{(2)}$

$c_{k}^{(3)}$

To summarize, for a given $k_{a}$ and $k_{d}$ , (3.15) is first solved for a series of $p_{k}$ , which are substituted into (3.21) to obtain the coefficients $a_{k}$ , $b_{k}$ and $c_{k}$ . The normalized longitudinal moments $M_{0}$ , $M_{1}$ and $M_{2}$ , the transport velocity $v$ and the dispersion coefficient $D_{L}$ are then determined by definitions (3.1) and (3.2).

3.3 Reduction to previous results

In the long-time limit, when the zeroth-order terms dominate, the transport velocity and dispersion coefficient are

(3.22a ) $$\begin{eqnarray}\displaystyle & \displaystyle v_{0}=\frac{b_{0}^{(2)}}{a_{0}}=\mathit{Pe}\frac{k_{d}}{k_{a}+k_{d}}=\mathit{Pe}\frac{1}{k+1}, & \displaystyle\end{eqnarray}$$

(3.22b ) $$\begin{eqnarray}\displaystyle & \displaystyle D_{0}=\frac{1}{2}\left(\frac{c_{0}^{(2)}}{a_{0}}-\frac{2b_{0}^{(1)}b_{0}^{(2)}}{{a_{0}}^{2}}\right)=\frac{1}{1+k}+\mathit{Pe}^{2}\frac{2}{105}\frac{1+9k+25.5k^{2}}{(1+k)^{3}}+\frac{\mathit{Pe}^{2}}{k_{d}}\frac{k}{(1+k)^{3}}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

$k>0$

At early but finite time, the first-order terms dominate and lead to an asymptotic velocity $v_{1}$ and dispersion coefficient $D_{1}$ given as

(3.23a,b ) $$\begin{eqnarray}v_{1}=\frac{b_{1}^{(2)}}{a_{1}}\quad \text{and}\quad D_{1}=\frac{1}{2}\left(\frac{c_{1}^{(2)}}{a_{1}}-\frac{2b_{1}^{(1)}b_{1}^{(2)}}{a_{1}^{2}}\right).\end{eqnarray}$$

In the limiting case of $k_{d}=0$ analysed by Lungu & Moffatt (Reference Lungu and Moffatt1982), the zeroth-order coefficients of the moments vanish, i.e. $a_{0}=b_{0}^{(2)}=b_{0}^{(1)}=c_{0}^{(3)}=c_{0}^{(2)}=c_{0}^{(1)}=0$ . Therefore, the first-order terms dominate and lead to the following asymptotic transport velocity and dispersion coefficient,

(3.24a ) $$\begin{eqnarray}v_{LM}=\frac{\mathit{Pe}(4k_{a}^{2}p_{1}^{2}+3k_{a}^{2}+3k_{a}+4p_{1}^{4}-3p_{1}^{2})}{4p_{1}^{2}(k_{a}^{2}+k_{a}+p_{1}^{2})},\end{eqnarray}$$

(3.24b ) $$\begin{eqnarray}\displaystyle D_{LM} & = & \displaystyle 1+ (\!\mathit{Pe}^{2}\! (\!-8k_{a}^{6}p_{1}^{4}+150k_{a}^{6}p_{1}^{2}-315k_{a}^{6}-56k_{a}^{5}p_{1}^{4}+750k_{a}^{5}p_{1}^{2}\nonumber\\ \displaystyle & & \displaystyle -\,945k_{a}^{5}-24k_{a}^{4}p_{1}^{6}+282k_{a}^{4}p_{1}^{4}+555k_{a}^{4}p_{1}^{2}-945k_{a}^{4}-192k_{a}^{3}p_{1}^{6}\nonumber\\ \displaystyle & & \displaystyle +\,1560k_{a}^{3}p_{1}^{4}-540k_{a}^{3}p_{1}^{2}-315k_{a}^{3}+24k_{a}^{2}p_{1}^{8}-78k_{a}^{2}p_{1}^{6}+1455k_{a}^{2}p_{1}^{4}\nonumber\\ \displaystyle & & \displaystyle -\,495k_{a}^{2}p_{1}^{2}-136k_{a}p_{1}^{8}+810k_{a}p_{1}^{6}+225k_{a}p_{1}^{4}-8p_{1}^{10}-210p_{1}^{8}+585p_{1}^{6}\! )\!)\nonumber\\ \displaystyle & & \displaystyle /(160p_{1}^{6}(k_{a}^{2}+k_{a}+p_{1}^{2})^{3}),\end{eqnarray}$$

$p_{1}$

$k_{a}>0$

$v_{1}$

$v_{0}$

$\mathit{Pe}$

$D_{1}$

$D_{0}$

$\mathit{Pe}$

$\mathit{Pe}$

3.4 Equilibrium sorption model

If the kinetics of the reactions are fast enough that local chemical equilibrium is valid, the linear kinetic sorption model reduces to the linear isotherm (i.e. equilibrium sorption model) $\unicode[STIX]{x1D6FE}=kc$ , with $k=k_{a}/k_{d}$ . For the equilibrium sorption model, (3.15) becomes

(3.25) $$\begin{eqnarray}\tan (p)=-kp,\end{eqnarray}$$

which can be solved for a series of $p_{k}$ . Taking the limit $k_{a}\rightarrow \infty ,~k_{d}\rightarrow \infty$ of (3.21) while keeping $k_{a}/k_{d}=k$ , the coefficients become functions of only the partition coefficient $k$ , as expected.

3.5 First-order approximation of the series solution

For the general case when $k_{d}$ is not zero, the fast transport described by (3.24) may survive at early times before desorption has come into play. In this case, the general series solution of the moments (3.19) allows us to study the transition from fast transport at early times described by first-order terms to slow transport at late times described by zeroth-order terms.

The zeroth- and first-order terms correspond to the residues $R_{0}$ and $R_{1}$ in (3.17). Figure 2(a,b) shows the comparison of the zeroth-order approximation $R_{0}$ and first-order approximation $R_{0}+R_{1}$ with the numerical inversion of the Laplace transform using Talbot’s method (Abate & Whitt Reference Abate and Whitt2006; McClure Reference McClure2013). As expected, the first-order approximation $R_{0}+R_{1}$ captures the solution at both the early and late times, while the zeroth-order approximation $R_{0}$ only describes the late-time behaviour. Additional tests show that the first-order approximation is sufficient to describe the solution for a large range of $k_{a}$ and $k_{d}$ after a short initial time. Therefore, we truncate the series solution (3.19) by retaining only the zeroth- and first-order terms,

(3.26a ) $$\begin{eqnarray}\displaystyle & \displaystyle m_{0}=a_{0}+a_{1}\exp (-p_{1}^{2}t), & \displaystyle\end{eqnarray}$$

(3.26b ) $$\begin{eqnarray}\displaystyle & \displaystyle m_{1}=(b_{0}^{(1)}+b_{0}^{(2)}t)+(b_{1}^{(1)}+b_{1}^{(2)}t)\exp (-p_{1}^{2}t), & \displaystyle\end{eqnarray}$$

(3.26c ) $$\begin{eqnarray}\displaystyle & \displaystyle m_{2}=(c_{0}^{(1)}+c_{0}^{(2)}t+c_{0}^{(3)}t^{2})+(c_{1}^{(1)}+c_{1}^{(2)}t+c_{1}^{(3)}t^{2})\exp (-p_{1}^{2}t), & \displaystyle\end{eqnarray}$$

Figure 2. Comparison of the zeroth-order approximation ( $R_{0}$ ) and the first-order approximation ( $R_{0}+R_{1}$ ) with the numerical inversion of the Laplace transform (squares) using Talbot’s method. $R_{0}$ and $R_{1}$ are the zeroth and first residues of the moments defined in (3.17) and (3.18). Results are shown for $\mathit{Pe}=10$ , $k_{a}=10$ and $k_{d}=1$ . (a) The transport velocity $v$ , where the mean flow velocity $u_{0}=U_{0}H/D=\mathit{Pe}$ . (b) The dispersion coefficient $D_{L}$ , where the dashed line labelled as Taylor denotes the Taylor dispersion $2/105\mathit{Pe}^{2}$ .

4.1 Two-dimensional simulations

Figure 3 shows two-dimensional simulations of the solute concentration at different times for $\mathit{Pe}=10$ , $k_{a}=50$ and $k_{d}=1$ . The full problem is numerically solved by the lattice Boltzmann method (LBM) (Chen & Doolen Reference Chen and Doolen1998; Wang & Kang Reference Wang and Kang2010; Zhang & Wang Reference Zhang and Wang2015). The $\unicode[STIX]{x1D6FF}$ -function initial condition is approximated by a piecewise-constant function that is non-zero in a small interval around the origin. This approximation of the initial condition only affects the results in a short diffusive transient and the results agree well with the analytical solution (figure 3 f–h).

Initially, the strong adsorption removes the solute from the slow-moving fluid near the wall. The remaining solute in the centre of the channel forms a fast-moving pulse (figure 3 b,c), particularly evident in the transversely averaged concentration shown in figure 3(e). This corresponds to the increased solute transport velocity in the irreversible sorption case (Lungu & Moffatt Reference Lungu and Moffatt1982). This regime persists as long as adsorption dominates.

However, the fast-moving pulse decays rapidly and eventually desorption releases solute in its wake (figure 3 d). As the amount of desorbed solute in the slow-moving fluid near the wall increases, the solute transport velocity declines. This process continues until desorption at the back balances adsorption at the front. The transport velocity and dispersion coefficient will approach the slow transport described by the one-dimensional model of the transversely averaged concentration in the reversible sorption case (Khan Reference Khan and Van Swaay1962).

Figure 3. Two-dimensional simulation of solute transport with sorption in Poiseuille flow with $\mathit{Pe}=10$ , $k_{a}=50$ and $k_{d}=1$ . (a–d) Concentration distribution and (e) transversely averaged concentration profile at $t=0.12$ , 0.73, 2.7, 12. (f–h) The evolution of the zeroth-, first- and second-order moments and comparison of the numerical simulation (LBM) with the first-order approximation of the analytical solution (Ana).

4.2 Equilibrium sorption model

Following the solution procedure in § 3.4, this section presents results and analysis for equilibrium sorption model, $\unicode[STIX]{x1D6FE}=kc$ . To demonstrate the different transport behaviours, we define a normalized transport velocity ${\mathcal{V}}$ and a normalized dispersion coefficient ${\mathcal{D}}$ as

(4.1a,b ) $$\begin{eqnarray}{\mathcal{V}}=\frac{v}{Pe}\quad \text{and}\quad {\mathcal{D}}=\frac{D_{L}-1}{D_{t}},\end{eqnarray}$$

where $D_{t}=2/105\mathit{Pe}^{2}$ is the Taylor dispersion coefficient for a tracer in Poiseuille flow and the unit contribution of diffusion has been subtracted in the numerator of (4.1b ). In this way, ${\mathcal{V}}>1$ ( ${\mathcal{D}}>1$ ) means increased velocity (dispersion) relative to a non-reactive tracer while ${\mathcal{V}}<1$ ( ${\mathcal{D}}<1$ ) means decreased velocity (dispersion).

Figure 4 shows the evolution of the position of the centre of mass $M_{1}$ and the normalized transport velocity ${\mathcal{V}}$ for different partition coefficients. For $k>10$ , a linear region emerges at early times in figure 4(a), corresponding to an initial plateau in figure 4(b). This corresponds to the well-developed early regime characterized by fast transport, approaching an asymptotic velocity $1+3/\unicode[STIX]{x03C0}^{2}\approx 1.3$ . This is consistent with the results in an adsorption-only case with $k_{a}\rightarrow \infty$ (Lungu & Moffatt Reference Lungu and Moffatt1982). After a transition period, a second linear region at late times appears, corresponding to the decreased transport velocity $1/(1+k)$ .

Figure 4. Evolution of (a) centre of mass $M_{1}$ and (b) normalized transport velocity ${\mathcal{V}}$ for different partition coefficients  $k$ . Dashed lines with labels LM and Taylor stand for the asymptotic regime of an adsorption-only case (Lungu & Moffatt Reference Lungu and Moffatt1982) and the asymptotic regime of a non-reactive tracer, respectively.

Figure 5 shows similar behaviours of the variance of solute mass in the fluid $M_{2}$ and the normalized dispersion coefficient ${\mathcal{D}}$ . In the early regime, the dispersion coefficient is reduced relative to a tracer with ${\mathcal{D}}\sim 0.14$ , which also agrees with the adsorption-only case with $k_{a}\rightarrow \infty$ . In the late regime, the dispersion coefficient is given by the first two terms of (3.22b ), first obtained by Golay (Reference Golay and Desty1958).

Figure 5. Evolution of (a) variance of solute mass distribution $M_{2}$ and (b) the normalized dispersion coefficient ${\mathcal{D}}$ for different partition coefficients  $k$ . Dashed lines with labels LM and Taylor stand for the asymptotic regime of an adsorption-only case (Lungu & Moffatt Reference Lungu and Moffatt1982) and the asymptotic regime of a non-reactive tracer, respectively.

Between the early and late times, there is a drastic transition of the transport behaviour. Especially when $k$ is large, both $M_{1}$ and $M_{2}$ decrease after reaching maxima during the transition, and this leads to negative velocity and dispersion coefficient. Physically, it means that desorption near the origin dominates over the fast-moving pulse in figure 3 so that the centre of mass shifts backwards and the variance reduces because the transversely averaged concentration distribution changes from a bimodal type (one peak near the origin and the other at the pulse front) to a unimodal type (single peak near the origin).

To compare the early and late behaviours as a function of $k$ , we define the normalized early and late velocities as

(4.2a,b ) $$\begin{eqnarray}{\mathcal{V}}_{e}=\frac{v_{1}}{\mathit{Pe}}\quad \text{and}\quad {\mathcal{V}}_{l}=\frac{v_{0}}{\mathit{Pe}},\end{eqnarray}$$

and similarly we define the normalized early and late dispersion coefficients as

(4.2c,d ) $$\begin{eqnarray}{\mathcal{D}}_{e}=\frac{D_{1}-1}{D_{t}}\quad \text{and}\quad {\mathcal{D}}_{l}=\frac{D_{0}-1/(1+k)}{D_{t}},\end{eqnarray}$$

where $v_{0}$ , $D_{0}$ and $v_{1}$ , $D_{1}$ are obtained from the zeroth- and first-order terms of the solution, the equilibrium limits of (3.22) and (3.23). Note that at early times the effective diffusion is not affected by sorption while it is reduced by a factor of $1/(1+k)$ at late times.

As shown in figure 6, for large $k$ the difference of normalized transport velocity between ${\mathcal{V}}_{e}$ and ${\mathcal{V}}_{l}$ is the largest and the normalized early-time dispersion coefficient ${\mathcal{D}}_{e}$ asymptotes to 0.14. The normalized late-time dispersion coefficient ${\mathcal{D}}_{l}$ first increases with $k$ , then reduces towards 0. For small $k$ , the early velocity ${\mathcal{V}}_{e}$ and dispersion coefficient ${\mathcal{D}}_{e}$ do not reach the asymptotic values 1.3 and 0.14. In this case, the early regime is not well developed and the first-order terms of the solution are not dominant. Therefore, ${\mathcal{V}}_{e}$ and ${\mathcal{D}}_{e}$ do not represent the transport behaviour in this case.

Figure 6. (a) Normalized early velocity ${\mathcal{V}}_{e}=v_{1}/\mathit{Pe}$ and late velocity ${\mathcal{V}}_{l}=v_{0}/\mathit{Pe}$ and (b) normalized early dispersion coefficient ${\mathcal{D}}_{e}=(D_{1}-1)/D_{t}$ and late dispersion coefficient ${\mathcal{D}}_{l}=(D_{0}-1/(1+k))/D_{t}$ as functions of partition coefficient  $k$ .

For a tracer in Poiseuille flow, the pre-asymptotic transport before equilibrium has been studied extensively (e.g. Gill & Sankarasubramanian Reference Gill and Sankarasubramanian1970; Haber & Mauri Reference Haber and Mauri1988; Mercer & Roberts Reference Mercer and Roberts1990; Latini & Bernoff Reference Latini and Bernoff2001; Dentz & Carrera Reference Dentz and Carrera2007; Bolster et al. Reference Bolster, Valdés-Parada, Leborgne, Dentz and Carrera2011; Wang et al. Reference Wang, Cardenas, Deng and Bennett2012). Typically, diffusion dominates when $t\ll t_{d}=\mathit{Pe}^{-2/3}$ ; and after the characteristic equilibrium time scale, $t=1$ , solute transport can be described by the transversely averaged model with the mean flow velocity and the dispersion coefficient $D_{t}$ . However, for a reactive case considered here, the time scale to reach equilibrium can be quite different from the tracer case because surface reactions introduce additional characteristic time scales.

In the first-order approximation (3.26), a series of time scales can be defined by comparing the zeroth-order and first-order terms. Taking $m_{0}$ as an example, by comparing $a_{0}$ and $a_{1}\exp (-p_{1}^{2}t)$ , we can define a time scale as

(4.3) $$\begin{eqnarray}t_{1}=\frac{1}{p_{1}^{2}}\ln \frac{a_{1}}{a_{0}}=\frac{1}{p_{1}^{2}}\ln \frac{2k^{2}(k+1)}{k^{2}p_{1}^{2}+k+1}.\end{eqnarray}$$

This time scale indicates the transition from the early regime when the transport is dominated by the first-order terms to the late regime dominated by the zeroth-order terms. Similarly, other time scales can be determined by comparing $b_{k}$ and $c_{k}$ . We notice that the coefficients in the series solution have the property

(4.4) $$\begin{eqnarray}\frac{a_{1}}{a_{0}}<\frac{b_{1}^{(2)}}{b_{0}^{(2)}}\sim \frac{c_{1}^{(2)}}{c_{0}^{(2)}}<\frac{c_{1}^{(3)}}{c_{0}^{(3)}},\end{eqnarray}$$

which means that $t_{1}$ will give the smallest time scale. At the same time, since the time scales for the velocity determined by $b_{1}^{(2)}/b_{0}^{(2)}$ and the dispersion coefficient determined by $c_{1}^{(2)}/c_{0}^{(2)}$ have the same scaling, we can choose

(4.5) $$\begin{eqnarray}t_{2}=\frac{1}{p_{1}^{2}}\ln \frac{b_{1}^{(2)}}{b_{0}^{(2)}}=\frac{1}{p_{1}^{2}}\ln \frac{k^{2}(k+1)^{2}(4k^{2}p_{1}^{4}+3k^{2}p_{1}^{2}-3k+4p_{1}^{2}-3)}{2p_{1}^{2}(k^{2}p_{1}^{2}+k+1)^{2}}\end{eqnarray}$$

as a critical time scale, after which the zeroth-order terms dominate and the late-time behaviour emerges. Note that the two time scales are not dependent on  $\mathit{Pe}$ . Consequently, the transient solute transport with sorption can be divided into the following three regimes:

1. (I) $0<t<t_{1}$ – early regime with fast transport;

2. (II) $t_{1}<t<t_{2}$ – transition period;

3. (III) $t_{2}<t$ – late regime with slow transport.

As shown in figures 4 and 5, the duration of the early regime, as well as the transition period, increase with increasing $k$ . In the limit of large $k$ , we have

(4.6a,b ) $$\begin{eqnarray}\lim _{k\rightarrow \infty }t_{1}=\frac{4}{\unicode[STIX]{x03C0}^{2}}\ln k\quad \text{and}\quad \lim _{k\rightarrow \infty }t_{2}=\frac{8}{\unicode[STIX]{x03C0}^{2}}\ln k,\end{eqnarray}$$

where we have used $\lim _{k\rightarrow \infty }p_{1}=\unicode[STIX]{x03C0}/2$ . Equations (4.6) predict a linear relationship between $t_{1}$ , $t_{2}$ and $\ln k$ when $k$ is large and $t_{2}\sim 2\,t_{1}$ . Figure 7 compares results from the numerical inverse Laplace transform with these analytically determined time scales as a function of  $k$ . When $k\gg 1$ , $t$ scales with $\ln k$ , as predicted by (4.6). For small $k$ , the early regime is so short that it is generally not observed.

Figure 7. Variation of the early time scale $t_{1}$ and the late time scale $t_{2}$ as functions of partition coefficient $k$ . The solid lines show $t_{1}$ and $t_{2}$ from (4.3) and (4.5), respectively, and the dashed lines indicate the scalings for large $k$ given by (4.6), which have been shifted to match the symbols determined by the numerical inverse Laplace transform. White symbols are determined by velocity and grey symbols are determined by dispersion coefficient. The transient solute transport with sorption is divided into three regimes: (I) early regime with fast transport; (II) transition period; and (III) late regime with slow transport.

4.3 Kinetic sorption model

In the kinetic sorption model, there are two additional governing parameters, namely, the dimensionless adsorption rate constant $k_{a}$ and the dimensionless desorption rate constant $k_{d}$ . Generally, a similar transition from early to late behaviour can be observed and the equilibrium results are recovered when kinetics are fast, i.e. $k_{a}\gg 1$ and $k_{d}\gg 1$ .

Similar to the way that the time scales are determined for the equilibrium model, we can obtain $t_{1}$ and $t_{2}$ for the kinetic model using (4.3) and (4.5),

(4.7a ) $$\begin{eqnarray}\displaystyle & \displaystyle t_{1}=\frac{1}{p_{1}^{2}}\ln \frac{2k_{a}^{2}(k_{a}+k_{d})}{k_{d}(p_{1}^{4}+(k_{a}^{2}+k_{a}-2k_{d})p_{1}^{2}+k_{a}k_{d}+k_{d}^{2})}, & \displaystyle\end{eqnarray}$$

(4.7b ) $$\begin{eqnarray}\displaystyle & \displaystyle t_{2}=\frac{1}{p_{1}^{2}}\ln \frac{k_{a}^{2}(k_{a}+k_{d})^{2}(4p_{1}^{6}+(4k_{a}^{2}-8k_{d}-3)p_{1}^{4}+(3k_{a}^{2}+3k_{a}+4k_{d}^{2}+6k_{d})p_{1}^{2}-3k_{a}k_{d}-3k_{d}^{2})}{2k_{d}^{2}p_{1}^{2}(p_{1}^{4}+(k_{a}^{2}+k_{a}-2k_{d})p_{1}^{2}+k_{a}k_{d}+k_{d}^{2})^{2}},\hspace{-6.0pt} & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The early regime is only observed when $k_{a}$ exceeds $k_{d}$ . In all other cases, transition occurs from the very beginning followed by a dominated late regime. When both $k_{a}$ and $k_{d}$ are large, the time scales of the kinetic model recover those of the equilibrium model.

However, if the rates decrease, the kinetic time scales become longer. In this case, the root $p_{1}$ of (3.15) can be approximated by $p_{1}^{2}\approx k_{a}+k_{d}$ using Taylor expansion for $\tan (p)$ . Then the ratios $a_{1}/a_{0}$ and $b_{1}^{(2)}/b_{0}^{(2)}$ used to obtain the time scales simplify to

(4.8a,b ) $$\begin{eqnarray}\frac{a_{1}}{a_{0}}=\frac{2k_{a}}{k_{d}(k_{a}+2)}\quad \text{and}\quad \frac{b_{1}^{(2)}}{b_{0}^{(2)}}=\frac{k_{a}^{2}(4k_{a}+4k_{d}+7)}{2k_{d}^{2}(k_{a}+2)^{2}}.\end{eqnarray}$$

This analysis shows that both time scales increase dramatically in the lower left region where $k_{a}$ and $k_{d}$ are small in figure 8. In this region, the duration of the early regime is long, but the deviations of velocity and dispersion coefficient from the tracer case are minor, as the limiting values given by the adsorption-only case approach unity with small $k_{a}$ . Physically, this region corresponds to a kinetically slow-sorbing ( $k_{a},k_{d}\ll 1$ ) solute with a large partition coefficient $k\gg 1$ .

The analytical solution presented in § 3.3 recovers the previous analysis in the limit of $k_{d}=0$ (Lungu & Moffatt Reference Lungu and Moffatt1982). This limiting solution puts an upper bound on the transport velocity and a lower bound on the dispersion coefficient in the early regime. If the early regime is well developed, the limiting solution given by Lungu & Moffatt (Reference Lungu and Moffatt1982) provides a good approximation for finite $k_{d}$ (see figure 8 c). The well-developed early regime is indicated by grey shadings in figure 8(a), where the early-time asymptotic transport velocity and dispersion coefficient given by (3.23) are within 10 % of the limiting values given by Lungu & Moffatt (Reference Lungu and Moffatt1982). Cases A, B and C give examples of the well-developed early regime, while the velocity does not reach the asymptotic value in case D. Generally, if $k>10$ ( $k>1000$ ), the early-time velocity (dispersion coefficient) is well developed. Note that the first-order analytical solution for transport velocity ( $R_{0}+R_{1}$ ) shown in figure 8(c) is computed by (3.2) and (3.26).

Figure 8. Contour lines of (a) the early time scale $t_{1}$ and (b) the late time scale $t_{2}$ for the kinetic sorption model. Dashed lines represent the time scales obtained from the equilibrium sorption model for large $k_{a}$ and $k_{d}$ and dotted lines represent approximations for small $k_{a}$ and $k_{d}$ , given by (4.8a,b ). The shaded area in (a) represents the region where the early regime is well developed, i.e. the velocity (light shaded area) and dispersion coefficient (dark shaded area) are close to the limiting values given by the adsorption-only case. (c) The evolution of velocity at early time for the conditions labelled as A, B, C and D in (a). Symbols are the results from full numerical inverse Laplace transform, solid lines are first-order approximation of the analytical solution and dashed line is the asymptotic value for $k_{d}=0$ at $k_{a}=1$ , given by Lungu & Moffatt (Reference Lungu and Moffatt1982).