1. Introduction
The hydrodynamic dispersion phenomenon has significant importance in various areas in science and engineering processes such as homogenization (Matsunaga, Lee & Nishino Reference Matsunaga, Lee and Nishino2013; Cheng et al. Reference Cheng, Feng, Cheng, Yang and Mao2015; Sadeghi Reference Sadeghi2016), separation (Garcia et al. Reference Garcia, Ista, Petsev, O'Brien, Bisong, Mammoli and López2005) and DNA amplification (Tripathi, Bozkurt & Chauhan Reference Tripathi, Bozkurt and Chauhan2005). It has also attracted the attention of researchers working on drug delivery systems (Siepmann & Siepmann Reference Siepmann and Siepmann2013; Thomas et al. Reference Thomas, Nair, Latha and Thomas2019) since the dependence of a drug release profile on the thickness of the hydrodynamic diffusion layer has been proven (Patel et al. Reference Patel, Panchal, Patel, Brahmbhatt and Suthar2011). Furthermore, closely related to the present topic, as the scattering of samples being transported by microflows must be minimized for efficient functionality of microchannels, the efficacy of microfluidic devices strongly depends on the control of the hydrodynamic dispersion (Ajdari, Bontoux & Stone Reference Ajdari, Bontoux and Stone2006).
Pioneering research works on hydrodynamic dispersion were done by Taylor (Reference Taylor1953) and Aris (Reference Aris1956) assuming a steady and laminar flow in a circular channel, who introduced the effective dispersion coefficient and presented a generalized solution including axial diffusion, respectively. Their seminal works were followed by that of Gill & Sankarasubramanian (Reference Gill and Sankarasubramanian1970) who were able to obtain a time-dependent dispersion coefficient and an exact solution for the local concentration via a generalized dispersion model; a steady and laminar flow in a circular geometry was considered by those authors while adopting a uniform solute injection. They soon extended their work to situations that involve non-uniform initial concentration distribution (Gill, Sankarasubramanian & Taylor Reference Gill, Sankarasubramanian and Taylor1971), time-dependent flow (Gill & Sankarasubramanian Reference Gill and Sankarasubramanian1972; Sankarasubramanian & Gill Reference Sankarasubramanian and Gill1972) and interphase reaction (Sankarasubramanian & Gill Reference Sankarasubramanian and Gill1973).
The works listed above constitute the most important progress made in the theoretical modelling of hydrodynamic dispersion. There are, however, more recent studies that have significantly contributed to the broadening of our knowledge of this phenomenon. By numerically simulating the dispersion of a pulse of solute in a tube with irreversible, heterogeneous reaction or reversible adsorption at the wall, Paine, Carbonell & Whitaker (Reference Paine, Carbonell and Whitaker1983) established the foundations required for further study of mass transport in porous media. In an attempt to extend the classical dispersion theory to multiphase flows, Phillips & Kaye (Reference Phillips and Kaye1998) developed asymptotic approximations for the transport of solutes in Poiseuille flow through a pipe having a permeable wall layer. More recently, Ajdari et al. (Reference Ajdari, Bontoux and Stone2006) theoretically showed that the hydrodynamic dispersion in shallow microchannels is controlled mostly by the width of the cross-section rather than by the channel height. Taghizadeh, Valdés-Parada & Wood (Reference Taghizadeh, Valdés-Parada and Wood2020) developed an upscaled formulation for the description of solute spreading at early times that was shown to be capable of truly modelling cases with multiple solute injections. Furthermore, Jiang and Chen extended the concept of hydrodynamic dispersion to the situation where particles are active, a condition that occurs in applications involving micro-organisms (Jiang & Chen Reference Jiang and Chen2019, Reference Jiang and Chen2020). Concurrently, Dehe, Rehm & Hardt (Reference Dehe, Rehm and Hardt2021) derived a two-dimensional dispersion model for inhomogeneous flow fields inside a Hele-Shaw cell. There is also a series of recently published papers on the hydrodynamic dispersion caused by new flow generation mechanisms that are particularly useful at the microscale such as electroosmosis (Dutta Reference Dutta2007; Paul & Ng Reference Paul and Ng2012; Arcos et al. Reference Arcos, Méndez, Bautista and Bautista2018; Hoshyargar et al. Reference Hoshyargar, Khorami, Ashrafizadeh and Sadeghi2018; Muñoz et al. Reference Muñoz, Arcos, Bautista and Méndez2018; Sadeghi et al. Reference Sadeghi, Saidi, Moosavi and Sadeghi2020; Talebi, Ashrafizadeh & Sadeghi Reference Talebi, Ashrafizadeh and Sadeghi2021), diffusioosmosis (Hoshyargar, Ashrafizadeh & Sadeghi Reference Hoshyargar, Ashrafizadeh and Sadeghi2017) and capillarity (Fridjonsson, Seymour & Codd Reference Fridjonsson, Seymour and Codd2014).
Although hundreds of papers have been published during the last seven decades, from the time Taylor published his groundbreaking work (Taylor Reference Taylor1953), the semicircular geometry has received no attention. A reasonable reason for such neglect might be the fact that this geometry does not have significant importance at the macroscale. However, the situation is totally different when it comes to the microscale. Before getting to the point, it is worth noting that there are different microfabrication techniques such as thin-film deposition, lithography and etching, among others, and nearly all of them are only capable of fabricating channels of specific cross-sectional geometries (Ziaie, Baldi & Atashbar Reference Ziaie, Baldi and Atashbar2010). In isotropic etching, since a material is under attack in all directions, a geometry is created that is very similar to a semicircle (Ziaie et al. Reference Ziaie, Baldi, Lei, Gu and Siegel2004; Sadeghi, Sadeghi & Saidi Reference Sadeghi, Sadeghi and Saidi2016). Accordingly, semicircular microchannels are used not only by the scientific community (Hisamoto et al. Reference Hisamoto, Shimizu, Uchiyama, Tokeshi, Kikutani, Hibara and Kitamori2003; Asadi-Saghandi et al. Reference Asadi-Saghandi, Karimi-Sabet, Ghorbanian and Moosavian2022) but also in commercial products such as the Dolomite Y-junction chips designed for the observation of the dispersion process (www.dolomite-microfluidics.com). The datasheet of the Dolomite product, sketched schematically in figure 1(a) by exaggerating the differences it has from a semicircular geometry, indicates a channel radius-to-depth ratio of 
$1.025$, which is very close to unity. Consequently, approximating such a geometry by a semicircular one, for better mathematical tractability, looks quite reasonable. This, together with the fact that geometry plays a very important role in hydrodynamic dispersion (Lee et al. Reference Lee, Luner, Marzuola and Harris2021), persuaded the authors to focus on the problem of hydrodynamic dispersion within semicircular microchannels in the present study. To this end, the generalized dispersion model (Sankarasubramanian & Gill Reference Sankarasubramanian and Gill1973) is invoked to track an injected solute band of arbitrary axial distribution and general shape in the cross-sectional area from the time of injection considering a steady and fully developed flow. It is assumed that an irreversible first-order reaction occurs at the curved channel wall while considering a no-flux boundary condition at the flat wall. Analytical solutions are obtained for the exchange, convection and dispersion coefficients while presenting formulas for calculating higher-order coefficients as well as solutions for the local and cross-sectionally averaged concentrations. Special solutions are also presented for uniform velocity and no-reaction cases. It should be pointed out that although the present modelling is motivated by the importance of the semicircular geometry in microfluidics, the results are general enough to be used at both the microscale and macroscale.
Figure 1. (a) Comparison between the cross-section of the Dolomite microfluidic channel under consideration (www.dolomite-microfluidics.com) and the semicircular geometry adopted here as well as the cross-sectional view of the CID that is of radius 
${R_D}$ and is centred at 
$R = {C_D}$ where 
$\phi = {\rm \pi}/2$; note that the difference between the actual and modelled geometries is exaggerated here to facilitate making a distinction between the two. (b) A three-dimensional sketch of the microfluidic channel wherein the coordinate system, the SID with radius 
$a{R_{SC}}$ and the broadening of a typical semicircular solute band over time are shown.
2. Problem formulation
2.1. Problem definition
The dispersion of a solute band injected into a straight semicircular microchannel with radius 
${R_{sc}}$ is theoretically studied. Because of its practical importance in modelling of short-term injections, it is first assumed that the extent of the solute band in the axial direction is vanishingly small and then the solutions are extended to arbitrary axial distributions. The band is supposed to be transported downstream by the flow of a Newtonian fluid with constant physical properties. The flow is considered to be steady, laminar and fully developed. Moreover, it is assumed that an irreversible first-order reaction occurs at the curved wall of the channel. As the flat wall is only a cover to the etched curved wall to create a closed fluidic passage, it is usually not used for reaction purposes; hence, the no-flux boundary condition is assumed to be valid at this wall. It is also assumed that the channel radius is at least two orders of magnitude larger than the Debye length, which is usually between 1 and 
$100\;\textrm{nm}$ (Karniadakis et al. Reference Karniadakis, Beskok and Aluru2005), to allow for the neglect of the surface effects. Two types of initial distributions are considered for the solute band in the cross-sectional area, namely circular and semicircular distributions, with the uniform distribution being a special case of the latter. The schematics of the circular initial distribution (CID) with radius 
${R_D}$ centred at 
$R = {C_D}$ where 
$\phi = {\rm \pi}/2$ and the semicircular initial distribution (SID) with radius 
$a{R_{sc}}$ are given in figures 1(a) and 1(b), respectively. Figure 1(b) also contains the details of the coordinate system along with three different concentration distributions that show the broadening of a typical semicircular solute band over time.
2.2. Velocity distribution
Considering a steady and fully developed flow, the basic law of momentum conservation is mathematically expressed as (Azari, Sadeghi & Chakraborty Reference Azari, Sadeghi and Chakraborty2020)
\begin{equation}\frac{1}{R}\frac{\partial }{{\partial R}}\left( {R\frac{{\partial U}}{{\partial R}}} \right) + \frac{1}{{{R^2}}}\frac{{{\partial ^2}U}}{{\partial {\phi ^2}}} = \frac{1}{\mu }\frac{{\textrm{d}P}}{{\textrm{d}X}}.\end{equation}
Here, U is the axial velocity, P is the pressure and 
$\mu $ stands for the dynamic viscosity. In this study, the no-slip boundary condition is considered at the walls. Therefore, the dimensionless form of (2.1) and the pertinent boundary conditions may be written as
\begin{gather}\frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial u}}{{\partial r}}} \right) + \frac{1}{{{r^2}}}\frac{{{\partial ^2}u}}{{\partial {\phi ^2}}} ={-} 4,\end{gather}
\begin{gather}u{|_{r = 1}} = u{|_{\phi = 0}} = u{|_{\phi = {\rm \pi}}} = 0,\quad u{|_{r = 0}} = \textrm{finite},\end{gather}
where 
$r = R/{R_{sc}}$ and 
$u = U/{U_R}$, with 
${U_R} ={-} R_{sc}^2(\textrm{d}P/\textrm{d}X)/4\mu $ representing the reference velocity, which is considered to be the same as the maximum fluid velocity in a circular channel of radius 
${R_{sc}}$. The solution of (2.2a) subject to the boundary conditions (2.2b) has been shown to be (Alassar & Abushoshah Reference Alassar and Abushoshah2012; Alassar Reference Alassar2014)
\begin{equation}u ={-} 2{r^2}{\sin ^2}\phi + \sum\limits_{w = 0}^\infty {{D_w}{r^{2w + 1}}} \sin [(2w + 1)\phi ],\quad {D_w} = \frac{{ - 16}}{{{\rm \pi} (8{w^3} + 12{w^2} - 2w - 3)}}.\end{equation}Considering the velocity distribution (2.3), it is straightforward to show that the dimensionless mean velocity over the channel cross-sectional area is
\begin{equation}{u_m} = \frac{{\mathop \smallint \nolimits_0^{\rm \pi} \mathop \smallint \nolimits_0^1 ur\,\textrm{d}r\,\textrm{d}\phi }}{{\mathop \smallint \nolimits_0^{\rm \pi} \mathop \smallint \nolimits_0^1 r\,\textrm{d}r\,\textrm{d}\phi }} ={-} \frac{1}{2} + \sum\limits_{w = 0}^\infty {\frac{{4{D_w}}}{{{\rm \pi} (4{w^2} + 8w + 3)}}} .\end{equation}2.3. Concentration distribution
The equation governing the solute concentration field can be considered as a time-dependent form of the convection–diffusion equation that takes the following form for the present problem:
\begin{equation}\frac{{\partial \varTheta }}{{\partial T}} + U\frac{{\partial \varTheta }}{{\partial X}} = D\left[ {\frac{{{\partial^2}\varTheta }}{{\partial {X^2}}} + \frac{1}{R}\frac{\partial }{{\partial R}}\left( {R\frac{{\partial \varTheta }}{{\partial R}}} \right) + \frac{1}{{{R^2}}}\frac{{{\partial^2}\varTheta }}{{\partial {\phi^2}}}} \right],\end{equation}
where 
$\varTheta $ is the number concentration of the solutes, T is the time and D represents the molecular diffusivity. Recalling the assumptions made in the problem formulation section, the solute concentration equation is constrained by the following initial and boundary conditions:
\begin{gather}\varTheta {|_{T = 0}} = {\varTheta _0}{R_{sc}}{\delta _1}(X){\psi _1}(R,\phi ),\end{gather}
\begin{gather}\varTheta {|_{X \to \infty }} = {\left. {\frac{{\partial \varTheta }}{{\partial X}}} \right|_{X \to \infty }} = {\left. {\frac{{\partial \varTheta }}{{\partial \phi }}} \right|_{\phi = 0}} = {\left. {\frac{{\partial \varTheta }}{{\partial \phi }}} \right|_{\phi = {\rm \pi}}} = 0,\end{gather}
\begin{gather}\varTheta {|_{R = 0}} = \textrm{finite},\end{gather}
\begin{gather}{\left. {\frac{{\partial \varTheta }}{{\partial R}}} \right|_{R = {R_{sc}}}} ={-} \frac{{{k_R}}}{D}\varTheta {|_{R = {R_{sc}}}}.\end{gather}
Equation (2.6a) is a manifestation of the facts that at time zero the solute band concentrates longitudinally near the channel inlet, represented mathematically via the Dirac delta function 
${\delta _1}$, and possesses a prescribed distribution over the channel cross-sectional area, represented mathematically via the dimensionless function 
${\psi _1}$. The parameter 
${\varTheta _0}$ in (2.6a) is an arbitrary reference concentration. In addition, (2.6b) denotes that the concentration and its axial derivative are both zero sufficiently downstream and there is no solute flux at the flat wall of the channel. Finally, (2.6c) and (2.6d) reflect the facts that the solute concentration is finite at the origin of the coordinate system and there is a first-order reaction taking place at the curved wall with rate constant 
${k_R}$, respectively. Considering the dimensionless parameters 
$\theta = \varTheta /{\varTheta _0}$, 
$t = TD/R_{sc}^2$, 
$Pe = {R_{sc}}{U_R}/D$ and 
$x = X/{R_{sc}}Pe$, the scaled forms of the convection–diffusion equation and the associated initial and boundary conditions may be written as
\begin{gather}\frac{{\partial \theta }}{{\partial t}} + u\frac{{\partial \theta }}{{\partial x}} = \frac{1}{{P{e^2}}}\frac{{{\partial ^2}\theta }}{{\partial {x^2}}} + \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial \theta }}{{\partial r}}} \right) + \frac{1}{{{r^2}}}\frac{{{\partial ^2}\theta }}{{\partial {\phi ^2}}},\end{gather}
\begin{gather}\theta {|_{t = 0}} = \frac{{\delta (x)}}{{Pe}}\psi (r,\phi ),\end{gather}
\begin{gather}\theta {|_{x \to \infty }} = {\left. {\frac{{\partial \theta }}{{\partial x}}} \right|_{x \to \infty }} = {\left. {\frac{{\partial \theta }}{{\partial \phi }}} \right|_{\phi = 0}} = {\left. {\frac{{\partial \theta }}{{\partial \phi }}} \right|_{\phi = {\rm \pi}}} = 0, \end{gather}
\begin{gather}\theta {|_{r = 0}} = \textrm{finite},\end{gather}
\begin{gather}{\left. {\frac{{\partial \theta }}{{\partial r}}} \right|_{r = 1}} ={-} Da\,\theta {|_{r = 1}},\end{gather}
in which 
$\delta (x) \equiv {R_{sc}}Pe{\delta _1}(X)$, 
$\psi (r,\phi ) \equiv {\psi _1}(R,\phi )$ and 
$Da = {R_{sc}}{k_R}/D$ is the Damköhler number, showing the ratio of the diffusion time scale 
$R_{sc}^2/D$ to the reaction time scale 
${R_{sc}}/{k_R}$. The procedure chosen for solving (2.7) subject to the initial and boundary conditions (2.8) is the powerful method proposed by Sankarasubramanian & Gill (Reference Sankarasubramanian and Gill1973) that starts with the consideration of following infinite series for the concentration distribution:
\begin{equation}\theta (t,r,\phi ,x) = \sum\limits_{k = 0}^\infty {{f_k}(t,r,\phi )} \frac{{{\partial ^k}{\theta _m}(t,x)}}{{\partial {x^k}}},\end{equation}
where 
${\theta _m}$ is the mean dimensionless concentration over the channel cross-section, given mathematically as
\begin{equation}{\theta _m}(t,x) = \frac{{\mathop \smallint \nolimits_0^1 \mathop \smallint \nolimits_0^{\rm \pi} r\theta \,\textrm{d}\phi \,\textrm{d}r}}{{\mathop \smallint \nolimits_0^1 \mathop \smallint \nolimits_0^{\rm \pi} r\,\textrm{d}\phi \,\textrm{d}r}} = \frac{2}{{\rm \pi} }\int_0^1 {\int_0^{\rm \pi} {r\theta \,\textrm{d}\phi \,\textrm{d}r} } .\end{equation}Making use of the boundary conditions (2.8b) after integrating (2.7) over the channel cross-sectional area, the following equation is obtained:
\begin{equation}\frac{{\partial {\theta _m}}}{{\partial t}} + \frac{2}{{\rm \pi} }\frac{\partial }{{\partial x}}\left( {\int_0^1 {\int_0^{\rm \pi} {ru\theta \,\textrm{d}\phi \,\textrm{d}r} } } \right) = \frac{1}{{P{e^2}}}\frac{{{\partial ^2}{\theta _m}}}{{\partial {x^2}}} + \frac{2}{{\rm \pi} }\int_0^{\rm \pi} {{{\left. {\frac{{\partial \theta }}{{\partial r}}} \right|}_{r = 1}}\,} \textrm{d}\phi .\end{equation}
Substituting for 
$\theta $ in (2.11) from (2.9), we will have
\begin{equation}\frac{{\partial {\theta _m}}}{{\partial t}} + \frac{2}{{\rm \pi} }\sum\limits_{k = 0}^\infty {\left( {\frac{{{\partial^{k + 1}}{\theta_m}}}{{\partial {x^{k + 1}}}}\int_0^1 {\int_0^{\rm \pi} {ru{f_k}\,\textrm{d}\phi \,\textrm{d}r} } } \right)} = \frac{1}{{P{e^2}}}\frac{{{\partial ^2}{\theta _m}}}{{\partial {x^2}}} + \frac{2}{{\rm \pi} }\sum\limits_{k = 0}^\infty {\left( {\frac{{{\partial^k}{\theta_m}}}{{\partial {x^k}}}\mathop \smallint \limits_0^{\rm \pi} {{\left. {\frac{{\partial {f_k}}}{{\partial r}}} \right|}_{r = 1}}\textrm{d}\phi } \right)} .\end{equation}Equation (2.12) can be rearranged to yield
\begin{equation}\frac{{\partial {\theta _m}}}{{\partial t}} = \sum\limits_{k = 0}^\infty {{K_k}(t)} \frac{{{\partial ^k}{\theta _m}}}{{\partial {x^k}}},\end{equation}wherein
\begin{equation}{K_k} = \frac{{{\delta _{k2}}}}{{P{e^2}}} + \frac{2}{{\rm \pi} }\int_0^{\rm \pi} {{{\left. {\frac{{\partial {f_k}}}{{\partial r}}} \right|}_{r = 1}}\textrm{d}\phi } - \frac{2}{{\rm \pi} }\int_0^1 {\int_0^{\rm \pi} {ru{f_{k - 1}}\,\textrm{d}\phi \,\textrm{d}r} } ,\quad k = 0,1,2, \ldots ,\end{equation}
where 
${f_{ - 1}} = {f_{ - 2}} = \cdots = 0$ and Kronecker delta 
${\delta _{ij}}$ equals one when 
$i = j$ and is zero in other situations. Expanding the series in (2.13) and truncating the terms corresponding to 
$k = 3$ and larger, due to their negligible effect on the mean concentration, the well-known dispersion model of Sankarasubramanian & Gill (Reference Sankarasubramanian and Gill1973) is recovered:
\begin{equation}\frac{{\partial {\theta _m}}}{{\partial t}} = {K_0}(t){\theta _m} + {K_1}(t)\frac{{\partial {\theta _m}}}{{\partial x}} + {K_2}(t)\frac{{{\partial ^2}{\theta _m}}}{{\partial {x^2}}}.\end{equation}
Note that 
${K_0}(t)$ is termed the exchange coefficient since it characterizes the degree of mass exchange at the channel boundaries and 
${K_1}(t)$ is called the convection coefficient and is due to the convective transport of the solutes. In addition, 
${K_2}(t)$, known as the dispersion coefficient, is a measure of solute dispersion due to the combined action of molecular diffusion and non-uniformity of the velocity profile. This means that all the information regarding the overall interphase mass transport at the wall, streamwise transport of the solute band and the degree of band broadening can be obtained as soon as the transport coefficients are known. Hence, different aspects of the hydrodynamic dispersion can be well described by the transport coefficients 
${K_k}$. In fact, these coefficients constitute all the design parameters required for liquid-phase transportation with surface reaction in flow conduits.
The determination of the transport coefficients is crucial for solving (2.15); to this end, first the functions 
${f_k}$ have to be obtained. This is sought by the substitution of (2.9) into (2.7), which results in
\begin{equation}\mathop \sum \limits_{k = 0}^\infty \left\{ {{f_k}\frac{{{\partial^{k + 1}}{\theta_m}}}{{\partial t\partial {x^k}}} + u{f_k}\frac{{{\partial^{k + 1}}{\theta_m}}}{{\partial {x^{k + 1}}}} + \left[ {\frac{{\partial {f_k}}}{{\partial t}} - \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial {f_k}}}{{\partial r}}} \right) - \frac{1}{{{r^2}}}\frac{{{\partial^2}{f_k}}}{{\partial {\phi^2}}}} \right]\frac{{{\partial^k}{\theta_m}}}{{\partial {x^k}}} - \frac{{{f_k}}}{{P{e^2}}}\frac{{{\partial^{k + 2}}{\theta_m}}}{{\partial {x^{k + 2}}}}} \right\} = 0.\end{equation}Recalling (2.13), the first term in (2.16) can be written as
\begin{equation}{f_k}\frac{{{\partial ^{k + 1}}{\theta _m}}}{{\partial t\partial {x^k}}} = {f_k}\frac{{{\partial ^k}}}{{\partial {x^k}}}\left[ {\sum\limits_{i = 0}^\infty {{K_i}(t)} \frac{{{\partial^i}{\theta_m}}}{{\partial {x^i}}}} \right] = {f_k}\sum\limits_{i = 0}^\infty {{K_i}(t)} \frac{{{\partial ^{i + k}}{\theta _m}}}{{\partial {x^{i + k}}}}.\end{equation}
Combining equations (2.16) and (2.17) and equating the coefficients of 
${\partial ^k}{\theta _m}/\partial {x^k}$ to zero, the following set of equations can be generated:
\begin{equation}\frac{{\partial {f_k}}}{{\partial t}} = \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial {f_k}}}{{\partial r}}} \right) + \frac{1}{{{r^2}}}\frac{{{\partial ^2}{f_k}}}{{\partial {\phi ^2}}} + \frac{{{f_{k - 2}}}}{{P{e^2}}} - u{f_{k - 1}} - \sum\limits_{i = 0}^k {{K_i}{f_{k - i}}} \quad \textrm{for}\;k = 0,1, 2, \ldots \end{equation}
Substituting equation (2.9) into the initial and boundary conditions (2.8) with the consideration of the fact that 
${\theta _m}{|_{t = 0}} = \delta (x){\psi _m}/Pe$ with 
${\psi _m} = \frac{2}{{\rm \pi} }\int_0^1 {\int_0^{\rm \pi} {r\psi \,\textrm{d}\phi \,\textrm{d}r} } $, the initial and boundary conditions associated with the functions 
${f_k}$ are derived to be
\begin{gather}{f_k}(0,r,\phi ) = {\delta _{k0}}\frac{{\psi (r,\phi )}}{{{\psi _m}}},\end{gather}
\begin{gather}{\left. {\frac{{\partial {f_k}}}{{\partial \phi }}} \right|_{\phi = 0}} = {\left. {\frac{{\partial {f_k}}}{{\partial \phi }}} \right|_{\phi = {\rm \pi}}} = 0,\quad {f_k}(t,0,\phi ) = \textrm{finite},\quad {\left. {\frac{{\partial {f_k}}}{{\partial r}}} \right|_{r = 1}} ={-} Da\,{f_k}{|_{r = 1}}.\end{gather}
Moreover, integrating equation (2.9) over the channel cross-sectional area provides the following constraints on 
${f_k}$:
\begin{equation}\int_0^1 {\int_0^{\rm \pi} {r{f_k}\textrm{(}t,r,\phi \textrm{)}\,\textrm{d}\phi \,\textrm{d}r} } = \frac{{\rm \pi} }{2}{\delta _{k0}}.\end{equation}
The governing equation and the associated initial and boundary conditions for 
${f_k}$ are now complete and we are ready to tackle them. The next step in the analysis is, therefore, to obtain the functions 
${f_k}$, which are then used for evaluating the transport coefficients 
${K_i}$ via (2.14). Afterward, the transport equations can be used to solve (2.15) for 
${\theta _m}$, which, along with 
${f_k}$, is the key factor in determining the concentration distribution within the channel through (2.9).
2.3.1. Exchange coefficient
The equation governing 
${f_0}$ can be recovered from (2.18) by setting 
$k = 0$ as
\begin{equation}\frac{{\partial {f_0}}}{{\partial t}} = \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial {f_0}}}{{\partial r}}} \right) + \frac{1}{{{r^2}}}\frac{{{\partial ^2}{f_0}}}{{\partial {\phi ^2}}} - {K_0}{f_0}.\end{equation}Equation (2.21) together with the initial and boundary conditions (2.19) forms a homogeneous partial differential equation that can be readily solved via the method of separation of variables (Constanda Reference Constanda2016), and the result is
\begin{equation}{f_0}(t,r,\phi ) = \mathrm{\mathbb{A}}(t)\sum\limits_{m = 0}^\infty {\sum\limits_{n = 0}^\infty {{A_{mn}}{J_m}({\lambda _{mn}}r)} } \cos (m\phi ){\rm exp} ( - \lambda _{mn}^2t),\end{equation}where
\begin{equation}\mathrm{\mathbb{A}}(t) = {\rm exp} ({\mathrm{\mathbb{C}}_0}),\quad {\mathrm{\mathbb{C}}_0}(t) ={-} \int_0^t {{K_0}\textrm{(}\tau \textrm{)}\,\textrm{d}\tau } ,\end{equation}
and 
${J_m}$ is the Bessel function of the first kind and the order m. Moreover, the eigenvalues 
${\lambda _{mn}}$ are the roots of the following nonlinear equation:
\begin{equation}m{J_m}({\lambda _{mn}}) - {\lambda _{mn}}{J_{m + 1}}({\lambda _{mn}}) ={-} Da\,{J_m}({\lambda _{mn}}).\end{equation}
To obtain the unknown coefficients 
${A_{mn}}$, the initial condition (2.19a) must be invoked. Following the orthogonality conditions after integrating the resultant equation over the dimensionless cross-sectional area, 
${A_{mn}}$ is found to be
\begin{equation}{A_{mn}} = \frac{1}{{{\psi _m}{\mathrm{\mathbb{N}}_{mn}}}}\int_0^1 {\int_0^{\rm \pi} {r\psi (r,\phi )} } {J_m}({\lambda _{mn}}r)\cos (m\phi )\,\textrm{d}\phi \,\textrm{d}r,\end{equation}where
\begin{align}{{\mathbb N}_{mn}} &= \int_0^1 {\int_0^{\rm \pi} {rJ_m^2({\lambda _{mn}}r){{\cos }^2}(m\phi )\,\textrm{d}\phi \,\textrm{d}r} }\nonumber\\ &= \frac{{\rm \pi} }{2}(1 + {\delta _{m0}})\left[ {\frac{1}{2}J_m^2({\lambda_{mn}}) - \frac{{m{J_m}({\lambda_{mn}}){J_{m + 1}}({\lambda_{mn}})}}{{{\lambda_{m,n}}}} + \frac{1}{2}J_{m + 1}^2({\lambda_{mn}})} \right].\end{align}
Now, the time-dependent coefficient 
$\mathrm{\mathbb{A}}$ may be evaluated utilizing the integral condition (2.20) as
\begin{align}\mathrm{\mathbb{A}}(t) &= \frac{{\rm \pi} }{2}{\left[ {\sum\limits_{m = 0}^\infty {\sum\limits_{n = 0}^\infty {{A_{mn}}{\rm exp} ( - \lambda_{mn}^2t)} } \int_0^1 {\int_0^{\rm \pi} {r{J_m}({\lambda_{mn}})\cos (m\phi )\,\textrm{d}\phi \,\textrm{d}r} } } \right]^{ - 1}}\nonumber\\ &= {\left[ {2\sum\limits_{n = 0}^\infty {{A_{0n}}} \frac{{{J_1}({\lambda_{0n}})}}{{{\lambda_{0n}}}}{\rm exp} ( - \lambda_{0n}^2t)} \right]^{ - 1}}.\end{align}
Having completed the solution of 
${f_0}$, we are now able to derive the expression of the exchange coefficient 
${K_0}$ by making use of (2.14):
\begin{equation}{K_0}(t) = \frac{2}{{\rm \pi} }\int_0^{\rm \pi} {{{\left. {\frac{{\partial {f_0}}}{{\partial r}}} \right|}_{r = 1}}\,\textrm{d}\phi } ={-} 2Da\,\mathrm{\mathbb{A}}(t)\sum\limits_{n = 0}^\infty {{A_{0n}}{J_0}({\lambda _{0n}}){\rm exp} ( - \lambda _{0n}^2t)} .\end{equation}
It should be pointed out that the exchange coefficient does not depend on the velocity distribution. However, as evidenced by (2.19), 
${K_0}$ is a function of both the initial form of the solute band and the rate constant of the reaction.
2.3.2. Convection coefficient
As soon as 
${f_0}$ and 
${K_0}$ are determined, the convection coefficient can be obtained via solving the equation governing 
${f_1}$, which is, according to (2.18), given as
\begin{equation}\frac{{\partial {f_1}}}{{\partial t}} = \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial {f_1}}}{{\partial r}}} \right) + \frac{1}{{{r^2}}}\frac{{{\partial ^2}{f_1}}}{{\partial {\phi ^2}}} - (u + {K_1}){f_0} - {K_0}{f_1}.\end{equation}Even though the associated boundary conditions, given by (2.19b), are all homogeneous, (2.29) itself is non-homogeneous and, hence, cannot be solved directly by the method of separation of variables. Here, the solution is sought via the eigenfunction expansion method (Constanda Reference Constanda2016), which provides the following solution:
\begin{equation}{f_1}(t,r,\phi ) = \mathrm{\mathbb{A}}(t)\sum\limits_{i = 0}^\infty {\sum\limits_{j = 0}^\infty {[{\mathrm{\mathbb{B}}_{1,ij}}(t) + {A_{ij}}{\mathrm{\mathbb{C}}_1}(t)]} } {J_i}({\lambda _{ij}}r)\cos (i\phi ){\rm exp} ( - \lambda _{ij}^2t),\end{equation}wherein
Here, 
represents the coefficients associated with the expansion of the term 
$- u{f_0}$ in (2.29), that is
Considering the orthogonality conditions, 
is obtained as
where
The time-dependent coefficient 
${\mathrm{\mathbb{B}}_{1,ij}}$ can now be evaluated as
The integral constraint (2.20) should be applied in order to obtain 
${\mathrm{\mathbb{C}}_1}(t )$. The result reads
\begin{equation}\begin{array}{ccccc} {\mathrm{\mathbb{C}}_1}(t) & = \dfrac{{\sum\limits_{i = 0}^\infty {\sum\limits_{j = 0}^\infty {\left[ {{\mathrm{\mathbb{B}}_{1,ij}}(t){\rm exp} ( - \lambda_{ij}^2t)\int_0^1 {\int_0^{\rm \pi} {r{J_i}({\lambda_{ij}}r)\cos (i\phi )\,\textrm{d}\phi \,\textrm{d}r} } } \right]} } }}{{\sum\limits_{i = 0}^\infty {\sum\limits_{j = 0}^\infty {\left[ {{A_{ij}}{\rm exp} ( - \lambda_{ij}^2t)\int_0^1 {\int_0^{\rm \pi} {r{J_i}({\lambda_{ij}}r)\cos (i\phi )\,\textrm{d}\phi \,\textrm{d}r} } } \right]} } }}\\ & = \dfrac{{\sum\limits_{j = 0}^\infty {\left[ {{\mathrm{\mathbb{B}}_{1,0j}}(t)\dfrac{{{J_1}({\lambda_{0j}})}}{{{\lambda_{0j}}}}{\rm exp} ( - \lambda_{0j}^2t)} \right]} }}{{\sum\limits_{j = 0}^\infty {\left[ {{A_{0j}}\dfrac{{{J_1}({\lambda_{0j}})}}{{{\lambda_{0j}}}}{\rm exp} ( - \lambda_{0j}^2t)} \right]} }}. \end{array}\end{equation}At this point, (2.14) can be again consulted, this time to evaluate the convection coefficient as
where
with 
${}_1\mathrm{\mathbb{F}}{\mathrm{\mathbb{R}}_2}$ given as
\begin{align}{}_1\mathrm{\mathbb{F}}{\mathrm{\mathbb{R}}_2}(s,m,n) &= \int_0^1 {{r^s}{J_0}({\lambda _{mn}}r\textrm{)}\,\textrm{d}r}\nonumber\\ &= \frac{{\lambda _{mn}^m}}{{{2^{1 + m}}}}\varGamma \left( {\frac{{1 + m + s}}{2}} \right){}_1F{R_2}\left( {\frac{{1 + m + s}}{2};\frac{{3 + m + s}}{2},1 + m; - \frac{{\lambda_{mn}^2}}{4}} \right),\end{align}
in which 
${}_1F{R_2}$ and 
$\varGamma $ are the regularized hypergeometric and gamma functions, respectively.
2.3.3. Dispersion coefficient
The governing equation for 
${f_2}$ can be recovered by setting 
$k = 2$ into (2.18) as
\begin{equation}\frac{{\partial {f_2}}}{{\partial t}} = \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial {f_2}}}{{\partial r}}} \right) + \frac{1}{{{r^2}}}\frac{{{\partial ^2}{f_2}}}{{\partial {\phi ^2}}} - {K_0}{f_2} - (u + {K_1}){f_1} + (P{e^{ - 2}} - {K_2}){f_0},\end{equation}
which can be solved utilizing much the same procedure adopted for solving 
${f_1}$ to yield
\begin{equation}{f_2}(t,r,\phi ) = \mathrm{\mathbb{A}}(t)\sum\limits_{s = 0}^\infty {\sum\limits_{q = 0}^\infty {[{\mathrm{\mathbb{B}}_{2,sq}}(t) + {A_{sq}}{{\mathbb C}_2}(t)]{J_s}({\lambda _{sq}}r)\cos (s\phi ){\rm exp} ( - \lambda _{sq}^2t)} } ,\end{equation}where
Here, 
stands for the coefficients associated with the expansion of the term 
$- (u + {K_1}){f_1}$ in (2.40), and is given as
Taking advantage of the integral condition (2.20), 
${\mathrm{\mathbb{C}}_2}$ is obtained as follows:
\begin{equation}{\mathrm{\mathbb{C}}_2}(t) ={-} \frac{{\sum\limits_{q = 0}^\infty {\left[ {{\mathrm{\mathbb{B}}_{2,0q}}(t)\frac{{{J_1}({\lambda_{0q}})}}{{{\lambda_{0q}}}}{\rm exp} ( - \lambda_{0q}^2t)} \right]} }}{{\sum\limits_{q = 0}^\infty {\left[ {{A_{0q}}\frac{{{J_1}({\lambda_{0j}})}}{{{\lambda_{0j}}}}{\rm exp} ( - \lambda_{0q}^2t)} \right]} }}.\end{equation}Now, the dispersion coefficient can be derived from (2.14) as follows:
2.3.4. General solutions for higher-order transport coefficients
Considering the infinite set of (2.18) along with the relevant initial and boundary conditions, given by (2.19), the following general solution can be proposed for the functions 
${f_k}$ when 
$3 \le k$:
in which
\begin{equation}{\mathrm{\mathbb{C}}_k}(t) ={-} {\delta _{k0}}\int_0^t {{K_0}\textrm{(}\tau \textrm{)}\,\textrm{d}\tau + (1 - {\delta _{k0}})} \int_0^t {\textrm{[}{\delta _{k2}}P{e^{ - 2}} - {K_k}\textrm{(}\tau \textrm{)]}\,\textrm{d}\tau } ,\end{equation}
where the coefficients 
are obtained via the following integral:
with
\begin{equation}{f^{\prime}_k} = \frac{{{f_{k - 2}}}}{{P{e^2}}} - u{f_{k - 1}} - \sum\limits_{i = 1}^{k - 1} {{K_i}{f_{k - i}}} \textrm{,}\quad k = 3,4,5, \ldots .\end{equation}
Making use of the integral constraint (2.20), 
${\mathrm{\mathbb{C}}_k}(t)$ is determined as
Consequently, one can obtain the following formula for the transport coefficient 
${K_k}$ with 
$3 \le k$:
2.3.5. Solute concentration distribution
It is now time to present a solution for the mean solute concentration since the exchange, convection and dispersion coefficients, appearing in (2.15), have all been obtained. The solution to (2.15) subject to the associated initial and boundary conditions, derived via the cross-sectional averaging of (2.8a) and (2.8b), can be obtained utilizing the infinite Fourier transform method. The pertinent procedure, along with details on how the solution can be extended to situations in which the axial variation of the initial solute concentration is represented by functions other than a Dirac delta function, is given in Appendix A. The final expression reads
\begin{equation}{\theta _m}(t,x) = \frac{{{\psi _m}}}{{2Pe\sqrt {{\rm \pi} \xi } }}{\rm exp} \left[ { - {\mathrm{\mathbb{C}}_0} - \frac{{{{(x - {\mathrm{\mathbb{C}}_1})}^2}}}{{4\xi }}} \right],\end{equation}
where 
$\xi (t) = P{e^{ - 2}}t - {\mathrm{\mathbb{C}}_2}$. Accordingly, the dimensionless local concentration can be readily determined utilizing (2.9) as
\begin{equation}\theta (t,r,\phi ,x) = \frac{{{\psi _m}}}{{2Pe\sqrt {{\rm \pi} \xi } }}\left\{ {{f_0} - \frac{{x - {\mathrm{\mathbb{C}}_1}}}{{2\xi }}{f_1} + \left[ { - \frac{1}{{2\xi }} + \frac{{{{(x - {\mathrm{\mathbb{C}}_1})}^2}}}{{4{\xi^2}}}} \right]{f_2}} \right\}{\rm exp} \left[ { - {\mathrm{\mathbb{C}}_0} - \frac{{{{(x - {\mathrm{\mathbb{C}}_1})}^2}}}{{4\xi }}} \right].\end{equation}2.3.6. Initial solute sources
As mentioned previously, two different shapes are considered for the initial form of the solute band: circular and semicircular distributions over the cross-sectional area, both having an infinitesimal length in the axial direction. The circular distribution considered, with radius 
${R_D}$, is centred at 
$R = {C_D},\; \phi = {\rm \pi}/2$. Using dimensionless coordinates, 
$\psi $ and 
${\psi _m}$ associated with the circular solute band are given as
\begin{equation}\psi (r,\phi ) = \left\{ {\begin{array}{*{20}{l}} 1&{\textrm{for}\;{\omega_1} \le \phi \le {\omega_2}\;\textrm{and}\;{{{\mathcal{R}}}_1} \le r \le {{{\mathcal{R}}}_2}}\\ 0&{\textrm{otherwise}} \end{array}} \right.,\quad {\psi _m} = 2r_D^2,\end{equation}
where 
${\omega _{1,2}} = ({\rm \pi} /2) \mp \textrm{asin(}{r_D}/{c_D}\textrm{)}$ and 
${{{\mathcal{R}}}_{1,2}} = {c_D}\sin (\phi ) \mp {[r_D^2 - c_D^2{\cos ^2}\phi ]^{1/2}}$ with 
${c_D}$ and 
${r_D}$ being the dimensionless forms of 
${C_D}$ and 
${R_D}$. The initial concentration distribution can only influence the coefficients 
${A_{mn}}$, given by (2.24), which, for this case, is modified as
\begin{equation}{A_{mn}} = \frac{1}{{{\psi _m}{\mathrm{\mathbb{N}}_{mn}}}}\int_{{\omega _1}}^{{\omega _2}} {\int_{{{{\mathcal{R}}}_1}}^{{{{\mathcal{R}}}_2}} {r{J_m}({\lambda _{mn}}r)\cos (m\phi )\,\textrm{d}r\,\textrm{d}\phi } } .\end{equation}
Because of the complex forms of the integral limits on both r and 
$\phi $, 
${A_{mn}}$ has to be calculated by standard numerical integration techniques for a circular initial concentration. This is the only place where numerical methods are used in the present work. For the case with a semicircular initial concentration distribution, the function 
$\psi $ and its cross-sectionally averaged value, 
${\psi _m}$, are given as
\begin{equation}\psi (r,\phi ) = \left\{ {\begin{array}{*{20}{l}} 1&{\textrm{for}\;0 \le r \le a}\\ 0&{\textrm{otherwise }} \end{array}} \right.,\quad {\psi _m} = {a^2},\end{equation}
where a is the dimensionless radius of the solute band. It can be shown that 
${A_{mn}} = 0$ when 
$m \ne 0$ for this case and
\begin{equation}{A_{0n}} = \frac{{\rm \pi} }{{{\psi _m}{\mathrm{\mathbb{N}}_{0n}}}}\int_0^a {r{J_0}\textrm{(}{\lambda _{0n}}r\textrm{)}\,\textrm{d}r} = \frac{{{\rm \pi} {J_1}({\lambda _{0n}}a)}}{{a{\lambda _{0n}}{\mathrm{\mathbb{N}}_{0n}}}} = \frac{{2{J_1}({\lambda _{0n}}a)}}{{a{\lambda _{0n}}J_0^2({\lambda _{0n}})}}.\end{equation}
It should be noted that the semicircular distribution reduces to a uniform distribution when 
$a = 1$.
The solutions developed in this section can be simplified under certain conditions. The two main special cases for which such simplified solutions are possible include when there is no wall reaction or when the velocity profile is uniform. The solutions to these special cases are given in Appendices B and C.
3. Results and discussion
Considering (2.7) and (2.8), it is clear that the parameters controlling the broadening of a solute band by fully developed flow in semicircular microchannels with surface reaction include the initial distribution function 
$\psi $, the Péclet number 
$Pe$ and the Damköhler number 
$Da$. Depending on the initial form of the injected band, the function 
$\psi $ may itself include several controlling parameters. Considering the two initial concentration distributions considered in the present study, these additional controlling parameters include 
${c_D}$, 
${r_D}$ and a. A parametric study is performed in this section to reveal the degree to which each of the above-mentioned parameters influences the mass transport characteristics. Unless otherwise stated, all the results are obtained by setting 
$Pe = 100$ and utilizing appropriate numbers of the series terms, given in table 1. This table also contains a convergence analysis that is performed by comparing the values of 
${u_m}$, 
${K_0}$, 
${K_1}$, 
${K_2}$ and 
${\theta _{m,max}}$ obtained utilizing different series terms. Note that the times considered for the last four parameters correspond to the situations where the highest sensitivity of these parameters to the number of series terms exists. The selection of a CID with 
${c_D} = \; {r_D} = 0.5$ is also for exactly the same reason. Hence, the required number of the series terms to obtain converged results, given in table 1, which is 100 for 
${K_0}$ and 10 for other parameters, holds for any circumstance since the worst cases have been considered in the convergence analysis. Moreover, to validate the solutions developed, numerical simulations were also performed using COMSOL Multiphysics software, version 5.5. The set of the governing equations, including (2.2a) and (2.7), was solved subject to the pertinent initial and boundary conditions, given by (2.2b) and (2.8), using PDEs branch of the Mathematics module of the software. The results of the numerical solutions are presented alongside the analytical solution results when discussing the mean concentration field.
Table 1. Dependence of the dimensionless mean velocity, 
${K_{0,1,2}}$, and the maximum dimensionless mean concentration 
${\theta _{m,max}}$ on the number of terms considered in the series solutions. A circular initial concentration distribution with 
${c_D} = {r_D} = 0.5$ along with 
$Da = 1$ is considered for obtaining the mass transport results.
3.1. Exchange coefficient
The discussion of results begins with investigating the effects of different influential parameters on the exchange coefficient, which is a measure of solute exchange between the wall and the flow due to surface adsorption or reaction. The influence of parameter a, which controls the extent of injection for a SID, on the exchange coefficient at different values of 
$Da$ is shown in figure 2. The first point drawing attention in figure 2 is that the exchange coefficient is an increasing function of 
$Da$ (and of the same order of magnitude as 
$Da$), which is anticipated recalling the definition of 
${K_0}$ and the fact that the rate of surface reaction depends directly on 
$Da$. It is also observed in figure 2 that the variation of 
${K_0}$ with time is a strong function of the solute injection. For instance, while the exchange coefficient is initially an increasing function of time for 
$a \ne 1$, quite the opposite is true for 
$a = 1$, that is, when solutes are initially distributed uniformly over the whole channel cross-section. For the former, solutes initially do not exist at the curved wall, rendering 
${K_0}$ zero at very small times. As time goes by, solutes start to spread everywhere including the near-wall area. Hence, surface reaction gradually starts to occur, rendering 
${K_0}$ an increasing function of time. At these conditions, an increase in a gives rise to larger exchange coefficients due to the fact that the higher a is, the sooner solutes reach the curved wall. For 
$a = 1$, however, there is an abundance of solutes adjacent to the wall at time zero; hence, 
${K_0}$ takes a large value. The situation is gradually changed because of solute depletion in the bulk due to the advection mechanism that hinders the diffusion of solutes to the curved wall, thereby reducing the exchange coefficient, as shown in figure 2. Irrespective of a, all the graphs collapse to one in the long term, where the influences of the initial conditions vanish.
Figure 2. (a–d) Variations of the exchange coefficient with the dimensionless time for the SID at different values of a and 
$Da$.
Considering different values for 
${c_D}$ and 
${r_D}$, variations of the exchange coefficient over time for a CID are plotted in figure 3. Comparing figures 2 and 3, it is deduced that, despite the long-term conditions for which no dependence of results on the form of injection is observed, the function 
$\psi (r,\phi )$ plays a very important role at short times. The location of injection, controlled via 
${c_D}$, is found to significantly affect the variations of 
${K_0}$ with time. At a fixed injection size, that is, for a fixed 
${r_D}$, higher values of 
${K_0}$ are achieved for a larger 
${c_D}$. According to figure 1, a larger 
${c_D}$ implies that the point of injection gets closer to the curved wall, where the reaction takes place. Accordingly, more solutes will be available to react, leading to larger solute exchange rates between the liquid and the wall. Exactly the same behaviour is observed when 
${r_D}$ is increased while keeping 
${c_D}$ fixed.
Figure 3. (a–d) Profiles of the exchange coefficient versus the dimensionless time for the circular initial concentration distribution considering different 
${r_D}$ and 
${c_D}$ for both 
$Da = 0.1$ and 
$Da = 10$.
As expected, 
${K_0}$ is zero at the beginning for the cases for which no solute is located at the wall at time zero, that is, when 
${r_D} + {c_D} < 1$. Under such circumstances, the value of 
${K_0}$ might be kept at zero up to the dimensionless times of the order of 
${10^{ - 2}}$ by injecting the solute band at an appropriate place and with a sufficiently small radius; this is important for applications where solute adsorption at the wall is to be avoided. Although for the cases with 
${r_D} + {c_D} = 1$, similar to SID with 
$a = 1$, non-zero values are observed for 
${K_0}$ at time zero; however, the values are now significantly smaller. This is due to the fact that, unlike the SID case with 
$a = 1$ for which all the curved wall is in contact with solutes at time zero, such contact exists only at one point for the CID case. Accordingly, solute exchange between the fluid and the wall will be much less significant for the CID case. An important phenomenon that is observed for both SID and CID cases, but is more pronounced for the latter, is the occurrence of non-monotonic trends in the plots of 
${K_0}$ versus time under some circumstances. Such a behaviour arises when a non-zero solute concentration exists in the close vicinity of the curved wall, but not over the wall, at time zero. Hence, the concentration front soon reaches the wall, leading to the sharp increase of 
${K_0}$ from zero to a large value. A large value of 
${K_0}$ implies that there is a high degree of solute consumption at the reactive wall, which leads to the depletion of the solution therein. Since, due to solute advection in the bulk, the solute concentration away from the wall is also low, the depletion of the solution near the wall cannot be compensated for; hence, a reduction occurs in 
${K_0}$. Since the depletion of the solution is more noticeable when 
$Da$ is high, because of higher reaction rates, the mentioned phenomenon is more pronounced for larger values of 
$Da$, as evidenced by both figures 2 and 3.
3.2. Convection coefficient
Now, we are ready to focus on 
${K_1}$. Plotted in figure 4 are the absolute values of 
${K_1}$ for SID versus the normalized time at different a and 
$Da$. When 
$a \ne 1$, the convection coefficient is found to be independent of 
$Da$ at small times, since there is no solute at the wall for the reaction to take place. For 
$a = 1$, however, the influence of 
$Da$ is felt even at very small times because solutes are present at the wall at the beginning of the process. For this case, a higher 
$Da$ gives rise to higher values of 
${K_1}$ (ignoring the negative sign) because the depletion of the solution within the low-velocity wall-adherent area due to the wall reaction alters the convection transport in favour of the high-velocity core region. Moreover, depending on a and 
$Da$, different time-dependence behaviours are observed. Whereas an increasing trend is generally observed for 
$a = 1$ as a result of the reaction-caused concentration alterations discussed above, quite the opposite is observed for 
$a = 0.7$ mainly because of the diffusion of solutes from the core to the wall-adherent area. For cases such as 
$a = 0.3$, non-uniform trends are observed: the convection coefficient increases first, and then it may either decrease or become nearly constant for a while followed by a subsequent increase, depending on 
$Da$. The reason that 
${K_1}$ increases first for 
$a = 0.3$ may be attributed to the fact that the majority of the sample is injected near the flat wall where fluid velocity is low. Hence, the diffusional transport of solutes towards the curved wall is initially accompanied by an increase in the convection rate. This trend is changed as soon as a considerable number of solutes reach adjacent to the curved wall. Now, the accumulation of solutes near the curved wall results in a reduction of 
${K_1}$ when 
$Da$ is not high enough, as is the case for 
$Da = 0.01$ and 
$Da = 1$. When 
$Da = 10$, the surface reaction is high enough to disrupt solute accumulation, rendering 
${K_1}$ a nearly constant function of time for a while after which an increasing trend is again observed as a result of further depletion of the solute band near the wall. As expected, all the results become independent of time in the long term where also there is no dependence on a due to the damping of the initial condition.
Figure 4. Time dependence of the convection coefficient for a semicircular initial injection at different values of a and 
$Da$. The top and the right-hand axes, drawn in green, belong to the dependence of the convection coefficient on a when 
$Da = 1$ and 
$t = {10^{ - 4}}$, the results of which are shown by a green line with circular symbols.
Besides the time dependence of 
${K_1}$, its dependence on a for 
$t = {10^{ - 4}}$ is also shown in figure 4, indicating an increasing trend at smaller values of a that is reversed at 
$a \cong 0.68$. This is because whereas at lower values of a increasing the radius of the solute band results in the addition of a solute strip that mostly lies in the bulk, the increase in a at its higher values is made via the addition of a solute strip that is located in the low-velocity wall-adherent area.
The impact of the Damköhler number on the time development of the convection coefficient for different circular injections is illustrated in figure 5. Shown in figure 5(a) are plots of the convection coefficient at different values of 
$Da$ and 
${r_D}$ while keeping 
${c_D} = 0.5$. While, as anticipated, the 
${r_D}$-dependence of the results vanishes in the long term, higher convection coefficients are achieved by decreasing 
${r_D}$ in the short term. This occurs because the centre of the injected band is very close to the point where the maximum velocity occurs 
$(r \cong 0.48)$; hence, the solute band is concentrated more around the maximum velocity point by decreasing 
${r_D}$, leading to enhanced downstream convection of mass. In the beginning, the Damköhler number is observed to have no effect on 
${K_1}$ because the solute band, or at least a considerable part of it, is still away from the curved wall. Since less time is required for a wider solute band to reach the curved wall and take part in the reaction, the Damköhler number effect is felt sooner for a larger 
${r_D}$. Unlike the SID case, the convection coefficient is now a monotonically decreasing function of time in response to the diffusion of solutes from the high-velocity core area to the low-velocity wall-adherent region. Such a monotonic behaviour is also observed in figure 5(b), which plots the graphs of the convection coefficient at different 
$Da$ and 
${c_D}$ while keeping 
${r_D} = 0.3$, except for the case with 
$Da = 10$ and 
${c_D} = 0.7$. For this special case, the solute band is sufficiently close to the curved wall and the reaction rate is so high that a considerable reduction occurs in solute concentration near the wall. Hence, the convection coefficient increases first before reducing to its long-term value. Inspecting the short-term results of the convection coefficient for the three values of 
${c_D}$ considered, namely 
$0.3$, 
$0.5$ and 
$0.7$, indicates that the maximum convection rate belongs to the case with 
${c_D} = 0.5$ for which the centroid of the solute band locates near the maximum velocity point.
Figure 5. Influence of the Damköhler number on the time dependence of the convection coefficient for the circular initial condition. (a) The graphs are plotted at three different 
${r_D}$ while keeping 
${c_D} = 0.5$ whereas (b) 
${c_D}$ is varied assuming 
${r_D} = 0.3$.
3.3. Dispersion coefficient
The influences of both 
$Da$ and a on the time development of the dispersion coefficient for the SID case are investigated in figure 6. It can be seen that 
${K_2} - P{e^{ - 2}}$ is initially zero and increases monotonically to its long-term value irrespective of 
$Da$ and a. The parameter 
${K_2} - P{e^{ - 2}}$ denotes the part of the dispersion coefficient that is solely due to the velocity variations. Hence, it takes a zero value at the beginning when the mass transport mechanism is solely due to molecular diffusion owing to the creation of large concentration gradients after the injection. As time passes, such large gradients are slowly damping; therefore, advection starts to contribute to the mass transport mechanism and 
${K_2} - P{e^{ - 2}}$ takes non-zero values until it reaches an asymptotic value at sufficiently long times. It is also observed in figure 6 that the transient dispersion coefficient is smaller for a smaller a due to the fact that the injected solute band gets farther from the curved wall where there is a region of high shear rate. In addition, a higher 
$Da$ results in a lower dispersion coefficient. This occurs because increasing the reaction rate gives rise to the reduction of the solute concentration near the curved wall, thereby weighting the concentration distribution in favour of the low-shear-rate bulk area.
Figure 6. Time development of the dispersion coefficient for the semicircular initial concentration distribution at different values of 
$Da$ and a.
An inspection of the effect of Damköhler number on the dispersion coefficient at different versions of the circular initial concentration distribution, illustrated in figure 7, confirms the results of figure 6 that 
$Da$ is the only factor governing the long-term values of the dispersion coefficient. Before reaching a quasi-steady state, however, the dispersion coefficient is a strong function of both position and the extent of injection. Comparing the results for 
${r_D} = 0.3$ and 
${r_D} = 0.5$, both with 
${c_D} = 0.5$, shows that the transient values of 
${K_2}$ are smaller for a smaller 
${r_D}$. This is not surprising recalling that when 
${r_D}$ decreases while having a fixed 
${c_D}$ of 
$0.5$ the injection is more concentrated in the low-shear-rate bulk area. On the other hand, when 
${c_D}$ increases from 
$0.3$ to 
$0.7$ while keeping 
${r_D} = 0.3$, a reduction occurs in the transient values of 
${K_2}$. The main reason is that, since the solute band is closer to the reactive wall for 
${c_D} = 0.7$, there is a higher consumption of the solutes for this case, leading to the depletion of the solution near the wall and, in turn, weighting the concentration field in favour of the low-shear-rate bulk region. This can be verified by the fact that the influence of 
${c_D}$ on the results is more pronounced for a higher 
$Da$. Note that scrutiny of the results for 
${r_D} = 0.3$ and 
${c_D} = 0.5$ indicates that the dispersion coefficient is practically zero up to 
$t \cong 0.002$ that corresponds to a time of 
$5\;\textrm{s}$ or a channel of radius 
$500\;\mathrm{\mu }\textrm{m}$ and a molecular diffusivity of 
${10^{ - 10}}\;{\textrm{m}^2}\;{\textrm{s}^{ - 1}}$. Then, considering a mean fluid velocity of 
$1\;\textrm{cm}\;{\textrm{s}^{ - 1}}$, it would be possible to transfer a solute band along a channel of length 
$5\;\textrm{cm}$ and shorter with practically no dispersion. This highlights the potentials of the semicircular geometry for liquid-phase transportation with little dispersion in typical microchannels by appropriately choosing the form and position of solute injection.
Figure 7. Time dependence of the dispersion coefficient for the circular initial concentration distribution at different 
$Da$, 
${r_D}$ and 
${c_D}$. The lower abscissa pertains to the graphs without a symbol while the upper abscissa corresponds to those with green symbols.
3.4. Long-term transport coefficients
The values of the transport coefficients, especially the dispersion coefficient, in the long-term limit, or as it is usually referred to, the Taylor–Aris regime, are of significant practical importance. The reason is that, for sufficiently long channels, the overall dispersion is mainly dominated by the long-term dispersion coefficients. Considering the very short length of most microfluidic devices, one might conclude that the Taylor–Aris dispersion regime is not observed in microchannels. However, this is not the whole story. By referring to the foundations of the Taylor–Aris theory, it can be inferred that the long-term dispersion coefficient is valid for time scales larger than the time required for the solutes to travel radially across the channel using the molecular diffusion mechanism, given as 
$R_{sc}^2/D$. This is verified by our presented results indicating that the values of the transport coefficients tend to asymptotic values when the dimensionless time, 
$TD/R_{sc}^2$, is of the order of 1. It is obvious that the molecular diffusion time scale 
$R_{sc}^2/D$ decreases rapidly with decreasing the channel radius. Hence, it takes very small values at the microscale. Moreover, flow velocities are much smaller at the microscale as compared with those in conventional channels, implying that the solute band may still need a significant time to reach the channel outlet. Hence, the Taylor–Aris dispersion regime might practically be important for microfluidic devices, especially when dealing with long microchannels like those with a serpentine geometry.
Based on the above discussion, the plots of the transport coefficients with respect to 
$Da$, given in figure 8, may be considered as one of the most important findings of this study. Figure 8 confirms that the long-term values of both the exchange and convection coefficients are monotonically increasing functions of 
$Da$, whereas the opposite is true for the dispersion coefficient. Note that the results should finally tend to asymptotic values at very high 
$Da$, but such values will be out of the practical range of the Damköhler number and, hence, are not considered here. Closely related to figure 8, one of the most important results for practitioners is the long-term value of the dispersion coefficient for the special case for which no reaction occurs at the walls. Our calculations show that the long-term value of 
${K_2} - P{e^{ - 2}}$ for 
$Da = 0$ is 
$5.49 \times {10^{ - 4}}$, which is one order of magnitude smaller than that for a circular channel (Sankarasubramanian & Gill Reference Sankarasubramanian and Gill1973). The main reason is that the mean velocity of a semicircular channel for a given pressure gradient is significantly smaller than that for a circular channel. In fact, considering the value of the dimensionless mean velocity, which is calculated to be 
$0.1894$, and the reference velocity, which is twice the mean velocity in a circular channel, the semicircular to circular channel mean velocity ratio is only 
$0.3788$. Note that, recalling that 
${K_2}$ scales inversely with the square of velocity, the long-term value of 
${K_2} - P{e^{ - 2}}$ based on the mean fluid velocity, which is calculated to be 
$5.49 \times {10^{ - 4}}/{0.1894^2} = 0.0153$, is still smaller than the well-known value for a circular geometry, given as 
$1/48 = 0.0208$ (Taylor Reference Taylor1953). However, this does not necessarily mean that the Taylor–Aris dispersion in a semicircular duct is lower than that in a circular channel because the dispersion coefficient is not the only determining factor. For a meaningful comparison, the same flow rates in the ducts should be considered. Since the cross-sectional area of a semicircular channel is half that of a circular channel with the same radius, the mean velocity in the former should be twice that in the latter to yield the same flow rate. Accordingly, the long-term hydrodynamic dispersion in a semicircular channel will be higher than that in a circular duct, as it should be.
Figure 8. Long-term values of the transport coefficients plotted as functions of the Damköhler number.
3.5. Solute concentration
It is now time to turn our attention to the dimensionless mean concentration profiles, which are plotted using (2.53) only. We start with studying the axial distributions of 
${\theta _m}$ at different 
$Da$ for a uniform initial distribution, illustrated in figure 9 for time 
$t = 1$. The first obvious point in figure 9 is the decreasing dependence of the maximum mean concentration on the reaction rate that is owing to the higher consumption of solutes at the wall. In addition, the location of the maximum 
${\theta _m}$ is shifted downstream when higher values of 
$Da$ are considered. This is a direct consequence of the fact that the convection coefficient is larger for a higher reaction rate, resulting in a faster overall movement of the solute band along the channel. Note that the connection between the location of the maximum 
${\theta _m}$ and the convection coefficient can also be inferred from a mathematical point of view. Considering the formula of 
${\theta _m}$, it can be deduced that the maximum mean concentration corresponds to the point 
$x = {\mathrm{\mathbb{C}}_1}(t) ={-} \int_0^t {{K_1}\,\textrm{d}\tau }$, which is a direct function of the convection coefficient. Note that for no reaction case 
${\mathrm{\mathbb{C}}_1}(t) = {u_m}t$, which means that the location of the maximum 
${\theta _m}$ moves with a speed equal to the average fluid velocity, as predicted by the classical Taylor dispersion theory (Taylor Reference Taylor1953).
Figure 9. Dimensionless mean concentration plotted versus the dimensionless axial coordinate at different 
$Da$. The results pertain to a uniform distribution of the initial concentration field and correspond to time 
$t = 1$.
In the following, we study the influence of changing the radius of the semicircular injection on the time development of the dimensionless mean concentration profiles when 
$Da = 1$. In figure 10, as time passes, the concentration profiles get wider as a result of the hydrodynamic dispersion. The concentration peak is decreased at the same time due to both widening of the profile and solute consumption at the wall because of surface reactions. As is clear, an increase in a leads to higher solute concentrations. The reason is that the amount of solute being injected into the channel increases with increasing a. Note that all of the graphs are very similar to a Gaussian distribution, implying that the Taylor–Aris dispersion regime is reached for 
$0.5\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel\prec\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }t$. This is consistent with the results presented in figure 6.
Figure 10. Axial dependence of the dimensionless mean concentration distribution at different times for the semicircular initial concentration distribution with 
$\; a = 0.5,\; 0.8$ considering 
$Da = 1$.
The influence of increasing the amount of solute being injected into the channel is much more visible in figure 11 that plots the dimensionless mean concentration graphs for the CID case at different 
${c_D}$ and 
${r_D}$ for fixed t and 
$Da$. A modified axial coordinate, given as 
$x - {\mathrm{\mathbb{C}}_1}(t)$, has been used here to fix the location of the maximum mean concentration at the centre of the coordinate system for a better comparison. The predictions of the numerical simulations are given alongside the analytical solution results in this figure to inspect the validity of the solutions developed. As observed, the results are in a good agreement; the small discrepancy between the results may be attributed to the limitations associated with accurately modelling the Dirac delta function in numerical simulations when applying the initial condition. It is clear in figure 11 that an increase in the radius of the concentration release at a fixed 
${c_D}$ gives rise to larger values of 
${\theta _m}$. On the other hand, the solute concentration reduces with increasing 
${c_D}$ at a given 
${r_D}$. The reason is that, by moving the centre of injection downward, the distance to the curved wall is reduced, implying that more solutes will be available to take part in surface reactions. The final outcome will be more consumption of solutes at the wall, leading to the depletion of the solute band. All in all, figure 11 indicates that the mean concentration can be effectively controlled by the initial solute distribution.
Figure 11. Plots of the dimensionless mean concentration versus the modified axial coordinate, 
$x - {\mathrm{\mathbb{C}}_1}(t)$, for the circular initial condition at different 
${c_D}$ and 
${r_D}$ while keeping 
$t = 0.5$ and 
$Da = 1$. Line graphs represent the results obtained utilizing the analytical solutions while the symbols show the numerical results.
The last illustration, figure 12, is devoted to studying the time evolution of the cross-sectional contours of 
$\theta $ at the channel inlet for two different Damköhler numbers. Along this line, a circular distribution with 
${c_D} = 0.5$ and 
${r_D} = 0.4$ is considered for the injected solute band and the Péclet number is fixed at 
$10$. As seen in figure 12, the contours are circular at the close vicinity of the centre of injection for 
$t = 0.01$ indicating that the predominant mass transfer mode is still molecular diffusion. Even though the situation is quite the same for both the Damköhler numbers considered at 
$t = 0.01$, significant differences are observed afterward. For example, while the contour lines intersect the curved wall perpendicularly for 
$Da = 0.1$ when 
$t = 0.05$ in response to low surface reaction rates, this is not true for 
$Da = 1$ for which the reaction rates are high. Moreover, although the high concentration area extends from the central region, where the injection is made, to the curved wall for 
$Da = 0.1$, it is only limited to the core for 
$Da = 1$ due to the depletion of the solution near the reactive wall. At 
$t = 0.15$, sufficient time has passed for the advection to take effect. Accordingly, the solutes located in the core are mostly conveyed downstream. Now, solute accumulation is observed in two different regions for 
$Da = 0.1$, one near the flat wall and the other near the curved wall where fluid velocities are low. The second region, however, is not observed for 
$Da = 1$ because of high solute consumption due to surface reactions. The high-concentration areas near the walls for 
$Da = 0.1$ still exist to a great extent at 
$t = 0.3$ while a minimum is established at the core due to high velocities therein that cause significant advection of solutes toward downstream. The situation is totally different for 
$Da = 1$: for this case, the minimum concentration occurs at the curved wall owing to the high reaction rates, and the area of maximum concentration is located near the middle of the flat wall that is at the farthest place from the reactive surface. In fact, the influences of surface reactions are so high as to overshadow the advection effects, leading to concentration contours that are almost parallel to the curved wall, especially in its close proximity, a situation that is expected to be observed only in purely diffusive transport.
Figure 12. Time evolution of the cross-sectional contours of the dimensionless concentration at the inlet for both 
$Da = 0.1$ (a,c,e,g) and 
$Da = 1$ (b,d,f,h), assuming 
$Pe = 10$. Here, a circular distribution with 
${c_D} = 0.5$ and 
${r_D} = 0.4$ is considered for the injected solute band.
4. Conclusions
The dispersion of a solute band by a steady-state and hydrodynamically fully developed laminar flow in a semicircular microchannel was theoretically studied. An irreversible first-order reaction was considered at the curved wall of the channel while considering no-flux boundary conditions at the flat wall so as to follow the physics observed in microfluidic applications. No constraint was imposed on the cross-sectional shape of the initial solute band nor on its axial distribution. The generalized dispersion model was used to obtain analytical solutions for the transport parameters including the exchange, convection and dispersion coefficients as well as for the local and the cross-sectionally averaged concentrations. Even though the solutions were obtained for a pressure-driven flow, special solutions were also presented by assuming a uniform fluid velocity that is typical of electroosmotic flow at high ionic concentrations. Special solutions were also obtained for the conditions where no reaction occurs at the channel walls, a special case that has applications in liquid-phase transportation. A parametric study of the transport properties was then performed by focusing on the influences of the Damköhler number, a measure of the surface reaction rate, and the initial concentration distribution function considering two different injection forms: circular and semicircular cross-sectional distributions with the uniform distribution being a special case of the latter.
It was found that the only factor determining the long-term values of the transport coefficients is the Damköhler number: the dispersion coefficient is a decreasing function of the Damköhler number whereas the opposite is true for the exchange and convection coefficients. In fact, our results showed that the exchange coefficient is of the same order as the Damköhler number. Besides dependence on the Damköhler number, the transport coefficients are strongly influenced by the form of injection in the short term. At such time scales, the dispersion coefficient always increases with time. The situation is more complex for the exchange and convection coefficients: depending on the injection form and the Damköhler number, each of increasing, decreasing and increasing–decreasing trends is observed. In addition, it was shown that by appropriately positioning and shaping the solute injection not only is it possible to avoid unfavourable solute adsorptions at the wall but it is also feasible to perform liquid-phase transportation of samples with little dispersion in semicircular microchannels of typical length.
The inspections of the mean concentration profiles revealed that the mean concentration is smaller for a higher Damköhler number but, at the same time, its location moves faster along the channel. It was shown that the mean concentration graphs become spread as time passes while the maximum concentration decreases. Finally, it was observed that the mean concentration can be effectively controlled by adjusting the initial solute concentration.
Funding
The corresponding author sincerely thanks the Alexander von Humboldt Foundation for its financial support during the course of this work.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Detailed solution of dimensionless mean concentration
The solution to (2.15) is here obtained using the Fourier transform technique. The Fourier transform of 
${\theta _m}$, given as 
$\vartheta $, and its inverse are defined as
\begin{equation}\vartheta (t,\alpha ) = \frac{1}{{\sqrt {2{\rm \pi} } }}\int_{ - \infty }^\infty {{\theta _m}\textrm{(}t,x\textrm{)}\,{\textrm{e}^{ - \textrm{i}\alpha x}}\,\textrm{d}\,x} \textrm{,}\quad {\theta _m} = \frac{1}{{\sqrt {2{\rm \pi} } }}\int_{ - \infty }^\infty {\vartheta \,{\textrm{e}^{\textrm{i}\alpha x}}\,\textrm{d}\alpha } .\end{equation}The application of the Fourier transform to (2.15) results in the following transformed equation:
\begin{equation}\frac{{\textrm{d}\vartheta }}{{\textrm{d}t}} = {K_0}\vartheta + \textrm{i}\alpha {K_1}\vartheta - {\alpha ^2}{K_2}\vartheta ,\end{equation}which is a separable ordinary differential equation and can be readily integrated to yield
\begin{equation}\vartheta = \vartheta (0,\alpha ){\rm exp} [ - {\mathrm{\mathbb{C}}_0} - \textrm{i}\alpha {\mathrm{\mathbb{C}}_1} - {\alpha ^2}\xi ],\end{equation}where
\begin{equation}\vartheta (0,\alpha ) = \frac{1}{{\sqrt {2{\rm \pi} } }}\int_{ - \infty }^\infty {{\theta _m}(0,x)\,{\textrm{e}^{ - \textrm{i}\alpha x}}\,\textrm{d}\,x} = \frac{1}{{\sqrt {2{\rm \pi} } }}\int_{ - \infty }^\infty {\frac{{\delta (x){\psi _m}}}{{Pe}}\,{\textrm{e}^{ - \textrm{i}\alpha x}}\,\textrm{d}\,x} = \frac{{{\psi _m}}}{{Pe\sqrt {2{\rm \pi} } }}.\end{equation}
Now that we have the transform 
$\vartheta (t,\alpha )$ in hand, 
${\theta _m}$ can be evaluated using the inverse formula as
\begin{equation}\begin{array}{ccccc} &{\theta _m} = \dfrac{1}{{\sqrt {2{\rm \pi} } }}\int_{ - \infty }^\infty {\dfrac{{{\psi _m}}}{{Pe\sqrt {2{\rm \pi} } }}\; {\rm exp} [ - {\mathrm{\mathbb{C}}_0} - \textrm{i}\alpha {\mathrm{\mathbb{C}}_1} - {\alpha ^2}\xi ]{\textrm{e}^{\textrm{i}\alpha x}}\,\textrm{d}\alpha } \\ & = \dfrac{{{\psi _m}}}{{2{\rm \pi} \,Pe}}\int_{ - \infty }^\infty {{\rm exp} [ - {\mathrm{\mathbb{C}}_0} + \textrm{i}\alpha (x - {\mathrm{\mathbb{C}}_1}) - {\alpha ^2}\xi ]\,\textrm{d}\alpha } \\ & = \dfrac{{{\psi _m}}}{{2{\rm \pi} \,Pe}}\int_{ - \infty }^\infty {{\rm exp} \left\{ { - {\mathrm{\mathbb{C}}_0} - {{\left[ {\dfrac{{i(x - {\mathrm{\mathbb{C}}_1})}}{{2{\xi^{1/2}}}} - \alpha {\xi^{1/2}}} \right]}^2} - \dfrac{{{{(x - {\mathrm{\mathbb{C}}_1})}^2}}}{{4\xi }}} \right\}\,\textrm{d}\alpha } \\ & = \dfrac{{{\psi _m}{\rm exp} \left[ { - {\mathrm{\mathbb{C}}_0} - \dfrac{{{{(x - {\mathrm{\mathbb{C}}_1})}^2}}}{{4\xi }}} \right]}}{{2{\rm \pi} \,Pe}}\int_{ - \infty }^\infty {{\rm exp} \left\{ { - {{\left[ {\dfrac{{i(x - {\mathrm{\mathbb{C}}_1})}}{{2{\xi^{1/2}}}} - \alpha {\xi^{1/2}}} \right]}^2}} \right\}\,\textrm{d}\alpha } \\ & = \dfrac{{{\psi _m}}}{{2Pe\sqrt {{\rm \pi} \xi } }}{\rm exp} \left[ { - {\mathrm{\mathbb{C}}_0} - \dfrac{{{{(x - {\mathrm{\mathbb{C}}_1})}^2}}}{{4\xi }}} \right]. \end{array}\end{equation}Due to the linearity of the mass transport equation, the solution given by (A5) can be easily extended to the situations in which the axial variation of the initial solute concentration is represented by functions other than a Dirac delta function utilizing the superposition approach. To this end, attention should be first given to the fact that the solution to a concentrated solute source of strength unity (the integral of the dimensionless mean solute concentration over the whole x axis being unity), located at the origin, is
\begin{equation}{\theta _m} = \frac{1}{{2\sqrt {{\rm \pi} \xi } }}{\rm exp} \left[ { - {\mathrm{\mathbb{C}}_0} - \frac{{{{(x - {\mathrm{\mathbb{C}}_1})}^2}}}{{4\xi }}} \right].\end{equation}
Hence, the solution corresponding to a concentrated solute source of strength 
$g(\eta )\,\textrm{d}\eta $, located at 
$x = \eta $, becomes
\begin{equation}\textrm{d}{\theta _m} = \frac{{g(\eta )\,\textrm{d}\eta }}{{2\sqrt {{\rm \pi} \xi } }}{\rm exp} \left[ { - {\mathrm{\mathbb{C}}_0} - \frac{{{{(x - \eta - {\mathrm{\mathbb{C}}_1})}^2}}}{{4\xi }}} \right].\end{equation}
Given the fact that the transport coefficients are not dependent upon the axial dependence of the initial concentration distribution, (A7) can be integrated to obtain 
${\theta _m}$ for any 
$g(\eta )$. Even though in general numerical integrations are required, however, solutions in terms of known functions are also possible. For example, the expression of 
${\theta _m}$ for the case for which the initial solute band extends uniformly in the axial direction from 
$x = 0$ to 
$x = {x_c}$ is obtained as
\begin{equation}\begin{array}{ccccc} &{\theta _m} = \int_0^{{x_c}} {\dfrac{{{\psi _m}}}{{2\sqrt {{\rm \pi} \xi } }}{\rm exp} \left[ { - {\mathrm{\mathbb{C}}_0} - \dfrac{{{{(x - \eta - {\mathrm{\mathbb{C}}_1})}^2}}}{{4\xi }}} \right]\,\textrm{d}\eta } = \dfrac{{{\psi _m}\,{\textrm{e}^{ - {\mathrm{\mathbb{C}}_0}}}}}{{2\sqrt {{\rm \pi} \xi } }}\int_0^{{x_c}} {{\rm exp} \left[ { - \dfrac{{{{(x - {\mathrm{\mathbb{C}}_1} - \eta )}^2}}}{{4\xi }}} \right]\,\textrm{d}\eta } \\ & = \dfrac{{{\psi _m}}}{2}\,{\textrm{e}^{ - {\mathrm{\mathbb{C}}_0}}}\left[ {\textrm{erf}\left( {\dfrac{{x - {\mathrm{\mathbb{C}}_1}}}{{2\sqrt \xi }}} \right) - \textrm{erf}\left( {\dfrac{{x - {\mathrm{\mathbb{C}}_1} - {x_c}}}{{2\sqrt \xi }}} \right)} \right]. \end{array}\end{equation}Appendix B. Special solutions for 
$Da = 0$
In this appendix, solutions are obtained for the situation where no reaction occurs at the channel walls, that is, when 
$Da = 0$. This situation arises when the channel is utilized for liquid-phase transportation only. In the absence of surface reaction, the boundary condition (2.8d) and its counterpart in (2.19b) are modified as
\begin{equation}{\left. {\frac{{\partial \theta }}{{\partial r}}} \right|_{r = 1}} = 0,\quad {\left. {\frac{{\partial {f_k}}}{{\partial r}}} \right|_{r = 1}} = 0.\end{equation}
From (2.28), it is obvious that the exchange coefficient 
${K_0}$ becomes zero when 
$Da = 0$. Moreover, under such circumstances, (2.14) and (2.18) are modified as
\begin{gather}{K_k} = \frac{{{\delta _{k2}}}}{{P{e^2}}} - \frac{2}{{\rm \pi} }\int_0^1 {\int_0^{\rm \pi} {ru{f_{k - 1}}\,\textrm{d}\phi \,\textrm{d}r} } ,\quad k = 1,2, \ldots ,\end{gather}
\begin{gather}\frac{{\partial {f_k}}}{{\partial t}} = \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial {f_k}}}{{\partial r}}} \right) + \frac{1}{{{r^2}}}\frac{{{\partial ^2}{f_k}}}{{\partial {\phi ^2}}} + \frac{{{f_{k - 2}}}}{{P{e^2}}} - u{f_{k - 1}} - \sum\limits_{i = 1}^k {{K_i}{f_{k - i}}} ,\quad k = 0,\; 1,\; 2,\; \ldots .\end{gather}
Note that even though 
${K_0}$ is zero for this case, we still have to determine 
${f_0}$ due to the fact that the determination of 
${K_k}$ is possible only when 
${f_{k - 1}}$ is available. It can be shown that the solution of 
${f_0}$ now is given as
\begin{align}{f_0}(t,r,\phi ) &= {A_{00}} + \mathop \sum \limits_{n = 1}^\infty {A_{0n}}{J_0}({\lambda _{0n}}r){\rm exp} ( - \lambda _{0n}^2t)\nonumber\\ &\quad + \sum\limits_{m = 1}^\infty {\sum\limits_{n = 0}^\infty {{A_{mn}}{J_m}({\lambda _{mn}}r)\cos (m\phi ){\rm exp} ( - \lambda _{mn}^2t)} } ,\end{align}
where the eigenvalues 
${\lambda _{mn}}$ satisfy the following equation:
\begin{equation}m{J_m}({\lambda _{mn}}) - {\lambda _{mn}}{J_{m + 1}}({\lambda _{mn}}) = 0.\end{equation}
Note that the reason for expanding the compact solution (2.22) in (B5) is to emphasize the fact that, unlike the general solution, 
${\lambda _{mn}} = 0$ is now one of the eigenvalues. This results in the appearance of the constant term 
${A_{00}}$ that, upon the application of the integral constraint (2.20), is determined to be 
$1$. For other 
${A_{mn}}$ coefficients, (2.25) may still be consulted. Having 
${f_0}$ in hand, 
${K_1}$ can be readily calculated:
where 
is still given by (2.38). For the SID, the solutions may simplify as
\begin{gather}{f_{0,SC}}(t,r,\phi ) = 1 + \sum\limits_{n = 1}^\infty {{A_{0n}}{J_0}({\lambda _{0n}}r){\rm exp} ( - \lambda _{0n}^2t)} ,\end{gather}
where
The solutions may be further simplified when the initial solute band is distributed uniformly over the channel cross-section. Under these circumstances, it can be shown that 
${A_{0n}} = 0$ for 
$n \ne 0$; hence, (B7) and (
) simplify to yield 
${f_{0,U}} = 1$ and 
${K_{1,U}} ={-} {u_m}$.
According to (B3), the equation governing 
${f_1}$ in the absence of surface reactions is given as
\begin{equation}\frac{{\partial {f_1}}}{{\partial t}} = \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial {f_1}}}{{\partial r}}} \right) + \frac{1}{{{r^2}}}\frac{{{\partial ^2}{f_1}}}{{\partial {\phi ^2}}} - (u + {K_1}){f_0},\end{equation}for which the following solution can be provided:
\begin{equation}{f_1}(t,r,\phi ) = \sum\limits_{i = 0}^\infty {\sum\limits_{j = 0}^\infty {[{\mathrm{\mathbb{B}}_{1,ij}}(t) + {A_{ij}}{{\mathbb C}_1}(t)]{J_i}({\lambda _{ij}}r)\cos (i\phi ){\rm exp} ( - \lambda _{ij}^2t)} } ,\end{equation}
with 
${\mathrm{\mathbb{B}}_{1,ij}}$ and 
${\mathrm{\mathbb{C}}_1}$ still being expressed by (2.31a,b). Now, the dispersion coefficient can be evaluated from (B2) as
The solutions to 
${f_1}$ for the case of a semicircular initial concentration distribution may be written as follows:
wherein
It should be noted that 
${\mathrm{\mathbb{C}}_{1,SC}}$ is obtained by using the integral constraint (2.20). Consequently, we can evaluate the associated dispersion coefficient in the following form:
The solutions are further simplified for a uniform solute band that corresponds to 
$a = 1$. For this case, one may show that
Lastly, according to (B3), the equation governing 
${f_2}$ is given as
\begin{equation}\frac{{\partial {f_2}}}{{\partial t}} = \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial {f_2}}}{{\partial r}}} \right) + \frac{1}{{{r^2}}}\frac{{{\partial ^2}{f_2}}}{{\partial {\phi ^2}}} +{-} (u + {K_1}){f_1} + (P{e^{ - 2}} - {K_2}){f_0},\end{equation}for which the following solution is obtained:
\begin{equation}{f_2}(t,r,\phi ) = \sum\limits_{s = 0}^\infty {\sum\limits_{q = 0}^\infty {[{\mathrm{\mathbb{B}}_{2,sq}}(t) + {A_{sq}}{{\mathbb C}_2}(t)]{J_s}({\lambda _{sq}}r)\cos (s\phi ){\rm exp} ( - \lambda _{sq}^2t)} } ,\end{equation}
where 
${\mathrm{\mathbb{B}}_{2,sq}}$ and 
${\mathrm{\mathbb{C}}_2}$ are given by (2.42). It is noteworthy that the expressions of the mean and local concentrations are still given by (2.53) and (2.54), respectively.
Appendix C. Special solutions for a uniform velocity field
The special solutions for a uniform velocity profile are relevant to the situations where the system enjoys an electroosmotic pumping mechanism. Electroosmotic flow fields are characterized by flat velocity profiles with sharp gradients near the wall. When the ionic concentration of the liquid is high, which is usually the case, such velocity gradients occur over extremely thin areas that constitute a very small portion of the channel cross-section (Sadeghi et al. Reference Sadeghi, Yavari, Saidi and Chakraborty2011). Hence, the velocity profile may be practically considered uniform over the whole fluidic area.
As mentioned previously, the exchange coefficient is not influenced by the velocity distribution. Accordingly, even in the presence of a uniform velocity field, 
${f_0}$ and 
${K_0}$ may still be obtained via (2.22) and (2.28). However, the velocity profile strongly influences the function 
${f_1}$ and the convection coefficient. For a uniform velocity profile, that is, 
$u = 1$, (2.29) is modified as
\begin{equation}\frac{{\partial {f_1}}}{{\partial t}} = \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial {f_1}}}{{\partial r}}} \right) + \frac{1}{{{r^2}}}\frac{{{\partial ^2}{f_1}}}{{\partial {\phi ^2}}} - (1 + {K_1}){f_0} - {K_0}{f_1}.\end{equation}It can be shown that the solution to (C1) is given as
\begin{equation}{f_1}(t,r,\phi ) = \mathrm{\mathbb{A}}(t)\sum\limits_{i = 0}^\infty {\sum\limits_{j = 0}^\infty {{A_{ij}}[ - t + {{\mathbb C}_1}(t)]{J_i}({\lambda _{ij}}r)\cos (i\phi ){\rm exp} ( - \lambda _{ij}^2t)} } .\end{equation}
The constraint (2.20) dictates the time-dependent function 
${\mathrm{\mathbb{C}}_1}$ to be t, resulting in 
${f_1} = 0$ and
\begin{equation}{K_1}(t) ={-} \frac{2}{{\rm \pi} }\int_0^1 {\int_0^{\rm \pi} {r{f_0}\,\textrm{d}\phi \,\textrm{d}r} } ={-} 1.\end{equation}
Accordingly, the equation governing 
${f_2}$ is simplified as
\begin{equation}\frac{{\partial {f_2}}}{{\partial t}} = \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial {f_2}}}{{\partial r}}} \right) + \frac{1}{{{r^2}}}\frac{{{\partial ^2}{f_2}}}{{\partial {\phi ^2}}} - {K_0}{f_2} + (P{e^{ - 2}} - {K_2}){f_0},\end{equation}to which the following solution may be obtained:
\begin{equation}{f_2}(t,r,\phi ) = \mathrm{\mathbb{A}}(t)\sum\limits_{s = 0}^\infty {\sum\limits_{q = 0}^\infty {{A_{sq}}{{\mathbb C}_2}(t){J_s}({\lambda _{sq}}r)\cos (s\phi ){\rm exp} ( - \lambda _{sq}^2t)} } .\end{equation}
Utilizing the constraint (2.20), it is straightforward to show that 
${\mathrm{\mathbb{C}}_2}(t) = 0$, leading to 
${f_2} = 0$ and 
${K_2} = P{e^{ - 2}}$. With 
${\mathrm{\mathbb{C}}_1} = t$ and 
${\mathrm{\mathbb{C}}_2} = 0$, the expression of the dimensionless mean concentration, given by (2.53), is simplified for a uniform velocity profile to yield
\begin{equation}{\theta _m}(t,x) = \frac{{{\psi _m}}}{{2\sqrt {{\rm \pi} t} }}{\rm exp} \left[ { - {\mathrm{\mathbb{C}}_0} - \frac{{{{(x - t)}^2}}}{{4P{e^{ - 2}}t}}} \right].\end{equation}Finally, the local dimensionless concentration becomes
\begin{equation}\theta (t,r,\phi ,x) = \frac{{{\psi _m}}}{{2\sqrt {{\rm \pi} t} }}{\rm exp} \left[ { - \frac{{{{(x - t)}^2}}}{{4P{e^{ - 2}}t}}} \right]\sum\limits_{m = 0}^\infty {\sum\limits_{n = 0}^\infty {{A_{mn}}{J_m}({\lambda _{mn}}r)\cos (m\phi ){\rm exp} ( - \lambda _{mn}^2t)} } .\end{equation}
