Closed-Form Solution of Radial Transport of Tracers in Porous Media Influenced by Linear Drift

Abstract

A new closed-form analytical solution to the radial transport of tracers in porous media under the influence of linear drift is presented. Specifically, the transport of tracers under convection–diffusion-dominated flow is considered. First, the radial transport equation was cast in the form of the Whittaker equation by defining a set of transformation relations. Then, linear drift was incorporated by considering a coordinate-independent scalar velocity field within the porous medium. A special case of low-intensity tracer injection where molecular diffusion controls tracer propagation but convection with linear velocity drift plays a significant role was presented and solved in Laplace space. Furthermore, a weak-form numerical solution of the nonlinear problem was obtained and used to analyse tracer concentration behaviour in a porous medium, where drift effects predominate and influence the flow pattern. Application in enhanced oil recovery (EOR) processes where linear drift may interfere with the flow path was also evaluated within the solution to obtain concentration profiles for different injection models. The results of the analyses indicated that the effect of linear drift on the tracer concentration profile is dependent on system heterogeneity and progressively becomes more pronounced at later times. This new solution demonstrates the necessity to consider the impact of drift on the transport of tracers, as arrival times may be significantly influenced by drift intensity.

1. Introduction

The study of the transport of tracers has become an essential technique for porous media characterisation, particularly in enhanced oil recovery (EOR) in hydrocarbon reservoirs (e.g., Baldwin [1]), hydrology (e.g., Rubin and James [2]), nuclear (e.g., Moreno et al. [3] and Herbert et al. [4]), drug transport in blood vessels (e.g., Mabuza et al. [5]) and geothermal engineering (e.g., Vetter and Zinnow [6]). Multiple processes and mechanisms are usually involved in the chemical interaction of the constituent components when the tracer is being transported through a porous medium. Two major processes involved in the transport phenomenon include convection and hydrodynamic dispersion. The convection process involves bulk movement of fluids, while hydrodynamic dispersion describes the dual actions of molecular dispersion and shear or mechanical mixing process. These complementary transport processes are adequately captured by the well-known convection–dispersion–diffusion equations with or without chemical reactions (e.g., Tomich et al. [7], Bear [8], Zhou and Zhan [9]).

Surfactant or biosurfactant partitioning and transport in the oil phase during enhanced oil recovery processes is usually neglected, because, it requires the solution of a system of nonlinear coupled partial differential equations whose solution is numerically challenging [10]. These diffusion equations are based on linear or one-dimensional geometry due to the relative ease with which such equations can be solved analytically. Recently, several authors have studied hydrodynamic transport in porous media using the random walk method (see a review paper by Noetinger et al. [11] and references cited therein). Approximate solutions have also been presented in modelling of radial geometry under conditions of shear mixing, albeit approximate in nature [12]. Exact analytical solutions have been obtained in cases where convective velocity and hydrodynamic dispersion functions were assumed constant (e.g., Carslaw and Jaeger [13]) and in porous media where tracer adsorption, non-uniform convection and variable dispersion manifest (e.g., Falade and Brigham [14]).

Attinger and Abdulle [15] studied the effective drift of transport problems in heterogeneous compressible flows. They discussed the impact of a mean drift and showed that static compressible flow with mean drift can produce a heterogeneity-driven large-scale drift or ballistic transport. A similar study was carried out by Vergassola and Avellaneda [16], where it was demonstrated that for static compressible flow without mean drift, there is no impact on the large-scale drift. The calculation of the effective ballistic velocity $V_b$ was reduced to the solution of one auxiliary equation. They derived an analytic expression for $V_b$ for some special instances where flow depends on a single coordinate, random with short correlation times and slightly compressible cellular flow. Transport will be depleted due of the trapping for arbitrary time-independent potential flow and for time-dependent potential flow or generic compressible flow, transport will be enhanced or depleted depending on the velocity field. Vergassola and Avellaneda [16] also discovered that trapping due to flow compressibility may enhance particle spreading, leading to ballistic transport that is very efficient.

In field applications, particularly during EOR involving chemical injections such as surfactants, alkali or polymer, fluid migration in an active or partially active aquifer formation may lead to displacement of the injected chemical during the shut-in period. Linear drift may also occur as a result of interference by the production/injection well, which is hydraulically connected to the formation of interest. Investigation conducted by Tomich et al. [7] indicated that in a single-well test involving the injection of ethyl acetate, fluid migration in the formation due to a reservoir water drive might lead to displacement of the tracer bank during the shut-in period. The injection of the ethyl acetate was followed by the injection of a water bank, allowing the system to hydrolyse during the chemical reaction to form ethanol, a secondary tracer. The difference in magnitude of the velocity of arrival of the two tracers was used in estimating the residual oil saturation.

The occurrence of linear drift may lead to the flow path being rerouted, leading to inaccurate and inconclusive tests with ultimate financial implications. Moench and Ogata [17] applied Laplace transform as described by Stehfest [18] to solve the dispersion in a radial flow in a porous medium. The resulting Airy function was computed using the series representation for $\left|z\right|>1$ [19]. Mashayekhizadeh et al. [20] applied Fourier series methods to numerically solve the Laplace transform of a pressure distribution equation for radial flow in porous media. Other authors (e.g., De-Hoog et al. [21], Dubner and Abate [22], Zakian [23] and Schapery [24] and Brzeziński and Ostalczyk [25]) have proposed improved techniques for numerical inversion of Laplace transforms, typically by accelerating the convergence of the Fourier series.

The occurrence of advection–dispersion with the influence of drift is vast and may occur in petroleum reservoirs with underlying aquifer, CO${}_2$–EOR processes and in contaminant hydrology. Understanding flow and transport behaviour in porous media where drift may occur is important in radioactive waste management due to the possible longevity of radionuclide materials and the possibility of being rerouted to the surface environment during transport processes. Despite the extensive research in this field, particularly in solving the advection–dispersion equation (ADE) both analytically and numerically, there is yet to be a consideration for the closed-form solution of the ADE in systems where the effect of linear drift may predominate.

A Fickian solution (Fick [26]) to the the convection–diffusion equation can be easily obtained for the simple cases where velocity and hydrodynamic dispersion are constant and the reaction term is either zero or first order in concentration (e.g., Falade and Brigham [14] and Skellam [27]). However, in cases where hydrodynamic dispersion is radially distributed and linear drift predominates, an exact analytical solution to the transport equation has not been reported in the literature. In this work, a closed-form solution of the transport of tracers in porous media under the influence of linear drift is presented. First, the radial transport equation is cast in the form of the Whittaker equation [28] by defining a set of transformation relations and a change of variables. Linear drift is incorporated by considering a coordinate-independent scalar velocity field within the porous medium. A special case of low intensity tracer injection where molecular diffusion controls tracer propagation but convection with linear velocity drift plays a significant role is presented and solved in Laplace space. Second, the concentration distribution around the source of tracer injection is solved analytically in radial coordinate and the obtained result transformed to the equivalent Cartesian coordinates system. A weak-form numerical solution is then obtained and used to analyse tracer concentration behaviour in enhanced oil recovery (EOR) processes where linear drift effect may interfere with the fluid flow path.

2. Radial Diffusion Models with Drift

Figure 1 shows a schematic representation of a chemical tracer injection in a single-well test, indicating (a) the injection of a chemical, (b) the reaction between the injected chemical and the injected water bank and (c) the production stage without the influence of drift. Figure 1d–f shows the same process as highlighted in Figure 1a–c, but, with underlying aquifer causing a noticeable drift effect during the production stage (f) (e.g., Tomich et al. [7]).

The transport of tracers in a constant flow of carrier fluid flowing in a porous medium governed by the convection–diffusion equation expressed in terms of resident concentration in radial coordinates can be written as (see, for instance, Falade and Brigham [14]):

$$\begin{array}{c}\hfill \frac{1}{r}\frac{\partial }{\partial r}\left(r\varphi D\frac{\partial C}{\partial r}\right)-\frac{1}{r}\frac{\partial }{\partial r}\left(r\varphi vC\right)-\gamma \left({\kappa }_r+s\right)C=\frac{\partial \varphi C}{\partial t},\end{array}$$

(1)

$$\begin{array}{ccc}\hfill v& =& \left(v_x\cos \theta +\frac{q_i}{2\pi r\varphi h}\right)\hfill \end{array}$$

(2)

$$\begin{array}{cc}=& \left(v_d+\frac{\alpha }{r}\right),\hfill \end{array}$$

(3)

where the composite velocity v (m/s) consists of the linear flow velocity $v_d$ (m/s) superimposed on a radial flow of the injected tracer of strength $q_i$ (m${}^3$/s). When the flow of the tracer is influenced by linear drift, the linear flow velocity will consist of both radial and tangential velocity components:

$$\begin{array}{ccc}\hfill v_d& =& v_r+v_t\hfill \end{array}$$

(4)

$$\begin{array}{cc}=& \frac{dr}{dt}{\widehat{e}}_r+r\frac{d\theta }{dt}{\widehat{e}}_{\theta }.\hfill \end{array}$$

(5)

For low intensity tracer injection, the tangential velocity is negligibly small compared to the radial velocity. Other variables $\alpha ,\gamma$ and D are defined as follows:

$$\alpha =\frac{q_i}{2\pi r\varphi h},$$

(6)

$$\gamma =\theta \left(\frac{S_{mp}+{\kappa }_l\left(1-S_{mp}\right)}{S_{mp}}\right),$$

(7)

and D is the flow hydrodynamic dispersion (m${}^2$/s), s is Laplace parameter, $S_m$ is the mobile fluid phase saturation, $q_i$ is the tracer injection rate (m${}^3$/s), r is radial distance (m), $\varphi$ is porosity and h is porous media thickness (m). Hydrodynamic dispersion D is generally believed to be made up of two components—molecular diffusion and shear mixing—which can be expressed as:

$$D=D_m+D_o\left|v_{m{r}_D}\right|.$$

(8)

In Equation (8), $D_m$ is the molecular diffusion constant (m${}^2$/s) and $D_o$ is the shear mixing constant (m).

The dimensionless form of the general convection–diffusion equation can be written in Laplace space as (see Appendix A for the derivation of the transport equation under the influence of linear drift):

$$\frac{d^2\mathsf{\Psi }}{d{r}_D^2}-\left\{\frac{{\alpha }^2{\left(\omega r_D+1\right)}^2-{\kappa }^2}{4{(\kappa r_D+\beta )}^2}+\frac{\alpha \omega }{2(\kappa r_D+{\beta }_w)}+\frac{\varphi ({\kappa }_r+s)r_D}{(\kappa r_D+\beta )}\right\}\mathsf{\Psi }=0,$$

(9)

where the variables are redefined as:

$$\begin{array}{ccc}\hfill D\left(r_{Dk}\right)& =& D_m+v_m{D}_o\left(\omega +\frac{1}{r_D}\right),\hfill \end{array}$$

(10)

$$\begin{array}{cc}=& \kappa +\frac{\beta }{r_D},\hfill \end{array}$$

(11)

and

$$\begin{array}{ccc}\hfill \kappa & =& D_m+\beta \omega \hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill v_w& =& \frac{\alpha }{r_w},\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill \beta & =& \frac{{\beta }_w\left(\alpha D_o\right)}{r_w},\hfill \end{array}$$

(14)

$$\begin{array}{cc}=& v_w{D}_o,\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill \gamma & =& \theta \left(\frac{S_m+{\kappa }_l\left(1-S_m\right)}{S_m}\right)r_w^2.\hfill \end{array}$$

(16)

3. Mathematical Formulation of the Radial Transport Equation with Linear Drift

In order to establish the radial transport equation where linear drift effect can be incorporated, the following transformation relations are defined:

$$\begin{array}{c}\hfill \eta =(\kappa r_D+\beta ),\end{array}$$

(17)

$$\begin{array}{c}\hfill r_D=(\eta -\beta )\frac{1}{\kappa },\end{array}$$

(18)

$$\begin{array}{c}\hfill \frac{d\eta }{d{r}_D}=\kappa ,\end{array}$$

(19)

$$\begin{array}{c}\hfill \frac{d^2\mathsf{\Psi }}{d{r}_D^2}={\kappa }^2\frac{d^2\mathsf{\Psi }}{d{\eta }^2},\end{array}$$

(20)

and applying the transformation to Equation (9):

$$\begin{array}{c}\hfill {\kappa }^2\frac{d^2\mathsf{\Psi }}{d{\eta }^2}-\left\{\stackrel{\stackrel{\phantom{(}i}{\downarrow }}{\frac{{\alpha }^2{\left[\frac{\omega }{\kappa }(\eta -\beta )+1\right]}^2-{\kappa }^2}{4{\eta }^2}}+\frac{\alpha \omega }{2\eta }+\stackrel{\stackrel{\phantom{(}ii}{\downarrow }}{\frac{\varphi ({\kappa }_r+s)\frac{1}{\kappa }(\eta -\beta )}{\eta }}\right\}\mathsf{\Psi }=0.\end{array}$$

(21)

Note:

$$\begin{array}{ccc}\hfill {\left[\frac{\omega }{\kappa }(\eta -\beta )+1\right]}^2& =& \frac{{\omega }^2}{{\kappa }^2}{(\eta -\beta )}^2+\frac{2\omega }{\kappa }(\eta -\beta )+1\hfill \\ & =& \frac{{\omega }^2}{{\kappa }^2}{\eta }^2-\frac{2\omega }{\kappa }\left(\frac{\omega \beta }{\kappa }-1\right)\eta +{\left(\frac{\omega \beta }{\kappa }-1\right)}^2.\hfill \end{array}$$

(22)

Therefore, the term highlighted as (i) in Equation (21) can be written out by considering the relational expression (Equation (22)), thus:

$$\begin{array}{c}\frac{{\alpha }^2{\left[\frac{\omega }{\kappa }(\eta -\beta )+1\right]}^2-{\kappa }^2}{4{\eta }^2}=\left[\frac{{\alpha }^2{\omega }^2}{{\kappa }^2}{\eta }^2-\frac{2{\alpha }^2\omega }{\kappa }\left(\frac{\omega \beta }{\kappa }-1\right)\eta +{\alpha }^2{\left(\frac{\omega \beta }{\kappa }-1\right)}^2-{\kappa }^2\right]\frac{1}{4{\eta }^2}\hfill \\ \hfill =\frac{{\alpha }^2{\omega }^2}{4{\kappa }^2}-\frac{{\alpha }^2\omega }{2\kappa }\left(\frac{\omega \beta }{\kappa }-1\right)\frac{1}{\eta }+\frac{1}{4{\eta }^2}\left[{\alpha }^2{\left(\frac{\omega \beta }{\kappa }-1\right)}^2-{\kappa }^2\right].\end{array}$$

(23)

Similarly, the term highlighted as (ii) in Equation (21) can be rewritten thus:

$$\begin{array}{c}\hfill \frac{\varphi ({\kappa }_r+s)(\eta -\beta )\frac{1}{\kappa }}{\eta }=\frac{\varphi ({\kappa }_r+s)}{\kappa }-\frac{\beta \varphi ({\kappa }_r+s)}{\kappa \eta }.\end{array}$$

(24)

Hence, Equation (21) can now be written as:

$$\begin{array}{c}{\kappa }^2\frac{d^2\mathsf{\Psi }}{d{\eta }^2}=0-\left\{\stackrel{\stackrel{\phantom{(}i}{\downarrow }}{\frac{1}{4{\eta }^2}\left[{\alpha }^2{\left(\frac{\omega \beta }{\kappa }-1\right)}^2-{\kappa }^2\right]}-\stackrel{\stackrel{\phantom{(}ii}{\downarrow }}{\frac{\omega }{2}\left[\frac{{\alpha }^2}{\kappa }\left(\frac{\omega \beta }{\kappa }-1\right)-\alpha +\frac{2\beta \varphi ({\kappa }_r+s)}{\omega \kappa }\right]}\frac{1}{\eta }\right\}\mathsf{\Psi }\hfill \\ \hfill +\left\{\frac{{\alpha }^2{\omega }^2}{4{\kappa }^2}+\frac{\varphi ({\kappa }_r+s)}{\kappa }\right\}\mathsf{\Psi }=0.\end{array}$$

(25)

Expanding the terms highlighted as i and $ii$ in Equation (25) and rearranging:

$$\begin{array}{c}\frac{d^2\mathsf{\Psi }}{d{\eta }^2}-\left[\frac{\frac{{\alpha }^2}{{\kappa }^2}{\left(\frac{\omega \beta }{\kappa }-1\right)}^2-1}{4{\eta }^2}-\frac{\omega }{2{\kappa }^2}\left[\frac{{\alpha }^2}{\kappa }\left(\frac{\omega \beta }{\kappa }-1\right)-\alpha +\frac{2\beta \varphi ({\kappa }_r+s)}{\omega \kappa }\right]\frac{1}{\eta }\right]\mathsf{\Psi }\hfill \\ \hfill +\left[\left(\frac{{\alpha }^2{\omega }^2}{4{\kappa }^2}+\frac{\varphi ({\kappa }_r+s)}{\kappa }\right)\frac{1}{{\kappa }^2}\right]\mathsf{\Psi }=0.\end{array}$$

(26)

Simplifying, Equation (26):

$$\begin{array}{c}\frac{d^2\mathsf{\Psi }}{d{\eta }^2}+\left[\frac{1-\frac{{\alpha }^2}{{\kappa }^2}{\left(\frac{\omega \beta }{\kappa }\right)}^2-1}{4{\eta }^2}+\frac{\omega }{2\kappa }\left[\frac{{\alpha }^2}{{\kappa }^2}\left(\frac{\omega \beta }{\kappa }-1\right)-\frac{\alpha }{\kappa }+\frac{2\beta \varphi ({\kappa }_r+s)}{\omega {\kappa }^2}\right]\frac{1}{\eta }\right]\mathsf{\Psi }\hfill \\ \hfill -\left[\left(\frac{{\alpha }^2{\omega }^2}{4{\kappa }^2}+\frac{\varphi ({\kappa }_r+s)}{\kappa }\right)\frac{1}{{\kappa }^2}\right]\mathsf{\Psi }=0.\end{array}$$

(27)

Equation (27) can be cast in the form of the Whittaker equation if a change of variable is defined thus:

$$\xi =2\eta \sqrt{a}.$$

(28)

where:

$$\begin{array}{ccc}\hfill a& =& \frac{1}{{\kappa }^2}\left(\frac{{\alpha }^2{\omega }^2}{4{\kappa }^2}+\frac{\varphi ({\kappa }_r+s)}{\kappa }\right)\hfill \end{array}$$

(29)

$$\begin{array}{cc}=& f{\left(s\right)}^2.\hfill \end{array}$$

(30)

Then:

$$\begin{array}{ccc}\hfill \frac{d\mathsf{\Psi }}{d\eta }& =& \frac{d\mathsf{\Psi }}{d\xi }\frac{d\xi }{d\eta }\hfill \end{array}$$

(31)

$$\begin{array}{cc}=& 2f\left(s\right)\frac{d\mathsf{\Psi }}{d\eta },\hfill \end{array}$$

(32)

and

$$\begin{array}{c}\hfill \frac{d^2\mathsf{\Psi }}{d{\eta }^2}=4f{\left(s\right)}^2\frac{d^2\mathsf{\Psi }}{d{\xi }^2}.\end{array}$$

(33)

Using Equations (32) and (33) in Equation (27) gives:

$$\begin{array}{c}\hfill \frac{d^2\mathsf{\Psi }}{d{\xi }^2}+\left[\frac{1-\frac{{\alpha }^2}{{\kappa }^2}{\left(\frac{\omega \beta }{\kappa }-1\right)}^2}{4{\xi }^2}+\frac{\omega }{4\sqrt{a\kappa }}\left[\frac{{\alpha }^2}{{\kappa }^2}\left(\frac{\omega \beta }{\kappa }-1\right)-\frac{\alpha }{\kappa }+\frac{2\beta \varphi ({\kappa }_r+s)}{\omega {\kappa }^2}\right]\frac{1}{\xi }-\frac{1}{4}\right]\mathsf{\Psi }=0\end{array}$$

(34)

$$\begin{array}{c}\hfill \frac{d^2\mathsf{\Psi }}{d{\xi }^2}+\left[\frac{1-\frac{{\alpha }^2}{{\kappa }^2}{\left(\frac{\omega \beta }{\kappa }-1\right)}^2}{4{\xi }^2}+\frac{\omega }{4\kappa f\left(s\right)}\left[\frac{{\alpha }^2}{{\kappa }^2}\left(\frac{\omega \beta }{\kappa }-1\right)-\frac{\alpha }{\kappa }+\frac{2\beta }{\omega \kappa }\left({\kappa }^{2}f{\left(s\right)}^2-\frac{{\alpha }^2{\omega }^2}{4{\kappa }^2}\right)\right]\frac{1}{\xi }-\frac{1}{4}\right]\mathsf{\Psi }=0\end{array}$$

(35)

$$\begin{array}{c}\hfill \frac{d^2\mathsf{\Psi }}{d{\xi }^2}+\left[\frac{1-\frac{{\alpha }^2}{{\kappa }^2}{\left(\frac{\omega \beta }{\kappa }-1\right)}^2}{4{\xi }^2}+\left[\frac{\frac{\omega {\alpha }^2}{4{\kappa }^3}\left(\frac{\omega \beta }{2\kappa }-1\right)-\frac{\alpha \omega }{4{\kappa }^2}}{f\left(s\right)}-\frac{\beta }{2}f\left(s\right)\right]\frac{1}{\xi }-\frac{1}{4}\right]\mathsf{\Psi }=0\end{array}$$

(36)

$$\frac{d^2\mathsf{\Psi }}{d{\xi }^2}+\left\{\frac{1-\frac{{\alpha }^2}{{\kappa }^2}{\left(\frac{\omega \beta }{\kappa }-1\right)}^2}{4{\xi }^2}+\left[\frac{\frac{\alpha \omega }{4{\kappa }^2}\left[\frac{\alpha }{\kappa }\left(\frac{\omega \beta }{\kappa }-1\right)-1\right]}{f\left(s\right)}-\frac{\beta }{2}f\left(s\right)\right]\frac{1}{\xi }-\frac{1}{4}\right\}\mathsf{\Psi }=0,$$

(37)

which can be expressed in the form of the Whittaker equation [28] as:

$$\frac{d^2\mathsf{\Psi }}{d{z}^2}+\left[-\frac{1}{4}+\frac{K}{z}+\frac{\left(\frac{1}{4}-{\mu }^2\right)}{z^2}\right]\mathsf{\Psi }=0,$$

(38)

where:

$$\mu =\frac{\alpha }{2\kappa }\left(\frac{\omega \beta }{\kappa }-1\right),$$

(39)

and:

$$K=\left[\frac{\frac{\alpha \omega }{4{\kappa }^2}\left[\frac{\alpha }{\kappa }\left(\frac{\omega \beta }{\kappa }-1\right)-1\right]}{f\left(s\right)}-\frac{\beta }{2}f\left(s\right)\right],$$

(40)

$$f\left(s\right)=\frac{1}{\kappa }\sqrt{\left(\frac{{\alpha }^2{\omega }^2}{4{\kappa }^2}+\frac{\gamma ({\kappa }_r+s)}{\kappa }\right)}$$

(41)

$$\begin{array}{ccc}\hfill \xi & =& 2(\kappa r_D+\beta )f\left(s\right),\hfill \end{array}$$

(42)

$$\begin{array}{cc}=& 2\eta f\left(s\right)\hfill \end{array}$$

(43)

3.1. Introducing the Linear Drift

Linear drift can be introduced to the convection–diffusion equation by applying it as a scalar velocity field $v_d$, since it is coordinate-independent, having a magnitude that acts on every point within the porous body. For the purpose of this analysis, linear drift is applied on the x-direction only. Using the general hydrodynamic description of the diffusivity coefficient:

$$\begin{array}{ccc}\hfill D\left(r\right)& =& D_m+D_o\left|v_d\right|,\hfill \end{array}$$

(44)

$$\begin{array}{cc}=& D_m+\alpha D_o\left(\overline{\omega }+\frac{1}{r}\right),\hfill \end{array}$$

(45)

$$\begin{array}{cc}=& D_m+v_m{D}_o\left(\omega +\frac{1}{r_D}\right),\hfill \end{array}$$

(46)

and defining the following relational variables:

$$\begin{array}{ccc}\hfill \alpha \overline{\omega }& =& v_m\omega ,\hfill \end{array}$$

(47)

$$\begin{array}{ccc}\hfill \frac{\alpha }{v_m}& =& \frac{\omega }{\overline{\omega }}=r_w,\hfill \end{array}$$

(48)

and dimensionless variables:

$$\begin{array}{ccc}\hfill \frac{\alpha }{r}& =& \frac{v_m}{r_D}\hfill \end{array}$$

(49)

$$\begin{array}{ccc}\hfill r_D& =& \frac{r}{r_w}\hfill \end{array}$$

(50)

$$\begin{array}{ccc}\hfill v_d& =& v-\frac{\alpha }{r}.\hfill \end{array}$$

(51)

Equation (46) can be written out after expanding and rearranging thus:

$$\begin{array}{ccc}\hfill D\left(r\right)& =& D_m+{\beta }_w\left(\overline{\omega }+\frac{1}{r}\right)\hfill \end{array}$$

(52)

$$\begin{array}{cc}=& \kappa +\frac{{\beta }_w}{r},\hfill \end{array}$$

(53)

where:

$$\begin{array}{ccc}\hfill \kappa & =& D_m+{\beta }_w\overline{\omega },\hfill \end{array}$$

(54)

$$\begin{array}{ccc}\hfill \overline{\omega }& =& \left(\frac{v_d\cos \theta }{\alpha }\right),\hfill \end{array}$$

(55)

$$\begin{array}{ccc}\hfill {\beta }_w& =& \alpha D_o,\hfill \end{array}$$

(56)

$$\begin{array}{ccc}\hfill \alpha & =& \left(\frac{Q}{2\pi h\theta }\right).\hfill \end{array}$$

(57)

3.2. Analytical Solution

The general solution of the Whittaker equation [28] is given as:

$$\begin{array}{ccc}\hfill \mathsf{\Psi }(\eta ,s)& =& A\left(s\right)M_{\kappa ,\mu }\left(\xi \right)+B\left(s\right)W_{\kappa ,\mu }\left(\xi \right).\hfill \end{array}$$

(58)

In Equation (58), $A\left(s\right)$ and $B\left(s\right)$ are arbitrary functions of ${}^{\prime }{s}^{\prime }$ to be determined by the requirements of the boundary conditions, while $M_{\kappa ,\mu }\left(\xi \right)$ and $W_{\kappa ,\mu }\left(\xi \right)$ are the Whittaker function, which can also be defined in terms of Kummer’s confluent hypergeometric functions as:

$$\begin{array}{c}\hfill M_{\kappa ,\mu }\left(\xi \right)=e^{-\frac{\xi }{2}}{\xi }^{\frac{1}{2}+\mu }M\left(\frac{1}{2}+\mu -\kappa ,1+2\mu ,\xi \right)\end{array}$$

(59)

$$\begin{array}{c}\hfill W_{\kappa ,\mu }\left(\xi \right)=e^{-\frac{\xi }{2}}{\xi }^{\frac{1}{2}+\mu }U\left(\frac{1}{2}+\mu -\kappa ,1+2\mu ,\xi \right).\end{array}$$

(60)

Therefore, Equation (58) can be presented in terms of the Kummer’s function as:

$$\psi (\xi ,s)=\left(A\left(s\right)M(a,b,\xi )+B\left(s\right)U(a,b,\xi )\right)e^{-\frac{\xi }{2}}{\xi }^{\frac{1}{2}+\mu }$$

(61)

In general, the Kummer’s function of the first kind:

$$M(a,b,\xi \to \infty )\to \frac{\mathsf{\Gamma }\left(b\right)}{\mathsf{\Gamma }\left(a\right)}{\xi }^{a-b}{e}^{\xi },$$

(62)

implies that the $M(a,b,\xi )$ function becomes unbounded when $\xi$ becomes large. Therefore, the arbitrary function coefficient $A\left(s\right)$ of Equation (61) must become zero for Equation (62) to satisfy the external boundary condition specified for the system. Equation (61) therefore reduces to:

$$\psi (\xi ,s)=e^{-\frac{\xi }{2}}{\xi }^{\frac{1}{2}+\mu }B\left(s\right)U\left(\frac{1}{2}+\mu -\kappa ,1+2\mu ,\xi \right).$$

(63)

The coefficient $B\left(s\right)$ is obtained from the application of the inner boundary condition. The function $U(a,b,\xi )$, variously referred to as the Kummer’s function of the second kind or the Tricomi function, decreases exponentially as $\xi$ increases and vanishes as $\xi$ becomes infinitely large as required by the inner boundary condition of this problem.

The tracer concentration $C(\eta ,s)$ can be now be written as:

$$C(\eta \left(r\right),s)=\mathsf{\Psi }{e}^{-\left\{\frac{1}{2}\int \left(\frac{\kappa -\alpha }{\kappa r+\beta }+\frac{\alpha \omega r}{\kappa r+\beta }\right)dr\right\}}.$$

(64)

Detailed mathematical transformation in dimensionless form and inverse Laplace transfom of the general solution are available in Appendix B. The solution to Equations (34)–(37) can therefore be expressed as:

$$\psi \left(\xi \right)=e^{-\frac{\xi }{2}}{\xi }^{\frac{1}{2}+\mu }U(\frac{1}{2}+\mu -\kappa ,1+2\mu ,\xi ).$$

(65)

However, the concentration C is defined in terms of $\psi$ as:

$$C=\psi e^{\left\{-\frac{1}{2}\int \left(\frac{D^{\prime }-\upsilon }{D}+\frac{1}{r_D}\right)d{r}_D\right\}},$$

(66)

or:

$$C=\psi e^{\left\{-\frac{1}{2}\int \left(\frac{\kappa -\alpha }{\kappa r_D+\beta }+\frac{\alpha \omega r_D}{\kappa r_D+\beta }\right)d{r}_D\right\}}$$

(67)

$$C=\psi {\left(\kappa r_D+\beta \right)}^{\left(\frac{{\kappa }^2-\alpha \kappa -\alpha \omega \beta }{{\kappa }^2}\right)}{e}^{\left\{-\frac{1}{2}\frac{\alpha \omega (\kappa r_D+\beta )}{{\kappa }^2}\right\}}$$

(68)

$$C={(\kappa r_D+\beta )}^{\left(\frac{{\kappa }^2-\alpha \kappa -\alpha \omega \beta }{{\kappa }^2}\right)}{e}^{\left\{-\frac{1}{2}\frac{\alpha \omega (\kappa r_D+\beta )}{{\kappa }^2}\right\}}{e}^{-\frac{\xi }{2}}{\xi }^{\frac{1}{2}+\mu }U\left(\frac{1}{2}+\mu -\kappa ,1+2\mu ,\xi \right),$$

(69)

so that:

$$C={\left(\eta \right)}^{\left(\frac{{\kappa }^2-\alpha \kappa -\alpha \omega \beta }{{\kappa }^2}\right)+(\frac{1}{2}+\mu )}{\left(2\sqrt{a}\right)}^{\frac{1}{2}+\mu }{e}^{\left\{-\left(\frac{1}{2}\frac{\alpha \omega }{{\kappa }^2}+\sqrt{a}\right)\eta \right\}}U(\frac{1}{2}+\mu -\kappa ,1+2\mu ,\left(2\sqrt{a}\right)\eta ),$$

(70)

or:

$$C={\left(\xi \right)}^{\left(\frac{k^2-\alpha k-\alpha \omega \beta }{k^2}\right)+(\frac{1}{2}+\mu )}{\left(2\sqrt{a}\right)}^{-\left(\frac{k^2-\alpha k-\alpha \omega \beta }{k^2}\right)}Exp\left\{-\left(\frac{1}{2}\frac{\alpha \omega }{\left(2\sqrt{a}\right)k^2}+1\right)\xi \right\}U(\frac{1}{2}+\mu -\kappa ,1+2\mu ,\xi ),$$

(71)

where ‘a’ is given as:

$$\begin{array}{ccc}\hfill a& =& \frac{1}{k^2}\left(\frac{{\alpha }^2{\omega }^2}{4k^2}+\frac{\varphi (k_r+s)}{k}\right)\hfill \\ & =& \frac{1}{4k^4}({\alpha }^2{\omega }^2+4k\varphi (k_r+s)),\hfill \end{array}$$

(72)

and:

$$\begin{array}{ccc}\hfill \xi & =& \left(2\sqrt{a}\right)\eta \hfill \\ & =& \left(2\sqrt{a}\right)(k{r}_D+\beta )\hfill \\ & =& \frac{1}{k^2}\sqrt{({\alpha }^2{\omega }^2+4k\varphi (k_r+s)}(k{r}_D+\beta )\hfill \end{array}$$

(73)

$${\xi }_{rD}=\frac{1}{k^2}\sqrt{({\alpha }^2{\omega }^2+4k\varphi (k_r+s)}(k+\beta ).$$

(74)

3.3. Weak-Form Numerical Solution of the Tracer Transport Equation

The weak-form solution of Equation (71) is now presented by considering (i) a ‘$pot$’ diffusion case where tracer flow is controlled by molecular diffusion with no hydrodynamic dispersion but with velocity drift and (ii) cases where convection dominated the molecular diffusion effect. In order to achieve this, the separation of variables is adopted with X-parameter $\left(X_p\right)$ and Y-parameter $\left(Y_p\right)$, defined thus:

$$C(x,y,s)=X_p\left(x\right)Y_p(y,s).$$

(75)

The X-parameter $\left(X_p\right)$ can be expressed as:

$$X_p\left(x\right)=\tau U(\frac{1}{2}+\mu -\kappa ,1+2\mu ,\sigma {\xi }_{X_p})e^{-\frac{1}{2}\left(\frac{1}{\sigma }-1\right){\xi }_{X_p}},$$

(76)

where the components and arguments of the Tricomi Kummer function $U(a,b,x)$ are defined thus:

$$\begin{array}{c}\hfill \mu =-\frac{1}{2},\end{array}$$

(77)

$$\begin{array}{c}\hfill {\xi }_{X_p}=\left[\frac{\beta v_0}{D_{0}d}+\frac{v_0}{D_0}x\right],\end{array}$$

(78)

$$\begin{array}{c}\hfill \sigma =\sqrt{1-\frac{4{\omega }^2{D}_0}{v_0^{2}d}},\end{array}$$

(79)

$$\begin{array}{c}\hfill \tau =\frac{\mathsf{\Gamma }(\frac{j}{k})}{\mathsf{\Gamma }\left(j\right)},\end{array}$$

(80)

$$\begin{array}{c}\hfill h=\frac{{\omega }^2\beta }{v_0{d}^2},\end{array}$$

(81)

$$\begin{array}{c}\hfill k=\sqrt{1-j},\end{array}$$

(82)

and Y-parameter $\left(Y_p\right)$:

$$Y_p(y,s)=\left(\frac{1}{3}{\xi }_{Y_p}^{\frac{1}{2}+\mu }\right)e^{-\frac{v_{0}y}{2D_o}}\left[I_{-\frac{1}{3}}\left(\frac{2}{3}L{\xi }_{Y_p}^{\frac{3}{2}+\mu }\right)+I_{+\frac{1}{3}}\left(\frac{2}{3}L{\xi }_{Y_p}^{\frac{3}{2}+\mu }\right)\right],$$

(83)

$$\begin{array}{ccc}\hfill \mu & =& 0,\hfill \end{array}$$

(84)

$$\begin{array}{ccc}\hfill {\xi }_{Y_p}& =& \left[\frac{1}{4}{\left(\frac{v_0}{D_0}\right)}^2+\frac{Rs+R\kappa +{\omega }^2}{D_0\lambda }y\right],\hfill \end{array}$$

(85)

$$\begin{array}{ccc}\hfill L& =& \frac{D_0\lambda }{Rs+R\kappa +{\omega }^2}.\hfill \end{array}$$

(86)

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, mirroring the positive half.

3.3.1. Separation Constant $\left({\omega }^2\right)$

The constant of separation $\left({\omega }^2\right)$ is obtained by rewriting the general advection–dispersion equation (ADE) for the flow of reactive tracers under the influence of linear drift thus:

$$D_0\left(\frac{\beta }{x}+d\right)C_{xx}-v_0\left(\frac{\beta }{x}+d\right)C_x+D_0\frac{\lambda }{y}{C}_{yy}-v_0\frac{\lambda }{y}{C}_y-R\kappa C=R{C}_t,$$

(87)

where the linear drift ratio is written as $d=\frac{v_d}{v_0}$. The linear drift ratio, d, is coordinate-independent; therefore, its magnitude can be applied equally to the y-axis or, in the case of a 3D system, the z-axis.

Substituting Equation (75) into Equation (87), dividing through by $X\left(x\right)Y(y,t)$ and rearranging (see Appendix C ) gives:

$$\begin{array}{c}\hfill D_0\left(\frac{\beta }{x}+d\right)X^{″}-v_0\left(\frac{\beta }{x}+d\right)X^{\prime }+{\omega }^{2}X\left(x\right)=0,\end{array}$$

(88)

$$\begin{array}{c}\hfill D_0\frac{\lambda }{y}{Y}^{″}-v_0\frac{\lambda }{y}{Y}^{\prime }-(Rs+R\kappa +{\omega }^2)Y(y,s)=0,\end{array}$$

(89)

where the component Equation (88) is time-independent, while Equation (89) is time-dependent and expressed in Laplace space with Laplace parameter s. Considering the following transformation parameters:

$$\begin{array}{c}\hfill \beta =r_w{\cos }^2\theta ,\end{array}$$

(90)

$$\begin{array}{c}\hfill \lambda =v_y{r}_w{\sin }^2\theta ,\end{array}$$

(91)

$$\begin{array}{c}\hfill x_w=r_w\cos \theta ,\end{array}$$

(92)

$$\begin{array}{c}\hfill y_w=r_w\sin \theta ,\end{array}$$

(93)

Equations (88) and (89) can be rewritten as:

$$\begin{array}{c}\hfill D_0\left(\cos \theta +d\right)X^{″}-v_0\left(\cos \theta +d\right)X^{\prime }+{\omega }^{2}X\left(x_w\right)=0,\end{array}$$

(94)

$$\begin{array}{c}\hfill {\mathcal{L}}^{-1}\left[D_0\left(v_y\sin \theta \right)Y^{″}-v_0\left(v_y\sin \theta \right)Y^{\prime }-(Rs+R\kappa +{\omega }^2)Y(y_w,s)\right]=0.\end{array}$$

(95)

Equation (94) is time-independent, while Equation (95) is time-dependent and expressed in Laplace space with Laplace parameter s.

3.3.2. Boundary and Initial Conditions

Equation (87) is governed by the following boundary conditions:

$$\begin{array}{cc}\hfill C(x,y,t=0)=C_i(x,y)& \mathrm{for}x=y=\mathbb{R}\end{array}$$

(96)

$$\begin{array}{cc}\hfill C(x=\pm \infty ,y,t)=0& \mathrm{for}y=\mathbb{R},t>0\end{array}$$

(97)

$$\begin{array}{cc}\hfill C(x,y=\pm \infty ,t)=0& \mathrm{for}x=\mathbb{R},t>0\end{array}$$

(98)

$$\begin{array}{cc}\hfill C(x_w,y_w,t)=C_0& \mathrm{for}t>0\end{array}$$

(99)

where, the transformation from the polar (radial) coordinate to Cartesian coordinate is given as $x_w=r_w\cos \theta$ and $y_w=r_w\sin \theta$. The concentration $C(x_w,y_w,t)$ of the tracer is known within the wellbore during a tracer test; thus, the solution for $t>0$ can be obtained by solving Equation (87) within the porous formation outside the wellbore.

The Y-parameter $Y_p$ is in Laplace space and requires a numerical inversion scheme, such as the Gaver–Stehfest algorithm [18,29], Talbot inversion algorithm [30] or Euler inversion [31,32] algorithm. The scripts for these algorithms are open source and are readily available for download from the Mathworks website [33]. Out of the two algorithms attempted for the numerical inverse Laplace operation, Gaver–Stehfest and Euler inversion, only the Euler inversion algorithm produced a stable result (see also Avdis and Whitt [34]). Hence, Euler’s Inversion Algorithm was used for the numerical inverse Laplace operation in the numerical code developed for this work.

Flow convection depends on the velocity of the system and is modelled by considering a heterogeneous system. Anisotropic porosity distribution was generated using the random probability density function (PDF) allowing for non-uniform velocity distribution to be computed. Typical porosity distribution is shown in Figure 2 with a mean value of $0.25$ and standard deviation of $0.64$. The corresponding computed velocity distribution is shown in Figure 3.

4. Analysis of Results

In this work, the effect of linear drift on the tracer propagation profile in a typical formation of thickness $h=9$ m and injection well of radius $r_w=0.127$ m was investigated. Injected particles are considered to be components of surfactants or polymers used in EOR processes, but in this case, the chemical reaction was neglected. The injection rate was fixed at $1.4\times {10}^{-3}$ m${}^3$/s and dispersion coefficient $D_0=1.4\times {10}^{-3}$ m${}^2$/s was applied. The tracer concentration distribution ratio $C_i$ was monitored under three continuous injection periods of ten (10) days, fifty (50) days and seventy (70) days respectively.

The developed solution was first tested by evaluating the error limit associated with the separation of variable parameter (${\omega }^2$) for different angle $\theta$. Simulation runs involved one hundred (100) values of ${\omega }^2$ as defined by Equation (81)—ranging from $0.01\times {\omega }_1^2$ corresponding to a value of ${\omega }^2=2.33\times {10}^{-8}$ to ${\omega }_{100}^2$ corresponding to a value of ${\omega }^2=2.59\times {10}^{-6}$—with an incremental value of $0.01\times {\omega }_{100}^2$. The minimum error corresponds to the value of separation variable ${\omega }^2=2.33\times {10}^{-8}$, as indicated by $X_p,Y_{p,t},X-Y_{p,t}$ plots (Figure 4, Figure 5 and Figure 6). The combined solution $X-Y_{p,t}$ indicated that there exist points of singularity at $0^{°},{180}^{°}$ and ${360}^{°}$ when taking the inverse of the coordinate point $y.$

A typical result of the concentration distribution as a function of angle $\theta$ obtained from the solution of Equation (75) is shown in Figure 6 for five 5 selected ranges of values of ${\omega }^2$. The concentration profile basically grows with increasing values of ${\omega }^2$ between ${\omega }^2=2.33\times {10}^{-8}$ and ${\omega }^2=1.06\times {10}^{-6}$ and between $\theta ={150}^{°}$ and $\theta ={300}^{°}$, after which the impact of drift sets in. With the onset of drift, the concentration profile shortens but with a much wider coverage for ${\omega }^2=1.58\times {10}^{-6}$. The region covered by ${\omega }^2=2.10\times {10}^{-6}$ can be seen to have drifted to $\theta ={240}^{°}$ to $-{300}^{°}$.

In order to further examine the impact of the drift parameter on separation constant ${\omega }^2$, $X_p$ and angle ${\theta }^{°}$, a full simulation run at different time intervals of $t=10$ d, $t=30$ d, 50 d and 70 d was carried out and results presented in Figure 7, Figure 8, Figure 9 and Figure 10. The magnitude of error due to ${\omega }^2$ generally reduces with increasing computational time. In Figure 10, the influence of drift is more pronounced. A typical concentration distribution within the porous media with linear drift after 10 days is shown in Figure 11.

Linear Drift Effect and Concentration Distribution Profile

In a single-well tracer or surfactant injection systems, a primary tracer bank is first injected into a formation containing oil at residual saturation. The bank is then followed by a bank of tracer-free water. The well is then shut in for a period of time, after which the well is produced and the concentration profiles monitored. Where there is fluid migration in the formation due to the movement of an underlying basal water displacing the injected tracer banks during shut-in, the production profile may be distorted and fluid pathway rerouted. In this situation, the effect of linear drift on fluid flow behaviour will have to be investigated.

Generally, in isotropic and homogeneous systems, where there is no linear drift or natural convection, the tracer propagation profile is expected to follow a cyclic pattern. The concentration distribution will be equal at an equidistant radial position from the injection well. In this work, a system with variable porosity distribution was modelled and the corresponding non-uniform velocity profile used in the computation of the drift ratio. In this case, the tracer propagation profile is expected to follow a natural pattern determined by the interplay of the forces associated with the system variables, such as the tracer injection rate.

In order to investigate the effect of the drift intensity on the tracer concentration profile, three (3) values of drift ratio $d=0.03,$ $0.09$ and $0.2$ were evaluated for time duration $t=10$ d, 50 d and 70 d. In the presence of linear drift in a heterogeneous system, however, there exists an unequal distribution of tracers along the principal x-axis and varies in a non-uniform manner along the positive and negative radial distance. Where the angle $\theta$ increases in the $+ve$ and $-ve$ y-axis, the system convection will lead to an increase or decrease in the tracer concentration distribution depending on the degree of system heterogeneity. It is important to note that an increase in the tracer concentration distribution will be observed in a homogeneous system.

The results of the tracer tests for different drift ratios at different time intervals are shown in Figure 12. At a fixed drift ratio (e.g., $d=0.2$), the effect of the linear drift on the tracer concentration profile is progressively more pronounced at later times; with the lowest tracer concentration ratio at later time indicating a high drift effect. Similarly, at any particular point in time (e.g., at time $t=70$ d), linear drift has a greater effect at a higher drift ratio (e.g., $d=0.2$).

5. Conclusions

The study of the transport of tracers with linear drift is an important aspect of porous media characterisation. Despite the extensive research in this field, particularly in solving the ADE both analytically and numerically, there is yet to be a consideration for the closed-form solution of the ADE in systems where the effect of linear drift may predominate and an exact analytical solution has not been reported in the literature. The following conclusions can be drawn from this work:

• A new closed-form analytical solution to the radial transport of tracers in porous media under the influence of linear drift and radial convection was developed. The radial transport equation was cast in the form of the Whittaker equation after adopting variable transformation and an exact solution for the tracer concentration derived therefrom.
• The weak-form solution was developed by splitting the transformed equation, adopting a common separation constant and invoking inverse Laplace transformation using the Euler inversion algorithm.
• Variable transformation from a radial to a Cartesian coordinate system was used to analyse the concentration distribution profiles in three-dimensional graphical plots.
• The obtained solutions are generally stable and dependent on the precision with which the separation constant $\left({\omega }^2\right)$ can be determined. This is important because the exponential term in the inversion formula may amplify the numerical error. The maximum error quantified by the separation constant is ${\omega }^2=2.10\times {10}^{-6}$.
• The influence of linear drift on the concentration profiles was evaluated in the x-direction for a system with nonhomogeneous porosity distribution and variable velocity profiles.
• The results of the analyses indicated that the effect of linear drift on the tracer concentration profile is dependent on system heterogeneity and progressively becomes more pronounced at later times.
• Practical application was demonstrated in a typical EOR process involving the injection of chemicals (e.g., surfactants or polymers), but without a chemical reaction. Another possible application is a single-well chemical tracer injection method for measuring residual oil saturation and fluid flow behaviour and characterisation in porous media.
• This work can be extended to the analysis of systems involving variation of tracer injection intensity, where spreading may occur in the r-$\theta$ or x-y plane. The developed solution can also be extended to systems where moderate-to-high intensity tracer flow with linear drift manifests. In this case, the tangential velocity component of the drift velocity becomes significant and will have to be included in the solution approach.

The new solution to the convective-diffusion equation developed and tested in this work demonstrates the need to study the impact of linear drift on transport of tracers in porous media. This is particularly important, since the arrival times of tracers may be significantly influenced by the drift intensity.