Analytical solution for a radial advection-dispersion equation including both mechanical dispersion and molecular diffusion for a steady-state flow field in a horizontal aquifer caused by a constant rate injection from a well

This study presents the analytical solution for a radial advection-dispersion equation for a 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. The derived analytical solutions were compared with finely-meshed finite difference solutions in steady-state and periodic steady-state problems with typical parameters. The results suggest that the analytical solution is exactly derived and ready for application. Comparisons with analytical solutions ignoring molecular diffusion suggest that the derived analytical solution should be used when the product of the decay constant and the retardation factor and the ratio of injection rate to diffusion coefficient are small. Comparisons with analytical solutions with Dirichlet-type boundary conditions confirmed that Robintype boundary conditions should be used to exactly evaluate the concentration profile.


INTRODUCTION
The radial advection-dispersion problem appears when analyzing contaminant transport in steady radial flow from an injection well that perfectly penetrates a homogeneous horizontal aquifer of uniform thickness ( Figure 1). The analytical solution of this problem is important, not only because it can be applied to the analysis of contaminant transport around wells, but also because it can provide benchmarks without numerical diffusion for the verification of numerical codes.
There have been several studies on the analytical solution of the radial advection-dispersion problem. Moench and Ogata (1981) derived the analytical solution using Laplace transform with respect to the time variable. Hsieh (1986) gave the alternative integral representation of an inverted version of Moench and Ogata (1981)'s solution that is much easier to compute. While these solutions con-Correspondence to: Masaatsu Aichi, Department of Environment Systems, Graduate School of Frontier Sciences, The University of Tokyo, Kashiwanoha 5-1-5, Kashiwa, Chiba 277-8563, Japan. E-mail: aichi@edu.k.utokyo.ac.jp sidered Dirichlet-type boundary condition at the well, Chen (1987) adopted Robin-type boundary condition based on the mass balance across the well screen and showed that the analytical solution with Dirichlet-type boundary conditions significantly overestimated the concentration. Phillip (1994) examined a case where the dispersion coefficient was proportional to the exponential of the Péclet number. Veling (2001) derived the integral representation of the generalized analytical solution in the time domain for all types of boundary conditions and arbitrary initial conditions based on the generalized Hankel transform. Hsieh and Yeh (2014) developed an analytical solution for a two-zone aquifer in the Laplace domain. Lai et al. (2016) showed that the combination of the generalized integral transform technique and the Laplace transform gives a robust integral representation of the analytical solution. However, these studies ignored the molecular diffusion term because they focused on cases with high injection rate.
On the other hand, tracer tests for low groundwater velocity field are important in studies related to radio-active waste disposal because molecular diffusion is not negligible (e.g. Gustafsson and Morosini, 2002). Also, the injection rate in the inter-aquifer or surface-aquifer contaminant transport through imperfectly packed wells might be small because the process is controlled by the natural flow rate of groundwater or rainfall through an imperfect packer, although it potentially continues for long periods (Doble  et al., 2018). Recently, Haddad et al. (2015) presented the analytical solution including molecular diffusion for the tracer injection test with Dirichlet-type boundary conditions at the well. Huang and Goltz (2017) also included the molecular diffusion term in the analytical solution for the extraction of contaminant in vadose zones. However, the number of studies on analytical solutions including molecular diffusion is far fewer than those ignoring molecular diffusion. Although Robin-type boundary conditions have been recognized as an appropriate expression to model mass transport at the well in injection problems according to Chen (1987), it has not yet studied for the case of finite molecular diffusion. Robin-type boundary conditions mean a boundary condition to specify a linear combination of the mass flux and the concentration. Note that Chen (1987) referred to this as a Cauchy boundary, although Cauchytype boundary conditions mean a set of boundary conditions specifying both the concentration and the mass flux at one side of boundary. In addition, previous studies adopted Laplace transforms with respect to time variables because they were interested in the constant concentration of contaminants in the injection fluid. Because the concentration of contaminants from agricultural activities such as pesticides, nitrogen, phosphorous, etc., can differ by season, periodic changes in concentration might be an important factor for problems related to contamination through imperfectly constructed wells. In addition, the inversion of the Fourier transform is much easier than that of the Laplace transform and it is not difficult to obtain the transformed boundary condition based on time series data of the measured concentration using a fast Fourier transform. Thus, Fourier transforms with respect to time variables are also considered worthy of focus.
In this study, the analytical solution for a radial advection-dispersion transport including molecular diffusion, a first-order decay, and a retardation, in a steady-state flow field caused by a continuous injection at a well under a Robin-type boundary condition is studied after Laplace or Fourier transform with respect to time variables.
The governing equation of the above mentioned problem is, where C is the concentration; t is the time; r is the distance from the well center; R is the retardation factor; D is the effective molecular diffusion coefficient; α is the dispersivity length; Q is the injection rate; H is the aquifer thickness; ϕ is the porosity; and λ is the decay constant. Applying Laplace or Fourier transform with respect to time variables, the governing equation becomes as follows: under the boundary conditions: where μ is the constant; r W is the well radius; C 0 is the concentration of the injected fluid; (ˉ) denotes the transformed variable. If the Laplace Transform with an initial condition of C(r, 0) = 0 is applied to Equation (1), where s is the Laplace operator. If the Fourier Transform is applied to Equation (1), where ω is the angular frequency; i is an imaginary unit. For a special case of the steady-state (ω → 0 in Equation (5)), The main objective of this study is providing the analytical solution to Equation (2) under Equation (3).

DERIVATION OF THE ANALYTICAL SOLUTION
By the following transformation of variable , Equation (2) becomes If the form of analytical solution is assumed to be the function f(r D ) must satisfy the following equation: so that Equation (9) satisfies Equation (8). Because Equation (10) is a Kummer's differential equation (Kummer, 1837), the general solution is known as follows (e.g. Olde Daalhuis, 2010): where C 1 and C 2 are integral constants; M and U are the Kummer's confluent hypergeometric function of the first and second kind, respectively; a and b are defined as follows: Based on Equation (3), the boundary conditions in terms of r D is Then, the final solution is obtained as follows: The numerical method to evaluate the function U has been described in the existing literature (e.g. Olde Daalhuis, 2010). For example, Mathematica software package, version 11.3 (Wolfram Research, Inc., 2018), is capable of computing Equation (15) and was used for the calculations in the following chapter.
On the other hand, in order to provide benchmarks which can be easily computed, it is worth mentioning the following two special cases where Equation (15) is expressed with a combination of elementary functions. Based on the known relations on the confluent hypergeometric function (e.g. Olde Daalhuis, 2010), one obtains when a = 0, that is, is satisfied. Another special case is when b = a + 1, that is, is satisfied.

COMPARISONS WITH FIRST-ORDER UPWIND FINITE DIFFERENCE SOLUTIONS
In order to demonstrate the validity and characteristics of the derived analytical solution, numerical examples with typical parameters were computed and compared with finely meshed implicit first-order upwind finite difference solutions, the analytical solution ignoring molecular diffusion after Equation (12) in Chen (1987), and the analytical solution considering molecular diffusion with Dirichlettype boundary condition at the well after Equation (13) in Haddad et al. (2015).

Steady-state solutions
The first example are steady-state solutions which were calculated with α = 1 m, r W = 0.2 m, D = 10 -9 m 2 s -1 , μ = 7.33 × 10 -9 s -1 , and νD = Q 2πHϕ was changed from 10 -10 to 10 -8 m 2 s -1 . Finite difference solutions were calculated with a spatial grid size of 10 -4 m to avoid numerical dispersion and the boundary condition for r → ∞ was approximately set at r = 10 m. These results of all analytical and numerical solutions are shown in Figure 2. Another example of a steady-state solution was computed with α = 1 m, r W = 0.2 m, D = 10 -9 m 2 s -1 , νD = 10 -9 m 2 s -1 and μ ranged from 7.33 × 10 -11 to 7.33 × 10 -9 s -1 . The boundary condition for r → ∞ was approximately set at r = 100 m for this case. The results are shown in Figure 3.
In these examples, the analytical solutions and the finite difference solutions satisfactory matched each other in all cases. As ν increases, that mean of the ratio of injection rate to diffusion coefficient becomes larger, the concentration profiles were transferred further away from the well and the concentration gradient becomes gentle due to increasing mechanical dispersion effects. As μ increases, that is the decay constant or retardation factor becomes larger, the concentration attenuates closer to the injection well. These results suggest that both analytical and numerical solutions give physically reasonable concentration profiles and the analytical solution is exactly derived.
Comparisons with the analytical solution ignoring molecular diffusions showed that ignoring molecular diffusions leads to significant overestimation when ν is small (Figure 2) or μ is small (Figure 3). However, the analytical solution of this study is asymptotic to the analytical solution ignoring molecular diffusion as ν increases ( Figure 2) and μ increases (Figure 3). The results also showed that the analytical solution with Dirichlet-type boundary condition provides significant overestimation of the concentration (Figures 2 and 3). These results suggest that diffusion and dispersion transport across the well screen are not negligible when the injection rate is small and it is strongly recommended to use Robin-type boundary condition.

Periodically steady-state solution
The analytical solution with a Fourier transform gives periodically steady-state solutions. Then, the analytical solution can provide benchmarks for a transient simulation. In order to illustrate this, a test problem with the following boundary condition for the concentration at the injection well, was considered, where T is the period of oscillating concentration. The period of oscillation was set to 1 year in this case study. The parameters for this example were α = 1 m, r W = 0.2 m, D = 10 -9 m 2 s -1 , μ = 7.33 × 10 -9 s -1 , and ν = 1. Based on the superposition of two analytical solutions, one with the boundary conditions C 0, t = C 0 2 and another C 0, t = C 0 2 cos 2π T t, the analytical solution for the boundary condition of Equation (20) was computed. The finite difference solutions were calculated with spatial grid size of 10 -2 m and the time step size of 1 day to avoid numerical dispersion effect and the boundary condition for r → ∞ was approximately set at r = 10 m. In order to reach a periodic steady-state, the transient simulation was conducted over 10 cycles of the period from the initial condition C(r, 0) = 0. This comparison, shown in Figure 4, demonstrates that the analytical solutions and finite difference solutions satisfactory match.

CONCLUSIONS
The analytical solution was obtained for a radial advection-dispersion equation for a steady-state flow field in a horizontal aquifer caused by a constant rate injection from a well, including both mechanical dispersion and molecular diffusion in addition to the retardation and the first-order attenuation under Robin-type boundary condition at the injection well. The derivation process was based on variable transform and reducing to Kummer's differential equation.
The derived analytical solution contains a confluent hypergeometric function, which is a famous special function. Although it can be calculated with existing mathematical software, for easier computing and in order to make benchmarks for numerical codes, two special cases in which the analytical solution can be expressed with a combination of elementary functions were also presented.
The derived analytical solutions were compared with finely-meshed finite difference solutions in steady-state and periodic steady-state problems. Because the analytical and numerical solutions satisfactory matched under several independent conditions, it is suggested that the analytical solution is exactly derived and ready for various applications, e.g. analysis of low injection rate tracer tests, and estimation of the contamination through leaky wells.
Comparisons with analytical solutions ignoring molecular diffusion showed that the derived analytical solution should be used when the product of the decay constant and the retardation factor and the ratio of injection rate to diffusion coefficient are small. Comparisons with analytical solutions with Dirichlet-type boundary condition confirmed that Robin-type boundary conditions should be used to exactly evaluate the concentration profile.
The obtained analytical solution is expected to converge to Chen (1987)'s solution in case D → 0. However, this limiting form of the confluent hypergeometric function has not been studied to the best of the authors' knowledge and remains a theoretical challenge of the derived formulation that is expected to be addressed in a future study.

ACKNOWLEDGMENTS
A part of this research was financially supported by Grants-in-Aid for Scientific Research 17H01360. The finite difference solution is only plotted as points every 0.02 m for easy observability. The first 1 m is selected to focus on a significant change profile. The solid line denotes the analytical solution of this study. See text for details