Environ. Sci. Technol. 2009, 43, 2412–2416
Unified Multilayer Diffusion Model and Application to Diffusion Experiment in Porous Media by Method of Chambers G A N G L I U , †,§ L E E B A R B O U R , ‡ A N D B I N G C . S I * ,† Departments of Soil Science and Civil and Geological Engineering, University of Saskatchewan, Saskatoon, Saskatchewan, Canada S7N5A8, and Department of Soil and Water, China Agricultural University, Beijing, China 100094
Received June 16, 2008. Revised manuscript received January 14, 2009. Accepted January 15, 2009.
Diffusion coefficient is an important parameter for examining contaminant transport in the environment. Chamber methods (with or without external mixing devices) are the most popular methods for measuring effective diffusion coefficients in porous media (Deff) through air or water. The objectives of this paper were to apply simplified and unified analytical methods for both perfectly mixed and nonmixed (one- or two-) chamber systems and to examine how mixing affects the estimation of Deff. An analytical solution for a multilayer transient diffusion model was applied to the chamber methods without external mixing. By increasing the diffusion coefficient in reservoirs (D1 and D3) more than 10 times from the value for air or water (D0), the model was sufficient to approximate the well-mixed condition and, consequently, can be used to model transient diffusion in chamber systems with external mixing devices. We demonstrated that at long time Deff was related to the first eigenvalue (β1) of a transcendental equation, which provided a quick method for determining Deff accurately from experimental data. The error caused by using the well-mixed approximation can be significant for a single-chamber system when there are no external mixing devices. This error increased rapidly with decreases in the experimental duration. A good fit for the concentration versus time curve could not be obtained for the wellmixed solution, especially when sampling ports were near the boundary (x ) 0) and interface (x ) l1). The proposed solutions are useful when the reservoir or chamber methods areusedformeasuringDeff andhavewideapplicationsinpredicting contaminates transport in porous media and groundwater.
1. Introduction Diffusion is an important gas transport process in the vadose zone (1-3). Results of field studies (4) indicated that aqueous diffusion is an important, if not the dominant, transport process in fine-grained soils. When the mechanical dispersion coefficient is small, diffusion dominates the hydrodynamic mixing transport processes (5, 6). * Corresponding author phone: (306) 966-6877; fax: (306)966-6881; e-mail:
[email protected]. † Department of Soil Science. ‡ Department of Civil and Geological Engineering. § Department of Soil and Water. 2412
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 43, NO. 7, 2009
FIGURE 1. Schematic graph of diffusion apparatus: (a) case 1, two-reservoir diffusion apparatus; (b) case 2, one-reservoir apparatus with impermeable boundary; (c) case 3, one-reservoir apparatus with one constant concentration boundary. The key parameter to determine in characterizing gaseous or aqueous diffusion in porous media is the effective diffusion coefficient Deff. The most common laboratory method for measuring Deff consists of either one or two reservoirs (Figure 1). Each apparatus has its own advantages and shortcomings (2, 4), but in terms of diffusion theory all of these approaches are based on a multilayer diffusion process. Because of the difficulty in obtaining an analytical solution for such a multilayer system, most existing models use either a numerical simulation (5, 7-9) or the well-mixed assumption for reservoirs (2, 9, 16). Although some researchers (10, 14, 16) used mechanical mixing devices (e.g., fans), there are experiments that make no use of mechanical mixing at all (7, 9, 13). Only a few papers (8, 15) indicated that significant error can be introduced when the well-mixed approximation in the single-chamber system is used. Because of the wellmixed approximation (15) and numerical methodology (8), the conclusion drawn by them should be reexamined. Some analytical models of two-layer diffusion (15, 17, 18) have been developed; however, the complexity of these solutions hinders their application in studying the diffusion in experimental apparatuses (Figure 1). For example, it is difficult to use these models to obtain Deff from experimental data by curve-fitting. Furthermore, for a two-chamber system, there is no analytical solution derived on the basis of these models that can also handle the other two cases (as shown in Figure 1). The objectives of this paper were to apply a multilayer transient diffusion model to estimate Deff for both perfectly mixed and nonmixed one- and two-chamber systems and to understand how the well-mixed assumption affects the estimation of Deff.
2. Theory Analytical approximation (19), Laplace transform (18, 20), or integral transform pair technique (17, 21) has been used to treat the diffusion or convection of gas and solute in layered porous media. Since the inverse Laplace transformation of an N-layer system and the integral transform pair technique with N > 2 are difficult (17, 18, 21), we use a special and simplified integral transform technique, the orthogonal expansion technique (22-25). Orthogonal expansion is more effective and direct in solving partial differential equations 10.1021/es801657x CCC: $40.75
2009 American Chemical Society
Published on Web 03/03/2009
when their eigenfunctions are known (25). The detailed procedure of obtaining analytical solutions to study the heat conduction problem can be found in ref 22. By using the similarity between the heat conduction equation and the diffusion equation (23, 24), we have developed an analytical method for the problem of diffusion in chamber systems, assuming the diffusion coefficient within each layer is constant. Consider a one-dimensional porous medium with N parallel layers: the mathematical formulations for diffusion in this system are
[
]
∂Ci(x, t) ∂Ci(x, t) ∂ D ) Ri ∂x i ∂x ∂t
for i ) 1, 2, 3 ,..., N
(1)
where Ci(x, t) denotes the volume-based (20, 23, 24) aqueous contaminant concentration or gaseous-phase concentration in layer i at location x and time t, Di is the diffusion coefficient, and Ri is the retardation factor (26) for the ith layer. For chamber systems in Figure 1, eq 1 is subject to following boundary conditions (23, 24): ∂C1 ∂x
|
)
x)0
∂CN ∂x
|
x)xN+1
) 0 or
∂C1 ∂x
|
x)0
) CN|x)xN+1 ) 0 (2)
Boundary conditions at the interface x ) xi+1 are εa,iDi
∂Ci ∂Ci+1 ) εa,i+1Di+1 ∂x ∂x
for x ) xi+1
(3)
where εa,i is the air-filled porosity of the porous medium in layer i. The initial conditions for eq 1 are Ci(x, 0) ) Fi(x) for xi < x < xi+1 and i ) 1, 2, 3 ,..., N (4) where xi is the position of the ith interface. The general solution for Ci(x, t), in the ith layer, is as follows (a detailed derivation can be found in Supporting Information): ∞
Ci(x, t) )
∑C
n
exp(-βn2t)ψi,n(x) + C∞ for i ) 1, 2, 3 ,..., N
n)1
(5) where the summation is over all the eigenvalues βn and eigenfunction ψn. The coefficients Cn can be determined by the initial conditions. The βn values are the roots of the transcendental equation that satisfies the boundary conditions. Note that the last term in the above equation was ¨ zisik (22). mistakenly omitted in the equations (6.50) in O Here C∞ is the concentration at infinite time. In the following sections, we used the above-introduced method to study the chamber systems and special cases of layered systems for measuring Deff (Figure 1). We selected data from three published experiments (9, 13, 27) for comparison with our solution. We chose Ri(x) to be unity, because Ri(x) ) 1 was used in the three published papers and is a good approximation for most gases (2, 27) and nonreactive solutes such as dissolved organic carbon in an aquitard (9). For diffusion with adsorption, we can set Ri(x) to the desired value in the solution and the procedure is exactly the same (23). The two-chamber system approximates a single-chamber system (24) (Figure 1b) if one of the end chambers is extremely small. Van Rees et al. (13) used a single-chamber system, while the McIntyre and Philip method (10, 28) for measuring Deff of gas in the field corresponds to a large D1 in the source reservoir (with a fan for mixing) and l2 . l1. When l2 ) 2l1 and the porosity in the source reservoir is equal to that of the soil sample, our model reduces to that of the half-cell method (4). In addition, the apparatus (Figure 1c) that has
FIGURE 2. (a) Comparison between gas-phase oxygen concentration curves obtained by solution (eq 5) and FEM simulation. (b) Comparison between liquid-phase CCl2F2 concentration vs time curves obtained by solution (eq 5) and other solutions. been used to measure Deff of gases in soil and other porous media (29), and carbon dioxide efflux (30) is similar if the boundary conditions of Figure 1b at the outer end, x ) l2, are changed to a constant concentration. Therefore, all solutions for diffusion in chambers (Figure 1) were unified by our model. For diffusion equations with more general transient boundary conditions, position-dependent diffusion coefficients and production terms, readers may consult refs 22-24.
3. Results and Discussion 3.1. Solution Verification. The infinite series in eq 5 can be truncated to a desirable accuracy. In our calculations, 200 terms were used for eq 5 and the accuracy of relative error of C(x, t) was controlled at a level 2000 s for this example), the infinite series in eq 5 can be truncated with only the first term being taken into account to a desirable accuracy (24). The magnitude of each term in the infinite series is controlled, at long times (which depends upon several factors, such as l1, l2, l3, εa, and Φ), by the exponential factor. Since the second eigenvalue was at least 2 times larger than the first eigenvalue (Table 1; Figure 6), the first term in the series (eq 5) dominates. Thus, at long times, it is possible to write Ci(x, t) ) C1 exp(-β12t)ψi,1(x) + C∞ for i ) 1, 2, 3 ,..., N (7) Therefore, a plot of ln (|1 - Ci/C∞|) versus time should yield a curve that approaches a straight line with a slope of -β12. Figure 5b (β2/β1 ) 2.302) showed that the slope of the curve generated with 200 terms in eq 5 was -1.6358 × 10-4 and corresponded to β1fit ) 0.012 79 when the curve was fitted by eq 7. The deviation of the fitted β1fit from the true β1 )
FIGURE 5. (a) Simulated gas-phase Freon-13 concentration evolution curves for case 1 with Nseries ) 200 and Nseries ) 1, respectively. Linearized analysis of simulated ln (|1 - Ci/C∞|) vs time curves with Nseries ) 200 for (b) case 1 (gas-phase Freon-13) and (c) case 2 (DOC). The long-time domain is labeled, and the straight-line portion of curves can be used to fit β1 and Deff.
FIGURE 6. Typical values of β2/β1 for gas-phase and liquid-phase tracer diffusion experiments with different length of soil column (l2 - l1), D1, D3, Φ, and εa. Point A corresponded to the smallest value of β2/β1 in this figure. 0.012 77 was less than 0.2%. Substituting β1fit into the determinant of matrix equation of the transcendental equation and solving for Deff, through Mathematica (32), we can obtain the value of Deff. By this method, the fitted Deff/D0 of Figure 5b was 0.396 85 while the theoretical value was 0.397 00 (Table 1), and the relative error was 0.03%. Similarly, in Figure 5c (β2/β1 ) 2.138) and for case 2 at long time, the ln (|1 Ci/C∞|) versus time curve also had a slope of -β12. Hence, Deff of both one- and two-chamber systems can be accurately determined by fitting the slope of experimental ln (|1 - Ci/ C∞|) versus time curve in the long-time domain. The fact that only the first term of the eigenvalues was involved in the fitting procedure greatly simplifies the calculation. As can be seen from the solid, dashed, and dotted lines of Figure 6, if all other parameters were fixed, then the larger the D1/D0 (or D3/D0) ratio, the larger the ratio of β2/β1. The solid, dashed-dotted, and dashed-double-dotted line of Figure 6 also showed that decreasing εa corresponded to increases in β2/β1; thus, decreasing the value of Deff/D0 will lead to increasing β2/β1. Here, the two-parameter Millington
FIGURE 7. (a) Measured (b) and calculated (s) concentration versus time for case 2. (b) Calculated concentration for case 2 and the well-mixed model. (c) Error associated with use of the well-mixed approximation as a function of x/l1 for case 2. Measured values are from Figure 3 of Van Rees et al. (13). Liquid-phase CCl2F2 was used as tracer. and Quirk (36) model Deff ) εa2/Φ2/3 was used, where Φ is the total porosity. The Millington and Quirk model (36) gives the upper bound for the estimated Deff (31). The solid line (Figure 6) containing point A had εa ) Φ and thus the largest Deff and Deff/D0 possible. Therefore, the solid line (Figure 6) gave the lower bound for estimated ratio of β2/β1. Furthermore, point A had the minimum β2/β1 values () 2.302). In Figure 5b, we used the smallest ratio of β2/β1 from point A (Figure 6). Because of the exponential factor in eq 5, the deviation of β1fit from β1 for other points in Figure 6 will be less than point A, and thus a more accurate fitted diffusion coefficient Deff can be expected. Here we chose the most popular geometric size of the soil columns and chambers that were used in experiments (2, 9, 13, 14, 27) for simulation. We assumed that the length of the source reservoirs was equal to that of the collective reservoir. On the basis of Figure 5b,c, we concluded that for diffusion experiments with an apparatus size similar to that of Figures 5 and 6 and Table 1, only the first eigenvalue (β1) is needed to calculate the concentration evolution. Therefore, we only need eq 7 instead of eq 5. For an apparatus with other sizes of reservoirs and soil columns, β2/β1 g 2.0 can be satisfied by adjusting the geometric size of the apparatus; consequently, only β1 is needed in estimating Deff. In addition, the values of β2/β1 were independent of the gases or solutes (Figure 6) because varying the value of D0 (changing from CCl2F2 to DOC) alone will not cause any change to β2/β1. Therefore, β2/β1 is determined by the intrinsic attributes of the chamber system (l1, l2, l3, εa, and Φ) rather than by specific gases or solutes. 3.5. Influence of Sampling Location and Well-Mixed Approximation for Case 2. For the one-reservoir system (Figure 1b), the concentration versus time curves, as calculated from eq 5 at three different x positions, deviated from each other (Figure 7a). The curve with x ) 0.6l1 fitted the experimental data (13) rather well. However, as in Figure 3b, without knowing the sampling location, Deff estimated by the nonlinear curve fitting method (32) will have a large VOL. 43, NO. 7, 2009 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2415
error. This was further corroborated by Figure 7b, which demonstrated that, without the value of x given (Figure 7b), Deff can vary continuously between 1.78 × 10-9 and 1.98 × 10-9 m2 s-1 depending on the choice of x. The strong dependence of concentration evolution on the sampling location in the source reservoir was a result of deviation from the well-mixed condition. 3.6. Error of Well-Mixed Approximation for Case 2. On the basis of the results shown in Figure 4, it was obvious that the difference between the well-mixed approximation and the multilayer diffusion model was not negligible, especially when the ratio of D0 to D1 and D3 was less than 10. Meanwhile, Figure 7b showed that the experimental data could be fitted well by either the well-mixed solution or the pure diffusion solution. Thus, the potential error in assuming a well-mixed condition was evaluated. Since error associated with using the well-mixed approximation can not be analyzed based on the Van Rees et al. (13) model, our result was compared with that of Van Rees et al. (13) to estimate the error. Using the parameters from Van Rees et al. (13) and different experimental durations, we generated three curves using our model and fitted the curves with eq 6 to obtain Deff as shown in Figure 7c, The magnitude of the error in the fitted Deff reached a positive maximum at the far end of the source reservoir (x ) l1) and a negative maximum at x ) 0. The positive and negative maxima were 20.3 and -27.5%, respectively (Figure 7c) for an experimental duration of 1.8 × 105 s. The error increased rapidly when the experimental duration declined, which agreed with El-Farhan et al. (15), who used a zero concentration boundary condition at l2. This is consistent with the large differences between concentration profiles at different locations (Figures 3b and 7a). As a result, the pure diffusion concentration versus time curve could not be fitted well by the well-mixed solution, especially for sampling ports near the boundary (x ) 0) and interface (x ) l1). Hence, if the method of reservoirs was used for studying the diffusion in a porous medium, our model will give more accurate parameters such as Deff, which is especially critical when there is no mixing device available to mix the gas or solute in the reservoirs during the experiment.
Acknowledgments This work was partially supported by the National Science and Engineering Research Council of Canada (NSERC) and by China National Natural Science Foundation (Grant G40671085 to G.L.).
(9)
(10) (11) (12) (13)
(14)
(15)
(16) (17)
(18)
(19)
(20)
(21)
(22) (23)
(24) (25) (26)
Supporting Information Available Detailed derivation of eq 5. This material is available free of charge via the Internet at http://pubs.acs.org.
(27)
Literature Cited
(28)
(1) Nazaroff, W. W. Radon transport from soil to air. Rev Geophys. 1992, 30, 137–160. (2) Rolston, D. E.; Moldrup, P. Gas diffusivity. In Methods of soil analysis. Part 4. Physical Methods; Dane, J. H., Topp, G. C., Eds.; Soil Science Society of America: Madison, WI, 2002; pp 11131139. (3) Lee, E. S.; Birkham, T. K.; Wassenaar, L. I.; Hendry, M. J. Microbial respiration and diffusive transport of O2, 16O2, and 18O16O in unsaturated soils and geologic sediments. Environ. Sci. Technol. 2003, 37, 2913–2919. (4) Shackelford, C. D. Laboratory diffusion testing for waste disposal - A review. J. Contam. Hydrol. 1991, 7, 177–217. (5) Ball, W. P.; Roberts, P. V. Long-term sorption of halogenated organic chemicals - part 2. Intraparticle diffusion. Environ. Sci. Technol. 1991, 25, 1237–1249. (6) Tokunaga, T. K.; Wan, J.; Pena, J.; Sutton, S. R.; Newville, M. Hexavalent uranium diffusion into soils from concentrated acidic and alkaline solutions. Environ. Sci. Technol. 2004, 38, 3056–3062. (7) Rowe, R. K.; Booker, J. R. POLLUTE Users Guide, 6th ed.; GAEA Environmental Engineering Ltd.:Whitby, Ontario, Canada, 1997. (8) Hutchinson, G. L.; Livingston, G. P.; Healy, R. W.; Striegl, R. G. Chamber measurement of surface-atmosphere trace gas exchange: Numerical evaluation of dependence on soil, interafacial 2416
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 43, NO. 7, 2009
(29) (30)
(31) (32) (33) (34) (35) (36)
layer, and source/sink properties. J. Geophys. Res. 2000, 105, 8865–8875. Hendry, M. J.; Ranville, J. R.; Boldt-Leppin, B. E. J.; Wassenaar, L. I. Geochemical and transport properties of dissolved organic carbon in a clay-rich aquitard. Water Resour. Res 2003, 39, 1194; doi 10.1029/2002WR001943. McIntyre, D. S.; Philip, J. R. A field method for measurement of gas diffusion into soil. Aust. J. Soil Res. 1964, 2, 133–145. Reible, D. D.; Shair, F. H. A technique for the measurement of gaseous diffusion in porous media. J. Soil Sci. 1982, 33, 165–174. Glauz, R. D.; Rolston, D. E. Optimal design of two-chamber, gas diffusion cells. Soil Sci. Soc. Am. J. 1989, 53, 1619–1624. Van Rees, K. C. J.; Sudicky, E. A.; Rao, P. S. C.; Reddy, K. R. Evaluation of laboratory techniques for measuring diffusion coefficients in sediments. Environ. Sci. Technol. 1991, 25, 1605–1611. Petersen, L. W.; Rolston, D. E.; Moldrup, P.; Yamaguchi, T. Volatile organic vapor diffusion and adsorption in soils. J. Environ. Qual. 1994, 23, 799–805. El-Farhan, Y. H.; Petersen, L. W.; Rolston, D. E.; Glauz, R. D. Analytical solution for two-region diffusion with two well-mixed end chambers. Soil Sci. Soc. Am. J. 1996, 60, 1697–1704. Jin, Y.; Jury, W. A. Characterizing the dependency of gas diffusion coefficient on soil properties. Soil Sci. Soc. Am. J. 1996, 60, 66–71. Liu, C.; Ball, W. P.; Ellis, J. H. An analytical solution to the onedimensional solute advection-dispersion equation in multi-layer porous media. Transp. Porous Media 1998, 30, 25–43. Yates, S. R.; Papiernik, S. K.; Gao, F.; Gan, J. Analytical solutions for the transport of volatile organic chemicals in unsaturated layered systems. Water Resour. Res. 2000, 36, 1993–2000. Wu, Y. S.; Kool, J. B.; Huyakorn, P. S. An analytical model for nonlinear adsorptive transport through layered soils. Water Resour. Res. 1997, 33, 21–29. Liu, C.; Ball, W. P. Analytical Modeling of Diffusion-Limited Contamination and Decontamination in a Two-Layer Porous Medium. Adv. Water Resour. 1998, 21, 297–313. Liu, C.; Szecsody, J. E.; Zachara, J. M.; Ball, W. P. Use of the generalized integral transform method for solving equations of solute transport in porous media. Adv. Water Resour. 2000, 23, 483–492. ¨ zisik, M. N. Heat Conduction, 2nd ed.; Wiley-Interscience: New O York, 1993. Liu, G.; Si, B. C. Analytical modeling of diffusion in layered systems with position-dependent diffusion coefficient under linear approximation. Adv. Water Resour. 2008, 31, 251–268. Liu, G.; Si, B. C. Multi-layer diffusion model and error analysis applied to chamber-based gas fluxes measurements. Agric. Forest Meteorol. 2009, 149, 169–178. Arfken, G.; Weber, H. Mathematical Methods for Physicists, 5th ed.; Harcourt Academic Press: San Diego, CA, 2001. Jury, W. A.; Spencer, W. F.; Farmer, W. J. Behavior assessment model for trace organics in soil, 1, Model description. J. Environ. Qual. 1983, 12, 558–564. Ball, B. C.; Harris, W.; Burford, J. R. A laboratory method to measure gas diffusion and flow in soil and other porous materials. J. Soil Sci. 1981, 32, 323–333. Elberling, B.; Nicholson, R. V. Field determination of sulphide oxidation rates in mine tailings. Water Resour. Res. 1996, 32, 1773–1784. Currie, J. A. Gaseous diffusion in porous media, Part 1 - A nonsteady state method. Br. J. Appl.Phys. 1960, 11, 314–317. Jassal, R. S.; Black, T. A.; Novak, M. D.; Morgenstern, K.; Nesic, Z.; Gaumont-Guay, D. Relationship between soil CO2 concentrations and forest-floor CO2 effluxes. Agric. For. Meteorol. 2005, 130, 176–192. Liu, G.; Li, B. G.; Hu, K. L.; van Genuchten, M. T. Simulating the gas diffusion coefficient in macropore network images: Influence of soil pore morphology. Soil Sci. Soc. Am. J. 2006, 70, 1252–1261. Wolfram, S. The Mathematica Book, 3rd ed.; Wolfram Media: Champaign, IL, 1996. Carslaw, H. S.; Jaeger, J. C. Conduction of heat in solids, 2nd ed.; Clarendon Press: Oxford, U.K., 1959. Campbell, G. S.; Norman, J. M. An Introduction to Environmental Biophysics, 2nd ed.;Springer-Verlag: New York, 1998. Detwiler, R. L.; Glass, R. J.; Bourcier, W. L. Experimental observations offracturedissolution:TheroleofPecletnumberonevolvingaperture variability. Geophys. Res. Lett. 2003, 30, 1648. Millington, R. J.; Quirk, J. M. Permeability of porous solids. Trans. Faraday Soc. 1961, 57, 1200–1207.
ES801657X
