# An Analytical Solution of the Advection Dispersion Equation in a Bounded Domain and Its Application to Laboratory Experiments

Authors: M. Massabó; R. Cianci; O. Paladino
Journal: Journal of Applied Mathematics
 (2011)

---

## Abstract

We study a uniform flow in a parallel plate geometry to model contaminant
transport through a saturated porous medium in a semi-infinite domain in
order to simulate an experimental apparatus mainly constituted by a chamber filled
with a glass beads bed. The general solution of the advection dispersion equation
in a porous medium was obtained by utilizing the Jacobiθ3 Function. The analytical
solution here presented has been provided when the inlet (Dirac) and the boundary
conditions (Dirichelet, Neumann, and mixed types) are fixed. The proposed solution
was used to study experimental data acquired by using a noninvasive technique.


---

## Body

## 1. Introduction

The contamination  of groundwater by substances of various kinds and the study of the behavior of compounds into natural or artificial porous media is a topic of growing interest. Laboratory flow experiments have been used in several works to study solute flow and transport phenomena at different spatial and temporal scales. Transmissive or reflective imaging techniques, in conjunction with dye tracer, allow to monitor the solute plume in a porous medium confined in a transparent box, satisfying the requirements of high spatial resolution and accuracy and achieving two or three orders of magnitude additional sampling points if compared to conventional analysis methods.A growing number of applications of imaging technique to investigate solute transport in porous media are present in the literature [1–8]. Experimental measurements are compared with simplified theoretical models, based upon advection-dispersion equation, and they show reasonable agreement [3, 4].The advection-dispersion equation is commonly used as governing equation for transport of contaminants, or more generally solutes, in saturated porous media [9]. Often the solution of this equation with particular boundary conditions requires the application of numerical methods. In other cases, where an analytical approach is possible, the solutions often deal with constant velocities. Many analytical solutions for constant velocity and dispersion coefficients and different boundary conditions are available. Lee in [10] offers exhaustive list of references and explanation of derivation techniques.Analytical solutions in semi-infinite domain with different initial and boundary conditions have been derived by several authors [11–14]. Extension to domain with finite thickness has been derived by Sim and Chrysikopoulos [15] and Park and Zhan [16]. Batu in [17, 18] proposes a two-dimensional analytical solution for solute transport in a bounded aquifer by adopting Fourier analysis and Laplace transform. Chen et al. [19] derive an analytical solution in finite thickness domain with zero  gradient boundary condition at the outlet of the domain and also consider a nonconstant dispersion coefficient. Analytical solutions for time-dependent dispersion coefficients have been derived by Aral and Liao [20] for a two-dimensional system and infinite domain. An analytical solution for the 2D steady-state convection-dispersion equation with nonconstant velocity and dispersion coefficients is given in [21]; the authors obtained the solution at steady state by the application of the Dupuit-Forchheimer approximation. Analytical solutions in cylindrical geometry, bounded domain and different source conditions are presented in [22]; these solutions, recovered by using a Bessel function expansion, have been used in order to estimate transport parameters [23, 24]. In this work an exact analytical solution of the advection-dispersion equation for the two-dimensional semi-infinite and laterally bounded domain is derived. The mathematical model wants to represent typical laboratory flow tank experiments where a tracer is injected in the porous medium and the domain is finite. The availability of an analytical solution taking into account the effects of the lateral borders and the influence of the upstream boundary condition in relation to the injection point can be necessary if one wants to use the whole set of experimental data acquired in the physical domain and if one wants to verify the experimental conditions for some fluctuations of the operative variables.The proposed analytical solution is discussed by means of comparison with the well-known two-dimensional infinite domain solution proposed by Bear [9] and with some limit analytical solutions. The influence of the boundary conditions on the solution is also discussed in terms of Péclet numbers, and some suggestions are given in order to choose the right analytical model depending on the experimental conditions. Finally, the analytical solution is compared with experimental data obtained from an experimental apparatus based on a noninvasive measuring technique, designed and constructed by some of the authors of [8].


## 2. Problem Formulation

Mass conservation of nonreactive solutes transported through porous media is described by a partial differential equation known as advection-dispersion equation. Here we consider the transport of a solute through a thin chamber filled with a homogeneous porous medium. Figure1 shows a graphical representation of the problem: the water flow is along the x coordinate, the length of the chamber is L, and the chambers  width is 2l. Fresh water is fed from the inlet of the chamber while a pulse injection of solute mass is initially provided from a source placed at  a distance q from the inlet. For thin chambers, where the thickness is much smaller than the other two dimensions in which the transport phenomenon occurs, the governing equation of solute concentration can be expressed by the two-dimensional advection-dispersion equation as follows [9]:(2.1)∂C∂t+u∂C∂x=DL∂2C∂x2+DT∂2C∂y2,
where u[L/T] is the pore scale velocity, C[M/L3]  is the resident fluid solute concentration, and DL[L2/T] and DT[L2/T]  are the longitudinal and transverse dispersion coefficients which are related to the pore scale velocity by [9]  (2.2)DL=dwmp+αLu,DT=dwmp+αTu,
where dwmp is the molecular diffusion coefficient in a porous medium, where dwmp≡ɛdwm is the molecular diffusion coefficient of the solute in water and dwmε is the porosity. DL and DT are the longitudinal and transverse dispersivity, respectively.Figure 1
A schematization of the experimental apparatus and mathematical domain considered.L represents the length of the domain.Initial and boundary conditions are to be set in order to solve (2.1). No solute flux across the lateral boundaries has to be imposed; the condition can be expressed as a second-type boundary condition as follows:  (2.3)∂C∂y|y=±L=0.
Concerning the longitudinal domain, two boundary conditions must be imposed. Exhaustive discussion of physical meaning of boundary condition is reported by Kreft and Zuber [25] and Parker and Van Genuchten [26]. The boundary condition for the inlet section (x=0) is represented by a mass balance equation, thus the solute flux across the left size of the section equals the flux from the right. The mass balance across x=0 is expressed in terms of resident concentration as [25](2.4)[uC]x=0-=[uC-DL∂C∂x]x=0+,
where Cx=0+ is the concentration of the inlet stream, which is zero in our case since the chamber is fed with pure water. In regard to the mass balance, the resulting boundary condition is then expressible as a mixed boundary value problem as follows:(2.5)[uC-DL∂C∂x]x=0+=0.It is important to notice that this boundary condition takes into account the fact that even if the source position is set inx=q, the solute can reach the inlet section by dispersion, so Cx=0- can be different from zero. Thus, the previous condition could represent a fresh groundwater that could be polluted by a spill located in a section relatively near to the inlet one, at a certain temporal instant. Similarly, mass balance should be imposed at x=l resulting in  (2.6)[uC-DL∂C∂x]x=l-=[uC]x=l+,
where Cx=l+ is the concentration at the outlet section, which cannot be known a priori, thus the problem cannot be solved unless simplifications are assumed. Two different ways exist in order to treat this issue, the most common in chemical engineering is to assume that concentrations are continuous across the outlet section [27], namely, Cx=l-=Cx=l+, resulting in  (2.7)∂C∂x|x=l-=0.
The second way is to assume an infinite chamber in the longitudinal direction (x axis for the problem considered). This assumption can be reasonably chosen when one is interested in the plume profile far from the outlet boundary condition. In this case we get a first type boundary condition expressed as(2.8)lim⁡x→+∞⁡C(x,y,t)=0.
Finally, the initial condition for the point pulse injection at y=0 and x=q yields(2.9)C(x,y,0)=Mδ(y)δ(x-q),
where δ represents the Dirac delta function. Note that M is the total mass per unit thickness of the chamber injected in the porous medium.


## 3. The Analytical Solution

The derivation of an analytical solution of (2.1) subject to boundary conditions (2.3), (2.5), and (2.8) and initial condition (2.9) is here illustrated. Due to the symmetry of both the transversal boundaries and the injection position with respect to the longitudinal axis, the solution of the PDE will be symmetric as well, resulting in C(x,y,t)=C(x,-y,t).Starting from this consideration, we consider a cosine Fourier series expansion for the solution, and we write(3.1)C(x,y,t)=12a0(x,t)+∑n=1+∞an(x,t)cos⁡⁡(nπyl),
where the coefficients of the expansion are  (3.2)an(x,t)=1l∫-l+lC(x,y,t)cos⁡⁡(nπyl)dy,forn=0,1,2,….It is important to notice that (3.1) satisfies the no-flux boundary condition (2.3). By substituting (3.1) into (2.1) we get a set of partial different equations for the coefficients:(3.3)∂2an(x,t)∂x2-n2π2DTl2DLan(x,t)-uDL∂an(x,t)∂x-1DL∂an(x,t)∂t=0.Similarly, boundaries and initial conditions (2.5), (2.8), and (2.9) may be written as  (3.4)[uan(x,t)-DL∂an(x,t)∂x]x=0+=0,(3.5)lim⁡x→+∞⁡an(x,t)=0,(3.6)an(x,0)=MLδ(x-q).The solution of (3.3) subject to (3.4), (3.5), and (3.6) is obtained by applying the Laplace transform. Details are reported in the appendix for completeness. By substituting (A.14) in (3.1) C(x,y,t) becomes  (3.7)C(x,y,t)=12a0(x,t)(1+2∑n=0+∞cos⁡⁡(nπyl)e-(n2π2DTt/l2)).The expansion present in this equation gives rise to the Jacobi functionθ3(z,r) whose definition can be found, for example, in the work of Abramowitz and Stegun [28] as  (3.8)θ3(z,r)=1+2∑n=1+∞rn2cos⁡⁡(2nz).
By using this definition, we obtain the final expression for C(x,y,t):  (3.9)C(x,y,t)=M4lθ3(πy2l,e(DTπ2/l2))(1DLπt[e((x-q-ut)2/4DTt)-(uq/DL)-e((x+q-ut)2/4DTt)]-uDLerfc⁡(x+q+ut2DLt)).


## 4. Influence of Boundary Conditions

The transport phenomenon into the chamber is mainly regulated by two time scales: the advection time scale related to the transport of the fluid from the inlet section to the outlet section and the dispersion time scale related to the dispersion of the solute from the source position in all the directions. The ratio between the longitudinaltLdisp(ch)=l2/DL dispersion time scale and the advection time scale of the chamber tadv(ch)=l/u is the longitudinal Péclet number, and it is defined as PeL(ch)=ul/DL (where the superscript ch stands for chamber).The ratio between the transverse dispersion time scaletTdisp(ch)=l2/DT and the advection time scale of the chamber tadv(ch)=l/u is the transversal Péclet number, and it is defined as PeT(ch)=tTdisp(ch)/tadv(ch)=ul2/LDT.The aim of this section is to evaluate the influence of lateral and upstream boundary conditions on the proposed solution. In order to illustrate this point the new solution has been compared with the solution of a 2D transport problem in a unbounded domain for a pulse injection point atx=q presented by [9](4.1)C(x,y,t)=M4πtDTDLe((x-q-ut)2/4DLt)-(y2/4DTt).
In principle, for PeT(ch)≫1 the lateral boundary condition has no influence on the solution since the transport of the solute by transversal dispersion is much slower than the transport by advection, and during the passage into the chamber the solute does not reach the borders. Figure 2 shows the concentration profiles of the two solutions for PeT(ch)=5 and PeT(ch)=10 at x/L=1 and t=((x-q)/u). Note that the series appearing  in the Jacobi θ3 function of Equation (18) converge very fast [28] for  (4.2)exp⁡(-DTπ2tl2)<1,
which is always satisfied for arbitrary positive t.Figure 2
Effects of the lateral boundary conditions on our solution and comparison with the Bear one in a transverse cross-section atx/L=1 and t=(x-q)/u (PeL(ch)=5 and PeT(ch)=10).The concentration profile of Figure2 roughly corresponds to the passage of the peak of the plume. The influence of the lateral boundary condition is, obviously, mainly evident close to the edge of the chamber and it is maximum at y/l=1. The difference between the two solutions becomes less than 1% at PeT(ch) greater than about 20. The influence of the lateral boundaries is much more evident at small PeT(ch). Small PeT(ch) means that the transverse dispersion coefficient is high or that the lateral dimension of the chamber is smaller than the longitudinal one. In this last case, the solution of our bounded two-dimensional analytical solution tends to the monodimensional solution.The influence of the upstream boundary condition is regulated by both advection and dispersion longitudinal time scales but defined taking into account the position of the source injection. The advection time scaletadv(s) (where the superscript s stands for source) is related to the transport of the fluid from the inlet section to the source position and the dispersion time scale tdisp(s) is related to the dispersion of the solute from the source position and so on, also back to the inlet section.We define another longitudinal Péclet number, related to the source position, as the dimensionless ratio betweentdisp(s) and tadv(s), namely,(4.3)PeL(s)=tdisp(s)tadv(s)=q2/DLq/u=uqDL.ThePeL(s) can be easily related to PeL(ch) through PeL(s)=PeL(ch)(q/L). In principle, for PeL(s)≫1 the time scale of advection is smaller than the time scale of back dispersion and the upstream boundary condition does not influence the solution.Figure3 shows the longitudinal profiles of the concentration solutions at y=0 for PeL(s)=0.1. The chamber longitudinal Péclet number is 100. It is possible to notice that, close to the inlet section, the profile of the new solution is different from the Bear one in both magnitude of the concentration (10% in the case proposed in Figure 3) and position of the center of mass. Both the former and the latter differences are induced by the presence of the upstream boundary condition which rebounds the solute back to the chamber. The influence of this boundary condition decreases with increasing the distance from the source injection point due to the longitudinal dispersion.Figure 3
Effects of the upstream boundary condition on our solution and comparison with the Bear one in a longitudinal cross-section aty=0 and two different times (PeL(ch)=.1 and PeT(ch)=10).Several simulations allowed to verify that the effect of the upstream boundary condition is still significant in proximity of the source forPeL(s) between 1 and 2, while negligible effects are present for PeL(s)>2 in any section.For experimental application, it is important to put in evidence that smallPeL(s) means either that the longitudinal dispersion coefficient is high or the velocity is small, or, equivalently, the injection point is very close to the inlet section (q/L small).


## 5. Application to Laboratory Experiments

The authors designed and constructed a laboratory scale flow chamber made of a Perspex box and filled it with transparent glass beads to study solute transport. The physical model represents a quasi-2-dimensional porous medium; indeed the thickness of the chamber is much smaller than the horizontal direction. The chamber is a Perspex box of internal dimension200×280×10 mm3 and is packed with glass beads of uniform diameter (1 mm).The solute is introduced by a pulse-like injection. The model is illuminated by a diffuse backlight UV source and uses the transmitted light technique to detect the dye emission (visible light) by a CCD camera; a scheme of the experimental apparatus is shown in Figure4. The acquired images are processed to estimate the 2D distribution of tracer concentrations by using a pixel-by-pixel calibration curve linking fluorescent intensity to concentration and taking into account the influence of the most common sources of errors due to the optical detection technique (photo-bleaching, nonuniform illumination, and vignetting).Figure 4
A schematization of the experimental apparatus and imaging set up.Sodium fluorescent (formula: C20H10Na2O5) has been chosen as tracer for experimental tests. It has a moderately high resistance to sorption [29]. Details of the experimental apparatus, illumination, imaging processing, and concentration estimation can be found in [8].The experimental campaign allows to obtain a complete set of fluorescein concentration data in the real domain. The dimension of the chamber and the solute mass are known, while pore-scale velocity (u) and the dispersion coefficients (DL and DT) are unknown.The proposed analytical solution has been used to simulate the experimental campaign. For the conducted experimental campaign the estimated pore-scale velocity (from flow-rate and porosity measurements) wasu=2.9×10-4 m/s and the ratio  between q and L was 0.1786. PeL(ch) and PeT(ch) have been estimated by residual minimization and least square best fit of Equation (25) and experimental data. PeL(ch) and PeT(ch) resulted to be, respectively, 87.3 and 37.5. Under these conditions, as stated in the previous paragraph, the lateral and the upstream borders should not influence the plume behavior at least at low times.In Figures5 and 6 the analytical simulations and the experimental data are compared, by considering the curves representative of the plume characteristics in the middle longitudinal cross-section (y=0) at different times (t1=375 s, t2=400 s, and t3=425 s) and at a fixed transverse cross-section in the middle part of the longitudinal domain (x=0.14 m) for the same three times previously chosen. Domains have been scaled.Figure 5
Comparison between analytical simulations and experimental data, for different times and fixed longitudinal cross-section aty=0.Figure 6
Comparison between analytical simulations and experimental data, for different times and fixed transverse cross-section atx=0.14 m.From Figure4 it is possible to notice that the real plume at longer times is more  advanced—in the longitudinal direction—than the analytical one. This happens in spite of using our analytical solution that gives, as discussed in paragraph 3, plume profiles whose center of mass is put forward with respect to the profile obtained with the Bear solution. This difference in position can be ascribed to the fact that the homogeneity of the porous medium is difficult to achieve in experiments, hence local heterogeneity of the hydraulic conductivity, which is not accounted for in the model, may produce local pore-scale velocity variation and hence the observed deviation of the center of mass.By observing Figure6 it is possible to notice again the effect of local medium heterogeneities (present in the physical domain of the experimental apparatus). Figure 6 is a transverse profile of plume concentration at different times; Figure 6 shows that the experimental plume has a slightly deviated shape with respect to the theoretical profile—in particular for t=425 s, the deviation may be ascribed to local heterogeneity that is responsible of plume meandering [24] and deviation of the transverse coordinate of the center of mass. However the fitting of experimental data and model solution is acceptable.


## 6. Conclusions

This study principally concerns the development of a new analytical solution of the advection-dispersion equation in 2D by taking into account a semi-infinite and laterally bounded domain and a point-like injection. The simulated domain is considered as constituted by a chamber filled with a porous medium, and the advective velocity is fixed as constant. The solution proposed, where both upstream and lateral boundaries are considered, shows a good applicability to real cases of typical flow tank experiments. In the second part of the present paper, a brief analysis of the new solution is given in terms of Péclet numbers, and some typical experimental situations are discussed.Finally, a comparison between the new solution and real experimental data is performed. Results show that the analytical solution correctly simulate the processes considered in all the physical domains chosen. With this analytical solution it will be possible to perform data analysis for parameters estimation using all the experimental data obtained by relatively long time experimental campaigns and also data taken near the injection point.



---

*Source: 493014-2011-01-12.xml*