Analytical solution for networks of irreversible first-order reactions

PII: S0043-1354(98)00273-5

Wat. Res. Vol. 33, No. 3, pp. 814±826, 1999 # 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0043-1354/98/$ - see front matter

ANALYTICAL SOLUTION FOR NETWORKS OF IRREVERSIBLE FIRST-ORDER REACTIONS GERALD R. EYKHOLT Department of Civil and Environmental Engineering, University of Wisconsin ± Madison, Madison, WI 53706, U.S.A. (First received August 1997; accepted in revised form June 1998) AbstractÐA new analytical solution technique for networks of irreversible chemical reactions in batch reactors with ®rst-order kinetics is presented. The solution technique involves the iterative calculation of a coecient matrix that is then applied for direct calculation of species concentrations at any time. The derivation, which is based on Laplace transforms of an arbitrary network, is presented in the appendix. Several examples and extensions to the solution are discussed. Extensions made to the solution include the consideration of phase equilibria, ®rst-order and zero-order sources and sinks, common reaction products outside the network, and combined networks of reactions and CSTRs. Rate-limited mass transfer e€ects and reversible reactions were not considered. The solution is general, compact, easy to use, and allows several extensions that will likely aid scientists and engineers in solving problems related to the dynamics of these networks. # 1998 Elsevier Science Ltd. All rights reserved Key wordsÐLaplace transforms, kinetics, parameter estimation, ®rst-order, zero-order, retardation

INTRODUCTION

Irreversible, sequential ®rst-order reactions are commonly used to model the degradation of chemicals. Networks of irreversible ®rst-order reactions are important for reductive dechlorination of chlorinated ethenes (Tondoi et al., 1994; Roberts et al., 1996), the production and depletion of radioisotopes during radioactive decay (Benedict et al., 1981), and for processing of polymers (Walas, 1991). Analytical solutions and derivations are widely available for simple systems, such as the production of a product C3 from C1 through its daughter product C2: k12

k23

C1 ÿ ÿ4 ÿ C2 ÿ ÿ4 ÿ C3 ,

…1†

where k12 and k23 are the ®rst-order degradation rate constants for C1 and C2. The analytical solution and its derivation are readily available (Holland and Anthony, 1989; Walas, 1991; Schnoor, 1996; Chapra, 1997). Because of the complexity of the integration problem, large systems of irreversible, parallel and consecutive ®rst-order reactions are commonly modeled numerically. In addition, numerical solutions are commonly invoked for systems with mass transport, second-order or mixed-order kinetics, or sources/ sinks from the system. One may also model complex processes with composite components and mechanisms, with a technique called kinetic lump814

ing. Li et al. (1996) demonstrated the use of kinetic lumping through several example problems relevant to wastewater treatment and chemical engineering processes. The method commonly involves ®rstorder networks and lumped parameter estimation (rather than providing full description of the dynamics of all components). An analytical solution for large networks of irreversible, ®rst-order reactions is practiced widely in nuclear engineering. The Bateman equation, ®rst published in 1910, is used to describe pure decay in straight decay chains (Benedict et al., 1981): Ni ˆ N 0m eÿli t ‡

iÿ1 X

N 0m mˆ1

iÿ1 Y jˆm

! lj

i X jˆm

eÿlj t i Y

…2†

…lk ÿ lj †,

kˆm k6ˆj

where Ni is the number of atoms of the daughter product i, and li is the ®rst-order decay constant. The index of the ®rst compound in the chain is one, and increases by one for each consecutive product. More complicated decay networks can be solved by decomposition of the network to linear chains and by summing the contributions to Ni from each chain. Solutions are also available for zero-order source/sinks (Benedict et al., 1981).

Solution for irreversible reactions

Since the set of ®rst-order di€erential equations is linear and homogeneous equations can be formed, the problem can also be solved as an eigenproblem (Walas, 1991; Chapra, 1997). This matrix method also applies to reversible systems and can be coded into a computer program equipped with matrix algebra routines. A set of linear equations yields solution roots, called eigenvalues. Calculation is generally faster than numerical integration for large systems. In addition, the general solution can be applied to a wide variety of initial conditions without much additional computational e€ort. The analytical solution for networks of ®rst-order systems can also be derived by Laplace transforms (Churchill, 1974; Walas, 1991). Resigno and Segre (1966) have provided a text on the use of Laplace transforms, with many examples of explicit, analytical solutions to pharmacokinetics problems. Using Laplace transforms, an explicit solution for a large network of irreversible ®rst-order reactions was derived. After several re®nements to the algebraic form of this particular solution, a general analytical solution technique for irreversible networks of ®rst-order equations was found. Like the Bateman equation, the solution is iterative. However, the solution strategy does not require network decomposition and superposition. This paper describes the new solution technique, gives several examples and extensions, and discusses the implications for reaction modeling and parameter estimation. It is important to note, however, that the new solution technique can only be used to approximate processes that have nonlinear rates or reversible reactions. Mixed-order and higher order processes are common in biological systems, and the applicability of the irreversible, ®rst-order network model will vary with each system.

DESCRIPTION OF GENERAL SOLUTION TECHNIQUE

Consider the following rate expression for a closed system of irreversible, ®rst-order reactions:  kRj‡N X i
…3†

… for all values of j †: The ®rst term on the right-hand-side represents the combined rate of production of nj moles of species j from ni moles of all possible parent species i. The second term represents the rate of depletion of species j toward the formation of other species k, where Nk is the number of daughter products from species j. Here, the index of any product species ( j or k) is always greater than the index of a parent species (i or j). For example, the above rate description would apply to the following network

815

of reactions:

…4†

The ®rst-order rate constants, kij, are considered positive for degradation reactions. The formulation also allows consideration of ®rst-order loss (or production) from outside the system, kjjCj, where kjj is positive for loss to outside the system. The rate constants applying directly to the depletion of species j can be summed into one rate constant, kj, kj ˆ

kRj‡N 1 Xk kjk : Rj kˆj

…5†

The retardation factor Rj will be described below for an extended solution, but can be set to one for this case. This step will simplify the derivation and form of the general solution. Given ®xed initial concentrations Bj for a wellmixed, constant-volume batch system, a general analytical solution for a network of ®rst-order reactions is Cj …t† ˆ

iRj X Dij eÿki t

… for all values of j †:

…6†

iˆ1

The coecient matrix Dij is described in three parts. The ®rst part describes diagonal elements, Djj ˆ Bj ÿ

i< j X Dij :

…7†

iˆ1

The ®rst element, D11, will be de®ned from the initial conditions and remaining elements will be de®ned upon calculation of the second and third parts. The second part is for direct links between all parent species i and daughter species j:   nj kij Dii Dij ˆ … for all i
i

kj ÿ km

… for all m
…9†

The summation in equation 9 applies only to multiple parents i of species j. The coecient equations (7±9) are closed form and are solved iteratively with increasing values of j. The Dmj coef®cients are zero for species m that do not have direct or indirect links with species j (for example,

816

Gerald R. Eykholt

species 2 and 3 from equation 4). However, no additional logic is needed to distinguish unrelated species m, since either kmj, kmi or Dmi will be zero for all i if m is not a parent or ancestor to j. Note also that equation 8 can be expressed by equation 9 if m is allowed to equal i. A more compact expression, for direct links, indirect links, and unrelated species is nj Dij ˆ

m
coecient matrix can be solved in the following order: D11 ˆ B1 ˆ C1 …0† k12 D11 k12 C1 …O† ˆ k2 ÿ k1 k23 ÿ k12

D12 ˆ

kj ÿ ki

…10†

D13 ˆ

…all possible links for i
Example 1 Consider the system described by equation 1, the consecutive ®rst-order reaction of C1 to C2 to C3. The general solution (equation 6) applies and the

…12†

k12 C1 …0† k23 ÿ k12

…13†

k23 D12 k23 D12 k23 C1 …0† ˆ ˆÿ k23 ÿ k12 k3 ÿ k1 0 ÿ k12

…14†

D22 ˆ B2 ÿ D12 ˆ C2 …0† ÿ

kmj Dim =nm

mˆi

…11†

D23 ˆ

k23 D22 k23 D22 ˆ ˆ ÿD22 k3 ÿ k2 0 ÿ k23

…15†

D33 ˆ B3 ÿ D23 ÿ D13 ˆ C3 …0† ‡ C2 …0† ÿ

k12 C1 …0† k23 C1 …0† ‡ k23 ÿ k12 k23 ÿ k12

ˆ C3 …0† ‡ C2 …0† ‡ C1 …0†:

…16†

The species concentrations are then calculated directly from equation 6. This solution is equivalent to earlier solutions (i.e., Schnoor, 1996; Chapra, 1997). An example solution and the corresponding rate constant and D matrices are shown in Fig. 1. Example 2 Consider the system described by equation 4. Since part of the network is similar to the network described in Section 3.1, part of the coecient matrix is applicable with minor changes (equations 11±13). However, k1 now equals the sum of k12 and k13 (see equation 5) and the stoichiometric coecients for terms including species 2 must be considered. The ®rst element (D11) of the coecient matrix is unchanged, and the remaining terms can be solved in the following order: D12 ˆ

n2 k12 D11 2k12 B1 ˆ n1 …k2 ÿ k1 † k24 ÿ k1

…17†

D22 ˆ B2 ÿ D12

…18†

D13 ˆ

k13 D11 k13 B1 ˆ k3 ÿ k1 k34 ÿ k1

D33 ˆ B3 ÿ D13 D14 ˆ

…19† …20†

…n4 =n2 †k24 D12 ‡ k34 D13 k4 ÿ k1 …k24 D12 =2† ‡ k34 D13 k1

…21†

n4 k24 D22 D22 ˆÿ n2 …k4 ÿ k2 † 2

…22†

ˆÿ D24 ˆ

Solution for irreversible reactions

817

Fig. 1. Kinetics solution example for sequential ®rst-order decay of species 1 to species 3 (see equations 1 and 11±16).

D34 ˆ

k34 D33 ˆ ÿD33 k4 ÿ k3

…23†

D44 ˆ B4 ÿ D34 ÿ D24 ÿ D14 ˆ B4 ‡ D22 ‡ D33 ÿ D14 :

…24†

The ®nal term, D44, can be reduced further to yield the sum of the initial concentrations, B1 through B4. An example solution and the corresponding rate constant and D matrices are shown in Fig. 2. Example 3 Roberts et al. (1996) have described the pathway for the reductive dechlorination and elimination of

chlorinated ethenes on zero-valent metals according to Fig. 3. Roberts et al. (1996) and Arnold (1997) have modeled the complex reaction kinetics of this system numerically using Scientist2 (MicroMath Scienti®c Software, Salt Lake City, Utah). For this example, Scientist2 was used with the integration settings of episode (sti€) integration, analytical Jacobian calculation, relative error of 1.0  10ÿ6, and a maximum step size of 5000. Preliminary estimates for rate constants were based on multiple laboratory experiments with aqueous solutions and zinc metal, and were expressed for a metal surface area of 1.0 m2 mLÿ1. Coecients for the analytical solution, shown in Fig. 4, were calculated using the rate constant estimates and an initial PCE concentration of 1.0 mM. Analytical and numerical sol-

Fig. 2. Kinetics solution example for parallel ®rst-order decay of species 1 to species 4 through species 2 and 3 (see equations 4 and 17±24).

818

Gerald R. Eykholt

Fig. 3. Mechanism for reductive dechlorination and elimination of chlorinated ethene system proposed by Roberts et al. (1996).

utions are shown in Fig. 5. The agreement between the solutions was excellent. Carbon mass balance errors were less than 1  10ÿ15 M for the analytical solution.

Extensions to the general analytical solution: Zero-order source/sinks and phase equilibria. Several extensions can be made to the general solution, including the consideration of zero-order sources

Fig. 4. First-order rate constants and D coecient matrices for the reductive dechlorination and elimination of chlorinated ethene mechanism proposed by Roberts et al. (1996). Initial concentration of PCE was 200 mg/L. Rate constants should be regarded as preliminary estimates.

Solution for irreversible reactions

819

Fig. 5. Agreement between analytical and numerical solutions for the reductive dechlorination and elimination pathway for the chlorinated ethene system contacted with zinc metal. Initial concentration of PCE was 200 mM. The reaction mechanism, rate constants and D coecient matrix are described in Figs 3 and 4. Numerical solution provided by Arnold (1997).

and sinks and multiphase systems in equilibrium with the reactant phase. The rate expression of daughter species j with a single parent i can be expressed: X dCj nj Rj kjk Cj ‡ pj , …25† ˆ kij Ci ÿ dt ni k where pj is the zero-order rate of production of species j per unit volume and Rj is a retardation factor describing instantaneous, linear equilibrium with other phases: Rj ˆ 1 ‡

mKd Vg K 0H ‡ : V V

…26†

The parameters m, V, and Vg refer to the solid mass, solution volume, gas volume, respectively. The equilibrium partition coecient for sorption to the solid is Kd, and K'H is the dimensionless Henry's Law constant for equilibrium between solution and the gas phase (which could be di€erent for each species j). The modi®ed rate expression is obtained if each side of equation 25 is divided by Rj: dCj nj kij pj Ci ÿ kj Cj ‡ : ˆ dt ni Rj Rj

…27†

The analytical solution for networks with these rate dynamics was derived in the same manner as the simpler solution. For consideration of a network of ®rst-order reactions with zero-order sources/sinks and multiphase systems in equilibrium with the reactant

phase, a new analytical solution is presented: Cj …t† ˆ Zj …Wj ‡ …1 ÿ Wj †t† ‡

iRj X Dij eÿki t iˆ1

…28†

… for all values of j †, where Wj is the maximum reaction order (0 or 1) for species j and ki is the lumped, ®rst-order rate constant de®ned by equation 5. Species j may be a terminal species, but all parent and ancestor species i are assumed to have a maximum reaction order of 1. The term Zi is a zero-order source/sink vector divided by kj and the retardation factor Rj: pj ‡ nj Zj ˆ

i< j X

kij Zi =ni

iˆ1

Rj …Wj kj ‡ …1 ÿ Wj ††

:

…29†

If the maximum order Wj is one, the solution de®ned by equation 28 expresses all time-dependency in the exponential terms and equation 29 requires a non-zero ®rst-order rate constant for each species j. The coecient matrix Dij requires only minor modi®cation, where the diagonal elements now re¯ect the contribution of the zeroorder terms: Djj ˆ Bj ÿ Wj Zj ÿ

i< j X Dij ,

…30†

iˆ1

where the initial aqueous concentration Bj can be calculated from the sum of mass of j in every phase

820

Gerald R. Eykholt

divided the solution volume V and Rj. The other matrix elements are also modi®ed in order to consider Rj:  Dij ˆ

nj ni



kij Dii Rj …kj ÿ ki †

works. Ignoring zero-order source/sinks for chloride, the general solution for network producing the common product chloride (such as the reaction network shown in Fig. 3) is: ‰Cl Št ˆ DCl,0 eÿkCl,0 t ‡

…direct links between parent i and daughter j †, …31† and nj Dmj ˆ

X

kij Dmi =ni

i

Rj …kj ÿ km †

…32†

iX
DCl,i eÿki t ,

where kCl,0 is the ®rst-order loss rate constant of chloride divided by the chloride retardation factor RCl (refer to equation 26), N is the number of species in the network (not counting the common product), and the common product coecient vector DCl,i is de®ned as:

…indirect links of ancestor m
DCl,i ˆ ÿ

jX
kCl,j Dij

jˆi‡1

Again, a more compact expression for direct links, indirect links, and unrelated species can be expressed nj Dij ˆ

m
kmj Dim =nm

mˆ1

Rj …kj ÿ ki †

…33†

…all possible links for i
…34†

iˆ1

kCl,i ÿ kCl,0

…for 1
…35†

and DCl,0 ˆ BCl ÿ

iX
DCl,i ,

…36†

iˆ1

where BCl is the total initial chloride mass in the system divided by the volume and RCl. The vector kCl,i relates to reaction rate constants from the network (kij), scaled by the number of moles of common product (nCl,ij) produced in the reaction of one mole of parent i to daughter j: jRN X

kCl,i ˆ

kij nCl,ij

jˆi‡1

RCl

…for 1
…37†

The concentrations of other common products and reactants can be calculated similarly. This common product balance may also be important in nuclear engineering for estimating the heat ¯ux in a reactor. It is important to note that diagonal terms of the network D coecients are not included in this analysis. Therefore, common product mass balances will re¯ect a mass balance error that can be attributed to ®rst-order or zero-order sources or sinks, or for ®rst-order reactions to other products outside the network. Combined networks for simulating chemical reaction networks in CSTR networks. Mixing e€ects in non-ideal reactors are commonly addressed by modeling networks of continuously stirred tank reactors, CSTRs (or completely mixed ¯ow reactors, CMFRs). CSTRs are also common compartments in water quality models (Chapra, 1997). As CSTR networks without reaction can be described by a set of ®rst-order and zero-order processes, the general analytical solution technique discussed above for reaction networks should also apply. Transformation in a CSTR with a step input of Cj,IN and a retention time of y can be described by

Solution for irreversible reactions

821

Fig. 6. Conceptual development of a combined network for modeling irreversible ®rst-order reactions in a CSTR network. Reaction networks are placed in each CSTR, components are connected between CSTRs to describe mass transfer, and a new numbering scheme is generated for the combined network.

the di€erential equation,

dCj 1 1 dCj : ˆ Cj,IN ÿ Cj ‡ Rj y y dt dt rxn

…38†

The step input can be modeled as a zero-order source per unit volume, pj, and the mass-transfer parameter yÿ1 can be treated as a ®rst-order sink kjj. However, as will be shown below, this formulation should be modi®ed for networks of CSTRs, since a portion of the ®rst-order sink may represent a ®rst-order source to another CSTR. On the hypothesis that irreversible networks of CSTRs and ®rst-order reactions may be expressed in the same mathematical form, it is possible to combine the networks such that the new analytical solution technique applies. The general concepts of the combined network model are (1) to reproduce the reaction network for as many compartments (CSTRs) in which the reactions apply, and (2) to provide the necessary mass-transfer linkages for each reaction network component, from one compartment to another. Conceptual diagrams for the development of a combined network for nitri®cation in a CSTR network are shown in Fig. 6. Here,

the same reaction network applies within CSTR compartments, and each component of the reaction network in one CSTR has a ®rst-order mass transfer linkage to the same reaction network position in another CSTR. Finally, the numbering scheme for species is revised such that general analytical solution technique can be applied over the combined network. As discussed above, step inputs to a CSTR can be modeled as a zero-order source. First-order mass transfer from one CSTR to another can be handled as a ®rst-order sink, as a ®rst-order linkage to the other reactor (analogous to the reaction mechanism rate constant kij), or as a combination of both. In order for the same network solution technique to work, the following conditions apply. First of all, the ®rst-order constant describing mass transfer out of the source CSTR i can be handled as the sum of a ®rst-order sink term, kii, and the mass transfer contribution to any other downstream CSTR j, kij. Secondly, the ®rst-order mass transfer contribution from CSTR i to the receiving CSTR j will be described by the parameter kij. For the source CSTR i with retention time yi and volume Vi, the

822

Gerald R. Eykholt

Table 1. Listing and description of constants used for a sample calculation of the combined nitri®cation reaction and CSTR ¯ow network shown in Fig. 6 Network model constants yA, yB, yC VA, VB, VC QA, QB, QC kv k12, k45, k78 k23, k56, k89 k17, k28, k39 k47, k58, k69 k11 k22, k33 k44 k55, k66 k77 k88, k99 p1 p4 Rj Bj

Description

Value(s) for example simulation

Retention times for CSTR A, CSTR B, and CSTR C 5.0, 8.33 and 10.0 d Reactor volumes for CSTR A, CSTR B, and CSTR C 20, 30 and 76 L Reactor ¯ow rates for CSTR A, CSTR B, and CSTR C 4.0, 3.6 and 7.6 L/d ®rst-order rate constant for volatilization loss of NH3±N 0.0100 dÿ1 ®rst-order rate constant for conversion of NH3±N to NOÿ 0.400 dÿ1 2 ±N ÿ ®rst-order rate constant for conversion of NOÿ 0.240 dÿ1 2 ±N to NO3 ±N ®rst-order coe€, for mass transfer from CSTR A to CSTR C, scaled to 0.05263 dÿ1 volume in CSTR C ®rst-order coe€, for mass transfer from CSTR B to CSTR C, scaled to 0.04737 dÿ1 volume in CSTR C kv (volatilization loss of NH3±N) added to ®rst-order sink compensation 0.1574 dÿ1 for changing volume from CSTR A to CSTR C ®rst-order sink compensation for changing volume from CSTR A to 0.1474 dÿ1 CSTR C kv (volatilization loss of NH3±N) added to ®rst-order sink compensation 0.08263 dÿ1 for changing volume from CSTR B to CSTR C ®rst-order sink compensation for changing volume from CSTR B to 0.07263 dÿ1 CSTR C combined ®rst-order coe€, for mass transfer from CSTR C to exit (1/yC) 0.110 dÿ1 and volatilization loss of NH3
N (kv) ®rst-order coe€. for mass transfer from CSTR C to exit 0.100 dÿ1 zero-order input rate (source) of NH3±N per unit volume for CSTR A 0.5 for t R15 d, 0.0 for t>15 d (mg dÿ1 Lÿ1) zero-order input rate (source) of NH3±N per unit volume for CSTR B 0.2 for t R15 d, 0.0 for t>15 d (mg dÿ1 Lÿ1) retardation coe€. for all species j 1.0 (no sorption) total initial concentration of species j expressed as aqueous B1=4 mg N/L, B4=2 mg N/L, all other Bj=0.

Fig. 7. First-order network constants and D coecient matrices for combined network problem described in Fig. 6 and Table 1. D coecient matrix applies to loading cycle only.

Solution for irreversible reactions

823

Fig. 8. Nitri®cation dynamics within a CSTR network described by the general analytical solution technique (see Figs 6 and 7, Table 1 for conditions).

mass transfer parameter for any component to another CSTR j with retention time yj and volume Vj, is kij ˆ

Vi : Vj yi

…39†

The ®rst-order sink term for CSTR i is then found by adding the normal ®rst-order sink kii,s/s (that may occur from processes such as volatilization) to the compensation for the possible imbal-

ance in reactor volumes, kii ˆ kii,s=s ‡

 X  Vj kij ÿ1 : Vi j

…40†

The ®rst-order mass transfer from the terminal CSTR in the network can be modeled with the normal ®rst-order sink parameter (kjj=1/yj). With these minor modi®cations and renumbering of the combined reaction network, the general analytical solution technique (equations 28±33) can be used.

824

Gerald R. Eykholt

An example calculation for the combined network shown in Fig. 6 was performed. Combined network parameters selected for the calculation are shown in Table 1. The CSTR reactors A and B have initial ammonia concentrations and receive step inputs of ammonia nitrogen over a period of 15 d. After 15 d, the ammonia input is ceased. Flow from CSTR A and CSTR B are combined in CSTR C (constant ¯ows, constant volumes). Volatile losses of ammonia are modeled by a ®rst-order sink term, and a sequential ®rst-order nitri®cation model from ammonia to nitrite to nitrate is used. Although di€erent reaction rate, sorption, volatilization rate parameters could be selected for each CSTR, these parameters were held constant over the mixing domain. The matrix of ®rst-order parameters kij for reaction and mass transfer in the combined network, and the D coecient matrix for the 15 d loading cycle are shown in Fig. 7. The analytical solution, calculated from a spreadsheet application, is shown in Fig. 8. The ammonia nitrogen concentration in each CSTR at 15 d was in agreement with the value calculated with the steady-state approximation. Implications for reaction modeling and parameter estimation. There are several important implications for reaction modeling and parameter estimation that may extend from broader use of the analytical solution. The analytical solution is easier to implement than existing solutions. Since the solution can be implemented on a spreadsheet, users will be able to consider complicated rate dynamics in systems with zero-order and ®rst-order sources/sinks and phase equilibria without use of more sophisticated numerical techniques. The analytical solution may also be implemented in more complicated reactor systems, numerical models, and for other simulations. For instance, Monte Carlo models for reactive barrier design (Eykholt, 1997) involving uncertainty in the chlorinated ethene reduction mechanism (Fig. 3) can now be developed. The analytical solution may also promote better estimates of reaction parameters from laboratory experiments. Analysis and design of experiments with volatile organic chemicals is challenging. Most careful researchers employ regular use of controls, blanks, and matrix spikes to improve quality assurance. The analytical solution allows the explicit consideration of rate dynamics observed from controls and blanks, so better estimates of the reaction dynamics will be obtained after proper adjustment of the zero- and ®rst-order source/sink terms. One may also be able to use multiple species for the prediction of mixed-rate kinetics of some species. If a particular species in the network does not follow ®rst-order kinetics, it may be possible to model the dynamics of that species and its products with imaginary combinations within the network. Future work is planned to develop functionally equivalent ®rst-order and zero-order networks to

model kinetics that are not ®rst-order or strictly irreversible. It may also be possible to model reversible systems with expanded network numerical approximations (Wilson and Henderson, 1997). CONCLUSION

An analytical solution for modeling the dynamics of networks with irreversible, ®rst-order reactions has been derived. The solution is general, compact, easy to use, and allows several extensions that will likely aid scientists and engineers in solving problems related to the dynamics of these networks. For instance, one working to describe the mechanisms for reductive dechlorination and elimination reactions with metallic zinc would be better able to consider alternate reaction mechanisms given di€erent starting conditions (di€erent network con®gurations), more rigorously model the impacts of sources and sinks and their e€ect on the experimental rate parameters, check mass balances, correct for sorption and gas±liquid equilibria, test the reactions in CSTRs, and perform a sensitivity analysis for better experimental design. Future work is being planned to extend the analytical solution in order to consider behavior of chemical networks in non-ideal reactors, species with mixed-order reaction kinetics, reversible systems, and systems governed by mass-transfer limitations. However, a multi-parameter model requires multiple experiments and careful controls in order to remain unique. The usefulness of the simulations will depend greatly on experimental design and the quality of the measurements. AcknowledgementsÐThis work was supported in part by a National Science Foundation CAREER award (No. 9625087) to the author. The author wishes to thank Professor A. Lynn Roberts and William Arnold of Johns Hopkins University for providing the kinetic parameter set from chlorinated ethene reactions on zinc metal, early review, and preliminary numerical solutions. Professor Thatcher Root and Professor Douglass Henderson of the University of Wisconsin also provided review of the analytical solution and discussed possible extensions. REFERENCES

Arnold W. A. (1997) Personal communication with G. R. Eykholt, 6/18/97. Benedict M., Pigford T. H. and Levi H. W. (1981) Nuclear Chemical Engineering, 2nd edn. McGraw-Hill, New York. Chapra S. C. (1997) Surface Water-Quality Modeling. McGraw-Hill, New York. Churchill R. V. (1974) Operational Mathematics, 3rd edn. McGraw-Hill, New York. Eykholt G. R. (1997) Uncertainty-based scaling of iron reactive barriers. In: In Situ Remediation of the Geoenvironment, Am. Soc. of Civil Engrs., Geotechnical Technical Publication No. 71, pp. 41±55. Holland C. D. and Anthony R. G. (1989) Fundamentals of Chemical Reaction Engineering, 2nd edn. Prentice-Hall, Englewood Cli€s, NJ.

Solution for irreversible reactions Li L., Crain N. and Gloyna E. F. (1996) Kinetic lumping applied to wastewater treatment. Water Environ. Res. 68(5), 841±854. Resigno A. and Segre G. (1966) Drug and Tracer Kinetics. Blaisdell Publishing Co., Waltham, MA. Roberts A. L., Totten L. A., Arnold W. A., Burris D. R. and Campbell T. J. (1996) Reductive elimination of chlorinated ethylenes by zero-valent metals. Environ. Sci. Technol. 30(8), 2654±2659. Schnoor J. L. (1996) Environmental Modeling. Wiley, New York. Tondoi V., DiStefano T. D., Bowser P. A., Gossett J. M. and Zinder S. H. (1994) Reductive dehalogenation of chlorinated ethenes and halogenated ethanes by a highrate anaerobic enrichment culture. Environ. Sci. Technol. 28(5), 973±979. Walas S. M. (1991) Modeling with Di€erential Equations in Chemical Engineering. Butterworth-Heinemann, Boston. Wilson P. P. H. and Henderson D. L. (1997) Qualitative analysis of physical and mathematical approximations necessary for induced radioactivity calculations in fusion devices. Fusion Engr. Design 36(1±2), 415±427. APPENDIX A

Derivation Of Analytical Solution Consider three species (i, i', and j) operating in a large network of ®rst-order reactions:

…A1†

Let kj=kjj. By taking the inverse Laplace transform of equations of the same form, rearranging terms, then applying the Laplace transform again, it can be shown that this equation is equivalent to: ! m Ri mRi X X0 ki 0 j Dmi 0 kij Dmi fj …s† ˆ Cj …0† ÿ ÿ =…s ‡ kj † kj ÿ km kj ÿ km

‡

mRi  X

m6ˆk

 mRi X0  ki 0 j Dmi 0 =…kj ÿ km †  kij Dmi =…kj ÿ km † ‡ …s ‡ km † …s ‡ km † m6ˆk

X kij Dki ‡ ki 0 j Dki 0 ‡ : …s ‡ kj †2

…A2†

…A3†

mˆ1

…A7†

kRi

kRi 0

While the last Laplace transform step is not necessary for the presentation of the solution, the step from equation A6 to equation A7 demonstrates an important factoring formula in the Laplace domain. The last term avoids singularities between species j and species k, where kk=k kj. However, the general characteristic solutions for the parent species i and i' need to be modi®ed (equations A3 and A4) and the derivation will become considerably more complicated. At this stage, singularities will be avoided in the general solution only by numerical approximation (adjustment) of the lumped rate constant kk such that no singularity occurs. The solution can then be obtained if an inverse Laplace transform is applied:  mRi mRi X X0 ki 0 j Dmi 0  kij Dmi Cj …t† ˆ Cj …0† ÿ eÿkj t ÿ k ÿ km mˆ1 kj ÿ km mˆ1 j mRi X kij Dmi eÿkm t mˆ1

Let the concentrations Ci and Ci' follow the characteristic equation (equation 6, Walas, 1991),

kj ÿ km

‡

mRi X0 mˆ1

ki 0 j Dmi 0 eÿkm t : kj ÿ km

…A8†

Since the species i and i' have common ancestors m, the summations can be decomposed to obtain the general solution: m
Cj …t† ˆ Dij eÿkj t ‡ Dij eÿki t ‡ Di 0 j eÿki 0 t ‡

m
X

Dmj eÿkm t …A9†

mˆ1

and

where: Ci 0 ˆ

mRi X0

Dmi 0 eÿkm t :

…A4†

mˆ1

The parent species i and i' have common ancestors m. If these expressions are substituted into equation A2 and a Laplace transform is applied, the resulting expression is: X  mRi mRi X0 Dmi 0  Dmi sfj …s† ÿ Cj …0† ˆ kij ‡ ki 0 j ÿ kjj fj …s†: s ‡ km s ‡ km mˆ1 mˆ1 …A5† This reduces to:

fj …s† ˆ

mˆ1

m6ˆk

m6ˆk

The rate expression for species j would be:

mRi X Ci ˆ Dmi eÿkm t ,

mˆ1

mˆ1

mˆ1

‡ dCj ˆ kij Ci ‡ ki 0 j Ci 0 ÿ kjj Cj : dt

825

Cj …0† ‡ s ‡ kjj

kij

Dij ˆ

kij Dii kj ÿ ki

… parent compound i, daughter product j †, …A10†

Di 0 j ˆ

k i 0 j Di 0 i 0 kj ÿ ki 0

… parent compound i 0 , daughter product j †, …A11†

Dmj ˆ

kij Dmi ‡ ki 0 j Dmi 0 kj ÿ km

… for ancestors m
mRi X mˆ1

Dmi =…s ‡ km † s ‡ kjj

ki 0 j ‡

mRi X0 mˆ1

Dmi 0 =…s ‡ km † s ‡ kjj

mRi

: …A6†

mRi 0

Djj ˆ Cj …0† ÿ

X mˆ1

Dmj :

…A13†

826

Gerald R. Eykholt

The solution is implemented iteratively, as the D coecients for daughter products depend on the D coecients obtained for ancestor species m. The solution can be re®ned further, by considering more than two parents i. equation A12 can be modi®ed to: i
Dmj ˆ

kij Dmi

iˆm‡1

kj ÿ km

:

…A14†

If there is no indirect link between species m and j, Dmj=0, and kij and/or Dmi will be zero for each value of i in equation A15. In addition, since no direct link between species i and j will result if kij=0 (Dij=0), equations A10 and A14 can be added without any e€ect on the solution. Therefore, distinguishing between parent, ancestor, and unrelated species becomes unnecessary, and all terms of the coecient matrix (with ki $ kj) above the diagonal can

be described by one equation: m
Dij ˆ

mˆ1

kmj Dim

kj ÿ ki

… for i
…A15†

where i is a potential parent or ancestor species to j, and m is a potential parent to species j. The order of iteration for equations A13 and A15 may vary. However a convenient strategy is to initialize j to one, then calculate Djj followed by Dij for all values of i less than j. The D coecients can be calculated in the same manner for increasing values of j. To reduce errors in calculation, the coecient matrix and any unde®ned rate constants should be initialized to zero. Once the D coecients are determined, the general solution for any species j can be found for any time t, Cj …t† ˆ

iˆj X Dij eÿki t : iˆ1

…A16†