Environ. Sci. Technol. 1995,29, 2725-2734
Solute Transpott in Porous Media with Sorption-Site Heterogeneity WENLIN CHEN' AND ROBERT J. WAGENET* Department of Soil, Crop, and Atmospheric Sciences, Cornell University, Ithaca, New York 14853
A model for describing solute nonequilibrium transport influen ced by so rption-site heterogeneity in porous media was proposed and analyzed in terms of its time moments. The sorption-site heterogeneity was conceptualized as variable sorbing fractions/ compartments that are grouped into classes according to a probability density function (pdf) for both equilibrium (linear partitioning, KD)and rate (first-order sorption time, &) parameters. Two typical pdfs, a simple two-site (TS) and a general y site (GS),were employed for comparison. The correlation between KD and T, was considered in t w o extreme cases: perfectly correlated as the linear free-energy relationship (LFER, linear log & log KO)holds or completely independent. As the LFER holds, the impact of sorption rate statistics presents on all moments (except the zeroth) of a solute breakthrough curve (btc), resulting in simultaneous enhanced peakedness and long time scale tailing on the btc. When KD and T, are uncorrelated, equivalency between the simple TS and the more general GS models can be established up t o the third btc moment. However, higher moments immediately become highly sensitive to the rate sorption statistics, rendering a substantial deviation between TS and GS in the description of the long time scale behavior. The model deviation can be enhanced in the existence of the LFER since the ultimate equivalence of the btc moments in this case can only be achieved up to the second. It is indicative in either case, therefore, that when sorption heterogeneity exists, simple distributions such as TS is inherently not sufficient to represent a more generally distributed sorption process (e.!., GS), which can even become more stringent as the time scale increases.
-
Introduction The movement of contaminants in natural porous media (soilslaquifers),subject to rate limitation of mass transfer between the mobile fluid and the regions of immobilization * Corresponding author e-mail address:
[email protected]. Present address: AgrEvo Research Center, 703 NOR-AM Rd., P.O. Box 538, Pikeville, NC 27863. t
0013-936X/95/0929-2725$09.00/0
0 1995 American Chemical Society
(sorbing solid phases and/or stagnant fluid in intraaggregate/particle pores), is increasingly the focus of research. Conventional approaches to deal with the ratelimited transport problem often assumethat natural porous media are homogeneous or can be represented as an equivalent homogeneous medium so that every ratelimiting process can be described in a deterministic manner. However, it is evident that this assumption is not always valid. First, physicallchemicalproperties of a soillaquifer material can vary at scales ranging from microscopic solid grainslinterstitial pores to macroscopic fieldslwatersheds ( 1 ) . Second,althoughan equivalenthomogeneous medium can sometimesbe defined as, for example, uniform spheres having a single radius averaged over a bulk sample (21, the impact of variations in the sorbent microscopic characteristics (e.g., sorbent structure and composition) on the rate-limiting process and transport can be significant, especially at both ends of the time scales experienced (early and prolonged stages) (3, 4). Heterogeneity at large scales (field and regional) has been investigated by correlating the equilibrium or rate sorption properties with the random hydraulic conductivity field. Jury (5) used a stochastic convective approach to investigate the influence of linear equilibriumsorption on transport by correlating the partition coefficient with the random velocity field. A similar approach was used by van der Zee and van Riemsdijk (a,who assumed sorption to be a nonlinear equilibrium process. More recently, this approach was further extended to rate-limiting cases ( 7 ) . Although the equilibriumor rate parameters in these largescale studies were correlated with the spatially variable hydraulic conductivity field, heterogeneity deriving from variable sorbent microscopic structure and composition was not considered. Heterogeneityat the sorbent grain scale has been amajor topic in most laboratory batch and column studies, primarily motivated by the observation that the entire concentration-time profile (sorptionldesorption or breakthrough curves, btc) cannot be sufficiently described by a model assuming a single uniform equilibriumor rate process (one site/compartment). A simple heterogeneous sorption model, the two-site approach, was subsequently proposed and found to be in better agreement with experimental observations (8-10). In this approach, two distinct sorbingfractions, equilibrium(instantaneous) and rate-limiting (time-dependent),were assumed to coexist in the soillaquifer material. While the instantaneous fraction commonly has been assumed to represent the equilibrium sorption onto the sorbing solid, quite different mechanisms have been proposed for the time-dependent fraction. Early versions of t h i s approach represented timedependentsorption as first-orderkinetics ( 111, The physical significanceofthe first-order kinetics assumptionwas questioned later by Wu and Gschwend (12),who proposed that the time-dependent process was controlled by pore diffusion (intrasorbent)within microporous particles retarded instantaneously by pore-wall sorption. The assumption of retarded intrasorbent diffusion was applied by Ball and Roberts (13)to the time-dependent fraction. An alternative description of the time-dependent fraction invoking the
VOL. 29, NO. 11. 1995 / ENVIRONMENTAL SCIENCE & TECHNOLOGY 12725
PSO Pore water (e) H
K Flow
having uniform local sorption properties. The sorbing mechanisms at each site are presumed to be either an intrasorbent/organic diffusion-controlled process or a combination of diffusion and sorption kinetics. It is also worthwhile to realize that the intrasorbent pore diffusion may still be present even though the immobile fluid within the sorbent micropores does not significantly contribute to the reduction of the mobile space (as assumed in Figure 1). Such is the case, for example, when micropore space is small but pore-wall retardation is large. Nevertheless,the equation of mass conservationfor the specified solute particle, neglecting local dispersion and molecular diffusion temporally, is given by ac pa% -+ - + u ( x ) - = ac o at e at ax
/
Soil surface (S) Ts (=I/K) FIGURE 1. Schematic illusuatiom d tLn solute flow and sorption in a porou8 medium witl, distributed sorpion hadions described br a continuous or a two-site discrete pdf. diffusion mechanism has been the diffusive partitioning of sorbate into organic matter (IO,14). Althoughthetwo-sitesorptionheterogeneitymodel with various mechanistic assumptions for the time-dependent fraction has resulted in improved description of measured concentration-time profiles obtained from sorption/desorption batch or packed column experiments. recent studies have furthershownthat thetwo-siteapproachoften fails to predict concentration changes at long time scales (4. IO). A general description of a population of sorption sites to extrapolate long time scale sorption behavior, therefore. is not only of theoretical concern but also of practical significance. The objective of this paper was to generalize the two-sitesorption heterogeneitymodel under transport conditions. with a probability distribution function describing both equilibrium and rate parameters in the first-order massa;lnsferprocessapproldmatedatalocal scale. The approach also employed a residence time distribution (RTD)to describe hydrodynamic transport. The methodology is general in the sense that the heterogeneity in hydraulic properties may also be incorporated using various RTDs of specific interest. With the presence of sorption-site heterogeneity. nonequilibrium conditionsfor a given dispersive flow system were examined. In the two subsequent papers of this framework. experimental data will be provided to further elucidate various model assumptionsand to show the flexibilityand limitations of the developed model by exuapolatingit toaseriesofherbicide leaching experiments accomplished at a variety of time scales.
Residence lime bistributioa for Hfllndpamic Transport Residence Time of a Single Solute Particle. A reactive solute particle moving through a porous medium in the meanflowdirectionxisillustrated in Figure 1. The porous medium is composed of a mobile fluid phase and a continuum of sorbing (solid)regions, sketched as a dashed arrow and different shaded patterns, respectively. The solute particlecan besorbed either instantaneouslyortimedependently onto a fraction of the sorbing sites during transport. Note that the term “site” here does not necessarily mean a distinct physical point but a sorbing domain or portion of the sorbent surface envisioned as a fraction 2726. ENVIRONMENTAL SCIENCE 8 TECHNOLOGY I VOL. 23. NO. 11.1995
(la)
where cis time: x is distance; cis concentrationin the mobile fluid u is pore velocity; sb is mass sorbed at a given sorbing site: p is medium bulk density; and 0 is porosity. All the symbols introduced here and later are summarized in the Glossary. The sorbed mass s b can generally be described by a convolution integral:
where h(t) is a density function describing the process of mass withdrawal and return between the mobile fluid and the sorbingregions and is also a function with flexibilityto be further specified as a first-order mass transfer or a diffusion process in later work. Note that, although local dispersion and molecular diffusion have been omitted in eq la, the macroscopic dispersiveeffect is retained through the ensemble approaches described later for the random velocity field (13. Introducingdimensionless variables, X= x / b ; T = t / & CI = c/G;S = psb/(eG);and V = vfolb,reduces eq 1 to the dimensionlessform:
as v(x)-= -+ -+ aT aT ax
0
and S(TX) = f i , ( T X ) h ( T - r) dT
(2b)
wherethereference p a e t e r s b , h,andGcanbechosen as the column length, mean residence time, and input concentration, respectively. Detailed mathematical derivations in the following sections are given elsewhere (16). Under zero initial conditions of c, and S and an intantaneous injection of unit mass at X = 0 cl(TJ = 0) = d+(n,where d+ is a Dirac function defined from the positivedirection;eq2canbesolvedintheLaplacedomain:
tl = exp(-st[l
+ h(s)l)
(3)
where s is the Laplace parameter: the symbol denotes that a variable has been Laplace transformed, and “1”
I d X “.LE
(4)
h(s) = 516,
(5)
Note that r in eq 4 is the residence time of the solute that was released at the inlet X = 0 and reached the outlet X =
1 with microscopic pore velocity Vunder no influence of any rate-limiting processes. The pore velocity V may be envisioned here, in the microscopicpoint ofview, as a mean projection of the three-dimensionally oriented velocity along the particle trajectory onto the Xcoordinate (71.h(s) in eq 5 is also called the transferfunction for the ratelimiting processes ( 1 7 ) . Detailed forms of h(s) will be discussed and specified in later sections. Statistical Ensemble Approach. Consider now a large number of solute particles that are uniformly and simultaneously applied to the inlet surface at X = 0. Due to the geometric irregularity of the pores composing the mobile phase, the microscopic velocitywith which each individual solute particle travels is variable. The velocity for a given solute particle may also change in a random manner during the trip from inlet to outlet due to the lateral shifting among different velocities on the plane normal to the mean flow direction X. As a result, the residence time of each solute in the porous medium, defined in eq 4, becomes a random variable. An extensive discussion of the residence time (or equivalently velocity) statistics and its mechanistic implications are not a subject of the current work, but are given by Beran (18). Here, we simply adopt the asymptotic behavior of the solute transport in the porous medium based on the central limit theorem. That is, all solute transport in the porous medium will eventually approach a Gaussian distribution in space whenever the transport scale is much longer than that of the velocity correlationlength (18).With this condition, the probability density function ( p a for the residence time t is inverse Gaussian ( 1 5 ) ; a pdf for positive-valued variates (19):
where 6, is the coefficient of variation (0 of the residence time, t,and p,(t;X = 1) d t is the probability for a particle to have residence time between z and d t before exiting the system at X = 1. Note that the mean residence time in dimensionless form is always unit. By knowing the RTD, p,(t; X), the expected output concentrationcan be calculated through eq3 in the Laplace domain:
where denotes the ensemble mean of CI in the Laplace domain. By substituting eqs 3 and 6 into eq 7, the final result is
Note that the inverse transform of eq 8 provides an expected concentration output across the exit plane at X = 1 under the correspondingDirac (impulse) input boundary conditions in a porous system, which is an equivalent btc in column studies. In particular, it can be shown that the relationship between the CV (eT) and the column Peclet number (PI is P = 2/€,2, rendering the inverse Gaussian RTD approach identical to the conventional convectiondispersion equation (15,167.The inverse Laplace transform in this workwas allperformed numericallywiththe Fortran
codes provided by Jury and Roth (20). The program is available upon request from the authors.
Distributed Rate-Limiting Process In the previous section, a convolution integral (eq 2b) was used to represent a general rate-limiting process, and the corresponding transfer function, hb), remained undetermined. Realizing that mechanismsgoverningrate-limiting can be diffusive (retarded intrasorbentIintra-organic)or a combinationof diffusive and kinetic for polarlweakly polar compounds, we here adopt the commonly used first-order mass transfer approximation at the local scale (a single sorption siteIregion). This simplifiesthe complicationsof sorbent geometric specification required in the diffusion formulation. The rate coefficient in the first-order formulation, therefore, inherently lumps all the local geometric influence as if a diffusive'mechanism prevails (21). The first-order equation is written as
where Kis the dimensionlessrate coefficient, defined as K = K ~ OK, is the rate coefficient; and KDis the dimensionless linear partition coefficient between fluid and solid at sorption equilibrium, defined as KD= &pIO. Taking the Laplace transform of eq 9 gives (10) where T, = 1IKis the dimensionlesssorption characteristic time. Defining T, as a reciprocal of K has the advantages of converting the large magnitudes of K to small values, giving a convenient representation of instantaneous reactions (i.e., Kat infinity corresponds to T, = 0). It should be noted that, although the transfer function (eq 10) of the rate-limitingprocess is derived from the local concentration cl, it remains the same at the ensemble concentration level (Le., eq 8) since the rate-limiting process is not treated as correlated to the random velocity field. While soillaquifer materials may be considered chemically homogeneous on the basis of conventional macroscopic properties, it is also true that heterogeneity in soil physical and chemical properties that are related to ratelimitingprocesses are variable at almost every microscopic scale (22).When intrasorbent/organic diffusion dominates as the rate-limiting process, any variations in the sorbent microgeometry (size and shape) can have a crucial impact on the sorption process over a wide range of time scales. Variations in sorbing components, the organic matter distributionamongvarious particle-size fractions, and even the variations within the soil aggregatedaquifer grains of the same size fraction may cause considerable fluctuations of the local rate process as the sorbate diffuses through the microporeslorganic matrices and/or sorbs onto the micropore walls. In the proposed first-order approach, these variations are lumped into the two parameters: equilibrium partitioning coefficient (KD)and sorption time constant (Ts). It is reasonable, therefore, to treat both KDand T, as randomly distributed variables. As a result, the transfer function of eq 10 is evaluated by taking a mathematical expectation as follows: VOL. 29, NO. 11,1995 /ENVIRONMENTAL SCIENCE &TECHNOLOGY
2727
ou=u
= 1 in eq 15 to retain the equilibrium site(s). Moreover,
fi and& serve as the fractions of both sites, andfi = 1 A. where H(s)is the mean transfer function for the distributed first-order mass transfer; and pks is the joint pdffor KDand Ts Since a linear correlation (linearfree-energyrelationship, LFER) is often found between the two log-transformed parameters, In KO and In K (231,the LFER in terms of Ts on linear scale is then likely
KD = UT,b
(12)
where a and b are two positive constants (because In KD and In K are often negatively correlated). Interestingly, two extreme cases of correlation are represented by eq 12 as b takes zero or nonzero values. That is, if b = 0, a completely independent relationship exists between KDand 7', (in this case KD = a, and a is then distributed as pk); whereas, any non-zero value of b gives a perfect correlation between KDand T, (Le.,given the value of T,, KDis completely determined by eq 12). Perfect correlation would simplify the bivariate pdf (pk,) to univariate (I, 24) and reduce eq 11 to a single integral. The resulting univariate pdf can be represented either by the pdf of KD (pk) or by the pdf of T, (p,), In the case b = 0, however, the bivariate pdf pks simply becomes equal to the product of the two marginal pdfs of KD and T,, denoted also as f i and p,, respectively. In this case, the first integral on KDin eq 11 would simply be equal to the mean of KD, thus rendering it identical to the assumption that KD is constant at all sorption sites. In either case, therefore, eq 11 is reduced to (13)
Continuous pdf. The discrete multisite model creates two parameters and T,i) for each individual type of sorption sites. When there are more than two types of sorption sites, eq 15 becomes less practical because of too many estimatedparameters. A pragmaticapproachwould be a continuous pdf of T, since continuous distributions usually need only two parameters (e.g., the mean and the variance). A candidate for such a continuous pdf for T, is the generalized y distribution:
where T(n) is the y function:
r(n)=
LEn-'exp(-E) d t
(17)
,9 and n are scale and shape parameters determining the mean and variance of the distribution, respectively, as
-
n
T, = -and
P
n 4 =-
P2
(18)
The y function is well-known for its flexibility in describing a wide range of skewed as well as very symmetric distributions. It has been used by Connaughton et al. (25) and Pedit and Miller (I) to represent sorption rate variability. It should be noted, however, that other pdfs are also available,but differentiation between different continuous pdfs was not the primary interest of this work. Given the y pdf of T, as eq 16, the mean transfer function H(s) is then derived by substituting eq 16 into eq 13 and integrating (16):
where ps is the pdf of Ts (or marginal pdf of Ts when b = 0) *
Discrete pdf. SoUaquifer materials are ofen assumed to consist of m different sized and shaped sorbent particle classes or of a number of sorbing components. In the lumped first-order approximation, we envision this as m distinct classes of sorption sites. The pdf of Ts (p,) is then discrete and is written as 1
m
(14)
wherefi is the probability of a solute molecule being sorbed onto site i with characteristic time T,i; and 6+ (5) is a Dirac function of 5, as defined before. Substituting eq 14 into eq 13 gives a simple algebraic sum for H(s):
(19)
where T(&zlis the incomplete y functionofzwith parameter 5 (26'). Equation 19 will be referred to as a y sorption (GS) model. Equations 15 (case m = 2) and 19 reveal that both TS and GS models include four parameters: two common constants, a and b two specific, fi and Ts2 for TS, 9, and n for GS. The four parameters, however, may be reduced to three after subjectingthe constraint of the sample-averaged KO,which is apotentially measurableparameter as pointed out by Pedit and Miller (1). That is, from eq 12, the expectation of KD is
E[KD]= RD = &[ (15)
Equation 15 is a multisite sorption model with two sorption parameters for each site, the fraction of the site 61,and the characteristic time (Tsj). The KD value at each site is correlated through eq 12. The well-known two-site (TS) sorption model, therfore, corresponds to m = 2 and b = 0 in eq 15, in which one site is T,I = 0 (equilibrium) and another site is Ts2 =- 0 (time dependent). It should be noted that in this case, we operationally represent Tslb= 2728
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 29, NO. 11.1995
r,b]
(20)
where E[TSb]is the bth moment (potentiallyfractional).For TS E[7',bI = (1 - fi)T,2b
(21)
and for GS (22)
Substituting eq 21 or eq 22 into eq 20, therefore, results in
four parameters that are all constrainedto the h o w n value of E[KDI. This means that whenever three of them are chosen, the fourth is automatically determined by eq 20.
EActs of Sorption-Si Heterogeneity on Transport
/ /
To quantitatively evaluate the effect of sorption heterogeneities in both equilibrium and rate properties on transport and to understand the implications, the time moment methodisutikedas describedindetailbyvalocchi (27).For example, time moments serve as descriptors of a concentration btc. The first four moments indicate the mean breakthrough time, the degree of spreading, the asymmetry, and the peakedness of the curve, respectively. It is expected therefore, that each higher time moment would reveal more information of the influence of ratelimiting processes that are embedded in the btc. The ith absolute moment of a btc generated by eq 8 is calculated as (27)
Given the absolute moments calculated as eq 23, the corresponding central moments can be obtainedaccording to Abramowitz and Stegun (26). The iirst four central moments of eq 8 with substitution of eq 13 are calculated as (16)
mbtc= 1
p4 = 3de:(1+
+ ~U(~:RE[T~”~I + E[T,”‘I)
(26)
56:)
+ 12a(R%:(1+
::
0.8
(25)
a(i + E:)(E[T,”~IY where mbtc and
(24)
+ 2aE[T,”’]
.‘b,= R%: p3 = 3 d 6 :
+ aE[T,q = R
FIGURE 2 y pdf m o m i n H t l , ofthe sorption characteristic time 1. a t a function of jand a WW the fl value fixed to 1.0. platted by Mathematica. version 22 Wolfram Research, Inc., Champaign, I1 61820.
-
.:i j.,: .:. :,....
- a=1.53. b = O a=2.36, b = O . 1
+
~E:)E[T~”~I
+ ~ E : R E [ T +~ ZE[T~”~I)(27)
are the btc mean and variance;p3 and
p4 are the btctbird and fourth central moments, respectively;
R is a modified retardation factor, which reduces to the regular value R as b = 0 in eq 24; and E[TA (i = b i, i = 0,1,2,...) is thejth absolute moment of Ts,determined by eq 21 for TS or by eq 22 for GS. Two interesting results are now demonstrated in the expressions of the first four moments for a given dispersive system he., eris given). First, when an LFER holds, i.e., KO and T, are perfectly correlated (b > 0). all btc moments becomedependentuponthe sorptiontimestatistics (EITA). Second, when an LFER is not valid, i.e., KO and T, are completelyindependent (b= 01, the first btcmoment (mean breakthrough time, m d becomes solely determined by the equilibrium properties (= retardation factor R). but higher moments describing btc spreading and tailing are all influenced by both equilibrium and rate sorption statistics. A general assessmentof the impact of the sorption time statistics on solute breakthrough can be difficult because of the complex nature of EITJ. E[TAcan be a function of the LFER constant b a s well as the two sorption time pdf parameters, for example, n and 9, in the y case (eq 22). Figure 2 depicts the change of EITJ as both j and n are independent variables while the scale factorb is fixed to 1. In general, EITJ is a concave “U-shaped function ofj at
+
VOL. 29. NO. 11.1995 I ENVIRONMENTAL SCIENCE &TECHNOLOGY. 272s
time scales than the noncorrelated case (b= 0). Increasingly sharpened and advanced peaks were also exhibited as the LFER constants increased. It is intuitively apparent that an LFER correlation would enhance the breakthrough behavior at both ends of the time scale because the fast sorption fractions controlling the peak (T, small) are assigned small KD (less retardation), while the opposite is true for the slower rate fraction that controls the tail. More accurately,an LFER correlation would completely exclude any real equilibrium process (T, = 0) since KD is zero at these sites (eq 12). An equivalent representation of a complex mechanistic model using simple equations such as first-order is commonly observed in the scientificliterature. The approach is usually accomplished by equating the first two btc moments between models (27) or by other means (21).As to the current effort, the question becomes whether sorption-site heterogeneity described by a generalpdf such as GS would be equivalently represented by a simple TS distribution when the first few btc moments are set equal between them. As we can see from the btc moments (eqs 24-27), it is possible in principle that the first two moments (the mean and variance) could be equal between TS and GS if E[Tsb]and E[TSb+’]can be equated between them by adjusting the two pdf parameters in either GS or TS. Such a manipulation, however, is not warranted because of the complex nature of the E[TJI functions (eqs 21 and 22). However this can be done when KDand Tsare not correlated ( b = 0). In this case, the first three btc moments become identicalbetween TS and GS so long as the first two sorption time moments (E[T,]and E[Ts21)are equal between them. The analysis for the uncorrelated case, therefore, can serve as a limiting case when the LFER does exist. Equating the first two sorption time moments (E[T,I and E[T?]) between TS and GS gives a simple relationship for the two pairs of parameters in each pdf (16):
fi = 1/(1 + n)
(28)
Ts2= (1 + n , / g
(29)
and
Given eqs 28 and 29, the first three btc moments obtained from the TS and GS models become identical. The discrepancy of the two models is then reflected in the difference of the fourth btc moment:
A4 = 24(R - 1)~:(1
+ E:)T:
(30)
where E , is the CV of Ts. Equation 30 demonstrates that the deviation of the fourth moment between TS and GS increases rapidly with the increase of variability of the sorption characteristic time T, (fourth power on E,) and the mean value ?, (third power). The high sensitivity is expected to hold for even higher moments. This means that deviations of the simple TS description of sorption-site heterogeneity from the more general GS approach would occur primarily at long-time transport behavior. Only when the variability of the sorption rate is small or when sorption is significantly faster than transport, as reflected by ?; (=?,/to, and ts = 1 / is~the nonreduced sorption characteristic time), would the two model predictions be similar at a particular time scale. Equation 30 also indicates that the differencebetween the fourth moments increases linearlywith the increase of the retardation coefficient R. 2730
ENVIRONMENTAL SCIENCE &TECHNOLOGY / VOL. 29, NO. 11, 1995
TABLE 1
Parameter Values for Simulation of Transport under Sorption Nonequilibrium GS model
cases
7,
la Ib
0.1 1.0 10.0 20.0 2.0 2.0 2.0 2.0 5.0 5.0 5.0 5.0
IC Id 2a 2b 2c
2d 3a 3b 3c 3d
cS
1
2.0 3.4 2.0 3.4 2.0 3.4 2.0 3.4 0.3 3.4 1.5 3.4 3.0 3.4 5.0 3.4 2.0 1.24 2.0 3.4 2.0 5.8 2.0 13.0
TS model
/3
n
fi
Ts~
2.5 0.25 0.025 0.0125 5.555 0.22 0.055 0.02 0.05 0.05 0.05 0.05
0.25 0.25 0.25 0.25 11.11 0.44 0.11 0.04 0.25 0.25 0.25 0.25
0.8 0.8 0.8 0.8 0.083 0.69 0.90 0.962 0.8 0.8 0.8 0.8
0.5 5.0 50.0 100.0 2.18 6.545 20.182 52.0 25.0 25.0 25.0 25.0
cr
0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
To further illustrate the effects of E,, T,, and R in eq 30 on solute transport and the discrepancy between the TS and GS model predictions, systematic simulation exercises were performed as three cases based on the parameter values given in Table 1. Simulation results are shown in Figures 4-6, where all concentrations in the figures are presented in the log scale to allow scrutiny of the btc tails. To facilitate the subsequent analysis of effects of R on the two model simulations,simulated btcs were plotted against the time variable (TIR)normalized by the retardation factor R. As expected from the proceeding discussions,increasing all three coefficientsincreases the difference between the two model predictions. Generally, predictions by the GS model are more spread but less peaked than the TS predictions in each case. Since in each case both models have the same mean and variance of the sorption characteristic time (T,) and all other transport-related parameters, the difference between model predictions is solely attributed to the differentassumptions regarding sorptionsite distribution: the TS model only considers two distinct regions, instantaneous and time-dependent; while the GS model consists of a continuum of sorption sites. The linear decline of all TS-simulated btc tails on the log scale is characteristicof the assumption of a first-orderrate process with a single rate constant assumed for the time-dependent site(s). Conversely, the extended tailing generated by the GS model signifies the effect of assuming a range of sorption sites, with a consequent time-dependent process leading to prolongation of leaching. This observation will be further elucidated with experimental data in the subsequent work on long-term leaching analysis currently in progress. Comparison between model-simulated btcs from the two pdfs of sorption-site heterogeneity under different ?, values is illustrated in Figure 4. The sharpness of btc peaks from the TS model developed rapidly as the Ts value increased from 0.1 to 20.0 (Figure 4a-d), indicating high dependence of the TS model upon the relative characteristic time of sorption and transport. When the value is high (fasterporewater velocity),sorption is generally slower than transport and nonequilibrium conditions present. For the TS model, this infers that most solutes bypass the singlevalued, time-dependent region(s),resultingin a sharp peak and a very small amount of sorbed solutes to produce a flattened btc tail. In the low value case (slower porewater velocity),nonequilibrium conditions are less severe, and most solutes have sufficienttime to diffuse to or react with
1oo lo-'
1o
- ~
1o
- ~
1 o-8 0
1
2
3
4
5
6
0
1
2
3
4
5
6
1
2
3
4
5
6
0
1
2
3
4
5
6
1 oo
lo-'
5;' F;
1o
- ~
L
c,
0
T/R
T/R
FIGURE 4. Effect of % on breakthrough curves of two-site (dashed lines) and y (solid lines) sorption models (case 1 in Table 1). Values of T, are indicated on each curve.
a . E, = 0.3
I
.I
91 0
2
4
6
8
10
0
2
4
6
e
io
10
0
2
4
6
8
10
1 oo lo-'
1o
- ~
I I
0
I
I
I
I
2
4
6
8
T/R
T/R
FIGURE 5. Effect of cS on breakthrough curves of two-site (dashed lines) and y (solid lines) sorption models (case 2 in Table 1). Values of cS are indicated on each plot.
the presumed single-valued, time-dependent region(s), giving a wider peak spreadingand a closer match to the GS
I
model predictions (Figure4a). Compared to the TS model, the GS model has been less affectedby ps,furtherindicating VOL. 29, NO. 11. 1995 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
2731
1
a. R = 1 . 2 4
b. R = 3.40
- GS
-----____ I I 0
loo
I
I
I
I
2
4
6
8
10
5.80
I
8
10
1
0
c. R =
2
4
6
T/R
0
2
4
6
8
10
d . R = 13.00
0
2
6
4
8
10
T/R
FIGURE 8. Effect of the retardation factor R on breakthrough curves of two-site (dashed line) and y (solid lines) sorption models (case 3 in Table 1). Values of R are indicated on each plot.
the complimentary effect of assuming a range of sorption sites. The high dependence of deviationsbetween different sorption-site distribution models on Fs was also found by Sardin et al. (28). They compared different distribution models having the same mean T, and attributed the model deviation to the variance of T,. As seen from the current simulations,where the same mean and variance were both used for T,in all comparison cases, substantial deviation between btcs still exists, indicating the intrinsic difference between TS and GS pdfs in describing sorption-site heterogeneity. The high dependence of the TS model on ?', may also imply that this model can be highly velocitydependent. The effects of different CVs of T, are shown in Figure 5. Theoretically, when the value of the CV tends to zero, both model predictions will converge to the results of a one-site model, representing homogeneous sorption behavior. This trend is shown in Figure 5, where a closer match between the TS and GS model simulations was achieved when the CV decreased from the highest 5.0 (Figure 5d) to the lowest 0.3 (Figure 5a). An interesting result shown in Figure 5 is that the deviationbetween these two models appears primarily on the tailing as the variability of the sorption rate increases. Although an apparently higher variability of sorption rate would spread the solute peak more significantly, peak spreading is limited in the TS model since the increase of the sorption rate variability means that the rate deviation increases between the two distinct sorption sites. That is, the single-valued timedependent site(s1has to be very slow, and most mass will bypass the slow site(s1,resulting in a sharpened instead of a more distributed peak (Figure 5a-d). 2732 a ENVIRONMENTAL SCIENCE & TECHNOLOGY /
VOL. 29, NO. 11. 1995
Figure 6 compares the influence of differentretardation values (R)on the two model simulation results. The time variable (TIR)was normalized against the first btc moment (eq 241, which allows the same relative time scale for the comparison among different R cases (Le., the mean breakthrough time is always unit). For the TS model (dashed lines), an increase in R increases the peak retardation as well as the concentration level in the btc tail. Comparatively, the GS model simulations are less influenced. Because the TS model assumes a single-valuedrate (desorption) constant ( K ) for the time-dependent sites, because no LFER between KO and Kwas assumed, and because KD is the ratio of the forward (sorption) to the reverse (desorption) rate constants, the increase of R (Le,, KDincreases)would only increasethe forward mass transfer rate, thus resulting in enhanced retardation and tailing for this model. By comparison, the GS model was hypothesized to contain a distribution of K. Thus, a change in R would not significantly change the model predictions. This indicates that the GS model, with no correlation between KD and T,,may well approximate the cases where KD and T,are not strictly correlated (as in eq 12). In summary, the deviations between TS and GS model simulations in all above cases tend to become more severe as time increases. An LFER between KDand T, would only increase the deviation between the models, as the ultimate model equivalence can be achieved onlyup to their second btc moments. This indicates that when sorption heterogeneity exists, simple assumptions such as TS are questionable to sufficiently represent a more generally distributed sorption process (e.g., GS). This is particularly true at both early and prolonged times. The attempt to lump
the effect of nonequilibrium processes into an “effective dispersion coefficient” or to fmd a simple equivalent representation for a more complex mechanisticmodel (27) may be effective for homogeneous systems having uniform equilibrium/rate sorptionproperties and for processes over a relatively short time scale. However, it becomes risky, as shown in the current work, when sorption-site heterogeneity and long-timebehavior are of concern. Comparisons between TS and GS model capability to extrapolate experimentalresults obtained over a variety of time scales is the next step in this framework.
Acknowledgments The authors would like to thank the two anonymous reviewers for their insight and valuable comments on this work. Thanks are also expressed to Drs. L. W. Lion and J. R. Stedinger of Cornell University for their reviews on the original version. This research was funded by the U.S. Geological Survey (Contract 14/08/0001/G1906).
Glossary empirical constant in the LFER between KD and T, empirical constant in the LFER between KD and T, local solute concentration, M/L3 = c/ Co, dimensionless local concentration input concentration, MIL3 ensemble mean concentration, dimensionless expectation operator fraction at sorption site i, dimensionless mean transfer function of an ensemble of ratelimiting process transfer function of a rate-limiting process = KtO, dimensionless first-order rate coefficient = & p l e , dimensionless sorption partitioning coefficient sorption partitioning coefficient, L3/M spatial reference length (column length), L ith absolute moment the first moment of a breakthrough curve shape factor in the y pdf = 2/62, Peclet number joint pdf of variates, KO and Ts marginal pdf of KD marginal pdf of Ts = 1+E[%],mean retardation factor, dimensionless = 1 uE[ modifled retardation factor, dimensionless = Sbp/ (Oca),dimensionlesssorbed concentration Laplace parameter sorbed concentration, MIM = 1/K, dimensionlesssorption characteristic time = 1/K, dimensionlesssorption characteristic time at sorption site i time, T mean residence time, T = 1 / ~sorption , characteristic time, T local porewater velocity, L/T = vto/Lo, dimensionless local porewater velocity spatial coordinate, L
+ GI,
= x/L,,, dimensionless spatial coordinate
dummy variable forward first-order mass transfer rate coefficient, T-1
scale factor in the y pdf the fourth btc moment difference between the TS and GS models y function the coefficient ofvariance of sorption characteristic time T, the coefficient of variance of residence time t reverse first-order mass transfer rate coefficient, T-1
the third central btc moment the fourth central btc moment bulk density, M/L3 btc variance variance of sorption characteristic time Ts dimensionless residence time porosity, L3/L3 dummy variable Laplace transform operator
Literature Cited (1) Pedit, J. A.; Miller, C. T. Enuiron. Sci. Technol. 1994, 28, 20942194. (2) Rao, P. S. C.; Jessup, R. E.; Addiscott, T. M. Soil Sci. 1982, 133, 342-394. (3) Ruthven, D. M.; Loughlin, K. F. Chem. Eng. Sci. 1971,26, 577584. (4) Rasmuson, A. Chem. Eng. Sci. 1985, 40, 621-629. (5) Jury, W. A. In Chemical mobility and reactivity in soil systems; Nelson, D. W., Elrick, D. E., Tanji, K. K., Eds.; Soil Science Society of America: Madison, WI, 1983; pp 49-64. (6) van der Zee, S. E. A. T. M.; van Riemsdijk, W. H. Water Resour. Res. 1987, 23, 2059-2069. (7) Cvetkovic, V. D.; Shapiro, A. M. Water Resour. Res. 1990, 26, 2057-2067. (8) Gamerdinger, A. P.; Wagenet, R. J.; van Genuchten, M. Th. Soil Sci. SOC. Am. J. 1990, 54, 957-963. (9) Brusseau, M. L.; Jessup, R. E.; Rao, P. S. C. Water Resour. Res. 1989, 25, 1971-1988. (10) Pignatello, J. J.; Ferrandino, F. J.; Huang, L. Q. Enuiron. Sci. Technol. 1993, 27, 1563-1571. (11) Cameron, D. A.; Klute, A. WaterResour. Res. 1977,13, 183-188. (12) Wu, S., Gschwend, P. M. Environ. Sci. Technol. 1986, 20, 717725. (13) Ball, W. P., Roberts, P. V. Enuiron. Sci. Technol. 1991,25, 12371249. (14) Brusseau, M. L.; Jessup,R. E.; Rao, P. S. C. Enuiron. Sci. Technol. 1991, 25, 134-142. (15) Simmons, C. S. Water Resour. Res. 1982, 18, 1191-1214. (16) Chen, W. Ph.D. Dissertation, Cornell University, 1994. (17) Villermaux, J. 1. Chromatogr. 1987, 406, 11-26. (18) Beran, M. J. Statistical Continuum Theories;Wiley Interscience: New York, 1968. (19) Johnson, N. L.; Kotz, S. Distribution In Statistics: Continuous Univariate distributions-1;John Wiley & Sons: New York, 1970; pp 137-153. (20) Jury, W. A.; Roth, K. Transfer functions and Solute Movement through Soil: TheoryandApplication;Birkhauser Verlag: Basel, 1990.
(21) van Genuchten, M. Th.; Dalton, F. N. Geoderma 1986,38, 165183. (22) Weber, W. J.; McGinley,P. M., Jr.; Katz, L. E. Enuiron. Sci. Technol. 1992, 26, 1955-1962. (23) Brusseau, M. L.; Rao, P. S. C. Geochemosphere 1989, 18, 16911706. (24) Meyer, P. L. IntroductoryProbabilityandStatisticalApplications; Addison-Wesley Publishing Company: Reading, 1970; pp 93114.
VOL. 29, NO. 11, 1995 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
2733
(25) Connauahton, D. F.; Stedinaer, J. R.; Lion, L. W.: Shuler, M. L.
Enuiron:ScL Technol. 1993r27,2397-2403. (26) Abramowitz, M.; Stenun, I. A. Handbook of Mathematical Functions with Form&, Graphs, and Mathematical Tables; Dover Publications, Inc.: New York, 1972; p 260. (27) Valocchi, A. J. Water Resour. Res. 1986,22, 399-407. (28) Sardin, M.; Schweich, D.; Leij, F. J.; van Genuchten, M. Th. Water Resour. Res. 1991,27, 2287-2307.
2734 1 ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 29, NO. 11,1995
I Received for review Sevtember 29, 1994. Revised manuscrivt received April 3, 1995.Accepted April 6, 1995.@
I
ES940611p
@
Abstract published in Aduance ACS Abstracts, October 15, 1995.
