Investigation of Drift Phenomena at the Pore Scale during Flow and Transport in Porous Media

This paper reports an analytical study conducted to investigate the behaviour of tracers undergoing creeping flow between two parallel plates in porous media. A new coupled model for the characterisation of fluid flow and transport of tracers at pore scale is formulated. Precisely, a weak-form solution of radial transport of tracers under convection–diffusion-dominated flow is established using hypergeometric functions. The velocity field associated with the radial transport is informed by the solution of the Stokes equations. Channel thickness as a function of velocities, maximum Reynolds number of each thickness as a function of maximum velocities and concentration profile for different drift and dispersion coefficients are computed and analysed. Analysis of the simulation results reveals that the dispersion coefficient appears to be a significant factor controlling the concentration distribution of the tracer at pore scale. Further analysis shows that the drift coefficient appears to influence tracer concentration distribution but only after a prolonged period. This indicates that even at pore scale, tracer drift characteristics can provide useful information about the flow and transport properties of individual pores in porous media.


Introduction
A thorough understanding of transport processes is significant for various applications in engineering, natural resources (for example waste management [1]) groundwater remediation [2] and porous media during tracer injection experiments in a core samples [3]. Tracer testing has become a highly vital tool in geothermal research, resource management and development. It is used for reinjection research and management, general hydrological studies of subsurface flow and flow rate measurement in pipeline [4]. It acts as the source of information for porous media characterisation, hydrocarbon recovery process and geothermal engineering [5].
Tracers are used to determine preferred flow direction, pore connectivity, flow velocity, dispersivity, heterogeneity of oil and gas reservoirs and residual oil saturation [6,7]. To enhance flow characterisation in porous media, an improved study of tracer transport at pore scale is important because pore-space geometry and topology are the key factors that influence flow and transport phenomena. Therefore, solute transport is significantly impacted by pore-scale heterogeneity found in natural porous media [8].
However, direct observation of pore-scale processes are not possible in the field scale [9]. Transport of tracers in porous media involves complex physical phenomena such as drift which often occur on widely varying scales for instance from a pore level to field level. Further, studies have shown that porous media in which transport takes place are heterogeneous and disordered on a microscopic scale, higher additional mixing is required by order of magnitude than the spreading due to pure molecular diffusion. Such an increase in the spreading of an initially narrow tracer pulse is due to mechanical dispersion.
Tracers introduced into a solvent, for example, water flowing through a porous system, will be exposed to dispersion in both longitudinal and transverse direction due to thermal, Fickian diffusion or drift due to variation of flow velocities within the pores [10]. Linear drift is a phenomenon experienced in subsurface flow and transport processes. It originates as a result of aquifer movement during production/injection in a distant reservoir with which there exist hydraulic connectivity and several factors. This drift of tracer particles in a liquid occurred at a pore scale due to microscopic variation of the fluid flow field [11]. Advancement in reservoir modelling which incorporates tracer partitioning between different phases and drift effect due to production as well as adsorption to grain surface requires adequate knowledge of tracer transport at pore scale [12].
Mass spreading of non-reactive tracers can be described by a combination of molecular diffusion, hydrodynamic dispersion and heterogeneous advection [13]. The hydrodynamic dispersion is usually quantified by the sum of the molecular diffusion and advection components. Many models have been developed to describe the space anomalous behaviour of dispersion coefficients (for example flow through rivulet [14], random capillaries [15]). Al Mukahal et al., 2017 provide a description of long-time Taylor-Aris dispersion of solute by using multiple scales to arrive at the expression of effective diffusivity. However, problems associated with Taylor dispersion in a capillary tube are momentous, i.e possibility of micro-macro duality and their interrelationship among scales manifest [16].
Taylor-Aris further examined a case where after a prolonged period, the concentration of solutes follows regular Gaussian distribution. To date, the majority of these studies have focused primarily on investigating flow and transport phenomena at Darcy's scale. Attention to a pore-scale investigation of flow and transport of solute under the influence of linear drift is yet to be reported in the literature.

Pore-Scale Radial Diffusion Models with Drift
Dispersion of solutes i.e tracers has wide applications in many fields of study (such as, for example, petroleum engineering, ecological studies and hydrology) and has been a huge subject of theoretical and experimental research. Figure 1 shows a schematic representation of a chemical tracer injection/extraction test in a porous media system. The drift effect fully manifested due to the microscopic variation of the fluid flow field at a later time during the extraction stage ( Figure 1f) [16]. However, a continuing challenge in mathematical and computational modelling is how to capture these mechanisms accurately and to handle them to relevant scales properly [17]. The mechanisms controlling the spreading differ in their behaviour and therefore require detailed nano-scale flow characterisation for accurate and reliable predictions of transport to be made.
To describe flow at pore scale, i.e smaller physical scales where pore and solid space are resolved separately, flow and transport are simulated using more detailed description such as Hagen-Poiseulle or Stokes flow [18]. The most basic approach to computing the spatiotemporal evolution of the concentration of tracers is provided by the advection-diffusion equation. The advection-diffusion equation is a standard approach based on Ficks law and the conservation of mass. It assumes Fickian processes and uses hydrodynamic dispersion coefficient in which the effect of solute mixing and spreading are embedded together [3]. It also uses an average linear velocity for a representative elementary volume and assuming perfectly mixed conditions within a representative elementary volume. However, at the pore scale when the velocity field is not fully sampled, incomplete mixing and spreading may exist at early times [8].
Recently, there has been increased interest in problems that involve moving interfaces at the pore scale [19]. These are problems that are characterised by sharp local gradients and mixing-controlled reactions. A large number of recent studies have explored these issues using pore-scale simulators, in which the geometry of pores-spaces and solid material is explicitly considered and processes are defined at the sub-pore [20]. Tsang and Neretnieks studied the flow and transport of solute through channels. They summarised channelling observations, concepts, and modelling of flow and transport in fractures and networks where preferential flow paths (channelling) take place [21].
Despite the numerous studies of flow and transport of tracers through geological porous media, our understanding still faces a significant challenge. One of the challenges is the observation of drift phenomena at a Darcy or field scale [5,16]. The signatures of drift phenomena affect tracer concentration distribution at field scale as reported in the literature [22]. Analytical model solutions that describe and predict tracer transports at field scale have been obtained in cases where tracer adsorption, non-uniform convection, and variable dispersion manifest [23].
Recently, investigations on closed-form solutions of radial transport of tracer in porous media at field scale have been reported in the literature (e.g., Akanji and Falade [16]). Analytical solutions at field scale have been obtained where variable dispersion coefficient and velocity are solved using Green's function method (GFM) [24] and in steady-state flow field in a horizontal aquifer caused by a constant rate injection from a well, including the mechanical dispersion and molecular diffusion terms in addition to the retardation and first-order attenuation under a Robin-type boundary condition at the well [25].
Akanji and Falade [16] studied the effect of drift on radial transport of tracers in porous media. We discussed the influence of linear drift and show that drift affects concentration profile for a typical system with non-homogeneous porosity distribution. Falade and Brigham [23] carried out a similar study where they used inverse Laplace transformation to solve radial transport equations in real-time-space.
Despite the recent advances in solving advection-diffusion equations both analytically and numerically, there is yet to be attention on the investigation of drift phenomena at the pore scale. All the previous studies mentioned above have focused on cases where convective velocity and hydrodynamic dispersion function were assumed to be constant, meaning the solute molecules i.e tracers are kinematically and dynamically indistinguishable from mobile phase molecules. However, it is important to note that this assumption is not appropriate for pore-scale flow and transport in porous media. Hence, the need for variable radial dispersion in nonuniform flows under the influence of drift is crucial for accurate description of pore-scale flow and transport Phenomena. To date, an exact solution or closed-form solution to the transport equation at pore scale has not been established for cases where hydrodynamic dispersion is radially distributed, and linear drift predominates. An investigation of drift phenomena at pore scale on fluid flow behaviour in porous media under the influence of linear drift is therefore presented in this work. The drift-velocity is encapsulated into the poiseuille flow velocity field. We establish and solve a novel pore-scale convection-dispersion equation and hydrodynamic dispersion function using the velocity obtained from Equation (1) as shown in Equations (4) and (5). This is the novelty of the work. The resulting radial transport equation is then cast in the form of the Whittaker equation using some set of transformation relations and change of variables. Concentration distribution around the source of the tracer or injected particles is computed from a weak-form numerical solution and used to analyse pore-scale tracer concentration behaviour in tracer injection/extraction experiment in a core sample during enhanced oil recovery processes where the influence of drift may obstruct the fluid flow path.

Pore-Scale Fluid Flow and Transport of Tracer Model
We consider the partial differential equations governing the flow of an incompressible Newtonian fluid at the pore scale (Akanji and Chidamoio [26,27]) thus: where the variables µ, ρ, ∇P and v are the viscosity, fluid density, reduced pressure gradient and fluid velocity vector respectively. The analytical solution to the Stokes equation is obtained for velocity fields as described in Appendix A.
To study the spatial, radial, and temporal non-uniformity of tracer-mixing flowing in a channel due to velocity variation of tracers flowing in a mobile phase water, time-dependent tracer concentration distribution is computed by solving the convection-diffusion equation expressed in terms of resident concentration in radial coordinates as (Akanji and Falade [16]): 1 r in Equation (2), D is the hydrodynamic dispersion function, v is the convective velocity at radial position r, C is the concentration of tracers, φ is the porosity of the system and γ(k r + s) is the source or sink term that accounts for phenomena such as tracer adsorption, physical loss or addition of tracer particles and reactions of tracer particles with the mobile fluid phase. The convective velocity v (m/s) in Equation (2) is decomposed into a drift velocity v d superimposed on a radial flow of injected tracer of strenght q i (m 3 /s).
A steady-state flow with a parabolic profile at every location in the channel is assumed. It implies that the transient acceleration term in Equation (1) goes to zero. The pressure gradient is only in the x-direction (unidirectional) and the plates are motionless. It is also assumed that the velocity vector and the velocity gradient tensor are off-diagonal, implying that the dot product is zero (i.e the second term in Equation (1) goes to zero), we also applied no-slip boundary conditions, since velocity is zero at the lower plate (i.e y = −h) of the channel. At the upper plate (i.e. y = +h), the velocity is zero as well, hence the velocity can be written as: This is a Poiseuille flow velocity, also called creeping flow velocity or weak inertial flow velocity (near zero Reynolds number) in a porous media. In Equation (3), dp dx , represents the average pressure gradient in the direction of the flow, µ represents the dynamic viscosity of the fluid, h is the channel thickness and v is the poiseuille fluid velocity. However, when the flow of tracer is motivated by linear drift, The convective flow velocity equals the poiseuille velocity which comprises linear drift and radial velocity components. Therefore Equation (3) represents the convective velocity that is incorporated into the transport equation. Drazin et al. 2006 [28] presented a classification and exact solutions to the Navier-Stokes equations for plane Poiseuille flow through some non-circular cross-sections.

Pore-Scale Convection-Diffusion Equation
In this section, we present only the key equations relevant to pore-scale convectiondiffusion and Taylor-Aris dispersion of solute undergoing creeping flow. Our interest is studying the influence of drift at creeping flow upon small scale transport. Thus, we simplified Equation (2) where the variables are defined as: Here, ψ is a variable which is a function of (r D ), r D is the dimensionless pore radius, k is Whittaker function, h is the channel thickness, y is the vertical distance from centre of the channel to top/bottom parallel plates shown in Figure 2, Φ is the Chromatographic response function, s is a Laplace transform parameter and 2 h is the distance between the parallel plates. The hydrodynamic dispersion D(h) in Equations (5) and (6) includes contribution from molecular diffusion and mechanical mixing. Further, the quadratic terms in pressure and pore size in Equations (4) and (5)

Mathematical Formulation of Pore-Scale Radial Transport Equation with Linear Drifts
To develop the radial transport equation where linear drift effect can be investigated, the following transformation relations are defined: However, applying the transformation to Equation (4) gives a pore-scale convectiondiffusion equation which can be expressed in the form of the Whittaker equation as: Equation (8) is a Kummers differential equation whose solution is known. Detailed transformation and simplification of Equation (4) into the Whittaker 8 can be found in Appendix C. Where: and:

Formulation of Linear Drift
For this investigation, linear drift is applied as a velocity field v d because it is coordinateindependent, carrying velocity magnitude that acts on every point within the porous system.
The tracers are carried along with the fluid in the channel while spreading due to diffusion. Consequently, the convective velocity of the fluid due to transport at pore scale is described by the Poiseuille velocity calculated from Equation (3) above which decomposes into drift velocity and a radial velocity part of an injected tracer of strength q i around the drift velocity. The convective fluid velocity in the channel system is therefore given by: where, v = p 2µ hy − y 2 , v p = α r p , and α = q i 2πhθ . Therefore, redefining the following variables from Equation (5) , β p = αD o which are described as effective dispersion constant, vorticity of the tracer particles in motion and parameter related to velocity along x-direction.

Pore-Scale Analytical Solution of Radial Transport Equation with Linear Drift
A detailed solution of the general Whittaker equation can be found in Akanji and Falade [16]. Here, we present the pore-scale analytical solution of the radial transport equation where linear drift component is accounted for. Tracer concentration can now be written as: However, constant of separation ω 2 by recasting the general advection-dispersion equation for the flow of tracers with embedded velocity from Navier-Stokes under the influence of linear drift as: where, the linear drift ratio is written as d = v d v . The linear drift ratio, d is coordinate independent, therefore, its magnitude can be used as either y-axis or, in the case of a 3D system z-axis.
Substituting C(x, y, s) = X p (x)Y p (y, s) into Equation (17) and dividing through by X(x)Y(y, s) and rearranging gives: The component in Equation (18) above is time independent, while Equation (19) is time dependent and expressed in Laplace space with Laplace parameters. Considering the following transformation parameters: β = r p cos 2 θ, λ = v p r p sin 2 θ, x p = r p cosθ and y p = r p sinθ.
Equations (18) and (19) can also be written as: Equation (20) is time independent, while Equation (21) is time dependent and expressed in Laplace space with Laplace parameter. However, Equation (17) is governed by the following boundary conditions: Considering the river-like nature of linear drift, radial flow as well as convection in every direction, the transformation from the polar ( radial) coordinate to Cartesian coordinates is given as x p = r p cosθ and y p = r p sinθ. However, for instantaneous injection, the concentration C x p , y, t = C o of the tracer in the pore is known during the tracer test study; therefore, the solution for t > 0 can be obtained by solving Equation (17) within the porous system. Several numerical inversion scheme are available for solving Equation (21) such as Euler inversion [29,30], Gaver-Steehfest algorithm [31,32] or Talbot inversion algorithm [33]. The Euler inversion algorithm was used for the numerical inverse Laplace operation in the numerical code used in this work. Flow velocity of the system is an important parameter. The flow of tracers within a carrier fluid depends on the velocity of the system and is modelled by considering a creeping flow in a heterogeneous system. Anisotropic porosity distribution was developed using the random probability density function (PDF) allowing for non-uniform velocity to be calculated. Typical normalised channel thickness as a function of velocities for various ranges of scales is shown in Figure 3 with different channel heights and Reynolds numbers.

Analysis of Results
Eleven channel thicknesses are selected to study the tracer-plate interaction effect on velocity profile. The displacement pressure gradient and viscosity are set to be 500 kg/s 2 and 0.000464 kg/ms. To illustrate the effect of linear drift on the tracer propagation profile at pore scale, flow in parallel plate of dimension 5.0 × 10 −6 m−1.3 × 10 −4 m was considered. Flow velocities in different directions as determined from the solution of Equation (3) as a function of channel thickness are computed. The maximum velocity patterns as a function of channel thickness were analysed based on Reynolds number in order to compare different flow regimes.
A graphical representation of the solution of Equation (3) for the parallel plate is shown in Figure 3. Good agreement is obtained in comparison with the asymptotic results of motion of particles between two parallel plane walls in low Reynold number Poiseuille flow [24]. Figure 3 immediately reveals the parabolic nature of the velocity fields of all the channels. The flow velocities with which the tracers enters the channel is not identical with the exit velocity of the channel. A more detailed description of the velocity field can be obtained by using 3D profiles of the velocity field in a single channel. Figure 3 shows the normalised channel thickness profile as a function of velocities computed in a spatially heterogeneous porous system. The velocity profiles indicate the magnitudes of velocity from fast flow, intermediate and slow flow. However, maximum velocity for each channel can be seen at the centre of Figure 3. In most of the previous investigations of flow and transport of tracers in porous media, Stokes velocity was not considered because of the difficulties imposed when solving analytically, coupled stokes velocity with advection-dispersion equation.
To investigate the effects of flow regime on the transport behaviour of tracers in porous media, the simulation was conducted under creeping, transition, and non-creeping flow regimes as shown in Figure 4. Velocity plots of the simulation data describing the maximum velocities of the selected channel thickness of the parallel plates versus Reynolds number are presented in Figure 4 indicating the existence of the three flow regimes. At higher Reynolds numbers, the trend of the velocity are similar to those obtained at a lower Reynolds numbers, entrance effects of tracers slightly appear at Reynolds number below 0.2 which is shown in Figure 4. However, for higher Reynolds numbers, a significant difference in magnitude of maximum velocity is observed. A well developed creeping flow regimes is established at 0 < Re < 1, transitional flow regime at Re = 1 and non-creeping flow regime at Re > 1. A set of 11 points were placed to uniformly cover each channel, at each point, the values of velocities were obtained but maximum velocity at each point were shown in Figure 4. The Reynolds number defines the flow regime and represents a ratio of inertial forces to viscous forces thus: where, ρ is the fluid density in kg/m 3 , v is the Poiseuille velocity in m/s, µ is the dynamic viscosity in kg/ms and h is the thickness of the channel in m. It can be seen from Figure 4 that, under the same flow condition, a higher Reynold number is observed for channel samples with higher velocities. Hence, the pressure gradient across the sample plate with a higher Reynolds number is more sensitive to the change of the flow condition. Therefore, flow greater than Reynolds number 1 (i.e., Re > 1) is above the creeping flow regime which is beyond the scope of this investigation. Sensitivity analysis was carried out on the developed solution by testing and evaluating the error limit associated with the separation of variable parameter (ω 2 ) for different angle θ. Simulation studies involved one hundred Further investigation on the impact of drift parameter at pore scale on separation constant, X p , Y p and angle θ , a simulation run was carried out at time interval of t = 10 d. The results is presented in Figure 7. Figure 7 also presents the combined threedimensional plots of all the plots in Figure 5a-f and the remaining data points of the other half-wing of Appendix F. Such peaks of Figure 6 would not be expected in a normal dispersion processes such as diffusion, mechanical mixing etc. But, as discussed in the literature such anomalous behaviour may be due to variation of flow velocity field inside the channel which subsequently leads to drift of tracer particles. This observation is key to understanding of drift phenomena at pore scale. By examining Figure 7 carefully, it is seen that the velocity profile increases from small data points to large data points. An analysis of Figure 7 reveals the general features of the tracer concentration distributions flowing in a heterogeneous porous media. All the data points in Figure 7 are pointing towards the upstream direction of ω 2 . Usually, in isotropic and homogeneous systems, where there is no linear drift or natural convection, the concentration of tracer ordinarily follows a cyclic pattern in nature. In this work, a system with variable porosity distribution was modelled and the corresponding Poiseuille velocity profile was used in the computation of the drift ratio. The tracer propagation profile is expected to follow a natural pattern resolved by the interaction of forces associated with the system's variables, such as the tracer injection rate.
In order to examine the effect of the drift magnitude on the tracer concentration profile at pore scale, three (3) values of the drift ratio d = 0.006, 0.03 and 0.2 were evaluated for variable dispersion coefficients of Do = 0.001389, 0.002778 and 0.004167 m 2 /s and at three selected time duration t of (t = 10, 50 and 70).
At time t = 10 days, the concentration shows a similar trend with exception of the contour plots at the left. i.e., plots (a, d and g). Previous studies did not take into account the effect of drift at the boundary. However, it is evident from contour plots a, d and g of Figure 8 that drift has manifest at the boundary which could be due to physical interaction between the channel and the fluid flowing through the channel. In general, 3D plots of all the contour plots distribution at time t = 10, 50 and 70 days can be seen in Appendix G which shows a good agreement with the tracer concentration having undergone large or field scale investigation. In order to find an appropriate mechanisms in predicting the tracer behaviour, we increase the observation time of the simulation to 50 days, drift effect was found to slightly affect the tracer concentration shown in contour plots 9. The drift effect may be as a result of velocity difference where by tracer particles drifts away from a regions of high diffusivity. 3D plots of concentration distribution of the contour plots of Figure 9 at time 50 days can be found in Appendix G.3 Figure A5.
Further investigation of drift phenomena was undertaken at time t = 70 days for d = 0.006, 0.03 and 0.2 and a variable dispersion coefficients of Do = 0.001389, 0.002778 and 0.004167 m 2 /s. The aim is to see if increase/decrease in dispersion coefficient will change the tracer distribution behaviour with time, it was found that drift effect appears more in the region of high diffusivity as shown in Figure 10. The results show that the effect of drift on concentration distribution appears at a later time for a variable dispersion coefficient.
This result is compared to tracer concentration distribution in porous media of Akanji and Falade (2019) for field scale investigation shown in Figure 11, indicating that drift effect become significant at a later time during flow and transport at pore scale.
This analysis shows that drift ratio effect manifested fully on the concentration distribution over time. However, as time increases from 10 to 70 days, the drift effect gradually manifests; suggesting that drift effect may arise over time. Thus, it is evident that concentration drifts away from the source with time [34].

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Derivation of the Analytical Solution to the Navier-Stokes Equation
Applying the following assumptions: • Steady state, the transient acceleration term i.e ρ ∂v ∂t = 0 • Unidirectional, pressure gradient is only in the X-direction and the plates are motionless. • Velocity vector and the velocity gradient tensor are off-diagonal; therefore, the dot product is zero. Hence Equation (A1) reduces to: From Equation (A2), dp dx , means the pressure gradient is only in the x-direction and is uniform. i.e., it is constant. Also µ ∂ 2 v ∂y 2 means there is a velocity vector that only has an x-component and gradient in the y-direction. Integrating Equation (A2) with respect to y we obtain: 1 2µ dp dx Applying the following boundary conditions:and y = h, v = 0. From the first boundary conditions,⇒ c 2 = 0 . From Equation (A4). therefore Equation (A4) reduces to: Applying the second boundary condition, Equation (A5) becomes: Therefore: Substitute for c 1 and c 2 into Equation (A4) to obtain:

Appendix B. Derivation of the Transport Flow Equation with Linear Drift
The Transport Equation (2) can be expanded by considering a steady-state condition thus: But γ φ = Φ is the chromatographic response function which measures the quality of separation of tracers from solvent during flow and transport.
where all variables are defined in the main text. Let: Then: and: Let: With the first order differential component dC dr defined by Equation (A12), then: Substituting Equations (A11) and (A17) in (A9) and rearranging gives: Expressing the effective fluid velocity v in dimensionless form: Thus: Let: so that: Therefore, Equation (A18) becomes: An expression of the form: is therefore obtained for the general pore-scale convection-diffusion equation in Laplace space. And it can be written in dimensionless form as: Recall that:

Appendix C. Transformation and Simplification of Pore-Scale Convection-Diffusion Equation to the General Whittaker Equation
However, applying the transformation of Section 3 to Equation (4) gives: Simplifying Equation (A48) above: Further, simplification of Equation (A52) gives: Equation (A53) can be written in the form of Whittaker equation by changing some variables defined as: where: Then: and Substituting Equations (A58) and (A60) into Equation (A53) gives: Hence Equation (A62) can now be written as: Equation (A63) is the equation that is written in the form of Whittaker Equation (8) in Section 3 of the main text.

Appendix D. Analytical Solution of Tracer Concentration and Its Weak-Form Solution
Detailed solution of the general whittaker equation can be found in Akanji and Falade [16]. Tracer concentration can now be written as: The solution to Equations (A62)-(A64) can therefore be expressed as: However, the concentration C is defined in terms of ψ as: or: (A68) so that: ) or: where a is given as: and: ξr D = 1 2µk 2 p 2 (hy − y 2 ) + 16µ 2 kΦ(k r + s)(kr + β) Having developed the transport equation for flow of tracers at pore scale, the weakform solution of Equation (17) is now presented in Laplace space by adopting the separation of variables with X-parameter (Xp) and Y-parameter (Yp), thus: The X-parameter (Xp) can be expressed as: where the components and arguments of the Tricomi Kummer function U (a, b, x) are defined thus: and Y-parameter (Yp): I is modified Bessel functions of the first kind and decays to zero rapidly with the concentration distribution of the Y (y,s) component in the negative half, depicting the positive half.