Closed-form expression for the dynamic dispersion coefficient in Hagen-Poiseuille flow

We present an exact expression for the upscaled dynamic dispersion coefficient (D) for one-dimensional transport by Hagen-Poiseuille flow which is the basis for modeling transport in porous media idealized as capillary tubes. The theoretical model is validated by comparing the breakthrough curves (BTCs) from a 1D advection-dispersion model with dynamic D to that from direct numerical solutions utilizing a 2D advection-diffusion model. Both Taylor dispersion theory and our new theory are good predictors of D at lower Peclet Number (Pe) regime, but gradually fail to capture most parts of BTCs as Pe increases. However, our model generally predicts the mixing and spreading of solutes better than Taylor’s theory since it covers all transport regimes from molecular diffusion, through anomalous transport, and to Taylor dispersion. The model accurately predicts D based on the early part of BTCs even at relatively high Pe regime (~62) where the Taylor’s theory fails. Furthermore, the model allows for calculation of the time scale that separates Fickian from non-Fickian transport. Therefore, our model can readily be used to calculate dispersion through short tubes of arbitrary radii such as the pore throats in a pore network model.


Introduction
Characterization of mass transport in Hagen-Poiseuille flow (i.e., tube model) is of biological, chemical and environmental importance.In geophysical and hydrological applications, it is typical to assume a tube or capillary model for pore throats.Thus, the tube model is the backbone of pore network models [1,2].In tube models, further mass mixing and spreading on top of molecular diffusion due to radial variation of velocity in a tube is quantified through a diffusion-type equation.In this case, a constant longitudinal dispersion coefficient (D) is used instead of the molecular diffusion coefficient after some asymptotic time scale [3].Following Taylor's [3] pioneering work, Aris [4] further extended Taylor's theoretical model by resorting to spatial moment methods.Since then, a series of papers have focused on the determination of dispersion coefficient under different initial and boundary conditions as well as transient flow fields (e.g., [5]).In most previous models, after some asymptotic time scale, dispersion resembles molecular diffusion where the classical advection-dispersion equation (ADE) governs mass transport.However, none of the previous results captures the non-Fickian transport behavior occurring at pre-asymptotic time scales.
Non-Fickian transport is present at multiple scales in geological settings [6] which could be possibly attributed to the failure of reaching the time or length threshold [7,8] necessary for enough mixing to have occurred such that transport becomes Fickian.Several methods (e.g., continuous time random walk, fractional advection-dispersion, multi-rate mass transfer models) have been applied to mathematically describe non-Fickian transport.However, one simple alternative approach is to quantify a dynamic dispersion coefficient for the ADE.For mass transport through a straight tube, previous studies using series expansion methods showed that dynamic D increases monotonically from the molecular diffusion coefficient to the value specified by Taylor's theory at either steady or transient flow [9,10].Here we applied the volume averaging method (VAM) with Green's function to derive the closed-form of dynamic D. The VAM was extensively used for upscaling properties of media from pore-scale to continuum scale [11][12][13][14].
This paper presents the theoretical development for the dynamic D with Hagen-Poiseuille flow to characterize non-Fickian and Fickian transport without time (or spatial) scale limitations when Peclet number (Pe) is reasonably high.We have recently shown the complete theoretical development for reaching a close-formed expression of dynamic D for Pouiseuille flow in between parallel plates by directly solving the 2D transport equation [7].Here, we utilized the same approach and assumptions as we did for parallel plates.The verification of dynamic D was conducted by comparing breakthrough curves from 1D macroscopic model to that from 2D direct advection-diffusion simulation with increasing Pe.The limitation of assumptions underlying the development of dynamic D was examined.

Two-dimensional advection-diffusion model
Taylor [3] investigated the D for Hagen-Poiseuille flow at asymptotic time scales.In this case, the advection-diffusion equation with appropriate boundary conditions is described as: C=C0, x=0, t>0 (4) ∂C/∂x=0, x=L, t>0 (5) where C is the concentration, x is the longitudinal direction, r is the radial direction, u is the xdirection velocity that is only a function of r (uniform in x), Dm is the molecular diffusion coefficient, t is time, R is the tube radius, L is the length of tube, and C0 (=1) is the step inlet concentration.

Macroscopic advection-dispersion model
Since we are mainly concerned on the evolution of areal mean of concentration, we decompose concentration and velocity in (1) into cross-sectional areal mean and its corresponding fluctuation components: where C and u are areal mean, which are defined as: and C′ and u′ are fluctuations about the mean.According to previous studies [3, 4], the C should also satisfy 1D advection-dispersion model (ADE): The exact expression for D in ( 8) is the goal of this paper.To solve D directly from the 2D advectiondiffusion equation, we first substitute ( 6) into (1) and note that C does not vary with r: Taking a coordinate transformation (10) and applying the chain rule ( 11) changes ( 9) into (12).
ξ=x-u t and τ=t (10) Instead of neglecting the diffusion term in the longitudinal direction as in [3], we retain it and apply the cross-sectional averaging operation to (12): . The cross-sectional averaging and noting that: reduces ( 12) to: By comparing ( 14) and ( 8) in the ordinary coordinate system, one can relate the second unknown term on the left-hand side of ( 14) to so as to quantify D in (8).Subtracting ( 14) from ( 12) leads to the transport equation for C′ : As we assumed previously [7], the longitudinal diffusion term is much smaller than the radial diffusion term when flow velocity is relative small (i.e.; pure mechanical dispersion is low).In addition, both  are much smaller than other terms in (15) and the difference between them is therefore negligible.Thus (15) simplifies to: With the given initial and no-flux boundary conditions in ( 2) and (3), the solution to C′ is obtained via a unique Green function [15]: where the Green function is defined by: ( ) J0 is the Bessel function of the first kind: Γ(z) is the gamma function, μn are positive zeros of the first-order Bessel function J1(μ)=0 in which the values of μn are shown in detail by Polyanin [15].
The second unknown term on the left-hand side of ( 14) is therefore given by:  20) and then substitute the resulting expression back into (14).With results expressed back in the ordinary coordinate system, this leads to: Therefore, the dynamic D in ( 8) is the parenthetical term in (21):

Hagen-Poiseuille flow model
A priori knowledge of u′ , i.e., the fluctuation of flow field, is crucial in evaluating dynamic D.
For a pressure-driven incompressible fluid flow through a tube under a fully-developed laminar regime, the velocity profile follows the well-known Hagen-Poiseuille flow model: with cross-sectional areal mean and corresponding fluctuation component: where p x

∂ ∂
is the pressure gradient.Substituting (24) into (22) gives the complete analytical expression for the dynamic D:

Validation of the theoretical model
To validate the theoretical model, we ran numerical experiments for a tube with length L=0.1 m and radius R=0.001 m with different pressure gradients and therefore different Pe.The Pe is defined by: The 2D advection-diffusion equation in ( 1) with initial and boundary conditions ( 2)-( 5) is solved numerically, where the steady-state flow field is known based on the Hagen-Poiseuille flow model.
The numerical experiments are conducted with COMSOL Multiphysics, a generic finite-element software.A very fine mesh with node spacing of ~ 2 × 10 -5 m is used and the mesh sensitivity has been checked.
The macroscopic models utilizing ADE with dynamic D and Taylor dispersion coefficient (i.e., appropriate forms of equation ( 8)) were also solved numerically with COMSOL Multiphysics.Breakthrough curves (BTCs) were generated both from the numerical 2D advection-diffusion model and the macroscopic models.We normalized time and concentration through the following: where Q [m 3 /s] is the outlet discharge, V [m 3 ] is the total volume of tube, and t' is the so-called "pore volume" as a dimensionless time.

Results and discussions
The model with the closed-form for dynamic D accurately predicts transport across all transport regimes, from molecular diffusion, anomalous transport and to Taylor dispersion (Figure 1).Our theory based on the 2D transport equation with few assumptions agrees well with results by Gill and Sankarasubramanian [9] acquired by series expansion methods (Figure 1).The dynamic D increases asymptotically towards the value predicted by Taylor's theory not surprisingly following the same behavior we showed earlier for parallel plates [7].Similar to [9], we also found that the dynamic D at pre-asymptotic time scale increases rapidly at τ<0.2 (τ=tDm/R 2 ), increases more gradually at 0.2≤τ<0.5, and finally approximates the maximum value predicted by Taylor's theory around τ=0.5, after which the variation of the initial concentration induced by the parabolic velocity profile has evened out.We define Τ as the asymptotic time scale delineating Fickian and non-Fickian transport behavior, where T is written in accordance toτ=0.5: In addition to quantifying the asymptotic time scale, the length scale can be determined by multiplying the asymptotic time scale by the mean velocity as was done for transport in parallel plates [7].Expectedly, our theory better predicts the pre-asymptotic transport behavior in Hagen-Poiseuille flow than Taylor's theory (Figure 2).BTCs generated by the 1D ADE either with dynamic D or Taylor D [4] agree well with BTCs from the direct 2D flow and solute transport simulations at lower Pe regimes (Figures 2a-2b), and little difference is shown between our theory and Taylor's theory in predicting D. As Pe increases (Figures 2c-2d), Taylor dispersion theory fails to capture most parts of BTCs at all longitudinal locations and expectedly overestimates the dispersion process due to its inability to capture the smaller D during preasymptotic times.At relatively higher Pe regime, Taylor dispersion theory leads to much earlier arrival of the solute due to the overestimated spreading (Figures 2c-2d).On the other hand, even at relatively higher Pe (~62), our theoretical model with the dynamic D is able to capture the first arrival of solute as shown in BTCs at all longitudinal locations, but fairly misses the middle and late parts of BTCs (Figure 2c).The errors between the 1D model with dynamic D and the 2D model are likely due to the few assumptions underlying the theoretical development; these apparently still lead to an overestimation of the dispersion process for the middle and late parts of BTCs.When Pe=~308, the assumption of constant spatial gradient of mean concentration at moving coordinate is invalid; it thus leads to large fitting errors (Figure 2d).This illustrates that the pure mechanical dispersion is dominant over molecular diffusion process, where unsteady closure treatment maybe more appropriate in the context of volume averaging method [11].This is beyond the scope of this study.The analytical model for dynamic D in Hagen-Poiseuille flow presented here is derived from solving the general transport equations directly.Although it is based on few assumptions, the model can readily replace Taylor's transport theory assumed in many tube models such as that of pore throats in pore network models [2].Due to the short lengths and times involved in transport through pore throats, it is difficult if not invalid to assume that mixing and spreading follows Taylor's theory for asymptotic dispersion in pore throats.Our model would therefore allow for more accurate modeling of dispersion through short tubes with arbitrary length and radii.

Concluding Remarks
This paper presents the analytical expression for the dynamic D for solute transport in Hagen-Poiseuille flow by directly solving the 2D transport equation with few assumptions.The model is validated by comparing breakthrough curves generated by 2D direct solute transport simulation.Both of Taylor's theory and our theory are good predictors of D at lower Pe regime, but gradually fail to capture most parts of BTCs as Pe increases.Unlike Taylor's theory which is valid only after some asymptotic time scale, the new theory predicts mixing and spreading from molecular diffusion to anomalous and then Taylor dispersion, which allows for capturing the first arrival of solute even at relatively high Pe regime where Taylor' theory fails.However, one should be cautious of using dynamic D when Pe is large (~308).Asymptotic time and length scales which separate Fickian from non-Fickian transport can also be determined using the theory.The method used here can be readily extended to other initial and boundary conditions.The new theory is useful in either verifying numerical models or describing scalar transport using tube models when Pe is reasonably high and should be used in lieu of Taylor's theory when pre-asymptotic transport is expected or prevalent.

Figure 2 .
Figure 2. Breakthrough curves (BTCs) for different Peclet Number (Pe).BTCs in a-d collected at different locations are represented by different colors: blue x/R=20, green x/R=60, fuchsia x/R =80, and red x/R=100.Dots represent direct 2D advection-diffusion simulation results, dash lines represent 1D ADE results with Taylor dispersion coefficient (D), and solid lines represent 1D ADE results with dynamic dispersion coefficient.

Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 30 May 2017 doi:10.20944/preprints201705.0211.v1
When pure mechanical dispersion is not yet dominant over molecular diffusion that the mean concentration gradient with respect to the moving coordinate (