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

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 [18]. 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 [1114]. 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. Figure 1 shows a graphical representation of the problem: the water flow is along the coordinate, the length of the chamber is , and the chambers width is 2. 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 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]: where is the pore scale velocity, is the resident fluid solute concentration, and and are the longitudinal and transverse dispersion coefficients which are related to the pore scale velocity by [9] where is the molecular diffusion coefficient in a porous medium, where is the molecular diffusion coefficient of the solute in water and is the porosity. and 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 second-type 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]. The boundary condition for the inlet section 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 is expressed in terms of resident concentration as [25] where 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 , the solute can reach the inlet section by dispersion, so 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 resulting in where 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, , resulting in The second way is to assume an infinite chamber in the longitudinal direction ( 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 and yields where represents the Dirac delta function. Note that 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 .

Starting from this consideration, we consider a cosine Fourier series expansion for the solution, and we write where the coefficients of the expansion are

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) becomes

The expansion present in this equation gives rise to the Jacobi function whose definition can be found, for example, in the work of Abramowitz and Stegun [28] as By using this definition, we obtain the final expression for ():

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 longitudinal dispersion time scale and the advection time scale of the chamber is the longitudinal Péclet number, and it is defined as (where the superscript ch stands for chamber).

The ratio between the transverse dispersion time scale and the advection time scale of the chamber is the transversal Péclet number, and it is defined as .

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 presented by [9] In principle, for 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 and at and . Note that the series appearing in the Jacobi function of Equation (18) converge very fast [28] for which is always satisfied for arbitrary positive .

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 . The difference between the two solutions becomes less than 1% at greater than about 20. The influence of the lateral boundaries is much more evident at small . Small 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 (where the superscript stands for source) is related to the transport of the fluid from the inlet section to the source position and the dispersion time scale 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 and , namely,

The can be easily related to through . In principle, for the time scale of advection is smaller than the time scale of back dispersion and the upstream boundary condition does not influence the solution.

Figure 3 shows the longitudinal profiles of the concentration solutions at for . 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 between 1 and 2, while negligible effects are present for in any section.

For experimental application, it is important to put in evidence that small 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 ( 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 dimension  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 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: 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 () and the dispersion coefficients ( and ) 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  m/s and the ratio between and was 0.1786. and have been estimated by residual minimization and least square best fit of Equation (25) and experimental data. and 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 () at different times ( s,  s, and  s) and at a fixed transverse cross-section in the middle part of the longitudinal domain  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  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.

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 is defined as where is the frequency variable of the Laplace transform. By applying the Laplace transform to (3.3) and rearranging, we get where is expressed by (3.6). Initial and boundary conditions then become Equation (A.2) is a nonhomogeneous ordinary differential equation that can be solved with classical methods. The general solution of (A.2) is for each where is the Heaviside function and is defined as follows: The functions and have to be determined by using boundaries conditions of (A.3) and (A.4). can be computed by using the right boundary condition for the longitudinal domain (A.4). If we compute the limits for , we get By imposing in (A.7) that the limit of is zero for , the terms can be found: The terms can be computed by using the upstream boundary condition (3.4). By substituting in (A.3) and rearranging, we get The expressions for can be found from (A.9): or, explicitly, The final expressions for , , can be written as follows: where we have set It is possible now to obtain the final solution by computing the inverse Laplace transform of ; so doing, we get