Analytical Solution of 1D Advection-Dispersion Equation In Finite Groundwater Reservoir With Spatially Dependent Dispersivity And Two Inputs Sources With Source-Sink Term Impact

This study derives an analytical solution of a one-dimensional (1D) advection-dispersion equation (ADE) for solute transport with two contaminant sources that takes into account the source term. For a heterogeneous medium, groundwater velocity is considered as a linear function while the dispersion as a nth-power of linear function of space and analytical solutions are obtained for n = 1 . 0 and 2 . 0. The solution in a heterogeneous ﬁnite domain with unsteady coeﬃcients is obtained using the Generalized Integral Transform Technique (GITT) with a new regular Sturm-Liouville Problem (SLP). The solutions are validated with the numerical solutions obtained using MATLAB pedpe solver and the existing solution from the proposed solutions. We exanimated the inﬂuence of the source term, the heterogeneity parameters and the unsteady coeﬃcient on the solute concentration distribution. The results show that the source term produces a solute build-up while the heterogeneity level decreases the concentration level in the medium. As an illustration, model predictions are used to estimate the time histories of the radiological doses of uranium at diﬀerent distances from the sources boundary in order to understand the potential radiological impact on the general public.


Keyworks
Two continuous sources · heterogeneous medium · ADE · Source-sink term · GITT List of symbols: The following symbols are used in this paper a Constant a 1 Constant, of dimension (L −1 ) a 2 Heterogeneity parameter, of dimension (L −1 ) C Solute concentration in groundwater flow at (x, t), of dimension (M L −1 ) C i Uniform initial concentration in liquid phase, of dimension (M L −1 ) C 1 Reference uniform concentration at the origin, of dimension (M L −1 ) C 2 Reference uniform concentration at the end oh the domain, of dimension (M L −1 ) DF Ingestion dose coefficient of the nuclide for the adult age group (Sv/Bq). D 0 Constant mechanical dispersion coefficient, of dimension (L 2 T −1 ) D * Constant molecular diffusion coefficient, of dimension (L 2 T −1 ) D x0 Effective dispersion coefficient in homogeneous medium, of dimension (L 2 T −1 ) f The fraction of sorbing site Eigenvalues of the SLP, m = 1, 2, 3, 4, and 5 τ The porosity of the porous media µ Effective first-order decay coefficient, of dimension (T −1 ) µ l First order decay rate of the liquid phase, of dimension (T −1 ) µ s First order decay rate of the solid phase, of dimension (T −1 ) ϕ m Eigenfunction of the SLP belonging to mth eigenvalue Ψ m Orthogonal function corresponding to mth eigenfunction φ Unsteady function ρ b Porous media bulk density, of dimension (M L −3 ) τ The tortuosity θ Volumetric water content θ m Time-dependent coefficients in Fourier series

Introduction
The pollution of groundwater has seriously increased this last years due to the socioeconomic development and, interest many Scientists such as chemists, geologists, engineers, physicians. The groundwater can be contam-inated through the infiltration leachate process causing a reduction of groundwater flow rate, which results in difficult rehabilitation of groundwater. The contaminants can be radionuclides, chemical fertilization, industrial waste storage or spills, hospital wastes, leakages from petrol pumps, septic systems, organics matter, and wells. This pollution has a serious impact on human lives, human activities as well as productivity. Thus, it is very important to study the solute transport characteristics in porous media for human and ecological health. The transport of a solute in porous media is governed by the processes of advection and dispersion. The mediums through which the solute transport occurs are homogeneous (when the dispersivity does not depend to the position) or heterogeneous (when the dispersivity is a distance dependent function). This characteristic of a medium as an important role to the solute transport trough them (Kumar and Yadav, 2014;Sanskrityayn et al., 2018;Yadav and Kumar, 2019). The advection-dispersion equation can be solved numerically or analytically. Although numerical solutions provide great flexibility, analytical solutions are still pursued by many scientists because they are relatively transparent with respect to model inputs and outputs, and they can provide better physical insights into the problems (Park and Zhan, 2001). Several methods are used to solved analytically the advection dispersion transport. Many analytical solutions of conservative and non-conservative contaminants transport in subsurface of water considering adsorption are present in the literature. Those models are carried out to investigate the effect of solute and medium properties to the solute transport, as well as the medium rehabilitation. Using the Laplace Transform Technique (LTT), Van Genuchten (1981) obtained analytical solutions for ADE with zero-order production, simultaneous adsorption, and first-order decay for chemical transport. Van Genuchten and Alves (1982) proposed several solutions that quantitatively describe the behaviour of solute in surface/subsurface in finite and semi-infinite medium considering adsorption, first-order decay and zero-order production in homogeneous medium. Considering the dispersion coefficient as a linear and exponential increasing space function, Yates (1992) developed an analytical solutions of advection dispersion equation in one-dimension. Zoppou and Knight (1997) obtained an analytical solution for the spatially dependent dispersion. Su et al. (2005) developed an analytical solution to ADE with a spatially and temporally varying dispersion coefficient for predicting solute transport in a steady and saturated subsurface flow through heterogeneous porous media. Zhan et al. (2009) obtained an analytical solution for two-dimensional solute transport using first and third type boundary conditions. Mazarheti et al. (2013) were first derived analytical solution of one-dimensional solute transport with several point sources and arbitrary time-dependent emission rate. The obtained analytical analytical was compared with the numerical solution. Using the LTT, . Kumar and Yadav (2014) developed a one-dimensional analytical solution for conservative solute transport in heterogeneous porous media for uniform and varying pulse type input point. Kumar et al. (2019) studied the effect of the source/sink term on the solute transport.
One of the most used methods to solve the ADE in finite domain with spatially dependent coefficients is the Generalize Integral Transform Technique (GITT). This technic has the particularity to apply to all ADES whatever the form of the coefficients. The GITT was applied by Liu et al. (2000) to solve the one-dimensional ADE in heterogeneous porous media with source/sink term, coupled with either linear or nonlinear sorption and decay. The GITT coupled with the LTT were used by .  to solve analytically an 1-D ADE in a finite spatial domain with an arbitrary time-dependent inlet boundary condition. For a finite spatial domain, the 1-D ADE considering the sorption and desorption of solute, with arbitrary space dependent coefficients was salve analytically using the GITT (Skaggs et al., 2007). . Pérez Guerrero et al. (2009) presented a new analytical method to solve a 3-D ADE in a finite domain with time varying boundary condition for both transient and steady-state regimes using change of variables in combination with the Classic Integral Transform technique (CITT).  presented an analytical solution of two-dimensional ADE in cylindrical coordinates using a combination of the second kind finite transform method and the GITT. Recently, Bharati et al. (2017), Bharati et al. (2018) and Bharati et al. (2019) presented an analytical solution of solute transport with distance depending coefficients without source term using the GITT with a new regular Sturm-Liouville problem (SLP) with a self-adjoint operator to derive analytical solutions in a finite domain. Despite those study proposed novel methods to solve the ADE or/and news SLP with self-adjoint operator, they do not consider many parameters on the solute transport such as the source term, the presence of many solute sources. Additionally, in the most of these studies the coefficients do not depend on time. Including many parameters in pollutant transport equation is helpful to cover many aspects of contaminant transport in groundwater in more natural way (Chaudhary et al., 2020).
In the present study, an analytical solution of one-dimensional advection-dispersion equation (ADE) originating from two point sources along the longitudinal direction of groundwater reservoir horizontal flow domain is derived. The coefficients of the ADE are considered spatio-temporal dependent. Temporal dependence is considered through different functions for velocity and dispersion coefficients, which may be considered by different combinations of linear, exponential forms. The analytical solution is obtained with the help of the GITT with associate advection-dispersion SLP using a self-adjoint operator. We demonstrate the influence of the heterogeneity parameter, source term, the nth-power parameter in the solute strength. The analytical solution are validated with the help of the numerical solution. The developed analytical model is applied to illustrate the time variation of uranium radioactivity concentration and the potential radiological impact through the time histories of the radiological doses at different distances from the origin.

Mathematical formulation and analytical solution
In this study, we consider one-dimensional solute transport in a horizontal aquifer. The one-dimensional advection-dispersion equation (ADE) in general form, describing the solute transport in groundwater domain from an injected source in a finite heterogeneous medium can be written as where C is the solute concentration of the dispersing pollutant mass in the liquid phase, S is the solid phase concentration, representing the sorbed solute per unit mass of solid, θ is the volumetric water content, f represents the fraction of sorption site, ρ b is the bulk density of porous media, η is the porosity of the porous media, D is the dispersion coefficient, u is the seepage velocity, µ l and µ s are the first order decay rate of the liquid an solid phases concentration respectively, the non-homogeneous term q(x, t) represents an arbitrary space and time variable zero-order production, x is the spatial coordinate and t is the time. The left-hand side of Eq.
(1) represents change in solute concentration in liquid and solid phases with time respectively for the first and the second terms. The right-hand side of the Eq. (1) represents the influence of the dispersion on the solute concentration distribution by the first term and the change of the solute concentration due to advective solute transport by second term. The third and the fourth terms of the right-hand represent the first-order decay of solute in the liquid and solid phases in the medium respectively. The fifth term represent the zero-order production (q(x, t) > 0) or sink (q(x, t) < 0) for solute which represents internal/external production or sink of the solute in the medium. The effective first-order decay coefficient may be constants or functions of time or space or both [28]. When study the contamination of groundwater flow for miscible solute with instantaneous sorption, it exists relationship which defines the sorbed concentration at all times. This relationship relates the sorbed S concentration to the liquid concentration through the sorption isotherm Γ Even though for more solute or geological formation Γ is not linear, the linear isotherm may be a reasonable approximation at low solute concentrations (Skaggs and Leij, 2002). Therefore, we considered the expression for a time-dependent linear equilibrium given by (Sim and Chrysikopoulos, 1996;Singh and Das, 2015) where F represents the mass fraction of sorption particles where sorption is instantaneous, k d is referred to as the distribution coefficient. It should be noted that for instantaneous sorption in all the mass F = 1 . The effective diffusion-dispersion coefficient in one direction is the sum of the soil molecular diffusion coefficient, taking into account the tortuous nature of the solute pathways, and the mechanical dispersion coefficient (Thorne, 2014). In this study, the groundwater is assumed to be governed by the Darcy equation. Due to the steady recharge, groundwater velocity will be increasing linearly with position (Serrano, 1992). Hence dispersion coefficient and velocity both are considered spatially dependant in general form. Thus, in general form, dispersion coefficient and velocity both are considered spatially and temporally dependent. For this study, the expression of groundwater velocity and the dispersion coefficient are considered to follow the dispersion theory according to which the space dependent dispersion is proportional to the nth-power of the space dependent velocity (ICRP, 2012). Thus, the expression of velocity and hydrodynamic dispersion coefficients are considered as follow where τ is the tortuosity, D 0 , D * and u 0 are respectively constant mechanical dispersion coefficient, molecular diffusion coefficient and velocity in a steady flow domain through a homogeneous porous medium. s is a coefficient whose dimension is inverse of the time variable, it represents the unsteady parameter that characterize the flow resistance quantity. Thus, φ(st) is an expression in non-dimensional variable known as the unsteady function. s = 0 corresponds to the temporally independent parameters. Generally, the expressions of the temporal function are chosen so that φ(st) = 1 for t = 0 or s = 0 . The former case represents the steady flow and the latter case represents the initial state. a 1 is a spatially non dimensional parameter and a 2 represents the spatially independent parameters known as heterogeneity parameter.
In this study, the source/sink is expressed by single function on space and time-dependent as that proposed by  q where q 0 is the uniform zero-first order production coefficient. Substituting Eqs.
(2), (3), (4) and (5) into Eq. (1) we get where R = θ + f ρ b η F k d is the retardation factor and µ = µ l + µ s ρ b η F k d is the effective first order decay constant in the medium.
The formulation of our problem to be completed by assuming a set of initial and boundary conditions. Groundwater flow is considered along the x-axis in a finite media, means the direction of the flow of water is from x = 0 to x = l. Before the pollutant's source is injected in the reservoir, it is supposed to contain a certain quantity of the solute concentration in the liquid phase with arbitrary space variable distribution in the flow direction and in the solid phase with uniform distribution as proposed by Thakur et al. (2019). For short special domain, Chaudhary et al. (2020) proposed to use exp(−sec(λx)) because its decreasing rate is much slower than exponential and could well represent slow movement of groundwater. The initial condition may be written in general form as follows where C i is the input liquid phase concentration and K i s the input solid phase concentration. The geological formation of the medium through whom the dispersion occurs is assumed to be bounded by two parallel planes, e.g. the planes at x = 0 , x = l with two pollutant sources located at the two planes so that the pollutant enters through the planes. Because the groundwater flows from x = 0 to x = l , this situation can referred to advection-dispersion with simultaneous input in the flow direction (case of input on the plane x = 0 ) and against the direction of the flow (case of input on the plane x = l). In this study, the surface concentrations variables are considered and may be expressed in general form as follows where f 1 (t) and f 2 (t) are arbitrary functions that quantify the amount of pollutant concentration entering through the planes and respectively.
The initial condition Eq. (7) and the boundary conditions Eqs. (8) and (9) are similar to the situation of the diffusion into a plane sheet of material whose surfaces, x = 0 , x = l , are maintained at constant concentrations C 1 and C 2 respectively with a general initial heat distribution studied by . Crank (1975). One of the most interesting cases of such problem is the situation where the two surfaces are maintained at the same temperature.
Usually, the analytical solution of the ADE defined in Eqs. (6)-(9) is obtained with the help of the Generalized Integral Transform Technique (GITT) which is based on the eigen-value expansion and the evaluation of the transition matrix. Due to the general form of the coefficients in Eq. (6) it is difficult to get the analytical solution of the system of Eqs. (6)-(9) using the GITT. So, the ADE needs to be reduced to a solvable form. Therefore, the following new time variable is introduced (Crank, 1975) Using Eq. (10), Eqs. (6)-(9) may be reduced in terms of a new dependent variable C(x, T ) as The GITT can be easily applied to obtain the analytical solution. To do so, it is necessary to homogenize the boundary conditions because solutions of non-homogeneous problems based on eigen-function expansions may converge slowly or even exhibit anomalous behavior, especially in the vicinity of the boundaries (Ozisik, 1980;Almeida and Cotta, 1995;Cotta and Mikhailov, 1997;Liu et al., 1998). The boundary conditions are homogenized by using a filter function that satisfies the original boundary conditions. Hence, the following new independent variable is defined where Q(x, T ) is the solution of the following problem Q(x, t) has the same forms of boundary conditions as (13) and (14) but with the right side now set equal to zero for both equations. The initial condition of Q(x, t) becomes When solving problem using the GITT, a pair of transforms, namely an integral transform and an inverse transform, has to be established (Almeida and Cotta, 1995;Cotta, 1993;Liu et al., 1998;Suk,2013). It is therefore important to choose an auxiliary problem so that constructing the pair of transforms be simplified and that the solution converges for much less number of term. As detailed by Bharati et al. (2017) and Pérez Guerrero and Skaggs, (2010), an eigenvalue problem of Standard Liuville Problem (SLP) with self-adjoint second order operator must be used. The regular SLP chosen for this study is the same as that proposed by Bharati et al., (2017), Bharati et al. (2018) and Bharati et al. (2019), which has the particularity of making concentration converge for the first five number of summation. The selected ODE of the SLP is written as: which could be expressed in explicit form with a self-adjoint operator as: with the following associated two homogeneous first type boundary conditions considered to be the same as those associated with the ADE in Eq. (16) expressed as follow: The trivial solution of this problem is ϕ = 0. The nontrivial solutions are called the eigenfunctions belonging to each eigenvalues β m , and may be express as ( Bharati et al., 2017Bharati et al., 2019) where the eigen values β m are given by β m = m 2 π 2 l 2 − 3 4 and β m = 1, 2, 3, .... The orthogonality property for the set of linearly independent eigenfunctions ϕ m (x) reference to the weight function ρ(x) = e x , associated with Eq. (19) is given by: where N m is the norm and δ mk is the Kronecker delta. The norm and the normalized eigenfunctions Ψ m (x) are respectively given by: Now, the unknown function Q(x, T ) is represent as a series expansion in terms of the eigenfunctions Ψ m (x, T ) as: where θ m (T ) is the transformed "potential". Eq. (26) is the inverse transform rule. The corresponding transform rule is obtained by following the procedure of Cotta (1993) and Ozisik(1993), i.e. applying the operator l 0 ρ(x)ψ m (x)(•)dx to both sides of Eq. (16) and using Eq. (24) (the orthogonality property) and Eq.
(27) to obtain The pair of Eqs. (26) and (27) is called the integral transform pair. Substituting this solution in the ADE in Eq. (16), multiplying by e x Ψ m (x), integrating over the given domain, and using the orthogonality properties in Eq. (27), the result is a system of first order ordinary differential equations, a system of IVP, in matrix form, as: with the transform initial condition where elements of the M × M matrices A and B and the M-length vectors G are given by: The functions Ψ m (x) being orthogonal to each other then, A is a diagonal matrix with elements R. The Fourier coefficients in Eq. (28) is given by where Φ m (T, τ ) is the transition matrix (Liu et al., 2000) 8 The temporal variation of µ leading to that of the matrix B, it is therefore difficult to analytically evaluate the transition matrix Φ m (T, τ ) because it depends on B. This evaluation can be done numerically but in this work, the analytical solutions are needs. To get around this difficulty, Liu et al. (2000) have proposed to work with broad classes of problems for which the transition matrix can be significantly simplified so that the analytical solution, can be of practical value. Therefore, in the rest of the section we will discuss two cases of problem, the first one concern the transitional flow, dispersion coefficient and source term in heterogeneous porous media and the second one the steady state flow dispersion coefficient.

Transitional flow, dispersion coefficient and source term
In this sub-section, we consider a conservative solute (µ l = µ s = 0) with constant input concentration means f 1 (t) = C 1 and f 2 (t) = C 2 . For such case of problem, neglecting the decay coefficients, the matrix B does not depend any more on time, then the analytical evaluation of transition matrix is made possible. Under these assumptions, the coefficients in matrix Eq. (28) becomes, where: sin mπ l x e x/2 dx It should be noted that the expression of the Fourier coefficients in Eq. (34) is a particular form of that in Eq. (17) of Liu et al. (2000). In fact, the introduction of new time dependent variable in Eq. (10) for the case of transitional flow reduces the transition matrix into a simple exponential matrix and then simplified the form of the coefficients.

Steady state flow, dispersion coefficient, and decay
Another class of transport problems studied by Liu et al. (2000) for which the evaluation of the transitional matrix Φ m (T, τ ) can be significantly simplified concern the steady flow in heterogeneous medium. We consider a non-conservative solute with exponential decreasing input boundary condition means f 1 (t) = C 1 exp(−λt) and f 2 (t) = C 2 exp(−λt). Using Eq. (10), the new time variable in Eqs. (28), (32) and (33) will become the old time variable that is T = t. As mentioned by Liu et al. (2000), it is clear that due to the non-time dependence of coefficients, the matrix B(t) becomes constant and the transition matrix be reduced into a matrix exponential. In these cases, the expression of coefficients in Eq. (28) becomes where elements of the M × M matrices A and B and the M-length vectors G are given by: sin mπ l x e x/2 dx For each broad classes of problem studied, analytical solutions are obtained in Euclidean framework i.e for n = 1 and 2 , representing an index of the spatial dependence of the dispersion coefficient in the dispersion theory Dαu n as proposed by Freeze and Cherry (1979). According to this theory, the dispersion coefficient increases linearly with the distance for n = 1 hence, it increases as a square of a linear dependence space function for n = 2. For n = 1 , the dispersivity (ratio of dispersion coefficient and velocity) is constant delineating that the medium is homogeneous, and for n = 2 it increases linearly with distance, in this case the medium is heterogeneous. In the proposed solutions, for the transitional flow, the unsteady function is considered as exponential decreasing form (φ(st) = exp(−st)) and exponential increasing form (φ(st) = exp(st)) representing the seasonal time variation in South-East Asia. As Bharati et al. (2017), Bharati et al. (2018) and Bharati et al. (2019), we found the analytical solution convergent to the designer pattern with the first five terms (M = 5) of the Fourier series. Thus, the solution in Fourier series with first five terms may be written as The coefficients of the Fourier series are obtained for specific values of the input parameters.

Numerical results and discussion
The solute concentration attenuation pattern due to two continuous point sources is obtained from the proposed analytical solutions in the longitudinal direction of groundwater finite flow domain and finite temporal domain for two types of solute as mentioned in the previous section.
• In the case of transitional flow, the diffusion coefficient is negligible in comparison to the hydrodynamic dispersion coefficient that means the dispersion coefficient is directly proportional to the nth-power of the seepage velocity. For this case the length of domain and the time scale study are considered to be 0 ≤ x(km) ≤ 1 and 0 ≤ t(day) ≤ 1000. The solute concentration strength is evaluated from the analytical solutions given in Eq. (36) and the five time-dependent Fourier coefficients are given by Eq. (34). Analytical solution of the solute concentration is interpreted for the following set of input parameters (Bharati et al., 2017): heterogeneity parameter a − 2 = 0.5km −1 , uniform velocity u 0 = 0.0005km/day, uniform dispersion coefficient D 0 = 0.001km 2 /day, length of domain l = 1km, a − 1 = 1 , uniform input concentrations at the origin of the domain C 1 = 1kg/m 3 . The uniform input concentration in the end of the domain (x = 1) is set to be C 2 = 1kg/m 3 . The space variation function of the source term is considered as the same of Kumar et al. (2019) in which the parameter a is set to be equal to 0.2km −1 . Additionally, the initial background concentration in the liquid phase is supposed to follow an exponential decreasing distribution with distance i.e h(x) = C i exp(−αx) with C i = 0.1kg/m 3 and α = 0.01day −1 . The other input parameters are as follow: uniform source term q 0 = 0.02kg/m 3 day, initial uniform concentration in the solid phase K i s = 0.1kg/m 3 , unsteady coefficient s = 0.01day −1 . The distribution coefficients and the volumetric water content are chosen so that the retardation factor must be equal to 1.
• In the case of steady coefficients, a non-conservative solute with exponential decreasing input boundary condition is considered, meaning that f 1 (t) = C 1 exp(−λt) and f 2 (t) = C 2 exp(−λt) where λ represents the flow resistance coefficient. In this case, we considered a finite spatial and temporal domain defined as 0 ≤ x(m) ≤ 1 and 0 ≤ t(year) ≤ 7. The solute concentration strength is evaluated from the analytical solutions given in Eq. (36) and the five time-dependent Fourier coefficients are given by Eq. (35). The initial concentration distribution in the liquid phase is considered to follow the one proposed by [6] i.e h(x) = C i exp(−sec(α 1 x)). If the solute distribution coefficient is negligible (i.e. K d → 0), there is no interaction between the solute and soil, the retardation factor R becomes equal to water content θ. In this study, water content (θ) is assumed to be less than 1, so the retardation factor R becomes less than 1 for the case mentioned above, this indicates that only a fraction of liquid phase concentration participates in the transport mechanism. The input parameters values used are given by Singh and Kumari (2014) : C 1 = 1.0mg/L,C i = 0.01mg/L, K i s = 0.01mg/L, u 0 = 0.01m/year, D 0 = 0.01m 2 /year, D * = 0.002m 2 /year, f = 0.8, K d = 0.01, m = 0.01year −1 , µ l = 0.0027year −1 , µ s = 0.13year −1 , a 2 = 0.8m −1 . The uniform input concentration in the end of the domain (x = l) is set to be C 2 = 0.5mg/L. The space variation function of the source term is also considered as the same of  with a = 0.2m −1 but the uniform source term q 0 = 0.33mg/Lyear. Three geological formations are considered here with average porosity η and bulk density ρ b as follows (Manger, 1963;Freeze, and Cherry, 1979): η = 0.3 (sandstone),0.1 (shale), 0.5 (gravel) ; ρ b = 2.49 (sandstone), 2.39 (shale), 2.68 (gravel). Figures 1-4 illustrate the solute concentration distribution in the case of transitional flow with two uniform continuous sources introduced at x = 0km and x = 1km. Figure 1 depicts the concentration distribution pattern at different times (t = 20 , 100, 250 and 500 days from the introduction of the two uniform continuous sources respectively) along non-uniform flow through a homogeneous medium (n = 1) and heterogeneous medium (n = 2) for exponential decreasing and increasing form of velocity. For each curves, whatever the expression of φ(st), the concentration at the origin (x = 0) is 1 and that at the end (x = 1) is 0.5 for the two values of n at each time, that means for those particular positions, the concentration is in conformity with the uniform continuous source described by Eq. (8), and the zero concentration condition at the far end of the finite domain in Eq. (9). The curves depict that for each value of n at early times (t = 20days, in this example), the solute profile is qualitatively similar to that predicted with standard advection-dispersion models. The similarity ends quickly, in fact, after a certain time, the solute concentration increases with traveling distance from 1 to peak concentration, then decreases back to 0.5. The distance concerned by the solute build-up is very important and increases with increasing time and the peak concentration value gradually increases with increasing time. The solute build-up in a large distance is due to the concentration effect of the additional source term in the aquifer that decreases with distance and also with the presence of two contaminant sources. Furthermore, the build-up depends to the value of n and the expression of the unsteady function. In fact, a comparison of Figures 1(a)-1(d) show that the peak of concentration is higher for n = 1 than for n = 2 for the same unsteady function, however, for the same value of n the exponential increasing form of unsteady function has the higher value of concentration than the exponential decreasing form. It is also observed that the increasing rate of solute concentration is higher for the exponential increasing form of unsteady function than that for exponential decreasing form, the same feature is observe for the solute decreasing rate however those processes are important for n = 1 than for n = 2. It depicts also that for the two values of n, the rate of increase of the pollutant concentration is higher for the long-time period while the rate of decrease of the pollutant concentration is slower for the long-time period along the flow direction. However, for each time the concentration level for n = 1 is higher than that for n = 2 at any intermediate position, that means solute transport is faster with n. This result should be obvious because the expressions of the spatially dependent velocity and the dispersion coefficient given by Eq. (4) that the values of these parameters increase with the position additionally, the dispersion coefficient increases with the value of n, hence solute transport will be faster in a heterogeneous medium (n = 2) than that in a homogeneous medium (n = 1). The comparison between numerical and analytical solution shows that there is not great difference between them. Therefore, the analytical solution obtained with the GITT can be validated. Figures 2(a), 2(b) (2c), and (2d) represent the concentration profiles for three heterogeneous parameters a 2 (km −1 ) = 0.8, 0.5, and 0.05 at three different times t = 20, 100 and 500 days. The higher values of a 2 represents spatially varying dispersion coefficients and velocity components in a heterogeneous medium, whereas for the much lower value (0.005)the coefficients of the ADE may be treated as constant and hence the medium through which solute transport is occurring may be treated as a homogeneous medium. The curves are plotted with the same input values used in Fig. (1). The figures have the same behaviors as those of Figure (1) whatever the value of a, however, the concentration level in the medium at a particular position increases with the decreasing of the heterogeneity parameter. The difference in concentration level is more significant for long times period in comparison to short times but also for n = 2 in comparison to n = 1 . Let us remind for n = 1 the medium is homogenous whatever the value of a and this does not affect the dispersivity of the medium but just the effective velocity of the solute which depends on a . Whereas for n = 2, the dispersivity and the effective velocity of the solute both increase with a . On the other hand, the decrease in the concentration with the increase in the heterogeneity parameter is due to the fact that the velocity and the dispersion coefficient vary significantly which thus causes a significant attenuation of the concentration at a position than that in a medium of lesser heterogeneity. It is reveals also that the rate of increasing of the concentration is higher for small value of heterogeneity parameter while the rate of decreasing of concentration is slower for small value of heterogeneity parameter thus, the solute transport is faster as heterogeneity increases. The influence of heterogeneity parameter on the concentration distribution may be easily accessed by varying the value of parameter a . Figure 3(a), 3(b) depict the solute concentration values for different values of dispersion coefficient and velocity in homogeneous and heterogeneous media at t = 250days by using the spatial dependence parameters a 1 = 1 and a 2 = 0.5km −1 , the other parameters remaining the same as in Figure 1. It can be observed that for the same value of the velocity the concentration values at different positions decreases with the increasing value of the dispersion coefficient and it decreases with the increasing value of the velocity for the same value of dispersion coefficient due to the additional source term and the presence of two contaminants sources. The curves also show for the same dispersion coefficient, the rate of increasing of concentration is higher for slow velocity than for high velocity and the rate of decreasing of concentration is higher for high velocity in comparison to the slow velocity. Even though it is well known that in the presence of a source term, the solute transport process is slow for higher value of velocity in comparison to slow velocity, this is not explicitly observed here. Seeing the shape of the curves, we were tempted to say that the solute transport is faster for slow velocity in comparison to slower velocity at small distances and that this feature is inverse for the long distances, but, this is not the case. This behavior in solute concentration is due to the presence of the two sources of contaminant located at the origin and at the end of domain moving in the opposite directions. Here, we see once again another impact of two sources on concentration.
To better undestand the effect of the two sources on the solute concentration in the reservoir with transitional flow, we considered the two continuous sources localized at x = 0 and x = 1 with the same value of the injected concentration i.e. C 1 = C 2 = 1.0kg/m 3 . Figure 4 contains the concentration curves obtained for the two values of n and the two forms of steady velocity at two times t = 20days and t = 250days. The other inputs parameters remain the same as in Figure 1. An analysis of those curves shows that at lower times the concentration starts at 1.0 at the origin, then decreases with the traveling distance up to a certain value then increases with the distance until reaching the value of 1.0 at the end of the domain. For long times period, we observe the rather a contrary tendency, i.e. the concentration increases from the value 1.0 at the origin until a certain then decreases until reaching the value 1.0 at the end of the domain. Those features are in agreement with the boundary conditions. Additionally, it can be observed that the maximum of the concentration tends to be reached in the middle of the domain length, especially for n = 2 and for the exponential increasing form of velocity and therefore, the concentration curve tend to be symmetrical whatever the time. The tendency for the concentration curves to be symmetrical is here is due to the presence of two identical solute sources localized at the beginning and at the end of the domain. Figures 5-9 contain plots that illustrate the solute concentration pattern in the case of steady flow with two uniform continuous sources introduced at x = 0m and x = 1m in three different geological formations. Figures 5 and 6 show the solute distribution for the three types of geological formations at t = 1, 2, 4 and 6 years in homogeneous (n = 1) and heterogeneous (n = 2) sites, respectively. Contaminant concentration pattern is investigated for sandstone (ρ b = 2.49, η = 0.3), gravel (ρ b = 2.68, η = 0.5) and shale (ρ b = 2.39, η = 0.1) formations. We found that in each geological formation in both sites, at earlier time (t = 1year for this example) the contaminant concentration strength decreases with traveling distance but after certain period the strength of contaminant increases with distance up to a maximum value and reduces to 0.5 concentration level. As mentioned in the case of transient coefficients this build-up of concentration is due to the presence of the source term. Additionally, the contaminant concentration strength increases with the increasing value of the bulk density for increasing time period at each of the position. The curves also show that the contaminant distribution depends on the properties and the characteristics of different geological formation. In fact, it can be observed that the concentration level is higher in gravel with bulk density (ρ b = 2.68) compared to sandstone (ρ b = 2.49) and shale (ρ b = 2.39) at each of the position and time. In the same way, the concentration increasing rate is higher in gravel than the two other formations but the concentration level decreased more rapidly in shale. This does not mean that the solute transport is faster in gravel than in the two other geological formations at small distances and that it becomes slower in gravel than in other formation at large distances. Keeping in mind that there are two contaminant sources localized at the bounding of the domain, there are the ones that give this impression. However, the contaminant concentration level is higher in homogeneous site in compare to heterogeneous site in all three geological formations and each time. The analytical solutions in each case are also compared with the respective numerical solutions obtained using MATLAB pdepe solver. Even though certain curve analytical curves move away a little from the numerical ones, most of the figures show an agreement between analytical and numerical solutions at each time. In addition, the shapes of the curves of the analytical solutions are the same as those of the numerical solutions. The root Means square between the analytical and numerical solution is between 0.01 and 0.2. This agreement validates the analytical solutions obtained using GITT. Figures 7.a and 7.b elucidate the concentration pattern due to the source uniform term (q 0 ) in homogeneous (n = 1) and heterogeneous (n = 2) sites respectively, for different geological formations for a fixed time t = 4 years. The solute concentration is obtained with the same parameters values as in Figs (5) and (6). The dotted curves are drawn for the gravel formation, the dashed are drawn for the sandstone formation and the solid curves are drawn for the shale formation. We observed that the concentration strength at different positions increases with the increasing value of the uniform source term (q 0 ). But the concentration level increases much faster in homogeneous sites as compare in heterogeneous site with the same increment in the source term. It is also observed that for the smaller value of uniform source term, the concentration strength decreases with space and reduces to 0.5 concentration level at x = 1 m in both sites and geological formations. The solute build-up in both sites is observed from a certain value of uniform source term and increases with the increasing of the source term. Figures 8.a and 8.b provide the concentration values evaluated from the analytical solution in Eqs. (33) and (35) in the case of steady flow in homogeneous and heterogeneous media respectively for different geological formations by using the following distribution coefficient K d = 0.1, 0.05 and 0.01. The curves are plotted at fixed time t = 4 years.It should be noted that a variation of the distribution coefficient in this case leads to a variation of the retardation factor (R) and the effective first order decay coefficient (µ). Consequently, looking at the effect of the distribution coefficient on the concentration amounts to looking at the simultaneous effect of the retardation factor and the effective decay coefficient. The values of those parameters in the different geological formations are: for K d = 0.1 sandstone (R = 0.732, µ = 0.0567), gravel (R = 0.5912, µ = 0.0338) and shale (R = 1.4721,µ = 0.1769), for K d = 0.05 sandstone (R = 0.566,µ = 0.0297), gravel (R = 0.4956,µ = 0.0182) and shale (R = 0.936,µ = 0.0898). The curves show that the distribution (K d ) coefficient influences the concentration pattern in both sites. For examples, as observed, the strength of the contaminant and the contaminant level at different positions increase with the decreasing value of the distribution coefficient. But the concentration level increases faster in homogeneous site as compared to heterogeneous site with the same increment in the distribution coefficient. The solute concentration increases more rapidly for low value of the distribution coefficient in comparison to the high values and decreases more rapidly for high values of distribution coefficient in comparison to the low value. However, the difference on the contaminant strength for each geological formation is more important in the homogeneous sites in comparison to the heterogeneous sites. But the shale geological formation presents the more evident difference on the contaminant strength in comparison to the two other geological formations. This is due to the fact that the variation of the distribution coefficient causes a considerable variation of the retardation coefficient and the decay constant for the shale geological formation in comparison to the two other formations. The low porosity of the shale formation is the cause of that considerable 13 variation.

Example of application
Many contaminants sources are present in the environment and some of them are subject of studies due to their importance or the damage they can cause to living organisms. These contaminants sources are of various nature and origin, among which the radionuclides. Several radionuclides are members of a radionuclide decay chain. A lot of researchers have developed various models involving sets of advective-dispersive transport equations coupled by first-order decay (Van Genuchten, 1985 ;Suk, 2013 ;Chen et al., 2019;Yu et al., 2019) but analytical solutions in closed form that consider both production and decay at the source are scarcely available Moranda et al., 2018;Sanskrityayn et al., 2017). In order to illustrate the model, the following example of application has been considered: a transport of a second member of radionuclide decay chain during its movement into groundwater with two sources In the chosen example, the radionuclide source undergoes at the same time a decay and a production, so that the boundary functions in Eqs. (8) and (9) are written as where λ s and λ p are the decay and the production constant respectively, K is the global kinetic rate. When studying the distribution of radionuclides in a medium, we are most often interest in their environmental impact, especially in the effective dose that human being can absorb. On the basis of the study model, one can determine the radionuclide absorbing dose from ingestion of drinking water. The time histories of concentrations at different distances from the sources boundary must be used to assess the committed effective dose rates consumed by members of the general public through ingestion of drinking water. The committed effective dose per person from a given radionuclide decay chain through groundwater can be calculated by where IR is the rate of intake (m 3 /day), C is the radioactivity concentration in groundwater of the nuclide (Bq/m 3 ), and DF is the ingestion dose coefficient of the nuclide for the adult age group (Sv/Bq).
Input parameters, used here, are kept from Carntrell et al. (2003), ICRP (2012) Table 1. Figure 10.a depicts the radionuclide concentration as a function of time at different distances from the origin (x = 50, x = 100 and x = 100 m) for different geological formations. The concentration of 234 U at each position increases in the early time period and start decaying after a certain year depending on the geological formation and the position. However, the concentration of 234 U maintains the increasing trends up to t = 10000 years. The gravel formation has the higher radioactive concentration than the two other geological formation and the highest is observed for x = 200 m at t = 1000 years with the value of 2.1 × 10 9 Bq/m 3 . It is seen that the radioactivity concentration level decreases with position and then increases due to the second radioactive source localized at the end of the domain. Figure 10.b depicts the time history of radionuclide dose at different distances from the origin (x = 50, x = 100 and x = 100 m) for different geological formations. A comparison of Figs. 10.a and 10.b shows that the amplitude of the doses at different position and for each geological formation follows the same sequence of the magnitude of the corresponding concentration. As the radioactive concentration, the higher dose values are obtained for gravel formation for an adult located at x = 200 m from the origin. At that location, the highest 14 value of the amplitude of the dose is about 0.021Sv/year . The smallest value of the dose is obtained for the shale formation at 100 m of the origin. However, for each geological formation and at each position studied here, the annual doses are all above the WHO guideline of 0.1 mSv/year for the drinking water pathway. The generalized analytical solutions can quickly and accurately predict the one-dimensioal solute (as radionuclide) migration and assess the radiological impact posed by radionuclides in the environment as a result of leakage from a nuclear waste repository or accidental discharge from a nuclear facility.

Conclusion
An analytical solution of 1-D ADE solute transport with time and distance dependent coefficients in a sorbing finite groundwater reservoir with an additional source-sink term was derived. Two continuous solute sources were assumed localized at the origin and at the end of the domain. The velocity was considered as a linear function of space function while the dispersion coefficient was considered as a nth power of the velocity as stated in the dispersion theory. The ADE was solved with the help of the GITT using an advection-dispersion SLP with a new self-adjoint operator. Two types of problems were considered concerning the transient flow and the steady flow. The effect of some parameters on the equation was investigated. The results reveal that: • The source term produces a solute build-up with a large peak of concentration • The presence of two sources contaminants sources increases the concentration level.
• The concentration level becomes stable much faster in homogeneous medium in comparison with heterogeneous medium whatever the type of problem.
• For the two broad classes of the studied problems, the solute concentration distribution increases with time in various geological formations for homogeneous and heterogeneous sites. Also, the concentration decreases with distance at earlier time and at the elther time it increases with distance up to a maximum value and then starts to decrease until reaching the concentration value of the source localized at the end of the domain.
• The accuracy of calculated analytical contaminant strength is analysed with their corresponding numerical results obtained by MATLAB pdepe Solver, which were found in acceptable compliance with each other for the two broad classes of studied problems in homogeneous and heterogeneous sites, medium.
•The obtained analytical solution should be useful for estimating the transport of contaminant in heterogeneous and homogeneous, sorbing groundwater reservoir with two sources of contamination. Furthermore, it can be recommended as a tool for human risk assessment by drinking water as illustrated by the example of application studied.

Author Declaration
• Funding: None • Conflicts of interest: None • Availability of data and materials: No materials are used in this study. The data will be accessible always, some of them are available in the literature. •

Acknowledgment
The first author acknowledges support by the Center of Excellence of Information and Communication Technologies (C.E.T.I.C) and the Department of Physic of the University of Yaoundé 1. The content of this manuscript does not necessarily reflect the views of the agencies and no official endorsement should be inferred. Figure 1: Comparison solute concentration in presence of two sources localized at x = 0 and x = 1km for unsteady coefficients obtained from homogeneous (n = 1) and heterogeneous site (n = 2) for exponential increasing and decreasing form of velocity Figure 2: Effect of heterogeneity parameter on solute concentration in presence of two sources localized at x = 0 and x = 1km for unsteady coefficients obtained from homogeneous (n = 1) and heterogeneous sites (n = 2) for exponential increasing and decreasing form of velocity Figure 3: Effect of dispersion and velocity on solute concentration in presence of two sources localized at x = 0 and x = 1km for unsteady coefficients obtained for homogeneous (n = 1) and heterogeneous sites (n = 2) for exponential increasing and decreasing form of velocity Figure 4: Solute concentration pattern with the same value of the two sources localized at x = 0 and x = 1km for unsteady coefficients obtained for homogeneous (n = 1) and heterogeneous sites (n = 2) for exponential increasing and decreasing form of velocity Figure 5:Comparison of contaminant concentration distribution on different geological formation in homogenous medium (n = 1)in presence of two sources localized at x = 0 and x = 1m for steady coefficients Figure 6: Comparison of contaminant concentration distribution on different geological formation in heterogenous medium (n = 2)in presence of two sources localized at x = 0 and x = 1m for steady coefficients Figure 7: Effect of source term of contaminant concentration distribution on different geological formation in homogeneous and heterogenous medium in presence of two sources localized at x = 0 and x = 1m for steady coefficients Figure 8: Effect of distribution coefficient on contaminant concentration distribution on different geological formation in homogeneous and heterogenous medium in presence of two sources localized at x = 0 and x = 6m for steady coefficients Figure 9: Effect of parameter a 1 on contaminant concentration distribution on different geological formation in homogeneous and heterogenous medium in presence of two sources localized at x = 0 and x = 1m for steady coefficients Figure Figure 1 Comparison solute concentration in presence of two sources localized at x = 0 and x = 1km for unsteady coe cients obtained from homogeneous (n = 1) and heterogeneous site (n = 2) for exponential increasing and decreasing form of velocity Effect of heterogeneity parameter on solute concentration in presence of two sources localized at x = 0 and x = 1km for unsteady coe cients obtained from homogeneous (n = 1) and heterogeneous sites (n = 2) for exponential increasing and decreasing form of velocity Effect of dispersion and velocity on solute concentration in presence of two sources localized at x = 0 and x = 1km for unsteady coe cients obtained for homogeneous (n = 1) and heterogeneous sites (n = 2) for exponential increasing and decreasing form of velocity Solute concentration pattern with the same value of the two sources localized at x = 0 and x = 1km for unsteady coe cients obtained for homogeneous (n = 1) and heterogeneous sites (n = 2) for exponential increasing and decreasing form of velocity Comparison of contaminant concentration distribution on different geological formation in homogenous medium (n = 1)in presence of two sources localized at x = 0 and x = 1m for steady coe cients Comparison of contaminant concentration distribution on different geological formation in heterogenous medium (n = 2)in presence of two sources localized at x = 0 and x = 1m for steady coe cients Effect of source term of contaminant concentration distribution on different geological formation in homogeneous and heterogenous medium in presence of two sources localized at x = 0 and x = 1m for steady coe cients Effect of distribution coe cient on contaminant concentration distribution on different geological formation in homogeneous and heterogenous medium in presence of two sources localized at x = 0 and x = 6m for steady coe cients Effect of parameter a1 on contaminant concentration distribution on different geological formation in homogeneous and heterogenous medium in presence of two sources localized at x = 0 and x = 1m for steady coe cients