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

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.


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 steadystate 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 wellknown 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 .

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.Figure 1 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 : 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 where C x 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: It is important to notice that this boundary condition takes into account the fact that even if the source position is set in x q, the solute can reach the inlet section by dispersion, so C x 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 where C x 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, 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 Finally, the initial condition for the point pulse injection at y 0 and x q yields C x, y, 0 Mδ y δ x − q , 2.9 where δ represents the Dirac delta function.Note that M is the total mass per unit thickness of the chamber injected in the porous medium.

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 where the coefficients of the expansion are C x, y, t cos nπy l dy, for n 0, 1, 2, . . . .

3.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: Similarly, boundaries and initial conditions 2.5 , 2.8 , and 2.9 may be written as 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 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 r n 2 cos 2nz .

3.8
By using this definition, we obtain the final expression for C x, y, t : 3.9

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 longitudinal t 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 at x q presented by 9 In principle, for Pe ch T 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.x − q /u .Note that the series appearing in the Jacobi θ 3 function of Equation 18 converge very fast 28 for which is always satisfied for arbitrary positive t.
The concentration profile of Figure 2 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 Pe ch T greater than about 20.The influence of the lateral boundaries is much more evident at small Pe ch T .Small Pe ch T 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 scale t s adv 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 t s disp 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 between t s disp and t s adv , namely, Pe The Pe Figure 3 shows the longitudinal profiles of the concentration solutions at y 0 for Pe s L 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.
Several simulations allowed to verify that the effect of the upstream boundary condition is still significant in proximity of the source for Pe s L between 1 and 2, while negligible effects are present for Pe s L > 2 in any section.For experimental application, it is important to put in evidence that small Pe s L 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 .

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 dimension 200 × 280 × 10 mm 3 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 Figure 4.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 .Sodium fluorescent formula: C 20 H 10 Na 2 O 5 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 D L and D T 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 was u 2.9 × 10 −4 m/s and the ratio between q and L was 0.1786.Pe 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 Figures 5 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 t 1 375 s, t 2 400 s, and t 3 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.
From Figure 4 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 Figure 6 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.

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.

Appendix
In this appendix we give some details of the procedure necessary to get the solution of 3.3 subject to boundaries and initial conditions 3.4 , 3.5 , and 3.6 .The Laplace transform of a n is defined as where s is the frequency variable of the Laplace transform.By applying the Laplace transform to 3.3 and rearranging, we get for each n 0, 1, 2, . . ., where H x is the Heaviside function and W 2 n is defined as follows: , n 0, 1, 2, . . . .

A.6
The functions Q n s and P n s have to be determined by using boundaries conditions of A.3 and A.4 .Q n s can be computed by using the right boundary condition for the longitudinal domain A.4 .If we compute the limits for x → ∞, we get lim s → ∞ A n x, s lim By imposing in A.7 that the limit of A n is zero for x → ∞, the terms Q n s can be found: The terms P n s can be computed by using the upstream boundary condition 3.4 .By substituting A n x, s in A.3 and rearranging, we get A.9 The expressions for P n s can be found from A.9 : or, explicitly, The final expressions for A n x, s , n 0, 1, 2, . .., can be written as follows: A n x, s Z 0 x H q − x e x−q where we have set Z 0 x M l e u x−q /2D L .A.13 It is possible now to obtain the final solution by computing the inverse Laplace transform of A n x, s ; so doing, we get a n x, t a 0 x, t e −n 2 π 2 D l t/l 2 A.14 a 0 x, t M 2l 1 D L πt e − x−q−ut 2 /4D L t e − x q−ut 2 /4D L t− uq/D L − u D L e ux/D L erfc x q ut 2t D L . A.15

Figure 1 :
Figure 1: A schematization of the experimental apparatus and mathematical domain considered.L represents the length of the domain.

Figure 2 :
Figure 2: Effects of the lateral boundary conditions on our solution and comparison with the Bear one in a transverse cross-section at x/L 1 and t x − q /u Pe ch L 5 and Pe ch T 10 .

ch Ldisp l 2
/D L dispersion time scale and the advection time scale of the chamber t ch adv l/u is the longitudinal Péclet number, and it is defined as Pe ch L ul/D L where the superscript ch stands for chamber .The ratio between the transverse dispersion time scale t ch T disp l 2 /D T and the advection time scale of the chamber t ch adv l/u is the transversal Péclet number, and it is defined as Pe ch T t ch T disp /t ch adv ul 2 /LD T .

Figure 3 :.Figure 4 :
Figure 3: Effects of the upstream boundary condition on our solution and comparison with the Bear one in a longitudinal cross-section at y 0 and two different times Pe ch L .1 and Pe ch T 10 .

Figure 2
Figure2shows the concentration profiles of the two solutions for Pe ch T

Figure 5 :Figure 6 :
Figure 5: Comparison between analytical simulations and experimental data, for different times and fixed longitudinal cross-section at y 0.
. In principle, for Pe s L 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.
by residual minimization and least square best fit of Equation 25 and experimental data.Pe ch L and Pe ch T εd wm is the molecular diffusion coefficient of the solute in water and d wm ε is the porosity.D L and D T are the longitudinal and transverse dispersivity, respectively.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 secondtype boundary condition as follows: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 .