Equilibrium distributions under advection–diffusion in laminar channel flow with partially absorbing boundaries

Abstract Advective–diffusive transport in Poiseuille flow through a channel with partially absorbing walls is a classical problem with applications to a broad range of natural and engineered scenarios, ranging from solute and heat transport in porous and fractured media to absorption in biological systems and chromatography. We study this problem from the perspective of transverse distributions of surviving mass and velocity, which are a central ingredient of recent stochastic models of transport based on the sampling of local flow velocities along trajectories. We show that these distributions tend to asymptotic equilibria for large times and travel distances, and derive rigorous explicit expressions for arbitrary reaction rate. We find that the equality of flux-weighted and breakthrough distributions that holds for conservative transport breaks in the presence of reaction, and that the average velocity of the scalar plume is no longer fully characterized by the transverse distribution of flow velocities sampled at a given time.


Introduction
Transport in Poiseuille flow within a channel is a classical problem that has served to advance the understanding of the fundamental features of advection-dispersion of scalars such as solutes and temperature, going back to the classical works on asymptotic dispersion of Taylor (1953) and Aris (1956).Partially absorbing channel walls, where a transported scalar that diffuses into the walls is consumed proportionally to its local concentration, can represent a variety of physical processes, such as reactive solute consumption or sorption in porous or fractured media (Battiato et al. 2011), heat transport in pipes, rock fractures or other conduits (Lungu & Moffatt 1982), absorption in biological tissues (Davidson & Schroter 1983), and chromatography (Sankarasubramanian & Gill 1973;Balakotaiah, Chang & Smith 1995).In what follows, we adopt the language of reactive solute transport for concreteness and ease of exposition.
Traditionally, the conservative problem is treated in terms of the method of moments (Aris 1956), in which differential equations can be derived for successively higher moments along the longitudinal (mean flow) direction, averaged along the transverse direction(s).In addition to asymptotic mean velocities and dispersion coefficients, the evolution of the longitudinal plume towards normality can also be quantified using this type of approach (Chatwin 1970).The question of how surface reaction affects effective advective velocities and dispersion has received much attention (Sankarasubramanian & Gill 1973;De Gance & Johns 1978b;Lungu & Moffatt 1982;Shapiro & Brenner 1986;Balakotaiah et al. 1995;Mikelić, Devigne & Van Duijn 2006;Biswas & Sen 2007).Sankarasubramanian & Gill (1973), followed by De Gance & Johns (1978a,b), generalized the asymptotic theory of Taylor (1953) and Aris (1956) to the case of absorbing walls, and computed dispersive approximations to the longitudinal concentration field.A treatment in terms of series expansions of the longitudinal moments in Fourier space is given by Lungu & Moffatt (1982).These works derive differential equations for longitudinal moments at asymptotic times, and in particular give a detailed account of effective plume velocities and dispersion coefficients.Intuitively, solute removal preferentially depletes low-velocity areas, and these works quantify the corresponding increase in average plume velocity, as well as the reduction in dispersion due to the associated decrease in the variability of sampled velocities.Later, Barton (1984) derived a general framework for computing the time-dependent moments and asymptotic solutions based on eigenfunction expansions, extending his previous work (Barton 1983) and that of Chatwin (1970) for the conservative problem.More recently, a stochastic formulation was developed by Biswas & Sen (2007).Although we focus here on a linear, irreversible reaction, we note that when forward reaction is accompanied by the reverse reaction, such as in chromatographic applications where sorption is accompanied by desorption, the net effect at late times is a decrease in plume velocity characterized by a retardation factor associated with the average time that the solute remains sorbed to the channel walls (Paine, Carbonell & Whitaker 1983;Balakotaiah et al. 1995;Jiang et al. 2022).These two scenarios and the time scales governing the different regimes can be reconciled within a unified description (Zhang, Hesse & Wang 2017).The linear irreversible picture applies, for example, over time scales where the reverse reaction may be neglected and local surface reactant concentrations may be approximated as constant, or when thermally conducting walls may be approximated as a thermal reservoir.
Despite these significant advances, previous work has focused mainly on the transverse-averaged moments and resulting approximations for the longitudinal concentration field.The present work offers a fresh perspective on this classical problem from the point of view of transverse mass and velocity distributions.The work of Barton (1984) provides a formal means to compute the dependency of longitudinal moments on transverse position and has been used to analyse the evolution of transverse averages and breakthrough curves (Yasuda 1984;Guan & Chen 2024).Here, we provide an alternative approach to analyse the asymptotic transverse variability in plume concentrations and advective velocities directly and in detail.Transverse distributions are relevant in the current context due to their fundamental role in spatial Markov models of transport in heterogeneous systems, which describe transport in terms of a stochastic evolution of velocities along solute trajectories over fixed spatial increments (Le Borgne, Dentz & Carrera 2008;Dentz et al. 2016;Sherman et al. 2019Sherman et al. , 2021;;Puyguiraud, Gouze & Dentz 2021;Aquino & Le Borgne 2021b).Recent theoretical, numerical and experimental advances regarding flow and transport in biological systems such as lung tissue (Sznitman 2021) and blood vessel networks in the brain (Goirand, Le Borgne & Lorthois 2021) and tumours (Dewhirst & Secomb 2017) underline the importance of transport and surface exchange for processes ranging from drug delivery to oxygen uptake and metabolic waste disposal.The question of how surface reaction modifies velocity distributions across solute plumes is thus central for the development of upscaled models of reactive transport in natural subsurface systems and other heterogeneous media.
We first address the following question: does the transverse distribution of solute concentrations, normalized by surviving mass, admit an asymptotic form for large times and travel distances from a given initial condition?We find that the answer is affirmative, and we derive explicit expressions for equilibrium mass, velocity and breakthrough distributions for arbitrary reaction rate.To characterize breakthrough distributions, it turns out to be necessary to determine the asymptotic mean plume velocity.To this end, we provide an analysis that highlights the role of subleading transverse variability in the mean plume position, while reproducing previous results based on the method of moments (Lungu & Moffatt 1982;Barton 1984).We show in particular that the equality of flux-weighted (i.e.weighted by local velocities) and breakthrough (i.e.associated with flux over a control plane) distributions that has been established for conservative transport (Dentz et al. 2016;Puyguiraud et al. 2021) breaks down in the presence of reaction, and that the mean plume velocity as a function of time is no longer fully characterized by the transverse distribution of flow velocities sampled by the plume at a given time.The extent of these effects depends on medium geometry and is more pronounced for flow in a cylindrical channel than between parallel plates.
The paper is organized as follows.In § 2, we formalize the problem and present the governing equations.In § 3, we introduce the transverse distributions and related average quantities that we then compute and compare to simulations for large times and travel distances in § 4, in both two dimensions (flow between parallel plates) and three dimensions (flow in a cylindrical channel).Section 5 presents an overall discussion and conclusions.Ancillary details and derivations can be found in the appendices.

Transport in a channel with reactive walls
Consider Poiseuille flow through an infinite channel of fixed cross-section, and a solute undergoing advective-diffusive transport within the channel.We begin by introducing some notation.We denote the channel by Ω and an arbitrary cross-section by Ω ⊥ , and the channel radius by a (half-width in two dimensions).We consider a Cartesian coordinate system with the origin at the channel centre and such that x is oriented along the channel.We denote spatial positions by x = (x, x ⊥ ) ∈ Ω, with x ⊥ ∈ Ω ⊥ equal to y in two dimensions and to ( y, z) in three dimensions.We use the notation |Λ| applied to a d Λ -dimensional spatial domain Λ to denote its d Λ -dimensional measure (volume for d Λ = 3, area for d Λ = 2, and number of points for d Λ = 1).For example, the cross-sectional area is given by |Ω ⊥ | = πa 2 for a cylindrical channel in three dimensions (Ω ⊥ has dimension two, so |Ω ⊥ | is an area) and |Ω ⊥ | = 2a for two spatial dimensions (|Ω ⊥ | is a length).Finally, we denote partial derivatives with respect to a variable ξ by ∂ ξ .
Within the channel, x ∈ Ω, solute concentrations obey the advection-dispersion equation where t is time, D is the diffusion coefficient, and for stratified, fully-developed channel flow, the spatial (Eulerian) velocity profile where v M is the maximum velocity, attained at the centre of the channel (x ⊥ = 0).We define the Péclet number as where is the Eulerian mean velocity.Solute diffusion coefficients in both environmental and engineering applications are commonly in the range 10 −10 -10 −9 m 2 s −1 , whereas common flow rates and relevant length scales (such as the pore, pipe or blood vessel diameter, or the fracture aperture) vary substantially and are often not independent.For example, subsurface flows through porous and fractured media usually range from approximately Pe = 10 −2 to Pe = 10 4 , with advection-dominated conditions (Pe > 1) being common, in particular in fractures (de Marsily 1986).In chromatography, the Péclet number is usually between unity and a few hundred (Gritti & Guiochon 2013), although it should be noted that turbulent effects may be non-negligible in engineering applications such as chromatography and flow cell batteries.Typical flow rates in blood capillary networks (Ivanov, Kalinina & Levkovich 1981) are of the order of mm s −1 , with vessel diameters of the order of μm, corresponding to a typical Pe of approximately 1-10.For heat transport in fractures, the thermal diffusivity is commonly approximately 10 −7 m 2 s, and Pe ∼ 10-10 2 (de Marsily 1986).
The cross-section concentration, obtained by integrating out the longitudinal coordinate x, is (2.5) Within Ω ⊥ , the cross-section concentration obeys the diffusion equation which can be verified by integrating out x in (2.1), and applying natural boundary conditions at x = ±∞, i.e. setting the limit as |x| → ∞ of c and its derivative to zero, which is needed to ensure finite mass in an infinite domain.Here and throughout, ∇ ⊥ represents the gradient with respect to the transverse coordinates.We consider this problem subject to a given initial condition c ⊥ (x ⊥ ; 0), which is determined from the initial concentration field c(x; 0) according to (2.5).Surface reaction can be represented through the boundary conditions.Linear depletion at a constant surface rate k A [LT −1 ] corresponds to partially absorbing boundaries, i.e. the Robin boundary conditions where n ⊥ (x ⊥ ) is the outward normal at the point x ⊥ on the boundary ∂Ω ⊥ .For the conservative problem, k A = 0, this reduces to a reflecting boundary, while for instantaneous reaction, k A → ∞, we have a fully absorbing boundary, c ⊥ | ∂Ω ⊥ = 0. We define the Damköhler number Da as the ratio of the characteristic diffusion and reaction times τ D and τ R associated with the channel radius a: Because the reaction rate k A depends on both the thermodynamics of the reaction and the surface concentrations of reactants at the interface, the Damköhler number for solute transport problems can vary over a very broad range of orders of magnitude.For heat transport, we have Da = aH/(2K), where H is the surface conductance or heat transfer coefficient, and K is the thermal conductivity of the surface (Carslaw & Jaeger 1986).For example, for subsurface transport in fractured rock, we have K ∼ W m −1 K −1 (Carslaw & Jaeger 1986) and H ∼ 10-10 2 W m −2 K −1 (Heinze 2021), so that, using a ∼ 10 −6 -10 −2 m, we estimate typical values of Da in the range 10 −5 -1.

Transverse distributions
In this section, we introduce different probability density functions (p.d.f.s) of interest across the transverse direction, before proceeding to discuss the corresponding asymptotic equilibria in § 4. In Appendix A, we discuss how these quantities can be represented in terms of concentration moments in the classical Aris-Barton formulation (Barton 1984).

Surviving mass
The total mass at time t can be expressed as We use the notation Ω dx = dx dy dz for multi-dimensional integrals.For M(t) / = 1, c ⊥ (•; t) is a density, but not a p.d.f., because it is not normalized to unit integral.The p.d.f.p t (•; t) of surviving mass, which describes the transverse mass profile along the 985 A16-5 T. Aquino cross-section at time t, is given by If it exists, the corresponding equilibrium p.d.f. is 3.2.Velocity at fixed time A velocity p.d.f.describes the probability density of encountering a certain velocity in a spatial domain, sampled according to some prescribed distribution.For the flow profile (2.2), we may disregard the longitudinal direction and consider only the cross-section.We first introduce the Eulerian velocity p.d.f.p E , which is defined as the probability density associated with finding a certain velocity magnitude value at a uniformly random spatial location.We also define the flux-weighted Eulerian p.d.f.
where the Eulerian mean velocity, defined in (2.4), obeys We denote the flux-weighted mean velocity by The transverse velocity p.d.f. at fixed time, p v t (•; t), is obtained by sampling spatial locations according to the p.d.f.p t (•; t), which describes the probability density of surviving mass at a given time as a function of cross-section positions, rather than uniformly.If p ∞ t exists, then this velocity p.d.f. also admits an equilibrium given by where are the radial distances from the channel centre where the velocity magnitude equals v for the Poiseuille flow profile (2.2).Here and in similar contexts, we use the slight abuse of notation if mass is distributed uniformly along the cross-section.Further details on the Eulerian and flux-weighted Eulerian p.d.f.s, including table 2 summarizing useful quantities and a derivation of (3.6), are given in Appendix B.

The mean velocity associated with p
Finally, we introduce the fixed-time flux-weighted velocity p.d.f., corresponding to assigning a weight proportional to the local velocity to the sampling of spatial locations: where the spatial sampling p.d.f. is obtained by flux-weighting the surviving mass p.d.f.(3.3): (3.10) The relationship between these p.d.f.s and solute breakthrough over a control plane is discussed in § 3.3.The associated mean velocity is One might reasonably suppose the asymptotic mean of the fixed-time velocity p.d.f., v ∞ t (3.8), to correspond to the asymptotic mean velocity of the solute plume at large fixed times.While this is true for conservative transport (Dentz et al. 2016), it is in fact not the case for the reactive problem.To see why, first recall that the mean plume velocity v P (t) as a function of time is defined as the rate of change of the mean plume position m P (t), which is given by (3.12) Taking the time derivative gives (3.13)where α = M −1 |dM/dt| is the effective reaction rate across the solute plume (units of inverse time).Integrating the advection-dispersion equation (ADE) (2.1) in Ω and assuming that the equilibrium distribution p ∞ t exists, we find that the late-time effective rate is constant: Again using the ADE (2.1) to substitute for ∂ t c(x, t) in (3.13) and performing the appropriate integrations, we find for late times that where are the mean longitudinal plume positions given an arbitrary transverse position x ⊥ and a position x ⊥ at the channel walls, respectively.
Reactive consumption cannot lead to a maximum velocity larger than v M , the centre-channel velocity.However, according to (3.15), if m P (t) − m W (t) grew asymptotically as time t → ∞, so would the mean plume velocity v P (t).Thus m P (t) − m W (t) must be at most constant at late times, and we find from (3.15) that v P (t) is constant at late times.For the conservative problem, it is well known that the mean plume velocity approaches the Eulerian mean velocity v E .Since reaction depletes mass in low-velocity regions, the late-time plume velocity must not be smaller than v E and must therefore be non-zero, so that m P (t) must grow asymptotically.Then the fact that m P (t) − m W (t) cannot grow at late times implies m P (t) = m W (t) to leading order.On the other hand, by definition, v P (t) = dm P (t)/dt, and we conclude that to leading order v P (t) = dm W (t)/dt also.Therefore, m W (t) ∼ v ∞ P t asymptotically, with the equilibrium mean plume velocity given by (3.17) As we will show below, it turns out that v , where the equalities hold only for the conservative problem.This may also be seen using the more traditional method of moments, which has been used to compute v ∞ P (Lungu & Moffatt 1982).The present approach clarifies the physical origin of this effect: even though the global mean plume position equals the mean plume position at the wall to leading order at late times, the subleading difference of these quantities contributes at leading order to the asymptotic mean plume velocity in (3.15).
In order to compute v ∞ P below, we consider the auxiliary quantity From this definition, we have at late times Multiplying the ADE (2.1) by x/M(t), integrating out x, and rearranging terms, we obtain an evolution equation for f .At late times, it reads with the reactive boundary condition (2.7) applied to f , and an initial condition f (x ⊥ , 0) determined by the initial concentration field according to (3.18).This equation has the form of the transverse ADE (2.6) in terms of differential operators, but with an additional linear production term at rate α ∞ , and a source term given by the velocity field weighted by the surviving mass distribution.We will analyse the late-time behaviour of (3.20) in more detail separately for the 2-D and 3-D channels in what follows.In particular, we will find that where the constant-in-time mean-position discrepancy at late times m(x ⊥ ) obeys (3.22)where m W = m(x ⊥ ∈ ∂Ω ⊥ ), so that (3.15) for the mean plume velocity is satisfied.

Breakthrough over a fixed control plane
The breakthrough of mass as a function of time is a standard metric that represents the total advective mass flux passing a control plane at a fixed longitudinal position.Here, we consider breakthrough profiles, by which we mean the profile of the mass flux per area as a function of transverse position over the control plane.As before, we are interested in the existence of equilibrium profiles when normalized by surviving mass.Thus we define the breakthrough p.d.f., i.e. the breakthrough profile normalized by the total mass that crosses a control plane at position x over all times, and we seek its equilibrium form: The advective mass flux is locally proportional to the flow velocity.At large distances, the crossing times are also large, so that we may expect the equilibrium breakthrough p.d.f. to result from flux-weighting the surviving mass p.d.f., i.e. p ∞ s = p ∞ t,F ; see (3.10).This is indeed the case for the conservative problem (Dentz et al. 2016), but, surprisingly, the equality does not hold exactly for the reactive problem due to the subleading contributions m(x ⊥ ) to the mean plume position discussed in the previous subsection (see (3.21)).To see this, consider that, as already discussed, the plume is Gaussian along the longitudinal direction at late times (Chatwin 1970;Biswas & Sen 2007).We thus write the concentration field as where G(•; ξ, σ 2 ) is the Gaussian p.d.f. of mean ξ and variance σ 2 , the mean plume position is given by (3.21), and the Taylor dispersion coefficient D e is given by (2.9).We have established in (3.14) that the mass decays at a constant rate α ∞ at late times, which is equivalent to saying the late-time decay of mass is exponential at rate α ∞ .Thus the time integrals in (3.23a) are proportional to the Laplace transform over time of the Gaussian p.d.f. in (3.24) evaluated at α ∞ , which, for x ≥ m(x ⊥ ), is given by The requirement x ≥ m(x ⊥ ) is needed because the Gaussian approximation deteriorates at smaller x.For a given time t > 0, the peak of the Gaussian is located at This means that under this approximation, the peak at a position x < m would occur before t = 0, resulting in a change of sign of the square root term within the exponential in the Laplace transform.Using (3.25), we conclude that where Here, we have used the fact that μ 2 1 = 2τ D α ∞ is a constant that depends only on the geometry of the problem and on the Damköhler number, as we will see below.Note that γ ≈ Pe/μ 1 for Pe In the last approximation, we have used the fact that the term involving η turns out to be small compared to 1 for both the 2-D and 3-D problems.We note that γ ≈ v E /v ∞ P is recovered if the plume is approximated as a Dirac delta rather than a Gaussian along the longitudinal direction.As shown in what follows, the dimensionless quantity α ∞ m(x ⊥ )/v E does not depend on the Péclet number, and γ increases monotonically with Pe.The correction factor β is thus largest and, more importantly, most spatially variable as Pe → ∞.Formally, the factor γ is also largest as η → 0 for a given Pe, with η depending on Da and geometry, although we find that the value of η has negligible impact for the problems treated here, irrespective of the value of Pe.This means that while the mean plume velocity v ∞ P is required to estimate γ , a detailed treatment of dispersion is not.For the conservative problem, α ∞ = 0, so β(x ⊥ ) = 1 and p ∞ s = p ∞ t,F , as expected.We thus conclude that for the reactive problem, breakthrough is favoured where m > 0, and suppressed where m < 0, in addition to the natural flux-weighting caused by considering advective breakthrough at fixed distance.From the preceding discussion, based on (3.24), we see that this enhancement or suppression of breakthrough is due to the fact that transverse positions where the mean plume position is locally larger are associated with a larger mass at the time of crossing, because solute arrives earlier and is therefore less subject to reaction on average.We find in what follows that the differences between p ∞ s and p ∞ t,F are non-zero but small for all Da > 0 and Pe > 0 for the 2-D channel, and are more pronounced for the 3-D channel, even at moderate Da and Pe.In both the 2-D and 3-D cases, the differences between the mean v ∞ t of the fixed-time transverse velocity p.d.f. and the mean plume velocity v ∞ P , which also result from the same effect, are appreciable starting at moderate Da and Pe.It is interesting to note that as we will see, m is non-zero even for the conservative problem, although it no longer plays a role at late times because it does not interact with reaction to create the necessary mass differences.
Finally, we consider the velocity p.d.f.associated with breakthrough over a control plane at a given distance from injection.Analogously to (3.6), we have the relationship represents the mean solute velocity upon crossing a far control plane, i.e. at a longitudinal position far from the initial condition irrespective of the crossing time.Note that for the conservative problem, this velocity p.d.f. is obtained by flux-weighting the fixed-time velocity p.d.f.That is, in that case, ; see (3.9) and (3.11).

Equilibrium distributions
In this section, we determine the transverse equilibrium distributions for surviving mass, breakthrough and velocity introduced in § 3. The main asymptotic quantities discussed in this section are summarized in table 1 for ease of reference.Consider first the conservative problem.The asymptotic concentration field for large times is constant over the cross-section, so that p ∞ t (x ⊥ ) = 1/|Ω ⊥ |.Thus, in this case, the asymptotic p.d.f. of velocities at fixed time coincides with the Eulerian p.d.f., p ∞ v t = p E .Since, as already discussed, the correction factor β(x ⊥ ) in (3.26) is unity in this case, breakthrough p.d.f.s are obtained by flux-weighting.Thus we have p (Dentz et al. 2016).In what follows, we analyse the reactive problem separately for flow between flat plates (2-D channel) and flow in a cylindrical pipe (3-D channel).We compare our results to numerical simulations obtained using a standard finite-volume method in OpenFOAM (OpenCFD Ltd 2022).For the simulations, we fix a moderate Péclet number Pe = 2, and consider three values of the Damköhler number, Da = 10 −1 , 1 and 10, corresponding to slow, moderate and fast reaction compared to diffusion at the scale of the channel radius (see Appendix C for further details).

Two-dimensional channel
In this case, the diffusion equation (2.6) with boundary conditions (2.7) and initial condition c ⊥ ( y; 0) has the solution (Polyanin & Nazaikinskii 2016)  where and μ n are the positive solutions to (4.1d) Trigonometric manipulation yields two useful alternate forms of this relation: The μ n dictate the rate of decay of mode n in (4.1a) through the exponential factor.While they cannot in general be determined analytically, by ordering the solutions so that μ n < μ n+1 , it holds that the n = 1 mode is dominant for t 2τ D /(μ 2 2 − μ 2 1 ), so that for large times, The dependency of μ 1 on Da is shown in figure 1, along with the 3-D case discussed in § 4.2 for comparison.Analytical results for the low and high Damköhler limits are derived in § E.1.The onset of the asymptotic regime depends only weakly on the reaction time τ R , through the dependence of the decay mode constants μ n on the Damköhler number.Indeed, μ 2 2 − μ 2 1 grows monotonically from ≈10 to ≈20 for the 2-D case, so that the characteristic time for the onset of the asymptotic regime decreases with Da but is always of the order of a tenth of the diffusion time.It may seem surprising that the time τ R does not enter more directly, but this can be understood as follows.In the high-Da limit, the reaction is transport-limited, so that the transition to the asymptotic regime is controlled by diffusion.On the other hand, in the low Da limit, the equilibrium profile is close to homogeneous, so that little reaction is necessary to achieve it, and the diffusive homogenization time remains the most important control.Similar considerations are valid for the 3-D case, for which μ 2 2 − μ 2 1 grows monotonically from ≈15 to ≈25.

Surviving mass and breakthrough
From (3.2) and (4.2), we obtain the equilibrium surviving mass p.d.f.
3) where we have used (4.1e) and some trigonometric manipulation.
We conclude that for arbitrary reaction rate and including Da → ∞, the transverse distribution of surviving mass admits an asymptotic form that is independent of the initial condition.From (4.1a), we also conclude immediately that the late-time reaction rate is constant and given by the slowest decay mode, This result can be verified to agree with (3.14) by direct calculation using (4.3) and the trigonometric relations (4.1e).The fact that the asymptotic rate is constant is consistent with the existence of a transverse equilibrium distribution.The relationship to recent descriptions in terms of first passage times (Aquino & Le Borgne 2021a; Aquino et al.

2023) is discussed in Appendix D.
To obtain the flux-weighted surviving mass p.d.f.p ∞ t,F , we flux-weight (4.3) according to (3.10): see also (4.13) below for the mean of the fixed-time velocity p.d.f.The limiting forms of p ∞ t and p ∞ t,F for small and large Da are derived in § E.1.To obtain the true breakthrough p.d.f.p ∞ s associated with flux over a fixed control plane, we must take into account the correction factor β( y) = exp[γ α ∞ m( y)/v E ]; see (3.26).Because of the normalization, any differences between p ∞ s and p ∞ t,F result from variability in the mean-position discrepancy m( y) with respect to the cross-section position y.To compute it, we must solve (3.20) for the auxiliary quantity f ( y, t).Since the differential operators are the same as for the transverse ADE (2.6), this equation has the same propagator, except that the decay modes are modified by the production term at rate α ∞ .The solution for the 2-D channel, reflecting the presence of a source term where the φ n and μ n are the same as in (4.1a).We assume that the initial condition f ( y, 0) ∝ dx x c(x, y, 0) does not depend on y, which holds for any initial condition of constant width along the longitudinal direction.We then set f ( y, 0) = 0 without further loss of generality by choosing the longitudinal origin of the coordinate system to coincide with the initial mean plume position.Using (4.3) and (4.4), we obtain, for The first term is associated with the mean plume velocity as discussed below.From the second term and (3.19a,b) and (3.21), we conclude that the late-time mean-position discrepancy obeys where we have used the fact that p ∞ t ( y) ∝ φ 1 ( y) and the φ n are a set of orthogonal modes (i.e. the product φ n φ m integrates to zero for n / = m).We note that for an initial condition where f ( y, 0) is not constant, there is an additional contribution to m( y) given by φ 1 ( y)/(ab 1 ) dy φ 1 ( y ) f ( y , 0), which cannot be fully removed by a choice of coordinate origin.We are not aware of a closed-form solution of (4.10) for arbitrary Da.However, rather surprisingly, the remaining integration and sum can both be performed analytically in the limit of large Da to yield the simple form Unfortunately, the denominator in (3.26) must still be computed numerically.As mentioned in § 3.3, m is non-zero even for the reactive problem.In that case, proceeding similarly using the conservative solution (Polyanin & Nazaikinskii 2016), we obtain We find by numerical computation according to (4.10), using the values of μ n up to n = 6 reported in Carslaw & Jaeger (1986), that this result is approached continuously in the limit of small Da.under the approximation p ∞ s ≈ p ∞ t,F are shown in figure 3 for different Da and Pe = 2.As expected, since the effect of β is small at large Da and Pe, the agreement is good across all values of Da.

Velocity
Using (2.2), (3.8) and ( 4.3), we determine the asymptotic mean of the fixed-time velocity p.d.f. as Damköhler number, Da Damköhler number, As discussed in § 3.2, this is smaller than the mean plume velocity v ∞ P , which describes the rate of change of the plume centre of mass at late times.We find from (3.17) and (4.8) that This result is equivalent to that reported in Lungu & Moffatt (1982).Using (4.3) for p ∞ t , (4.11) for α ∞ m( y)/v E , (4.13) and (4.14) for v ∞ t /v E and v ∞ P /v E , and m W = m(±a), we find by direct calculation that (3.22) is satisfied as expected.
In line with the discussion for the breakthrough p.d.f., we approximate the true average breakthrough velocity v ∞ s by the flux-weighted mean v ∞ t,F .Proceeding similarly to the calculation of v ∞ t , using (3.11), we obtain It is interesting to note that for both the 2-D and 3-D channels, we have In this high-Da limit, there is a noticeable difference corresponding to an increase of approximately 7 % and 13 % between v ∞ t and v ∞ P in two and three dimensions, respectively.The plume velocity v ∞ P differs from the fixed-time flux-weighted velocity v ∞ t,F by only approximately 1 % and 2 %, respectively, and the fixed-time average velocity v ∞ t is similar to the Eulerian flux-weighted mean velocity v F rather than to the Eulerian mean velocity v E , differing by about 1 % in the 2-D case, and 4 % in the 3-D case.The average breakthrough velocity v ∞ s is approximately halfway 985 A16-16  vs as the flux-weighted p.d.f.p ∞ vt,F as discussed in the text.The solid lines in (a,c) show the analytical limit forms (E13) and (E15) for low Da, and (E14) and (E16) for high Da, respectively.The solid lines in (b) show the analytical solutions (4.16) and (4.17), with μ 1 computed numerically according to (4.1d).Dotted lines show the Eulerian and flux-weighted Eulerian velocity p.d.f.s p E (v) and p F (v), which equal p ∞ vt and p ∞ vs = p ∞ vt,F for the conservative problem.
as for the conservative case, irrespective of Da (see § E.2 for further details on the lowand high-Da limits).
In line with the results for the breakthrough p.d.f., we approximate the true velocity p.d.f.p ∞ v s associated with breakthrough at fixed distance by the flux-weighted p.d.f.p ∞ v t ,F .Flux-weighting (4.16) and using (4.13), we find (4.17) The low-velocity scaling is now linear for finite Da, and becomes quadratic as Da → ∞.
As before, for high velocities of Da.The flux-weighted approximation agrees well with numerical simulations, as can be seen in figure 5.

Three-dimensional channel
Similarly to the 2-D case, decay modes of order higher than the first can be neglected at late times for transport in a cylindrical channel.Furthermore, in this case, symmetry ensures that any angular variability in the initial concentration along the channel cross-section vanishes at late times.The late-time solution for transverse concentrations is thus associated with the slowest radial decay mode and is given by (Polyanin & Nazaikinskii 2016) where J ν is the Bessel function of the first kind of order ν, and μ 1 is the smallest positive root of Using the well-known recurrence relation 2α J α (x)/x = J α−1 (x) + J α+1 (x) (Abramowitz & Stegun 1964) for α = 1, 2, we find that the latter is equivalent to The behaviour of μ 1 as a function of Da is shown in figure 1(b); see § F.1 for the derivation of the asymptotic forms.Note that the nth radial mode decays exponentially at rate μ 2 n /(2τ D ) as for the 2-D case, although the values of μ n are different because they now solve (4.18b).

Surviving mass and breakthrough
From (4.18a), we obtain the surviving mass p.d.f.
Again, we conclude in particular that an equilibrium form independent of the initial condition exists for all Da, including Da → ∞.As before, the asymptotic reaction rate is given by Flux-weighting, we find see also the result for the fixed-time average velocity in (4.24) below.The limiting forms of these two p.d.f.s for Da → 0 and Da → ∞ are discussed in § F.1.To find the breakthrough p.d.f.p ∞ s , we must solve (3.20).As for the 2-D case, we assume that the initial condition f (x ⊥ , 0) ∝ dx x c(x, x ⊥ , 0) does not depend on x ⊥ , and we set f (x ⊥ , 0) = 0 without further loss of generality.Analogously to the 2-D channel, the solution for t  As before, the first term is associated with the mean plume velocity (see below), and from (3.19a,b) and (3.21) we conclude from the second term that We are not aware of a closed-form solution for either the integral or the sum in this case.
Recall that the mean-position discrepancy affects the breakthrough p.d.f.p ∞ s through the factor β(x ⊥ ) in (3.26).The β profiles for different values of the Péclet number at infinite Da are shown in figure 6(a), along with a comparison between the resulting p ∞ s and the flux-weighted p.d.f.p ∞ t,F in figure 6(b).Although the difference to the flux-weighted p.d.f. is not dramatic, it is substantially more pronounced than in the 2-D case, and clearly discernible at the channel centre.At Pe = 2, the results are similar to those for Pe → ∞, which explains the good agreement found between simulations at Pe = 2 and theory for Pe → ∞ for the mean breakthrough velocity in figure 4(b).The results for the surviving mass and breakthrough p.d.f.s are shown in figure 7 for different Da at Pe = 2.As expected, the differences between p ∞ s and p ∞ t,F become more pronounced with increasing Da, but are already visible at Da = 10 −1 .

Velocity
Proceeding as for the 2-D case, but integrating (3.8) in polar coordinates and consulting table 2 for the 3-D case, we obtain for the asymptotic mean of the fixed-time velocity p.d.f.:   4.21), with μ 1 computed numerically according to (4.18b).In all plots, the solid lines representing p ∞ s are computed using (3.26) and (4.23), using terms up to n = 6 in the sum according to the values of μ n reported in Carslaw & Jaeger (1986).Using (4.22), arguments similar to those for the 2-D case lead to the asymptotic plume velocity This result is equivalent to that reported in Lungu & Moffatt (1982).For the flux-weighted mean velocity, we find The different mean velocities are shown as a function of Da in figure 4(b); see also the associated discussion in § § 4.1.2and F.2 for the limiting forms at low and high Da.
We now turn to the asymptotic velocity p.d.f.s, again proceeding as for the one-dimensional case.First, according to (3.6), the fixed-time p.d.f. is given by Note that in this case, the Eulerian p.d.f.p Finally, for the flux-weighted p.d.f., we have

Discussion and conclusions
We have provided a detailed analysis of the asymptotic transverse distributions of surviving mass, velocity and breakthrough fluxes for advective-diffusive transport in Poiseuille flow between two flat plates and in a cylindrical channel.Just as for the conservative case, asymptotic forms exist, even in the limit of instantaneous reaction.In addition to qualitatively intuitive changes associated with the depletion of low-velocity regions due to reaction at the channel walls, we have also found that the equality between flux-weighted and breakthrough quantities, as well as between the mean plume velocity and the average of the transverse velocity distribution at fixed time, are modified by the presence of reaction.These effects result from variability in the mean plume position over the channel cross-section, which, despite being of subleading order in time, contributes at leading order to the mean plume velocity and results in non-vanishing asymptotic effects.
As we have seen, the effect of this mechanism on breakthrough distributions is noticeable but mild, especially in the two-dimensional case.The effect on the mean plume velocity is more substantial in both cases.The simplicity of the set-up used here in terms of the geometry, flow profile and chemical kinetics allows most of the calculations to be carried out analytically.This elucidates the underlying mechanisms and also provides analytical building blocks that can be used to develop models that are applicable to more complex scenarios.In this context, the results raise a number of questions that remain open.Does the ordering under fast reaction of the different mean velocities discussed here hold for arbitrary geometries?What geometric features lead to larger or smaller differences between the mean plume velocity and the average of the transverse velocity distribution, and between breakthrough and flux-weighted quantities?To what extent do these results affect the formal structure and predictions of stochastic models of transport based on velocity distributions, and how can they be adapted to account for more complex geometry, flow heterogeneity and rheology, as encountered, for example, in natural porous and fractured rock formations or in biological systems?Answering these questions, in particular for the single channel and for ensembles of connected channels, will be the subject of future work.

Appendix C. Numerical simulations
Numerical simulations of the reactive transport problem used in the text were performed using a standard second-order finite-volume solver in OpenFOAM (OpenCFD Ltd 2022).The mean Eulerian velocity was set to v E = 1, the domain radius (half-width in two dimensions) to a = 1, and the diffusion coefficient to D = 1/2.A flux-weighted pulse injection one discretization cell wide near x = x 0 = 10 was used in all simulations.The channel inlet was placed at x = 0, and the outlet at x = 100.The normalized equilibrium results presented in the text are not affected by these parameter and initial condition choices, except for the value of Péclet number Pe = 2 (defined according to (2.3)) in the case of breakthrough quantities.The boundary conditions at the inlet and outlet of the channel were absorbing, although care was taken to verify that no mass left the domain by the inlet or outlet within the times of interest.The Robin boundary condition (2.7) was imposed at the channel walls.The simulations were run up to time t = 100.Damköhler numbers Da = 10 −1 , 1 and 10, defined according to (2.8a-c), were simulated.
For the 2-D channel, we employed a regular grid discretization with 1000 cells in the x direction and 20 cells in the y direction.For the 3-D channel, we used a regular discretization in a square central region of side 1/2 to avoid numerical issues at the centre of the channel, as shown in figure 9.The discretization along the x direction was again regular using 1000 cells.Along each of the y and z directions, 30 cells were used: 10 in the square central region, and 10 on each side.The outer regions employed an expansion ratio of 25 %, such that cells decreased in size towards the channel wall from the centre.
The equilibrium profiles were evaluated at times τ and longitudinal positions x 0 + v E τ .For the 2-D case, τ = 30 for all values of Da.For the 3-D case, τ = 20 for Da = 10 −1 , and τ = 5 for Da = 1 and Da = 10.These smaller values of τ in the latter cases were used to avoid numerical issues due to very small concentration values, which affect the computation of breakthrough fluxes at large distances.In addition, in these two cases, the initial mass was multiplied by a large value to avoid errors associated with roundoff procedures employed by OpenFOAM to prevent overflow on division by small numbers.This does not affect the results after normalizing by mass because the problem is linear.
Surviving mass p.d.f.s p t are computed by integrating concentrations numerically over the longitudinal direction at the time of interest, using the underlying grid discretization.The breakthrough profiles p s are computed by integrating concentrations multiplied by on the first passage time picture provides good approximations and generalizes well to more complex geometries (Aquino & Le Borgne 2021a;Aquino et al. 2023).This also shows that since for Da → ∞ we have α ∞ independent of the initial condition, the actual late-time tail of the first passage time distribution is independent of the initial mass distribution and is instead fully determined by the diffusion coefficient and the geometry of the problem.

Appendix E. Damköhler expansions for the 2-D channel
In this appendix, we derive limiting forms for low and high Damöhler number Da of some of the equilibrium quantities discussed in the main text for the 2-D channel.

E.1. Surviving mass
First, consider the slowest concentration decay mode μ 1 .For Da 1, we can Taylor expand the second form in (4.1e) around μ = 0, to obtain In the opposite limit of fast reaction, Da 1, we expand (4.1d) around μ = π/2 and 1/Da = 0 to obtain In particular, we recover v ∞ t = v E for the conservative problem, as expected.Note that in order to obtain this result to second order in Da with the expansion (E1) for μ 1 , it is necessary to use the first form in (4.13); for the second form, a higher-order expansion of μ 1 in Da would be necessary.The second-order accuracy is necessary to achieve first-order accuracy in v ∞ t,F below.For Da 1, we have For the mean plume velocity (4.14), we find, for Da 1, ) and, for Da For Da 1, using (E1) and (E7), (4.15) for the flux-weighted average velocity reduces to Again as expected, we recover v ∞ t,F = v F in the limit.For Da 1, Proceeding similarly, the Da 1 limit of (4.16) for the equilibrium fixed-time velocity p.d.f. is and for Da 1 we have (E14) For the equilibrium flux-weighted velocity p.d.f.(4.17), the low-Da limit is and for the high-Da limit, we obtain Appendix F. Damköhler expansions for the 3-D channel Here, we derive limiting forms for low and high Damöhler number Da of some of the equilibrium quantities discussed in the main text for the 3-D channel.
F.1.Surviving mass Expanding (4.18b) to sixth order in μ, solving for μ, and expanding for small Da, we find for the slowest concentration decay mode: For Da → ∞, we have J 0 (μ) = 0, so that μ 1 = j 0,1 ≈ 2.4048, the smallest positive root of J 0 .Expanding μ to first order around this value in (4.18b), we obtain the high-Da expansion Using (4.19) and (F1), we find the low-Da expansion for the surviving mass p.d.f.: Note that, as expected, in the conservative limit we recover a uniform transverse p.d.f.p ∞ t (x ⊥ ) = ⊥ | = 1/(πa 2 ).Based on (4.19) and (F2), we obtain the high-Da expansion For the plume velocity (4.25), we find for Da 1, and for Da 1, Note that these approximations agree with those reported in Sankarasubramanian & Gill (1973) and Barton (1984) to the relevant order in Da.
For Da 1, we obtain the flux-weighted average velocity (4.26), and, for Da For the equilibrium fixed-time velocity p.d.f.(4.27), we have for low Da, and for high Da, (F14) For the equilibrium flux-weighted velocity p.d.f.(4.28), the low-Da limit is and the high-Da limit is 16 J 1 ( j 0,1 )

Figure 1 .
Figure 1.Slowest decay mode constant μ 1 as a function of Da for the 2-D and 3-D channels.The solid lines are obtained by finding the smallest solution to (4.1d) and (4.18b) numerically.The analytical limiting forms (E1) and (F1) for small Da, and (E2) and (F2) for large Da, are shown as dashed lines for Da ≤ 1 and Da ≥ 1, respectively.

Figure 2 Figure 2 .Figure 3 .
Figure 2. Effect of the mean-position discrepancy m( y) on the breakthrough p.d.f.p ∞ s for the 2-D channel.The mean-position discrepancy generates differences between the flux-weighted p.d.f.p ∞ t,F and the true breakthrough p.d.f.p ∞ s through the factor β( y) = exp[−γ α ∞ m( y)/v E ] in (3.26).(a) Plot of γ as a function of Péclet number Pe in the limit of infinite Da, with the low-and high-Pe behaviours shown as dashed lines.(b) Comparison of p ∞ s ((3.26) and (4.10)) and p ∞ t,F (E6) in the limit of large Da and Pe, showing that differences are minor for the 2-D channel; the factor β( y), which is below 1 near the channel walls and above 1 near the channel centre but always close to unity, is shown as a dashed line.
.15) These results are shown as a function of Da in figure 4(a) (see § E.2 for the derivation of the high-and low-Da forms).The 3-D case, discussed in detail in § 4.2, is shown in figure 4(b) for comparison.

Figure 5 .
Figure 5. Equilibrium velocity p.d.f. in time p ∞ vt and breakthrough velocity p.d.f.p ∞ vs in the 2-D channel, for (a) Da = 10 −1 , (b) Da = 1, and (c) Da = 10.Symbols show the results of numerical simulations.The theoretical predictions approximate p ∞vs as the flux-weighted p.d.f.p ∞ vt,F as discussed in the text.The solid lines in (a,c) show the analytical limit forms (E13) and (E15) for low Da, and (E14) and (E16) for high Da, respectively.The solid lines in (b) show the analytical solutions (4.16) and (4.17), with μ 1 computed numerically according to (4.1d).Dotted lines show the Eulerian and flux-weighted Eulerian velocity p.d.f.s p E (v) and p F (v), which equal p ∞ vt and p ∞ vs = p ∞ vt,F for the conservative problem.

Figure 6 .
Figure 6.Effect of the mean-position discrepancy m(x ⊥ ) on the breakthrough p.d.f.p ∞ s for the 3-D channel.The mean-position discrepancy generates differences between the flux-weighted p.d.f.p ∞ t,F and the true breakthrough p.d.f.p ∞ s through the factor β(x ⊥ ) in (3.26).(a) Behaviour of β in the limit of infinite Da, for different values of Péclet number Pe.(b) Comparison of the resulting p ∞ s ((3.26) and (4.23), using terms up to n = 6 in the sum according to the values of μ n reported in Carslaw & Jaeger 1986) to p ∞ t,F (F6).

Figure 7 .
Figure 7. Equilibrium surviving mass and breakthrough p.d.f.s p ∞ t and p ∞ s , and flux-weighted mass p.d.f.p ∞ t,F , in the 3-D channel, for (a) Da = 10 −1 , (b) Da = 1, and (c) Da = 10.Symbols show the results of numerical simulations.The solid lines in (a,c) for p ∞ t and p ∞ t,F show the analytical limiting forms (F3) and (F5) for low Da, and (F4) and (F6) for high Da.The corresponding solid lines in (b) show the analytical solutions (4.19) and (4.21), with μ 1 computed numerically according to (4.18b).In all plots, the solid lines representing p ∞ s are computed using (3.26) and (4.23), using terms up to n = 6 in the sum according to the values of μ n reported inCarslaw & Jaeger (1986).
.28) These velocity p.d.f.s are illustrated in figure8(see § F.2 for low-and high-Da forms) along with the breakthrough velocity p.d.f.p ∞ v s from (3.27) that takes into account the mean-position discrepancy (4.23).The differences between p ∞ v s and p ∞ v t ,F are in line with the discussion of the breakthrough p.d.f.p ∞ s ; see figure 7.

Figure 8 .
Figure 8. Equilibrium velocity p.d.f. at fixed time p ∞ vt , corresponding flux-weighted p.d.f.p ∞ vt,F , and breakthrough velocity p.d.f.p ∞ vs in the 3-D channel, for (a) Da = 10 −1 , (b) Da = 1, and (c) Da = 10.Symbols show the results of numerical simulations.The solid lines associated with p ∞ vt and p ∞ vt,F in (a,c) show the analytical limit forms (F13) and (F15) for low Da, and (F14) and (F16) for high Da.The corresponding solid lines in (b) show the analytical solutions (4.27) and (4.28), with μ 1 computed numerically according to (4.18b).The solid lines associated with p ∞ vs use (3.27) in all plots, with p s computed as in figure 7. Dotted lines show the Eulerian and flux-weighted Eulerian p.d.f.s, which equal p ∞ vt and p ∞ vs = p ∞ vt,F for the conservative problem.

Table 1 .
.28) Main equilibrium quantities discussed in the text, along with important equations.

)
Substituting (E1) in (4.3) and Taylor expanding for Da 1 yields for the surviving mass p.d.f.Velocity For Da 1, we obtain the average of the fixed-time velocity p.d.f.(4.13),