Envlron. Scl. Technol. 1994, 28, 2094-2104
Heterogeneous Sorption Processes in Subsurface Systems. 1 Model Formulations and Applications Joseph A. Pedit and Cass T. Miller'
Department of Environmental Sciences and Engineering, CB 7400, 104 Rosenau Hall, University of North Carolina, Chapel Hill, North Carolina 27599-7400 ~
The development of an appropriate model for sorption processes in natural systems has proven elusive. We hypothesize that standard deterministic models cannot capture the inherent variability in most natural systems. This work responds to deficiencies in current models by formulating two classes of sorption models: an extended set of deterministic models based upon a generalized multiple domain concept and a class of stochastic models. Stochastic sorption models are formulated to describe microscale variability, which may be interparticle or intraparticle in origin, in terms of both sorption equilibrium and rate properties. Stochastic model formulations treat sorption equilibrium and rate properties as continuously distributed random variables described by lognormal or r probability density functions. Long-term sorption equilibrium and rate data are shown for a herbicide on a subsurface material. Results from parameter estimation for 16 different rate models applied to the data set are summarized. Model-data comparisons show a range of agreement between individual model fits and experimental data, with several trends evident: diffusion models outperformed first-order models; the addition of an instantaneously sorbing fraction significantly improved model-data agreement; multiple-particle class models more closely fit the data than single-particle class models, without an increase in the number of fitted parameters; and stochastic first-ordermodels agreed well with the data.
Introduction
Contamination of sediments and subsurface environments by hydrophobic organic compounds has led to significant research in sorption processes. Many mathematical models have been formulated for modeling the sorption process; most use a deterministic framework and assume an idealized homogeneous sorbent structure. However,natural systems display heterogeneities at scales ranging from the megascale used to describe large natural systems to the microscale of an individual solid grain. Experimental investigations of sorption processes are often conducted in batch systems in an effort to isolate the sorption process from advective and dispersive transport processes. Modeling of such systems requires a description of sorption equilibrium and rate properties. Most sorption models assume the solid phase of a bulk sample to be homogeneous with respect to sorption equilibrium and rate properties, although many experimental investigationshave demonstrated otherwise. The first objective of this work was to extend currently used deterministic modeling frameworks to account for heterogeneities that have been observed within bulk samples; the second objective was to formulate a class of models that treat sorption equilibrium and rate properties as
* Corresponding author; e-mail address:
[email protected]. 2094
Environ. Sci. Technol., Vol. 28, No. 12, 1994
continuously distributed random variables; and the third objective was to examine the ability of the various models to describe the results of a sorption rate experiment and to discuss the appropriateness of these models for systems with heterogeneous sorption properties. Background
The extent to which a hydrophobic organic compound will sorb to anatural sorbent in aqueous-saturated systems is primarily dependent upon the organic carbon content of the sorbent. Adsorption to mineral surfaces and partitioning into structured fluid near mineral surfaces can be significant contributors to overall sorption when the organic carbon content is low ( I , 2). Variations in organiccarbon content and mineralogy have been observed at the field scale (3-7), core-sample scale (81, and bulksample scale (9-14). Bulk samples are defined as samples that are a composite, typically from auger cuttings or excavations. Given these variations, it follows that variations in the equilibrium distribution of a solute between the fluid and solid phases should be expected at many scales. Field-scale variations in sorption equilibrium were observed by Elabd and co-workers (31, who measured the distribution coefficient of napropamide on 36 soil samples from a 0.6-ha field site and found a coefficient of variation of 31 % Field-scale variations in sorption equilibriumwere also observed by MacIntyre and co-workers (6),who measured the distribution coefficient of naphthalene on 10 different composited core samples from an aquifer at a Columbus,MS, site and found the distribution coefficient to vary by a factor of 9. Variations in sorption equilibrium at the core-sample scale were demonstrated by Mackay and co-workers (8), who measured the distribution coefficient of tetrachloroethene on 10-cmsegments of a 2-m core sample from an aquifer at the Borden, Ontario, site and found the distribution coefficient to vary by nearly an order of magnitude. Allen-King and Mackay (15) found similar variations in tetrachloroethene sorption behavior over a 1-m core sample from the Borden aquifer. The variation of sorption equilibrium as a function of particle size within a bulk sample has been demonstrated in numerous studies (9-11,14, 16-18). Ball and Roberts (14) also found variations in the distribution coefficient within a particle size fraction that had been separated by magnetic properties. The sorption process is assumed to be diffusion limited. The time required for a solute-sorbent system to reach equilibrium has been observed to be related to the magnitude of the equilibrium distribution relationship (19, 20) and, to some extent, particle size (14, 21); both observations are consistent with a diffusion-limited process. Given the variations in equilibrium distribution relationships and particle size that are possible within a
.
0013-936X194/0928-2094$04.50/0
0 1994 American Chemical Society
bulk sample, it is not surprising that sorption modeling approaches that incorporate equilibrium distribution relationships that are measured on a bulk sample and incorporate a first-order mass transfer rate model have met with only limited success in batch systems (e.g., refs 22 and 23). An instantaneous equilibrium fraction is often incorporated into rate modeling of batch systems to improve model fits (e.g., refs 21, 22, and 24). The need to use an instantaneous equilibrium fraction in rate modeling of batch systems is likely to be a result of model misspecification caused by variations in sorption equilibrium and rate properties at the grain scale.
Model Formulations In the first part of this section, sorption rate models are formulated based upon multiple-particle class extensions of Fickian diffusion models and first-order mass transfer models. In the latter part of this section, aclass of sorption rate models that treat equilibrium and rate properties as continuously distributed random variables are formulated. Regardless of whether a deterministic or stochastic formulation is used, sorption rate models require a description of the equilibrium distribution relationship of a solute between the fluid and solid phases. Equilibrium Modeling. The equilibrium distribution of a solute between the fluid and solid phases is often assumed to be a linear relationship, particularly at low fluid-phase solute concentrations. Many systems display nonlinear equilibrium distribution relationships, particularly if fluid-phase solute concentrations span several orders of magnitude. In this work, the Toth model is used to describe nonlinear sorption because it reduces to the linear model as the fluid-phase solute concentration approaches zero. The Toth model is given by
where qe is the equilibrium solid-phase solute concentration, QOis a sorption maximum, b is an affinity parameter, Ceis the equilibrium fluid-phase solute concentration, and PT is a heterogeneity parameter (25). The Toth model reduces to the Langmuir model when PT is equal to 1 (homogeneity). The Toth model has great flexibility in its ability to fit experimental data, but it has three adjustable parameters. Fickian Diffusion Modeling. The most common sorption rate model formulations assume that the sorption process is limited by diffusion. The diffusion process can be modeled by Fick's law for a well-defined geometry or by a first-order mass transfer approximation. An extensive review of sorption rate modeling is not warranted here; reviews are available (26, 27). Diffusion models based on Fick's law often assume that the sorption process is controlled by pore or surface diffusion into a well-defined particle geometry, usually a sphere. Pore diffusion is envisioned to occur in the intraparticle pore spaces of minerals and aggregates. Severalinvestigators have used the pore diffusion approach tomodel batch systems (21,23,28-31). Wuand Gschwend (23,32) accounted for the effect of particle size variations within a bulk sample by allowing for multiple particle classes with different particle radii, using a pore diffusion approach. The pore diffusion modeling approach has been
extended to include an instantaneous equilibrium fraction in an effort to improve model fits (21, 31). Surface diffusion is envisioned to occur along intraparticle pore walls, the surface diffusion coefficient (DJ being dependent upon the solid-phase solute concentration (33). The surface diffusion coefficient is often assumed to be constant for the sake of simplicity or because of the lack of information about its concentration dependence (34). Several investigators (22, 24, 35) have used the surface diffusion model assuming a constant surface diffusion coefficient to model batch systems. A general model based on Fick's law should incorporate the major features that have been used by other researchers. The general model formulated for this work allows for both pore and surface diffusion in a fashion similar to the approach of Crittenden and co-workers (36). The general model allows for multiple particle classes with different physical properties and sorption parameters for each class, using an approach similar to previous works (e.g.,refs 23,32, and 37-39). The generalmodel also allows for an instantaneous equilibrium fraction for each particle class, following previous work (21, 31). The governing system of equations and numerical methods used to solve the governing equations are discussed elsewhere (40). First-Order Mass Transfer Modeling. First-order mass transfer model formulations for the sorption process are frequently used for modeling batch systems (22, 23, 29,41-46). Wu and Gschwend (32)accounted for the effect of particle size variations within a bulk sample by allowing for multiple particle classes with the first-order mass transfer coefficient for each particle class being inversely related to the square of particle radius. A fraction of equilibrium-type sorption sites are often included in first-order mass transfer models that have been used to model batch systems (e.g., refs 19 and 22). Firstorder mass transfer models that include a fraction of equilibrium-type sorption sites are often referred to as bicontinuum models. Physical process interpretations of bicontinuum models, often referred to as two-region or mobile-immobile fluid-phase models, assume that all of the sorption sites are at equilibrium with the adjacent fluid-phase solute concentration. However, transfer of solute between the mobile and immobile fluid phases is diffusion controlled. In the model formulation that follows, the mobile and immobile fluid phases are viewed as the bulk and the intraparticle fluid phases, respectively. A general first-order mass transfer model, which allows for multiple particle (or sorption site) classes with equilibrium and diffusion-controlled sorption sites, is given by
M , n~ V
i=l
I -
dqf(t) dCb(t)
l+--Cf'f'-
'ddCb(t)
--
dt
where M , is the mass of solids, V is the volume of bulk fluid, i is the particle-class index, np is the number of particle classes, f is the mass fraction of particle class i, f is the fraction of equilibrium-type sorption sites of particle class i, qf(t) is the solid-phase solute concentration of particle class i at equilibrium with the bulk fluidphase solute concentration, Cb(t) is the bulk fluid-phase
i
Environ. Sci. Technol., Vol. 28, No. 12, 1994 2095
solute concentration, t is time, k;o is the first-order mass transfer coefficient for class i, S:(t) is the sorbent-phase solute concentration for the class i diffusion-controlled fraction at equilibrium with the bulk fluid-phase solute concentration, Si(t)is the sorbent-phase solute concentration for the class i diffusion-controlled fraction, hi is the equilibrium-type solid-phase first-order reaction rate coefficient for particle class i, and hb is the bulk fluidphase first-order reaction rate coefficient. Sk(t) and Si(t) have units of mass of solute in the class i sorbent phase per mass of solid phase in the class i sorbent phase. Note that the diffusion-controlled sorbent-phase solute concentration includes solute in the intraparticle pore fluid as well as solute sorbed to the solid phase. If the linearlocal-equilibrium assumption can be invoked within the sorbent phase, then S l ( t ) and Si(t) for each class are given by
where 06 is the intraparticle porosity of particle class i, the apparent particle density of particle class i, K i is the distribution coefficient of particle class i, and Cc(t) is the intraparticle fluid-phase solute concentration ofparticle class i. The mass balance equation for each class when the linear model is used to describe equilibrium is p i , is
where A; is the intraparticle fluid-phase first-order reaction rate coefficient for particle class i, and h; is the intraparticle solid-phase first-order reaction rate coefficient for particle class i. If solute is initially only present in the bulk fluid phase, then the appropriate initial conditions are C,(t = 0) = CbO
(6)
~ ; ( =t 0) = o for i = 1, ..., np
(7)
where CbO is the initial bulk fluid-phase solute concentration. Stochastic Modeling. The two general model formulations outlined above attempt to address the variability that is possible in a bulk sample. The multiple-particle class formulation of the Fickian diffusion model can only address variability down to the individual particle scale. The formulation assumes an individual particle to be homogeneous and does not allow for intraparticle variations in sorption rate and equilibrium. The multipleparticle class formulation of the first-order mass transfer model allows for intraparticle variations in sorption rate and equilibrium, but experimentally such measurements are not possible. Each particle class of a bulk sample could be assigned physical properties and sorption parameters with these model formulations, but such an approach is impractical from an experimental and modeling perspective. A more viable approach is to describe the sorption 2096
Envlron. Sci. Technol., Vol. 28, No. 12, 1994
equilibrium and rate parameters of the bulk sample with a continuous multivariate distribution, which is applied to a first-order mass transfer framework below. A batch system at equilibrium does not have a uniform sorbent-phase solute concentration throughout the bulk sample. Instead, the sorbent-phase solute concentration varies due to local variations of the equilibrium distribution relationship. These variations can be viewed as occurring a t the interparticle scale and the intraparticle scale. Variations will also exist in the rate of solute transfer to and from the bulk fluid phase, even for particles or collections of sorption sites with the same equilibrium distribution relationship. At the interparticle scale, these variations may occur due to differences in particle size and shape. At the intraparticle scale, these variations may occur due to differences in depth or accessibility of collections of sorption sites within the particle. In the model formulation that follows, it is assumed that the linear equilibrium model and first-order mass transfer model apply throughout the bulk sample. The sorbentphase solute concentration, distribution coefficient, and first-order mass transfer coefficient are treated as random variables. Random variables are designated by bold italics with particular values of the random variables designated by ordinary italics. The model formulation ignores the presence of equilibrium-type sorption sites because these sites are represented by large values of the first-order mass transfer coefficient random variable. The model formulation also ignores direct mass transfer between adjacent collections of sorption sites within a particle (i.e., only mass transfer in parallel is allowed). With the assumptions outlined above, the bulk fluid-phase solute concentration mass balance may be written as
where fm(Kd,kfo) is the mass-fraction probability density function, S,(t) is the sorbent-phase solute concentration for a particle or collection of sorption sites at equilibrium with the bulk fluid-phase solute concentration, and S(t) is the sorbent-phase solute concentration for a particle or collection of sorption sites with locally defined sorption parameter values Kd and kf,. The form of the probability density function is unknown but is restricted to functions that equal zero for negative values of the continuously distributed random variables Kd and kf,. A variety of probability density functions satisfy this restriction, the most convenient being associated with the bivariate lognormal distribution. The probability density function for the bivariate lognormal distribution is (47)
- &I2
(In(Kd) - kK)(ln(kfo)- kk) + 'KOk
2 Ok
where p is the correlation coefficient; p~ and Pk are the means of ln(Kd) and ln(kfo),respectively; and u i and u: are the variances of ln(Kd) and ln(kf,), respectively.
There are five parameters ( p , p ~ pk, , &, and needed to describe the bivariate lognormal distribution model. However, the mean and variance of In(&) are restricted because the expected value of Kd must equal the experimentally determined bulk distribution coefficient, effectively reducing the number of parameters to four. This restriction is given by (47)
E(Kd) = exp(pK + 4 / 2 ) Also note that the arithmetic variance of
(11) Kd
v(Kd) = [E(Kd)12[exp(&)- 11
is given by (12)
Stochastic Model Simplifications. The bivariate lognormal distribution model can be simplified when the correlation coefficient is equal to f1. The model is reduced to a univariate lognormal distribution with the massfraction probability density function expressed in terms of one variable given by
The mass-fraction probability density function for ln(kf,) is of the same general form as eq 13. A linear relationship between In&,) and In(&) exists in this case and may be written ln(kfo)= cyo
+ a1ln(Kd)
(14)
where ‘YO is the intercept and 011 is the slope. From eq 14, it follows that the mean and variance of ln(kf,) are pk = ‘YO + CYIPK and ut = ala;, respectively. The expected value and arithmetic variance of kf, are of the same form as eqs 11 and 12. There are four parameters ( p ~u, i , ‘YO, and ‘ ~ 1 ) needed to describe the simplified bivariate lognormal distribution model. However, the mean and variance of In(&) are again restricted by eq 11,effectively reducing the number of parameters to three. The linear free energy relationship (or LFER) developed by Brusseau and Rao (20) is of the same general form as eq 14, the LFER being developedbased on bulk behavior. The model given by eqs 8, 11, 13, and 14 will be referred to as the lognormal linear free energy relationship (or LLFER) model. A model similar to the LLFER model can be formulated if it is assumed that Kd follows a r distribution instead of a lognormal distribution. The I’ mass-fraction probability density function is given by (48)
where ‘YK is a shape parameter, PK is a scale parameter, and r(‘YK) is the r function given by I’(aK)= K u a K - ’exp(-u) du
(16)
The I? distribution is versatile in the sense that it can take on a variety of shapes depending on the values of the shape and scale parameters. The I’ distribution reduces to the exponential distribution when the shape parameter is equal to one. If the linear relationship between ln(kf,,) and ln(Kd) given by eq 14 is assumed, then the mass-fraction
probability density function for kf, is given by
There are four parameters (CYK, PK, ‘YO, and a l ) needed to describe this r distribution model. However, the scale and shape parameters of Kd are restricted because the expected value of Kd must equal the experimentally determined bulk distribution coefficient, effectively reducing the number of parameters to three. This restriction is given by (48)
The arithmetic variance of
Kd
is given by
The model given by eqs 8 and 15-18 will be referred to as the r linear free energy relationship (or GLFER) model. The above formulations have assumed a linear equilibrium relationship where Kd is allowed to vary within the bulk sample. A simple extension to nonlinear equilibrium relationships can be developed if one assumes that all particle or sorption site classes exhibit the same nonlinearity and only vary in their sorption maximum. For the Toth model, this assumption would require that band pT be uniform and QObe allowed to vary in the same manner that Kd was allowed to vary. Such an approach for nonlinear equilibrium relationships is not without merit. Ball and Roberts (14),in estimating the Freundlich model parameters for sorption of 1,2,4,5-tetrachlorobenzene on two particle size fractions (0.42-0.85 and 0.1250.18 mm particle diameter ranges) of the Borden aquifer material found that while the affinity parameters varied by a factor of 9 between the two fractions, the Freundlich exponents were essentially identical (0.81 vs 0.82). When Allen-King and Mackay (15) estimated the Freundlich model parameters for the sorption of tetrachloroethene on 20 5-cm segments of a l-m core sample from the Borden site, they found that while the affinity parameters varied by an order of magnitude, only a small variation in Freundlich exponents was observed. However, not all studies have had this result. In investigating the sorption of 1,2,4-trichlorobenzene on shale-like materials isolated from three subsurface samples, Weber and co-workers (49) found that the Freundlich exponents for the shale-like materials and their parent samples compared favorably for one sample but varied by up to 0.1 for the other two samples. An alternate approach, albeit more difficult, for nonlinear equilibrium relationships would be to incorporate the underlying sorption energy distributions used to develop many of the nonlinear equilibrium models. The drawback to this approach is that a set of equilibrium data may often be equally well described by a number of different nonlinear models. The underlying sorption energy distribution functions for the competing nonlinear models may vary considerably in their shapes (50). The bivariate lognormal distribution model can be further simplified when Kd is constant throughout the bulk sample (i.e., v i is zero). The bivariate probability density function reduces to an univariate probability Environ. Sci. Technol., VoI. 28, No. 12, 1994
2097
density function given by
There are only two parameters (Pk and ai) needed to describe the univariate lognormal distribution model. The expected value and arithmetic variance of kfoare of the same form as eqs 11 and 12. This simplified model will be referred to as the lognormal rate (or LR) model. A model similar to the LR model can be formulated if it is assumed that kf, follows a r distribution instead of a lognormal distribution. Such an approach was developed by Connaughton and co-workers (511, who applied the I' distribution rate model to desorption rate experiments. The r mass-fraction probability density function for kfo is given by (48)
The expected value and arithmetic variance of kf, are of the same form as eqs 18 and 19. The model based on eqs 8 and 21 will be referred to as the r rate (or GR) model. The numerical methods used to solve all of the first-order mass transfer based models are discussed elsewhere (40).
Application of Models
Materials. A subsurface material (Wagner) was collected from a sand and gravel pit at a depth of about 25 m below ground surface, and it has been used in previous sorption studies (35, 49, 52). The material was air dried and sieved to remove particles larger than 2 mm in diameter. The material was then "prewashed" to remove nonsettling particles and easily dissolved organic matter by a procedure similar to that described by Gschwend and Wu (53) and again air dried. A portion of the material was pulverized in a tungsten-carbide ball mill (SPEX Industries, Inc.) to reduce particle size, because the use of pulverized material in sorption equilibrium experiments has been shown to substantially reduce the time required to reach equilibrium without significantly altering the final equilibrium (14,21). The material was prewashed because it was to be used later in desorption studies. Nonsettling particles and easily dissolved organic matter are possible sources of observed desorption equilibrium hysteresis (53). The unaltered material was characterized by particle size distribution (median particle size 0.49mm),uniformity (uniformity coefficient 3.0),and solid-phase density f2.673 f 0.006 (mean f one standard deviation) g . ~ mbased - ~ on four measurements] by the methods given by Black (54). The results of sieve analyses on the unaltered and pulverized samples are shown in Table 1. The majority of the pulverized material had a particle size less than 0.09 mm. Specific surface area (1.6f 0.3m2.g-l based on three measurements) and intraparticle porosity (0.012f 0.004 based on three measurements) were estimated by nitrogen adsorption-desorption on an Autosorb-1 (Quantachrome Corp.) for the unaltered material. Inorganic carbon content of the unaltered material was 3.78f 0.12% (based on four measurements) when determined by a method given by Ball and co-workers (12). Organic carbon content was determined by the persulfate oxidation ampule technique on an Oceanographic 2098
Environ. Scl. Technol., Vol. 28, No. 12, 1994
Table 1. Particle Size Distribution of Wagner Material
sieve range (mm)
geometric mean radius (mm)
1.40-2.00 1.00-1.40 0.71-1.00 0.50-0.71 0.355-0.50 0.25-0.355 0.18-0.25 0.125-0.18 0.09-0.125
