- Email: [email protected]
Evaluation of surface and film diffusion coefficients for carbon adsorption
Wat. Res. Vol. 23, No. 3, pp. 267-273, 1989 Printed in Great Britain. All rights reserved
0043-1354/89 $3.00+0.00 Copyright © 1989 Pergamon Press pie
EVALUATION OF SURFACE A N D FILM DIFFUSION COEFFICIENTS FOR CARBON ADSORPTION U. K. TRAEGNER and M. T. SUIDAN~) Department of Civil Engineering, University of Illinois, Urbana, IL 61801, U.S.A. (First received March 1988; accepted in revised form September 1988)
Abstract--The dynamics of adsorption for various adsorbate-adsorbent systems have been predicted successfully using the homogeneous surface diffusion model (HSDM). Kinetic parameters and equilibrium parameters have to be determined in order to solve the model. Batch adsorption experiments are usually used for this purpose. This paper presents a parameter search procedure which is capable of uniquely determining these parameters for the HSDM. The advantage of the proposed parameter search procedure is that it is not limited to cases of dominant transport mechanism, such as cases limited by either external film or surface diffusional resistances. Key words---adsorption, activated carbon, surface diffusion, liquid film diffusion, parameter search, closed batch reactors, homogeneous surface diffusion model
NOMENCLATURE
INTRODUCTION
Bi = Biot number based on surface diffusion (dimensionless) cb(t ) = bulk liquid phase adsorbate concentration (M/L 3) c * (t) = nondimensional bulk liquid concentration cs(t ) = liquid phase adsorbate concentration at solidliquid interface (M/L 3) c*(t)=dimensionless liquid phase concentration at interface co = initial liquid phase concentration (M/L 3) dp = adsorbent particle diameter (L) D s = solute distribution parameter (dimensionless) D, = surface diffusion coefficient (L2/T) f ( x ) = function of residual surface K = Freundlich isotherm capacity constant (L3/M) n kf = liquid film mass transfer coefficient (L/T) L --liquid volume in closed batch test (L3) m = total mass of carbon in closed batch test (M) n = Freundlich isotherm intensity constant (dimensionless) q(r, t)--absorbed phase adsorbate concentration (M/M) qavs= average adsorbed phase adsorbate concentration (M/M) q*vs= nondimensional average surface concentration qs(t) = adsorbed phase adsorbate concentration at solid-liquid interface (M/M) q*(t) = nondimensional adsorbed phase concentration at interface q0 ffi initial adsorbed phase adsorbate concentration (M/M) r = radial coordinate (L) ri(x ) ffi residual at i's observation R --- dimensionless radial coordinate sk ffi Newton step t = physical time (T) T--dimensionless time based on surface diffusion x = parameter vector X k = current parameter vector Xk+~ ffi updated parameter vector p/~= apparent particle density (M/L 3) = partial derivative 02ffi second partial derivative V = grad V2 ffi Laplacian.
The homogeneous surface diffusion model ( H S D M ) has been successfully used to model the dynamics of the adsorption process for various organic compounds on granular activated carbon (GAC). It can be used as a tool for the design engineer who is interested in planning the scope of pilot scale studies, interpreting pilot scale test results or estimating preliminary costs for G A C adsorber units. The mathematical formulation of the H S D M requires the numerical evaluation of differential equations which describe the adsorption process. The H S D M , however, contains physical parameters such as the liquid film mass transfer coefficient, kf, and the surface diffusion coefficient Ds, which cannot be measured directly by analytical means, but have to be known to get unique solutions to the H S D M . The usual way to determine these kinetic parameters is by performing closed batch tests and minimizing the difference between model calculations and actual batch data, while varying the kinetic parameters kf and D,. H a n d et al. 0983) developed a procedure for determining the surface diffusion coefficient by experimentally eliminating liquid film mass transfer resistance and comparing batch adsorption data to empirical solutions to the H S D M . The method presented in this paper is applicable to determining the kinetic parameters without any restriction on mass transport patterns. After presenting the differential equations which govern the dynamics of adsorption and presenting the basic concept of the unconstrained optimization technique, this technique will be applied to find the kinetic coefficients kc and Ds that best describe closed batch data collected for the adsorption of the dye R 6 G on 80317 activated carbon. Furthermore, the search will be applied to data collected on the adsorption o f p - n i t r o p h e n o l on
267
268
U.K. TRAEGNEgand M. T. SUIDAN
F400 activated carbon in closed batch experiments. In this case, the kinetic parameters and the equilibrium constants for the Freundlich adsorption isotherm will be determined, simultaneously. MODEL DEVELOPMENT AND GOVERNING EQUATIONS
The schematic diagram of an activated carbon granule, given in Fig. 1, illustrates the mechanism of mass transport into a spherical particle occurring in the HSDM (Weber and Chakravorti, 1974; Mathews and Weber, 1975; Crittenden and Weber, 1978). The porous particle is surrounded by a stagnant liquid layer which represents a mathematical expression for mass transfer resistances occurring outside the particle. Once the molecules diffuse through the stagnant liquid layer, instantaneous equilibrium is assumed to be established at the outer surface of the carbon particle between the adsorbate molecules in the liquid and adsorbed phases. The molecules are assumed to penetrate further into the particle by creeping along the abundant inner surfaces. The driving force for this movement is the local concentration gradient of the adsorbate along the inner surfaces of the particle. This movement is assumed to be diffusional only and is, therefore, modeled using Fick's first law. The governing equations of the HSDM as applied to the closed batch experiment are the subject of the next section.
surface load qavs. The average carbon load, which is only a function of time, is given by: 3 ffp/2 qavs = (dp/2)3 q(r, t)r 2 dr (2) with q(r, t) representing the concentration along the inner surface of the particle and dp representing the particle diameter. The variable r denotes the radial coordinate with an origin at the center of the particle. The partial differential equation from which q(r, t) is obtained is written in spherical coordinates as:
Oq(r,_~t) = D s ~~ + r { t 3 :t)q ( r '20q(r,~r t)t/
(3)
in which D,, the surface diffusion coefficient, represents a physical parameter. The magnitude of D~ is a measure of how fast the molecules diffuse inside the carbon particle and therefore sets a time scale for the adsorption process. Two boundary conditions and one initial condition have to be specified in order to obtain unique solutions to equation (3). Initially the particle is free of adsorbate, i.e.:
q(r, t) = 0.
(4)
The boundary condition at the center of the particle is:
aq(r = O, t) 0r
= 0
(5)
MODEL EQUATIONS
i.e. no substrate flux across the center. Finally, the continuity of flux at the solid-liquid interface has to be satisfied:
The overall mass balance for the closed batch experiment is:
= dp/2, t) ppDsaq(r Or = kf(cb- cs)
dCb L = - m dq~,,s
with cb and c~ denoting the bulk liquid and solid-liquid interface adsorbate concentrations, respectively. This boundary condition contains the second important kinetic parameter kf, which represents the liquid film mass transfer coefficient. The parameter kr is a measure for how fast the molecules diffuse across the stagnant liquid film layer. The instantaneous equilibrium between liquid phase and solid phase concentrations at the solid-liquid interface is usually expressed by an isotherm expression, such as that given by the Freundlich equation:
_
_
dt
(1)
dt
where Cb is the bulk liquid adsorbate concentration, qavg the average surface load on the carbon and m and L represent the total mass of carbon in the reactor and the total liquid volume, respectively. Since no adsorbate is allowed to enter or leave the reactor during the process of the experiment, equation (1) states that the change in the bulk liquid concentration causes a proportional change in the average carbon
qs = Kc~
0
r
~"
r=dp/2
Fig. 1. Schematic of mass transport.
(6)
(7)
with K and n being the Freundlich parameters. It is convenient to introduce dimensionless parameters for the HSDM since the absolute values of the two kinetic parameters, Ds and k c, do not give satisfactory information to characterize the adsorption process. The ratio of the kinetic parameters together with other parameters in the HSDM give information about whether surface diffusion or liquid film mass transfer limits the process of transport of adsorbate into the particle. Defining the dimensionless adsorbate liquid phase concentration c*, the dimensionless adsorbed phase concentration q*, the
Parameter estimation for adsorption dimensionless time T and the dimensionless radial coordinate R as: c*
c q; =--;c0 q*=--q0
4Ds
2
T f - ~ p t; R = z r .
(8)
Equations (2) and (3) can be rewritten in a nondimensional form as: q*vg = 3
Io'
q * R 2 dR
aq*(R, T) = 02q*(R, T) 2 0q*(R, T) q 0T OR z R OR
(9)
(10)
with initial and boundary conditions: q*(R, T = O) = 0
aq*(R
= O, T)
~R
aq*(R = 1, T) OR
(12)
= Bi (c~' - c~*)
(13)
(14)
where q0 represents a solid phase concentration which is in equilibrium with a liquid phase concentration co. Equation (13), the boundary condition at the solid-liquid interface, contains the nondimensional Biot number, Bi: B i = kfdpc° .
batch data. This method, however, is time consuming, uncertain and requires a certain feeling or experience relative to how variations in these parameters affect the model output. In the following section a parameter search procedure is presented which is based on sound mathematical foundations and eliminates all uncertainties associated with manual trial and error procedures. After a short explanation of the method, the computational procedure will be applied to a set of real batch adsorption data to find appropriate kinetic parameters kf and Ds. The program will then be extended to search for four variables, namely kf, Ds and the Freundlich parameters, K and n.
(II)
=0
q* = c~*n
269
(15)
2DsPpqo The Biot number represents the ratio of the rate of transport across the liquid layer to the rate of diffusion within the particle. For Bi << 1 external mass transport resistance is the controlling mass transfer step, while for Bi>> 100, surface diffusion is the controlling mass transfer mechanism. Bi numbers between 1 and 100 indicate that both mass transfer mechanisms are important for the particular process (Sontheimvr et al., 1984). The overall mass balance in nondimensional form, integrated to any time T and with appropriate initial condition, is:
SEARCH TECHNIQUE AND THE LEVENBERG-MARQUARDT ALGORITHM
The objective behind the development of a parameter search procedure is to automatically minimize the difference between the model solution, which strongly depends on the parameters and experimental data collected from batch adsorption experiments. Starting from some initial guesses for the values of the unknown parameters, the program should vary these values until the difference between the model solution and the experimental data is a minimum. The parameter values, in this case, are optimum values. This search method is known as a nonlinear least square fit, and can be programmed in using a number of different procedures (Dennis and Schnabel, 1983). The steps involved in obtaining the optimum parameter values are as followed. Starting from some initial estimates of the unknown parameters, the HSDM is solved to compute the closed batch adsorption curve. From this curve, residuals are calculated as: r i ----[C b (ti)experiment- ¢~'(ti)HSDM ].
Their sums squared: ri2
1 = c~ + Dsq,* 8
(16)
with D8 representing the nondimensional solute distribution parameter: D =re_q_0 g L Co"
(17)
All parameters necessary to solve the HSDM are now conveniently grouped in the solute distribution parameter D s and the Biot number Bi. The geometrical and volumetrical quantities (pp, m, L, alp) are directly measurable while K and n can be found from isotherm data. The kinetic parameters D, and kf cannot be measured directly and alternatives for determining these parameters have to be found. The usual method by which these parameters are evaluated is by varying both kinetic coefficients by trial and error until the model solution best mimics the actual experimental
(18)
(19)
i-l
is a measure of how well the model and the actual data agree for a given set of proposed parameters. The mathematical problem is then to minimize the residual function f ( x ) : minimize [f(x)] = ~ ri(x) 2
(20)
where m is the number of data points. The symbol x stands for the unknown parameter vector. The search procedure is iterative in nature since perturbations in the values of the unknown parameters give nonlinear responses of the model solution. The technique applied to solve the minimization problem [equation (20)] is a modified version of the Lavenberg-Marquardt algorithm and is provided by the International Mathematics and Statistical
270
U.K. TRAEGNERand M. T. SUIDAN
Library (IMSL, 1984). This algorithm is a modification of the Newtons method for solving: V 2 f ( X k ) Sk =
- - V f ( x k)
Xk+l = Xk + Sk
(21)
where Xk is the current parameter vector, Xk+~ the newly improved parameter values. Sk in equation (21) represents the Newton step. Vf(Xk) and V2f(xk) represent the gradient and Laplacian operators applied to the residual function defined in equation (20). These derivatives are not available analytically and have to be approximated by finite differences. Levenberg (1944) and Marquardt (1963) introduced several modifications to the full Newton step in order to avoid numerical evaluation of the Hessian matrix. The description of this technique is beyond the scope of this paper and can be found in the cited literature (IMSL, 1984). One major advantage of using the Levenberg-Marquardt algorithm lies in the fact that it is invariant to the solution technique used to solve the HSDM model. Numerical experiments with a finite difference algorithm and an orthogonal collocation technique show identical results for the unknown parameter values. This fact favors the Levenberg-Marquardt algorithm over the quasilinearization technique in conjunction with orthogonal collocation used by Kim et al. (1978). In his work, Kim estimated kinetic coefficients in the pore diffusion adsorption model as applied to data from closed batch experiments. Special features of the Levenberg-Marquardt parameter search procedure as pertaining to the search algorithm, stopping criteria and the acceptability of the iterated Parameter values are presented in the next paragraph. Starting from some initial guess, the procedure suggests improved parameter values according to the change in the magnitude of the residuals with respect to variations in the values of the unknown par-
ameters. These newly proposed parameters are then used to reevaluate the model and the residual funct i o n f ( x ) . This process is repeated until the unknown parameters converge to constant values. Figure 2 depicts the parameter search algorithm. The sum of the residuals is a function of the unknown parameters, x~, and can be seen as a surface spanned by the axes which represent the unknown parameters. The coordinates of the optimum parameter values are the locus of the minimum on the residual surface f ( x ) . The search pathway as suggested by the algorithm while searching for a local minimum on the residual surface is illustrated in Fig. 3. The contours in the plot represent lines of constant residuals. The straight lines on the plot represent the pathways of the iteration procedure. At each end point of the straight lines the algorithm computes the descent direction for a decrease in the magnitude of the residual function, f ( x ) . These steps are repeated until the location of the minimum on the surface along with the corresponding values of parameters are found. It should be pointed out that the number of iterations, i.e. the number of steps (solid lines) involved in finding the minimum on f ( x ) , strongly depends on how close the initial guesses were to their optimum values. Iteration is stopped when the computed parameter values in the new step do not change in the leading four significant digits. The iterated parameter values are then unique solutions to the minimization problem. Confidence region calculations as suggested by Hand et al. (1983) can then be used to provide a sensitivity analysis of the parameter values to the data fit. Final acceptance of the iterated parameter values, however, should be done by plotting the residuals versus time as shown in Fig. 4. The residuals should be scattered in a random fashion. Systematic error will emerge if residuals tend to lie blockwise under or over the time axis (Draper and Smith, 1981).
°
+
++
I READ DATA I" I INITIAL GUESS I OF PARAMETERS x + = Xguess
"I
LEVENBERG MARQUARDT ALGORITHM:
NEW ESTImTES k+l
" "
• • • o" • T
model evaluation . k ,or x
compute residueIs "
•
.
I,,-%(,+l
*
-%(,,)
(expj
jj q I
HSCKI 1' I ~ - - - - - ' ~ '
/ "o
/
[
(HSDM)
I
~_~'~-~-
'
Fig. 2. Flow diagram of the parameter search procedure.
-"
1
Parameter estimation for adsorption
f
~
" ''":
P
~
271
ments were conducted in 6 and 41. of solution for the R6G dye and the p-nitrophenoi, respectively. Mixing was provided by a stirrer, whose speed was kept constant for all the experiments. Duplicate samples were taken for the determination of the initial concentrations and then an accurately weighted mass of GAC was added to the reactor. Samples of 10 ml vol were then withdrawn periodically and filtered. The concentrations of the samples were analyzed by u.v. iight absorption on a spectrophotometer.
~0,9
RESULTS
0
Results for R6G
Fig. 3. Schematic of the search pathway for the optimum kf-Ds values. In this case the experiment should be repeated, or if the quality of the data is not questionable, the HSDM has to be rejected as a suitable model for describing the interaction between a given adsorbate and an activated carbon. The Levenberg-Marquardt algorithm used to minimize the residual function was that provided in the International Mathematics and Statistical Library (IMSL, 1984) and is available in most computer libraries. The calculations for the following examples were performed on a personal computer (IBM-AT) which was found to be powerful enough for performing this task. APPLICATION TO CLOSED EXPERIMENTS
The first solute-carbon combination (Thacker
et al., 1981) is R6G, which is a water soluble red dye, and 30 x 40 U.S. Mesh (dp -- 0.050 cm) 80317 granular activated carbon (Carborundum Co.). The initial concentration and the carbon concentration were co = 0.101 mmol/l and m / L = 0.293 g/l, respectively, while the Freundlich isotherm parameters n and K were 0.035 and 0.362, based on the units: mmol, 1, g-carbon. Applying the algorithm to the norrealized R6G data, the kinetic parameters kf and D, were determined to be 1.022x 10-3cm/s and 1.282 x 10-1°crn2/s. The Blot number in this experiment was 81.20, indicating that mass transfer is controlled by surface as well as external film transport resistance. For Blot numbers higher than 100 (Hand et al., 1983) mass transport is completely controlled by surface diffusion. Figure 5 shows the result of the curve fitted to the R6G adsorption data using the optimum parameter values.
BATCH
Results for PNP
Two applications of the proposed search algorithm are presented in this section. In the first application only the kinetic parameters kf and D, are determined, while the Freundlich isotherm parameters are known from isotherm data. This situation is probably the most common application of the search routine. For the second example, all of the parameters, the kinetic as well as the isotherm parameters, are assumed to be unknown. The simultaneous kinetic and equilibrium parameter estimation was possible since the last points of a set of closed batch experiments represents the isotherm relationship.
The second adsorbate-adsorbent pair (McKinnon, 1977) is p-nitrophenol (PNP) and 60 x 80 U.S. mesh (dp--0.0213cm) F400 granular activated carbon (Calgon Corp.). In this case all four important parameters were simultaneously determined from the computer search. Six adsorption experiments with initial concentrations ranging from 0.455 up to 2.800 mmol/I and a constant carbon concentration of 0.5 g/l were run. In contrast to the first experiment, the batch tests were monitored for up to 8 days in
1.0 MATERIALS
AND
'
i
,
o.a ~ 0 . 1 0 1 E - 3
METHODS
The carbon for both batch experiments was carefully prepared by first washing the sieved carbon in deionized water, in order to remove carbon fines, and then drying at 105°C for several days. The adsorbate solutions were prepared with deionized water and buffered to assure that the solutes were in their molecular form. The batch rate experi-
; mole/l
0.8
0.4 • Da~ 0.3
--
I ~ M ~lutlon
"7
÷ ~
"'t
X
X
~1
X X
~
x
I~
x l
X :d " :~x
-
time Fig. 4. T y p i c a l
residual b e h a v i o r
0.0
,
I O00
,
I 1600
, ?,400
"' plot.
Fig. 5. Closed batch adsorption experimental data and optimum model solution for R6G and 80317 GAC.
272
U.K. TgAEO~ER and M. T. SUIDAN LO
I
10"a
i
C, = 2 . 7 3 5 E - 3
mole/l
0.8
O.B lO.S o 0.t x adsorption resorption
o Da~ --
0.2
HSDM solution 10-4
0.0
I
'
{
EO0
100
Time, min.
order to get equilibrium isotherm data. Since all of the experiments on PNP were done under the same operational conditions, identical kinetic parameters could be expected for all experiments. This is due to the fact that kf, the liquid film mass transfer coefficient, is only a function of the stirrer agitation and D s is only a function of the adsorbate-adsorbent pair. The 8 day period for the batch experiment was assumed to be sufficient to equilibrate the system and the last points of the batch experiments, therefore, belonged to the adsorption isotherm. Additional isotherm points were obtained by running desorption tests on the equilibrated carbon taken from the previous adsorption batch experiments. The search for all four parameters was then carried out by using all adsorption experiments simultaneously as input for the program. The objective function, f ( x ) , was built from the residuals of all of the adsorption experiments together. Figures 6 and 7 show examples of the simultaneous curve fitting procedure, pointing out that the drop in bulk liquid concentration is moderate for the experiment with an initial adsorbate concentration of 2.735x 10-3mol/l, whereas the change in bulk liquid adsorbate concentration was more pronounced for the experiment with an
i
,
, ,,,,,,~
t0-¢
, --
......
t0-~
10"a
Fig. 8. Freundlich isotherm for PNP on F400 GAC, experimental data and best fit obtained from dynamic test. initial adsorbate concentration of 9.343 x I 0- 4 mol/l. The liquid film mass transfer coefficient and the surface diffusion coefficient were computed to be 6.563x 10-3cm/s and 1.915x 10-Scm2/s, respectively. The Freundlich exponent was determined by this procedure as 0.144 and the coefficient as 6.562 x 10 -3, based on the units: mol, 1, g-carbon. The Biot number of 5.23 indicates again that intraparticle and liquid layer mass transfer resistance are important. Only for Bi<< 1 (Sontheimer et al., 1984) could surface diffusion be neglected and the process assumed to be completely controlled by external mass transport. The isotherm, as predicted by the search program along with the final data points of the batch adsorption and desorption experiments, can be seen in Fig. 8. The isotherm data shown in Fig. 8 justify the assumption that the Freundlich equation can be used to model equilibrium. The result for the Freundlich equilibrium parameter values obtained by the simultaneous curve fit are identical to those obtained by separately fitting the Freundlich expression to the isotherm data and searching for only the kinetic constants. The newly determined kinetic and isotherm parameters were then used to predict the desorption experiments. Using the same values for the kinetic parameters kf and D s is justified since the experimental conditions responsible for the liquid
r
mole/l
qo =
0.8
0,8
0,8
0.8
0.4
,,,,,
t0 -~ C. mole/I
L0
i
Co = 9 . 3 4 8 E - 4
, ,,~,,,I
300
Fig. 6. Closed batch adsorption experimental data and optimum model solution for PNP and F400 GAC (co = 2.735 x 10-3 mol/1).
1.0
, 10--a
,
/
1.771E-3
mole/l
~
0.4 Data
0.2
-
,
0.0
0
I
B
100
-
Data
HSDM solution
I
200
--
0.2
SO0
Time, mln.
Fig. 7. Closed batch adsorption experimental data and optimum model solution for PNP and F 4 ~ CAC (co = 9.343 x 10-4 tool/l).
0.0
HSDMsolution
I
I
100
200
300
Time, mi~.
Fig. 9. Closed batch desorption experimental data and optimum model solution for PNP and F400 GAC (co was chosen to be that concentration in equilibrium with q0).
Parameter estimation for adsorption layer are not changed and the surface diffusion coefficient is independent of operation. Figure 9 shows agreement between the predicted and actual desorption experiment. SUMMARY
In cases where the HSDM is an appropriate model for describing the dynamics of adsorption, kinetic parameters must be found before the model may be used for predictive purposes. These kinetic parameters, and for certain cases even the isotherm parameters, can be uniquely determined using the Levenberg-Marquardt numerical algorithm which minimizes the differences between actual batch data and HSDM numerical solutions. Applied to two model components, namely R6G, a high molecular weight and water soluble dye, and p-nitrophenol, the algorithm was able to determine kinetic and equilibrium constants. Both components in this research were found to be limited both by intraparticle and liquid film mass transfer. Acknowledgement--Funding for this work was provided by the U.S. DOE under grant No. DE-FG22-86PC9051 and the U.S. EPA, Cooperative Agreement CR812582-02. REFERENCES
Crittenden J. C. and Weber W. J. (1978) A predictive model for design of fixed-bed adsorbers: single component model verification. J. envir. Engng Div., Am. Soc. cir. Engrs 104, 433-443. Dennis J. E. and Schnabel R. B. (1983) Numerical Methods
WR 23/3---B
273
for Unconstrained Optimization and Nonlinear Equations. Prentice-Hall, Englewood Cliffs, N.J. Draper N. R. and Smith H. (1981) Applied Regression Analysis. Wiley, New York. Hand D. W., Crittenden J. C. and Thacker W. E. (1983) User oriented batch reactor solutions to the homogeneous surface diffusion model. J. envir. Engng Div., Am. Soc. cir. Engrs 109~ 82. International Mathematics and Statistical Library (IMSL) (1984) Math~PC-Library. IMSL, Houston, Tex. Kim B. R., Schmitz R. A., Snoeyink V. L. and Tauxe G. W. (1978) Analysis of models for dichloramine removal by activated carbon in batch and packed-bed reactors using quasilinearization and orthogonal collocation methods. Wat. Res. 12, 317. Levenberg K. (1944) A method for the solution of certain problems in least squares. Quart. appl. Math. 2, 164. Marquardt D. (1983) An algorithm for least-squares estimation of nonlinear parameters. Siam J. appl. Math. 11, 431. Mathews A. and Weber W. J. (1975) Mathematical modeling of multicomponent adsorption kinetics. Presented at the November, 1975 68th Annual Meeting, American Institute of Chemical Engineers, Los Angeles, Calif. McKinnon J. T. (1977) An evaluation of equilibrium in activated carbon adsorption. School of Civil Engineering, Georgia Institute of Technology. Sontheimer H. et al. (1984) Adsorptionsverfahren zur Wasserreinigung. Engler-Bunte-Institut at the University of Karlsruhe, F.R.G. Thacker W. E., Snoeyink V. L. and Crittenden J. C. (1981) Modeling of activated carbon and coal gasification char adsorbents in single-solute and bisolute systems. UILUWRC810161, Research Report 161, University of Illinois at Urbana-Champaign, Ill. Weber W. J. and Chakravorti R. K. (1974) Pore and solid diffusion models for fixed bed adsorbers. Am. Inst. chem. Engrs J. 20, 228-238.
0043-1354/89 $3.00+0.00 Copyright © 1989 Pergamon Press pie
EVALUATION OF SURFACE A N D FILM DIFFUSION COEFFICIENTS FOR CARBON ADSORPTION U. K. TRAEGNER and M. T. SUIDAN~) Department of Civil Engineering, University of Illinois, Urbana, IL 61801, U.S.A. (First received March 1988; accepted in revised form September 1988)
Abstract--The dynamics of adsorption for various adsorbate-adsorbent systems have been predicted successfully using the homogeneous surface diffusion model (HSDM). Kinetic parameters and equilibrium parameters have to be determined in order to solve the model. Batch adsorption experiments are usually used for this purpose. This paper presents a parameter search procedure which is capable of uniquely determining these parameters for the HSDM. The advantage of the proposed parameter search procedure is that it is not limited to cases of dominant transport mechanism, such as cases limited by either external film or surface diffusional resistances. Key words---adsorption, activated carbon, surface diffusion, liquid film diffusion, parameter search, closed batch reactors, homogeneous surface diffusion model
NOMENCLATURE
INTRODUCTION
Bi = Biot number based on surface diffusion (dimensionless) cb(t ) = bulk liquid phase adsorbate concentration (M/L 3) c * (t) = nondimensional bulk liquid concentration cs(t ) = liquid phase adsorbate concentration at solidliquid interface (M/L 3) c*(t)=dimensionless liquid phase concentration at interface co = initial liquid phase concentration (M/L 3) dp = adsorbent particle diameter (L) D s = solute distribution parameter (dimensionless) D, = surface diffusion coefficient (L2/T) f ( x ) = function of residual surface K = Freundlich isotherm capacity constant (L3/M) n kf = liquid film mass transfer coefficient (L/T) L --liquid volume in closed batch test (L3) m = total mass of carbon in closed batch test (M) n = Freundlich isotherm intensity constant (dimensionless) q(r, t)--absorbed phase adsorbate concentration (M/M) qavs= average adsorbed phase adsorbate concentration (M/M) q*vs= nondimensional average surface concentration qs(t) = adsorbed phase adsorbate concentration at solid-liquid interface (M/M) q*(t) = nondimensional adsorbed phase concentration at interface q0 ffi initial adsorbed phase adsorbate concentration (M/M) r = radial coordinate (L) ri(x ) ffi residual at i's observation R --- dimensionless radial coordinate sk ffi Newton step t = physical time (T) T--dimensionless time based on surface diffusion x = parameter vector X k = current parameter vector Xk+~ ffi updated parameter vector p/~= apparent particle density (M/L 3) = partial derivative 02ffi second partial derivative V = grad V2 ffi Laplacian.
The homogeneous surface diffusion model ( H S D M ) has been successfully used to model the dynamics of the adsorption process for various organic compounds on granular activated carbon (GAC). It can be used as a tool for the design engineer who is interested in planning the scope of pilot scale studies, interpreting pilot scale test results or estimating preliminary costs for G A C adsorber units. The mathematical formulation of the H S D M requires the numerical evaluation of differential equations which describe the adsorption process. The H S D M , however, contains physical parameters such as the liquid film mass transfer coefficient, kf, and the surface diffusion coefficient Ds, which cannot be measured directly by analytical means, but have to be known to get unique solutions to the H S D M . The usual way to determine these kinetic parameters is by performing closed batch tests and minimizing the difference between model calculations and actual batch data, while varying the kinetic parameters kf and D,. H a n d et al. 0983) developed a procedure for determining the surface diffusion coefficient by experimentally eliminating liquid film mass transfer resistance and comparing batch adsorption data to empirical solutions to the H S D M . The method presented in this paper is applicable to determining the kinetic parameters without any restriction on mass transport patterns. After presenting the differential equations which govern the dynamics of adsorption and presenting the basic concept of the unconstrained optimization technique, this technique will be applied to find the kinetic coefficients kc and Ds that best describe closed batch data collected for the adsorption of the dye R 6 G on 80317 activated carbon. Furthermore, the search will be applied to data collected on the adsorption o f p - n i t r o p h e n o l on
267
268
U.K. TRAEGNEgand M. T. SUIDAN
F400 activated carbon in closed batch experiments. In this case, the kinetic parameters and the equilibrium constants for the Freundlich adsorption isotherm will be determined, simultaneously. MODEL DEVELOPMENT AND GOVERNING EQUATIONS
The schematic diagram of an activated carbon granule, given in Fig. 1, illustrates the mechanism of mass transport into a spherical particle occurring in the HSDM (Weber and Chakravorti, 1974; Mathews and Weber, 1975; Crittenden and Weber, 1978). The porous particle is surrounded by a stagnant liquid layer which represents a mathematical expression for mass transfer resistances occurring outside the particle. Once the molecules diffuse through the stagnant liquid layer, instantaneous equilibrium is assumed to be established at the outer surface of the carbon particle between the adsorbate molecules in the liquid and adsorbed phases. The molecules are assumed to penetrate further into the particle by creeping along the abundant inner surfaces. The driving force for this movement is the local concentration gradient of the adsorbate along the inner surfaces of the particle. This movement is assumed to be diffusional only and is, therefore, modeled using Fick's first law. The governing equations of the HSDM as applied to the closed batch experiment are the subject of the next section.
surface load qavs. The average carbon load, which is only a function of time, is given by: 3 ffp/2 qavs = (dp/2)3 q(r, t)r 2 dr (2) with q(r, t) representing the concentration along the inner surface of the particle and dp representing the particle diameter. The variable r denotes the radial coordinate with an origin at the center of the particle. The partial differential equation from which q(r, t) is obtained is written in spherical coordinates as:
Oq(r,_~t) = D s ~~ + r { t 3 :t)q ( r '20q(r,~r t)t/
(3)
in which D,, the surface diffusion coefficient, represents a physical parameter. The magnitude of D~ is a measure of how fast the molecules diffuse inside the carbon particle and therefore sets a time scale for the adsorption process. Two boundary conditions and one initial condition have to be specified in order to obtain unique solutions to equation (3). Initially the particle is free of adsorbate, i.e.:
q(r, t) = 0.
(4)
The boundary condition at the center of the particle is:
aq(r = O, t) 0r
= 0
(5)
MODEL EQUATIONS
i.e. no substrate flux across the center. Finally, the continuity of flux at the solid-liquid interface has to be satisfied:
The overall mass balance for the closed batch experiment is:
= dp/2, t) ppDsaq(r Or = kf(cb- cs)
dCb L = - m dq~,,s
with cb and c~ denoting the bulk liquid and solid-liquid interface adsorbate concentrations, respectively. This boundary condition contains the second important kinetic parameter kf, which represents the liquid film mass transfer coefficient. The parameter kr is a measure for how fast the molecules diffuse across the stagnant liquid film layer. The instantaneous equilibrium between liquid phase and solid phase concentrations at the solid-liquid interface is usually expressed by an isotherm expression, such as that given by the Freundlich equation:
_
_
dt
(1)
dt
where Cb is the bulk liquid adsorbate concentration, qavg the average surface load on the carbon and m and L represent the total mass of carbon in the reactor and the total liquid volume, respectively. Since no adsorbate is allowed to enter or leave the reactor during the process of the experiment, equation (1) states that the change in the bulk liquid concentration causes a proportional change in the average carbon
qs = Kc~
0
r
~"
r=dp/2
Fig. 1. Schematic of mass transport.
(6)
(7)
with K and n being the Freundlich parameters. It is convenient to introduce dimensionless parameters for the HSDM since the absolute values of the two kinetic parameters, Ds and k c, do not give satisfactory information to characterize the adsorption process. The ratio of the kinetic parameters together with other parameters in the HSDM give information about whether surface diffusion or liquid film mass transfer limits the process of transport of adsorbate into the particle. Defining the dimensionless adsorbate liquid phase concentration c*, the dimensionless adsorbed phase concentration q*, the
Parameter estimation for adsorption dimensionless time T and the dimensionless radial coordinate R as: c*
c q; =--;c0 q*=--q0
4Ds
2
T f - ~ p t; R = z r .
(8)
Equations (2) and (3) can be rewritten in a nondimensional form as: q*vg = 3
Io'
q * R 2 dR
aq*(R, T) = 02q*(R, T) 2 0q*(R, T) q 0T OR z R OR
(9)
(10)
with initial and boundary conditions: q*(R, T = O) = 0
aq*(R
= O, T)
~R
aq*(R = 1, T) OR
(12)
= Bi (c~' - c~*)
(13)
(14)
where q0 represents a solid phase concentration which is in equilibrium with a liquid phase concentration co. Equation (13), the boundary condition at the solid-liquid interface, contains the nondimensional Biot number, Bi: B i = kfdpc° .
batch data. This method, however, is time consuming, uncertain and requires a certain feeling or experience relative to how variations in these parameters affect the model output. In the following section a parameter search procedure is presented which is based on sound mathematical foundations and eliminates all uncertainties associated with manual trial and error procedures. After a short explanation of the method, the computational procedure will be applied to a set of real batch adsorption data to find appropriate kinetic parameters kf and Ds. The program will then be extended to search for four variables, namely kf, Ds and the Freundlich parameters, K and n.
(II)
=0
q* = c~*n
269
(15)
2DsPpqo The Biot number represents the ratio of the rate of transport across the liquid layer to the rate of diffusion within the particle. For Bi << 1 external mass transport resistance is the controlling mass transfer step, while for Bi>> 100, surface diffusion is the controlling mass transfer mechanism. Bi numbers between 1 and 100 indicate that both mass transfer mechanisms are important for the particular process (Sontheimvr et al., 1984). The overall mass balance in nondimensional form, integrated to any time T and with appropriate initial condition, is:
SEARCH TECHNIQUE AND THE LEVENBERG-MARQUARDT ALGORITHM
The objective behind the development of a parameter search procedure is to automatically minimize the difference between the model solution, which strongly depends on the parameters and experimental data collected from batch adsorption experiments. Starting from some initial guesses for the values of the unknown parameters, the program should vary these values until the difference between the model solution and the experimental data is a minimum. The parameter values, in this case, are optimum values. This search method is known as a nonlinear least square fit, and can be programmed in using a number of different procedures (Dennis and Schnabel, 1983). The steps involved in obtaining the optimum parameter values are as followed. Starting from some initial estimates of the unknown parameters, the HSDM is solved to compute the closed batch adsorption curve. From this curve, residuals are calculated as: r i ----[C b (ti)experiment- ¢~'(ti)HSDM ].
Their sums squared: ri2
1 = c~ + Dsq,* 8
(16)
with D8 representing the nondimensional solute distribution parameter: D =re_q_0 g L Co"
(17)
All parameters necessary to solve the HSDM are now conveniently grouped in the solute distribution parameter D s and the Biot number Bi. The geometrical and volumetrical quantities (pp, m, L, alp) are directly measurable while K and n can be found from isotherm data. The kinetic parameters D, and kf cannot be measured directly and alternatives for determining these parameters have to be found. The usual method by which these parameters are evaluated is by varying both kinetic coefficients by trial and error until the model solution best mimics the actual experimental
(18)
(19)
i-l
is a measure of how well the model and the actual data agree for a given set of proposed parameters. The mathematical problem is then to minimize the residual function f ( x ) : minimize [f(x)] = ~ ri(x) 2
(20)
where m is the number of data points. The symbol x stands for the unknown parameter vector. The search procedure is iterative in nature since perturbations in the values of the unknown parameters give nonlinear responses of the model solution. The technique applied to solve the minimization problem [equation (20)] is a modified version of the Lavenberg-Marquardt algorithm and is provided by the International Mathematics and Statistical
270
U.K. TRAEGNERand M. T. SUIDAN
Library (IMSL, 1984). This algorithm is a modification of the Newtons method for solving: V 2 f ( X k ) Sk =
- - V f ( x k)
Xk+l = Xk + Sk
(21)
where Xk is the current parameter vector, Xk+~ the newly improved parameter values. Sk in equation (21) represents the Newton step. Vf(Xk) and V2f(xk) represent the gradient and Laplacian operators applied to the residual function defined in equation (20). These derivatives are not available analytically and have to be approximated by finite differences. Levenberg (1944) and Marquardt (1963) introduced several modifications to the full Newton step in order to avoid numerical evaluation of the Hessian matrix. The description of this technique is beyond the scope of this paper and can be found in the cited literature (IMSL, 1984). One major advantage of using the Levenberg-Marquardt algorithm lies in the fact that it is invariant to the solution technique used to solve the HSDM model. Numerical experiments with a finite difference algorithm and an orthogonal collocation technique show identical results for the unknown parameter values. This fact favors the Levenberg-Marquardt algorithm over the quasilinearization technique in conjunction with orthogonal collocation used by Kim et al. (1978). In his work, Kim estimated kinetic coefficients in the pore diffusion adsorption model as applied to data from closed batch experiments. Special features of the Levenberg-Marquardt parameter search procedure as pertaining to the search algorithm, stopping criteria and the acceptability of the iterated Parameter values are presented in the next paragraph. Starting from some initial guess, the procedure suggests improved parameter values according to the change in the magnitude of the residuals with respect to variations in the values of the unknown par-
ameters. These newly proposed parameters are then used to reevaluate the model and the residual funct i o n f ( x ) . This process is repeated until the unknown parameters converge to constant values. Figure 2 depicts the parameter search algorithm. The sum of the residuals is a function of the unknown parameters, x~, and can be seen as a surface spanned by the axes which represent the unknown parameters. The coordinates of the optimum parameter values are the locus of the minimum on the residual surface f ( x ) . The search pathway as suggested by the algorithm while searching for a local minimum on the residual surface is illustrated in Fig. 3. The contours in the plot represent lines of constant residuals. The straight lines on the plot represent the pathways of the iteration procedure. At each end point of the straight lines the algorithm computes the descent direction for a decrease in the magnitude of the residual function, f ( x ) . These steps are repeated until the location of the minimum on the surface along with the corresponding values of parameters are found. It should be pointed out that the number of iterations, i.e. the number of steps (solid lines) involved in finding the minimum on f ( x ) , strongly depends on how close the initial guesses were to their optimum values. Iteration is stopped when the computed parameter values in the new step do not change in the leading four significant digits. The iterated parameter values are then unique solutions to the minimization problem. Confidence region calculations as suggested by Hand et al. (1983) can then be used to provide a sensitivity analysis of the parameter values to the data fit. Final acceptance of the iterated parameter values, however, should be done by plotting the residuals versus time as shown in Fig. 4. The residuals should be scattered in a random fashion. Systematic error will emerge if residuals tend to lie blockwise under or over the time axis (Draper and Smith, 1981).
°
+
++
I READ DATA I" I INITIAL GUESS I OF PARAMETERS x + = Xguess
"I
LEVENBERG MARQUARDT ALGORITHM:
NEW ESTImTES k+l
" "
• • • o" • T
model evaluation . k ,or x
compute residueIs "
•
.
I,,-%(,+l
*
-%(,,)
(expj
jj q I
HSCKI 1' I ~ - - - - - ' ~ '
/ "o
/
[
(HSDM)
I
~_~'~-~-
'
Fig. 2. Flow diagram of the parameter search procedure.
-"
1
Parameter estimation for adsorption
f
~
" ''":
P
~
271
ments were conducted in 6 and 41. of solution for the R6G dye and the p-nitrophenoi, respectively. Mixing was provided by a stirrer, whose speed was kept constant for all the experiments. Duplicate samples were taken for the determination of the initial concentrations and then an accurately weighted mass of GAC was added to the reactor. Samples of 10 ml vol were then withdrawn periodically and filtered. The concentrations of the samples were analyzed by u.v. iight absorption on a spectrophotometer.
~0,9
RESULTS
0
Results for R6G
Fig. 3. Schematic of the search pathway for the optimum kf-Ds values. In this case the experiment should be repeated, or if the quality of the data is not questionable, the HSDM has to be rejected as a suitable model for describing the interaction between a given adsorbate and an activated carbon. The Levenberg-Marquardt algorithm used to minimize the residual function was that provided in the International Mathematics and Statistical Library (IMSL, 1984) and is available in most computer libraries. The calculations for the following examples were performed on a personal computer (IBM-AT) which was found to be powerful enough for performing this task. APPLICATION TO CLOSED EXPERIMENTS
The first solute-carbon combination (Thacker
et al., 1981) is R6G, which is a water soluble red dye, and 30 x 40 U.S. Mesh (dp -- 0.050 cm) 80317 granular activated carbon (Carborundum Co.). The initial concentration and the carbon concentration were co = 0.101 mmol/l and m / L = 0.293 g/l, respectively, while the Freundlich isotherm parameters n and K were 0.035 and 0.362, based on the units: mmol, 1, g-carbon. Applying the algorithm to the norrealized R6G data, the kinetic parameters kf and D, were determined to be 1.022x 10-3cm/s and 1.282 x 10-1°crn2/s. The Blot number in this experiment was 81.20, indicating that mass transfer is controlled by surface as well as external film transport resistance. For Blot numbers higher than 100 (Hand et al., 1983) mass transport is completely controlled by surface diffusion. Figure 5 shows the result of the curve fitted to the R6G adsorption data using the optimum parameter values.
BATCH
Results for PNP
Two applications of the proposed search algorithm are presented in this section. In the first application only the kinetic parameters kf and D, are determined, while the Freundlich isotherm parameters are known from isotherm data. This situation is probably the most common application of the search routine. For the second example, all of the parameters, the kinetic as well as the isotherm parameters, are assumed to be unknown. The simultaneous kinetic and equilibrium parameter estimation was possible since the last points of a set of closed batch experiments represents the isotherm relationship.
The second adsorbate-adsorbent pair (McKinnon, 1977) is p-nitrophenol (PNP) and 60 x 80 U.S. mesh (dp--0.0213cm) F400 granular activated carbon (Calgon Corp.). In this case all four important parameters were simultaneously determined from the computer search. Six adsorption experiments with initial concentrations ranging from 0.455 up to 2.800 mmol/I and a constant carbon concentration of 0.5 g/l were run. In contrast to the first experiment, the batch tests were monitored for up to 8 days in
1.0 MATERIALS
AND
'
i
,
o.a ~ 0 . 1 0 1 E - 3
METHODS
The carbon for both batch experiments was carefully prepared by first washing the sieved carbon in deionized water, in order to remove carbon fines, and then drying at 105°C for several days. The adsorbate solutions were prepared with deionized water and buffered to assure that the solutes were in their molecular form. The batch rate experi-
; mole/l
0.8
0.4 • Da~ 0.3
--
I ~ M ~lutlon
"7
÷ ~
"'t
X
X
~1
X X
~
x
I~
x l
X :d " :~x
-
time Fig. 4. T y p i c a l
residual b e h a v i o r
0.0
,
I O00
,
I 1600
, ?,400
"' plot.
Fig. 5. Closed batch adsorption experimental data and optimum model solution for R6G and 80317 GAC.
272
U.K. TgAEO~ER and M. T. SUIDAN LO
I
10"a
i
C, = 2 . 7 3 5 E - 3
mole/l
0.8
O.B lO.S o 0.t x adsorption resorption
o Da~ --
0.2
HSDM solution 10-4
0.0
I
'
{
EO0
100
Time, min.
order to get equilibrium isotherm data. Since all of the experiments on PNP were done under the same operational conditions, identical kinetic parameters could be expected for all experiments. This is due to the fact that kf, the liquid film mass transfer coefficient, is only a function of the stirrer agitation and D s is only a function of the adsorbate-adsorbent pair. The 8 day period for the batch experiment was assumed to be sufficient to equilibrate the system and the last points of the batch experiments, therefore, belonged to the adsorption isotherm. Additional isotherm points were obtained by running desorption tests on the equilibrated carbon taken from the previous adsorption batch experiments. The search for all four parameters was then carried out by using all adsorption experiments simultaneously as input for the program. The objective function, f ( x ) , was built from the residuals of all of the adsorption experiments together. Figures 6 and 7 show examples of the simultaneous curve fitting procedure, pointing out that the drop in bulk liquid concentration is moderate for the experiment with an initial adsorbate concentration of 2.735x 10-3mol/l, whereas the change in bulk liquid adsorbate concentration was more pronounced for the experiment with an
i
,
, ,,,,,,~
t0-¢
, --
......
t0-~
10"a
Fig. 8. Freundlich isotherm for PNP on F400 GAC, experimental data and best fit obtained from dynamic test. initial adsorbate concentration of 9.343 x I 0- 4 mol/l. The liquid film mass transfer coefficient and the surface diffusion coefficient were computed to be 6.563x 10-3cm/s and 1.915x 10-Scm2/s, respectively. The Freundlich exponent was determined by this procedure as 0.144 and the coefficient as 6.562 x 10 -3, based on the units: mol, 1, g-carbon. The Biot number of 5.23 indicates again that intraparticle and liquid layer mass transfer resistance are important. Only for Bi<< 1 (Sontheimer et al., 1984) could surface diffusion be neglected and the process assumed to be completely controlled by external mass transport. The isotherm, as predicted by the search program along with the final data points of the batch adsorption and desorption experiments, can be seen in Fig. 8. The isotherm data shown in Fig. 8 justify the assumption that the Freundlich equation can be used to model equilibrium. The result for the Freundlich equilibrium parameter values obtained by the simultaneous curve fit are identical to those obtained by separately fitting the Freundlich expression to the isotherm data and searching for only the kinetic constants. The newly determined kinetic and isotherm parameters were then used to predict the desorption experiments. Using the same values for the kinetic parameters kf and D s is justified since the experimental conditions responsible for the liquid
r
mole/l
qo =
0.8
0,8
0,8
0.8
0.4
,,,,,
t0 -~ C. mole/I
L0
i
Co = 9 . 3 4 8 E - 4
, ,,~,,,I
300
Fig. 6. Closed batch adsorption experimental data and optimum model solution for PNP and F400 GAC (co = 2.735 x 10-3 mol/1).
1.0
, 10--a
,
/
1.771E-3
mole/l
~
0.4 Data
0.2
-
,
0.0
0
I
B
100
-
Data
HSDM solution
I
200
--
0.2
SO0
Time, mln.
Fig. 7. Closed batch adsorption experimental data and optimum model solution for PNP and F 4 ~ CAC (co = 9.343 x 10-4 tool/l).
0.0
HSDMsolution
I
I
100
200
300
Time, mi~.
Fig. 9. Closed batch desorption experimental data and optimum model solution for PNP and F400 GAC (co was chosen to be that concentration in equilibrium with q0).
Parameter estimation for adsorption layer are not changed and the surface diffusion coefficient is independent of operation. Figure 9 shows agreement between the predicted and actual desorption experiment. SUMMARY
In cases where the HSDM is an appropriate model for describing the dynamics of adsorption, kinetic parameters must be found before the model may be used for predictive purposes. These kinetic parameters, and for certain cases even the isotherm parameters, can be uniquely determined using the Levenberg-Marquardt numerical algorithm which minimizes the differences between actual batch data and HSDM numerical solutions. Applied to two model components, namely R6G, a high molecular weight and water soluble dye, and p-nitrophenol, the algorithm was able to determine kinetic and equilibrium constants. Both components in this research were found to be limited both by intraparticle and liquid film mass transfer. Acknowledgement--Funding for this work was provided by the U.S. DOE under grant No. DE-FG22-86PC9051 and the U.S. EPA, Cooperative Agreement CR812582-02. REFERENCES
Crittenden J. C. and Weber W. J. (1978) A predictive model for design of fixed-bed adsorbers: single component model verification. J. envir. Engng Div., Am. Soc. cir. Engrs 104, 433-443. Dennis J. E. and Schnabel R. B. (1983) Numerical Methods
WR 23/3---B
273
for Unconstrained Optimization and Nonlinear Equations. Prentice-Hall, Englewood Cliffs, N.J. Draper N. R. and Smith H. (1981) Applied Regression Analysis. Wiley, New York. Hand D. W., Crittenden J. C. and Thacker W. E. (1983) User oriented batch reactor solutions to the homogeneous surface diffusion model. J. envir. Engng Div., Am. Soc. cir. Engrs 109~ 82. International Mathematics and Statistical Library (IMSL) (1984) Math~PC-Library. IMSL, Houston, Tex. Kim B. R., Schmitz R. A., Snoeyink V. L. and Tauxe G. W. (1978) Analysis of models for dichloramine removal by activated carbon in batch and packed-bed reactors using quasilinearization and orthogonal collocation methods. Wat. Res. 12, 317. Levenberg K. (1944) A method for the solution of certain problems in least squares. Quart. appl. Math. 2, 164. Marquardt D. (1983) An algorithm for least-squares estimation of nonlinear parameters. Siam J. appl. Math. 11, 431. Mathews A. and Weber W. J. (1975) Mathematical modeling of multicomponent adsorption kinetics. Presented at the November, 1975 68th Annual Meeting, American Institute of Chemical Engineers, Los Angeles, Calif. McKinnon J. T. (1977) An evaluation of equilibrium in activated carbon adsorption. School of Civil Engineering, Georgia Institute of Technology. Sontheimer H. et al. (1984) Adsorptionsverfahren zur Wasserreinigung. Engler-Bunte-Institut at the University of Karlsruhe, F.R.G. Thacker W. E., Snoeyink V. L. and Crittenden J. C. (1981) Modeling of activated carbon and coal gasification char adsorbents in single-solute and bisolute systems. UILUWRC810161, Research Report 161, University of Illinois at Urbana-Champaign, Ill. Weber W. J. and Chakravorti R. K. (1974) Pore and solid diffusion models for fixed bed adsorbers. Am. Inst. chem. Engrs J. 20, 228-238.