UNCERTAINTY ANALYSIS OF THE
EMISSION INVENTORY OF AN INDUSTRIAL
WASTEWATER TREATMENT PLANT
A CASE STUDY AT THE IWTP OF
TINKER AIR FORCE BASE, TINKER, OK
BY
KAISHAN ZHANG
Bachelor of Engineering
Hefei University of Technology
Hefei, P.R. China
1995
Submitted to the Faculty of the
Graduate College of the
Oklahoma State University
In partial fulfillment of
The requirements for
The Degree of
MASTER OF SCIENCE
August, 2002
UNCERTAINTY ANALYSIS OF THE
EMISSION INVENTORY OF AN INDUSTRIAL
WASTEWATER TREATMENT PLANT
A CASE STUDY AT THE IWTP OF
TINKER AIR FORCE BASE, TINKER, OK
Thesis Approved:
Thesis Advisor
II
UNCERTAINTY ANALYSIS OF THE
EMISSION INVENTORY OF AN INDUSTRIAL
WASTEWATER TREATMENT PLANT
A CASE STUDY AT THE IWTP OF
TINKER AIR FORCE BASE, TINKER, OK
Thesis Approved:
Thesis Advisor
II
Dedicated to my grandmother
Mrs. Suzhen Yang
with reverence
III
ACKNOWLEDGEMENT
I wish to express my sincere appreciation to Dr. John N. Veenstra, member of
my advisory committee for his intelligent supervision, constructive guidance,
inspiration, encouragement and friendship throughout the completion of my work.
My sincere appreciation extends to Dr. Dee Ann Sanders, my principal advisor,
for her encouragement, assistance and guidance during my study in the United
States. Thanks to Dr. Gregory Wilber for his invaluable assistance during my
research and serve as my advisory committee member.
Moreover, I wish to express sincere thanks to Dr. Williams F. McTernan,
who leads me to the world of uncertainty study, for his intelligent guidance,
warmly assistance during my study in uncertainty analysis and completion of my
research work. Many thanks also to my friend, Manivannan Nagaiah, who had
provided a lot of data useful for my research.
Finally, I would thank my parents, Mr. Yongchu Zhang and Mrs. Yali Li, my
wife Haiying Chen and many other friends and relatives in my home country for
their continuous encouragement, inspiration and support.
This thesis is dedicated to my grandmother, Mrs. Suzhen Yang, for her
devoting efforts in her life to our family.
IV
Chapter
TABLE OF CONENTS
Page
I. INTRODUCTION.............................................................. 1
1.1 Background............................................................ \
1.2 The emission uncertainty of the IWTP................... 3
1.3 General Fate Models (GFMs)....................................... 4
1.4 Significance of VOCs emissions analysis from IWTPs.......... 6
II. LITERATURE REVIEW.. 7
2.1 Determination of Emission Inventory.... .. .. . .. .. . . . . .. . .. .. . ... 7
2.2 Principles of GFMs............................................... 8
2.3 Uncertainty and Variability in an Emission Inventory.. 11
2.4 Methods of Uncertainty Analysis.. 12
2.5 Distributions in Uncertainty Analysis........ 19
Ill. METHODOLOGy.... 25
3.1 Emission estimating models................................... 25
3.2 Monte Carlo simulation........................................ 25
3.3 Bootstrapping............................................ 28
3.4 Maximum Likelihood Estimators (MLEs).................. 30
3.5 Analysis tools............................................... ..... 32
3.,6 Model setup for IWTP with TOXCHEM+V3 and WATER9... 33
3.7 General procedures................................... 36
[V. RESUTS.... 39
4.1 Sensitivity study ,. . 39
4.2 Models inputs setups..... 41
4.3 Fitting distributions.............. 42
4.4 Bootstrapping simulations........................................ 56
4.5 Backsolver study.................. 63
4.6 Simulation adequacy study........ .. . 65
4.7 Models simulation results..................................... 70
V. DISCUSSIONS................................... 81
5. I Distribution fitting..... 81
5.2 Interpretation of probabilistic results.......................... 83
v
Chapter Page
VI. CONCLUSION AND SUGGESTIONS . 92
BIBLIOGRAPHy........ 94
APPENDIX.. .. .. .. 98
APPENDIX A1 The un it process diagram of the IWTP
used in case study...................... 99
APPENDIX A2 The units dimensi.ons of the lWTP
in case study.......................... ... 100
APPENDIX A3 Raw flow rate data for distribution fitting
in case study. .. .. . .. .. . ... .. . . .. .. . .. . .... 104
APPENDIX A4Raw compounds data for distribution fitting
in case study................. ...... ...... 105
APPENDIX B1 Units parameters setup
forTOXCHEM+V3..................... 107
APPENDIX B2Units parameters setups
forWATER9............................. JOg
APPENDIX C An example of Bootstrapping simulation..... 117
VI
LIST OF TABLES
Table Page
80
80
88
91
16
41
45
46
48
48
50
51
53
53
55
56
62
63
63
64
Summary of method lISed to estimate emissions uncertainty .
Inputs of variables for WATER9 and TOXCHEM+V3 Simulation .
Summary of statistics for phenol. , ..
Comparisons of test methods for phenol. , ,.,.
Summary of fitting statistics of acetone , ,., .
Comparisons of test methods for acetone. , ..
Summary of fitting statistics for methylene chloride , , ,
Comparisons of test methods for methylene chloride , .
Summary of fitting statistics for 2butanone , ..
Comparisons of test methods for 2butanone , .
Summary of fitting statistics for flow rate ..
Comparisons of test methods for flow rate " ..
Comparisons of fitted and bootstrapped means .
Comparisons of fitted and bootstrapped standard deviation .
Summary of fitted distributions for compounds and flow rate .
Backsol ving studies of all the compounds ..
Summary of probabilistic and deterministic runnings
(TOXCHEM+V3 and WATER9 using their own Henry's constants)... 79
XVIII. Summary of probabilistic and deterministic runnings
(TOXCHEM+V3 and WATER9 using same Henry's constants)....... 79
XIX. Cumulative frequency of deterministic results on the simulated cumulative
distribution (TOXCHEM+V3 and WATER9 using their own Henry's
constant ,. . . . .. 80
XX. Cumulative frequency of deterministic results on the simulated cumulative
distribution (TOXCHEM+V3 and WATER9 using same Henry's
constants) .
XXI. Henry's constants used in both models (dimensi.onless) ..
XXTI. Comparisons of unit emissions in TXOCHEM+V3 and WATER9 ..
XXIII. Emission rate of compounds based on field measurement. .
I.
II.
III.
IV.
V.
VI.
Vll.
VIII.
IX.
X.
XI.
XII.
Xlll.
XIV.
XV.
XVI.
XVII.
VIl
Figures
LIST OF FIGURES
Page
1. Standard normal distribution................................................... 20
2. Standard exponential distribution........... .. .. 21
3. Pearson5 (3,3) distribution............... 21
4. Monte Carlo Sampling..... .. 27
5. Latin Hypercube Sampling..... 28
6. IWTP Configurations with WATER9. 34
7. IWTP Configuration with TOXCHEM+V3. . . .. .. . ... . .. .. . .. . .. . . .. . .. . .. 35
8. Diagram of simulating procedures............................................. 38
9. Sensitivity of wind speed on emission rate.................................... 39
10. Sensitivity of chemical concentmtion on emission rate.................. .... 40
11. Sensitivity of flow rate on emission rate....................................... 40
12. Fitted distribution of phenol................. . 43
13. PP plot of fitted distribution for phenol................................... 44
14. QQ plot of fitted distribution for phenol.. . 44
IS. Fitted distribution of acetone.................................................... 46
16. PP plot of fitted distribut~on for acetone...................................... 47
17. QQ plot of fitted distribution for acetone..................................... 47
18. Fitted distribution of methylene chloride................... 49
19. PP plot of fitted distribution for methylene chloride... . . . 49
20. QQ plot of fitted distribution for methylene chloride........................ 50
21. Fitted distribution of 2butanone........... 51,
22. PP plot of fitted distribution for 2bulanone.................... 52
23. QQ plot of fitted distribution for 2butanone................................ 52
24. Fitted distribution of flow rate.......................... 54
25. PP plot of fitted distribution for flow rate.............................. 54
26. QQ plot of fitted distribution for flow rate.............. .. . 55
27. Cumulative distribution of flow rate with bootstrapping.............. 57
28. Cumulative distribution of standard deviation of flow rate with
bootstrapping........................................................................ 57
29. Cumulative distribution of acetone with bootstrapping..................... 58
30. Cumulative distribution of standard deviation of acetone
with bootstrapping.... . . . 58
31. Cumulative distribution of methylene chloride with bootstrapping........ 59
VIII
Figures Page
32. Cumulative distribution of standard deviation of methylene chloride with
bootstrapping................ .. 60
33. Cumulative distribution of 2butanone with bootstrapping................. 60
34. Cumulative distribution of standard deviation of 2butanonewith
bootstrapping..... .. 61
35. Cumulative distribution of phenol rate with bootstrapping................. 61
36. Cumulative distribution of standard deviation of phenol with
bootstrapping................ . 62
37. Linearity of backsolver for chemicals of concern....... . 64
38. Adequacy of simulation using TOXCHEM+V3 for benzyl alcohol, phenol,
acetone and 2butanone.. 66
39. Adequacy of simulations using TOXCHEM+V3 for
methylene chloride.. .. 66
40. Adequacy of simulations using WATER9
for benzyl alcohol... 67
4 J. Adequacy of simulations using WATER9 for phenol..... .. 67
42. Adequacy of simulations using WATER9 for methylene chloride, 2butanone
and acetone 68
43. Adequacy of simulations using TOXCHEM+V3 with standard deviation... 69
44. Adequacy of simulations using WATER9 with standard deviation.... ....... 69
45. Adequacy of simulation using WATER9 with standard deviation for phenol ...70
46. Simulation results of benzyl alcohol by using TOXCHEM+V3............ . 71
47. Simulation results of phenol by using TOXCHEM+V3.. 7\
48. Simulation results of methylene chloride by using TOXCHEM+V3....... 72
49. Simulation results of 2butanone by using TOXCHEM+V3.......... 72
50. Simulation results of acetone by using TOXCHEM+V3........ ..... 73
51. Simulation results of benzyl alcohol by using WATER9.................... 7J
52. Simulation results of phenol by using WATER9.............................. 74
53. Simulation results of methylene chloride by using WATER9............... 74
54. Simulation results of 2butanone by using WATER9........................ 75
55. Simulation results of acetone by using WATER9............................ 75
56. Simulation results of phenol with WATER9 Henry's constant using
TOXCHEM+V3.......... 76
57. Simulation results of methylene chloride with WATER9 Henry's constant using
TOXCHEM+V3.................................................................... 76
58. Simulation results of 2butanone with WATER9 Henry's constant using
TOXCHEM+V3 _. 77
59. Simulation results of benzyl alcohol with WATER9 Henry's constant using
TOXCHEM+V3. 78
60. Simulation results of acetone with WATER9 Henry's constant usmg
TOXCHEM+V3..... 79
61. Prediction of VOCs emissions with GFMs...... .. .. .. 89
IX
R~res ~~
62. Estimated VOCs emissions by process...................................... .... 90
x
Chapter I
Introduction
1.1 Background
The School of Civil and Environmental Engineering at Oklahoma State University is
engaged in a project to improve and validate tbe accuracy, reliability, and repeatability of
target pollutant emissions estimates through monitoring, process unit sampling, and
computer modeling of the OCALC (Oklahoma City Air Logistics Center) / IWTP
(Industrial Wastewater Treatment Plant) air emission sources for Tinker Air Force Base.
The sources of air emissions include: the primary paint chip clarifier, oilwater
separators, aerated equalization basins, storage / stabilization tanks, metals treatment
basins, solid contract clarifiers, lift stations, and gravity thickener. In addition, this
project will include development of an air emission sampling strategy to improve the
accuracy of the air emissions reporting data (Veenstra, et aI., 200 l).
In production and maintenance processes at the OCALe, industrial wastewater
streams are generated which contain organic and inorganic compounds. Most of the
wastewater is generated from electroplating, chemical cleaning, and chemical depainting
operations (Hall, 1999). Before being discharged from the treatment facihty to the city of
Oklahoma City, these wastewaters are conected and treated in a variety of ways.
However, since many of these collection and treatment systems are open to the
atmosphere, they allow organiccontaining, heavy metalladen wastewaters to contact
ambient air. Based on logs of chemi.cal consumption and field monitoring, benzyl
 1 
alcohol, phenol, acetone, 2butanone, and methylene chloride are the main pollutants
emitted via the air phase from this IWTP (Hall, 1999).
The Clean Air Act Amendments of 1990 requires wastewater treatment facilities to
identify sources and quantify emissions of volatile organic compounds and hazardous air
pollutants. Regulation under 40 CFR Paft 63, National Emission Standards for Hazardous
Air Pollutants: Publfucly Owned Treatment Works, require operations to quantify and
report chemical releases to the environment (CFR Title 40: Protection of Environment,
US EPA, 1995).
Since residential neighborhoods are located in close proximity to the IWTP,
residents are potentially exposed to the chemicals released from the IWTP and could be
at risk. To evaluate effects of these chemicals on the environment and residents,
measurements of the air contaminants are required.
The overall research project involves three phases. Phase 1 involves the acquisition
of facility data and application of various air emission models. Most work for this phase
bas been accomplished. Two air emission models have been chosen, i.e., WATER9
(USEPA) and TOXCHEM+V3 (Enviromega Inc, Canada). Phase 2 includes the analysis
of model output and comparison to field and pilot plant data, while Phase 3 will involve
establishing the emission factors from each of the individual industrial wastewater
treatment plant process units. The project now is in Phase 2.
 2
Based on the summary of Phase 1, it was found that variability of components in the
influent is significant (Veenstra, et a1, 2001). For example, benzyl alcohol is not actually
measured but estimated, based on the measured concentration of phenol since the
consumption of benzyl alcohol is proportional to that of phenol in the chemical
consumption logs. It is estimated to have a mean concentration of 15.4 mgll if the
diffusivity of benzyl alcohol and phenol in the liquid phase are considered, and an
estimated mean concentration of 52.4 mgll if the diffusivity of both chemicals are not
considered (Veenstra, et aL, 2001). The concentrations of the contaminants of concern
vary seasonaUy. Methlyene chloride averages range from 18.8 mgll in winter to 22.75
mgll in summer. 2butanone averages 7.97 mg/l in the winter and 8.35 mg/l in the
summer, while acetone averages range from 14.75 mg/I to 18.05 rngll in summer and
winter, respectively. Phenol is estimated to have a mean concentration of 14.3 mg/l over
19992001, in both summer and winter (Veenstra, et al., 2001).
1.2 The emission uncertainty of the IWTP
WATER9 and TOXCHEM+V3 were used to estimate air emissions from the IWTP
based on the concentrations of the five compounds mentioned above. But it is not known
which value is more representative and should be selected for use in determining the
emission factors. Thus, uncertainty inevitably exists. Generally, uncertainty has three
sources: modeling uncertainty, parameter uncertainty and input uncertainty (Monte,
1996). The more that is known about the model, parameter, or input value, the less
uncertainty remains. Therefore, exactly modeling the actual units for the treatment
facilities is important, and it helps reduce the uncertainty in the estimation of air
 3 
emissions from this system. Sensitivity studies can help determine which factor the
system is the most sensitive to and potentially has a high contribution to the uncertainty.
In this thesis, the uncertainty contributed by the variability of concentrations of
compounds in water is the primary focus and will be discussed in detail.
1.3 General Fate Models (GFMs)
Field measurements and general fate models (GFMs) are the two main means to
estimate the emissions from IWTPs (Curto and Daly, 1995). GFMs are computational
models that perform a mass balance around each specified wastewater unit operation and
certain solids handling facilities, as well as the whole wastewater treatment facility. The
mass balances are usually performed considering five mechanisms in the treatment
system: 1) volatilization across the exposed wastewateratmosphere interface, 2) stripping
to diffused air bubbles, 3) adsorption to solid particles or biomass, 4) absorption to
immiscible liquids, and 5) biodegradation (Corsi and Olson, 1998). The development of
GFMs for air emissions was prompted by the complexity and high cost of direct air
samphng and measurement from wastewater treatment facilities. Using GFMs, it is
possible to estimate emissions from complex treatment configurations while considering
split flows, liquid streams, quiescent surfaces, weirs, drops, as well as aerated, biological,
and covered processes or any single operation or process (Curto and Daly, 1995). Influent
wastewater characteristics, physical design characteristics, and operational data are
required to set up and run one of these models. Furthermore, GFMs are particularly
advantageous when projecting potential emissions under varying flow and design
conditions (Curto and Daly, 1995). The comparisons between field measurements and
 4
1.4 Significance of VOCs emissions analysis from IWTPs
Volatile Organic Compounds (VOCs), such as benzyl alcohol, phenol,
methylene chloride, 2butanone and acetone, emitted from Tinker's lWTPs, can cause
serious environmental problems due to their odor and toxicities. Some of them are
harmful and carcinogenic if they are directly contacted or inhalated. For example,
methylene chloride, is carcinogenic and might cause hepatocellular adenomas or
carcinomas, hepatocellular cancer and neoplastic nodules with oral exposure ORIS
data base, USEPA, 1986). Potential health problems created by these VOCs would
affect not only the workers of the treatment plants, but aJso the general public and
residents in the surrounding area. Before measures are taken to protect the
environment from being polluted and people from being at risk, the VOCs emission
inventory of the IWTPs and the emission rate should be determined. The emission
inventory and emission rate determine what protective actions should be taken.
However, the variation and uncertainty of VOCs emissions make it tougher for
decision makers, especially in environmental risk assessment. The analysis of
variation and uncertainty of VOCs emissions are very important and necessary for an
efficient action of protection; however, uncertainty analysis of VOCs emissions from
IWTPs was not discussed in literature.
This thesis presented a method for uncertainty analysis of VOCs emissions from
IWTPs with a case study.
 6 
Chapter II
Literature Review
2.1 Determination of Emission Inventory
The determination of an emission inventory for air phase pollutantreleasing
facilities is essential for overall environmental management, environmental risk
assessment, and environmental impact assessment. It is required and regulated by the
USEPA (United States Environmental Protection Agency) (Curto and Daly, 1995).
Methods that are available to estimate air emissions include stack/field testing, published
emission factors, engineering equations, and a new and innovative tool/methodgeneral
fate modeling (Curto and Daly, 1995). Some authors have also llsed the inputoutput
table method to predict the emissions from facilities (Jin, 1986; Ni, et aI, 2001). USEPA
AP42 (Compilation of Air Pollutant Emission Factors) is widely used in the prediction
of ail" pollutants from mobile sources, coalfired power plants and other energy facilities,
for which published emissi.on factors for the regulated facilities are available. For the
prediction of emissions of Volatile Organic Compounds (VOCs) from wastewater
treatment facilities, the general fate model is a new and innovative method to estimate the
emissions.
GFMs are computational models that perform a mass balance around each specified
wastewater unit operation and certain solids handling facilities, as well as the entire
wastewater treatment facility. Available GFMs include BASTE (Bay Area Sewage
Toxics Emissions) by R.L. Corsi of the University of Texas at Austin, WATER8
(USEPA), CINe} CUSEPA Cincinnati model) by Richard Dobbs at the University of
 7 
Cincinnati, CORAL (Collection System Organic Release Algorithm) by R.L. Corsi of the
University of Texas at Austin while at the University of CaliforniaDavis, EPA FATE
(Fate and Treatability Estimator) by ABB Environmental, NOCEPM (NCASI Organic
Compound Elimination Pathway Model) by D.A. Barton with the National Council of the
Paper Industry for Air and Stream Improvement (NCASI), PAVE (Programs to Assess
Volatile Emissions) by the U.S. Chemical Manufacture Association, SIMS (Surface
Impoundment Modeling Systems1990) by RADIAN for USEPA Office of Air Quality
Planning and Standards, TORONTO (specially for biological wastewater treatment
facility) by B. Clark and n. Mackay of the Institute of Environmental Studies, University
of TORONTO, and TOXCHEM+V3 (Enviromega, Canada). Among these models,
WATER8 and TOXCHEM+V3 were the only two models with temperature correction
for Henry's Constant, which is important for volatiles in a wastewater plant (Hall, 1999).
Furthermore, WATER8, BASTE, and TOXCHEM were selected by USEPA as
appropriate models for wastewater treatment systems (Card, 1995). WATER9 is the
newly version of WATER8 developed by USEPA for estimating air emissions from
water and wastewater sources.
2.2 Principles of GFMs
Corsi and Olson (1998) discussed the fundamental principles of these emission
estimating models or VOC fate models. Under steadystate condition, the mass balance in
a treatment system can be expressed as follows:
de
V =QCo  QC + Rv + R~ + Rad + Rab + Rb
d/
 8 
(J)
where:
V: reactor volume (m\
C, Co: dissolved contaminant concentration leaving and entering the reactor,
respectively (mg/m3
),
Q : volumetric flow rate entering into the reactor (m'/s).
Rv, Rs.Rad,Rab•Rb: contaminant removal rate by volatilization, stripping, adsorption,
absorption, and biodegradation, respectively (mg/s).
Since absorption is generally a complex process that is not well understood,
absorption is not considered for applications involving municipal wastewater, but may be
important in industrial wastewater systems (Corsi and Olson, 1998). The rate of
volatilization, stripping, adsorption and biodegradation can be modeled by the following
equations:
The rate of volatilization is typically modeled as:
where:
KL is an overall mass transfer coefficient (m/s),
A is the interfacial area over which mass transfer occurs (m\
(2)
CO" is the contaminant concentration in the gas phase adjacent Lo the wellmixed liquid
b
He is the Henry's law constant for chemical interested (mli//ffigas\
The rate of stripping by air bubbles can be modeled as:
where:
 9 
(3)
Qb is the bubble volumetric flow rate (m3/s),
Yis a variable which represents the degree of saturation.
The adsorption rate can be modeled as:
where:
Qw is the volumetric sludge wastage rate (m3/s),
Cs is the solids or biomass concentration (mg/m3
),
Kp is a linear liquidsolid partition coefficient (m3/mg).
(4)
The rate of biodegradation is assumed to follow Monod kinetics and can be modeled as:
where:
kb is the first order biodegradation rate constant (rn3/mg.s),
X is the active biomass concentration (mg/m\
Combining all the above from equations 1 through 5, for a steadystate condition and
open process (Cg = 0), yields:
1
= 
1+ K1.A/ Q + (Qh / Q)yf(. + (Q". / Q)KpC, + kbX / Q
(5)
(6)
Equation 6 serves to illustrate differences in existing models for estimating VOC
emissions from wastewater. For example, some models account for all the terms in the
denominator, while others do not. In addition, models differ in how they estimate or
prescribe parameters such as KL, He, y. Kp, and kb (Corsi and Olson, 1998).
. 10
2.3 Uncertainty and Variability in an Emission Inventory
Most emissions inventories are not obtained by field acquisition of data, but
predicted based on mass balance and computational models (McKay, et a1., 1998). The
GFMs are based on computational equations to simulate the actual scenarios occurring in
a wastewater plant, and as such, variability and uncertainty will inevitably exist.
Emission inventory uncertainty and variability analysis is not a new concept and many
studies on this topic have been conducted over several decades (McKay, et aI., 1998).
The U.S. Environmental Protection Agency (EPA) and the State and Territorial Air
Pollution Control Officers' Association and Association of Local Air Pollution Control
Officers (STAPPNALAPCO), published guidelines for evaluating the uncertainty of
emission estimates as part of the Emission Inventory Improvement Program in 1997 (Roe
and Reisman, 1998). Uncertainty and variability are confusing concepts. A dear
understanding of these two terms is required before using these two concepts to express
the analysis results when there is uncertainty or variability. or both, in an emission
inventory. USEPA (1997) has laid out a detailed definition and deseri ption of the two
concepts in "Guiding Principles for Monte Carlo Analysis" (USEPA, 1997).
"Uncertainty refers to lack of knowledge about specific factors, parameters,
or models. For example, we may be uncertain about the mean concentration of a
specific pollutant at a contaminated site or we may be uncertain about a speci fic
measurement of uptake. Uncertainty includes parameter uncertainty
(measurement errors, sampling errors, systemic errors), model uncertainty
(uncertainty due to simplification of realworld processes, misspecification of the
model structure, model misllse, use of inappropriate surrogate variables), and
 II 
scenario uncertainty (descriptive errors, aggregation errors, errors in professional
judgment, incomplete analysis), while variability refers to observed differences
attributed to true heterogeneity or diversity in a population or exposure parameter.
Sources of variability are the result of random processes and stem from
environmental, lifestyle, and generic differences among humans. Examples
include human physiological variation (e.g., natural variation in bodyweight,
height, breathing rates, and drinking water intake rates), weather variability, and
variation in soil types and differences in contaminant concentrations in the
environment. Variability is usually not reducible by further measurement or study
(but can be better characterized) (USEPA, 1997)."
2.4 Methods of Uncertainty Ana ysis
Many authors also provided commentary and examples to help the nonexpert reader
have a better understanding of these two concepts using actual applications (Frey, et aI.,
1995; 1998; 1999). Uncertainty refers to lack of knowledge and mea 'uremenl error while
variability refers to temporal variation. Usually, the uncertainty of the output of an
estimating model is the mixture of uncertainty and variability.
Simulation variability, input uncertainty, and structure uncertainty are three sources
of uncertainty or variability in the prediction of simulation models (McKay et aI., 1999).
Many authors have tried different ways to study uncertainty of models. Hanson (1999)
presented a framework for assessing uncertainties in simulation predictions. The methods
used in the framework included individual experiment analysis, many experiment
 12 
anajysis, mathematical approximation (i.e. Gaussian approximation), and Markov Chain
Monte Carlo simulations. Monte, et aJ. (1996) used:
where:
d2 is the error term,
Mj is experimental value,
OJ is the predicted value, and
n is the number of couples of experimental and predicted value,
(7)
as a measure to perform uncertainty studies and validate environmental models. This
method is also called EBUA (empirically based uncertainty analysis). McKay et.al (1999)
used:
11 2 =V(y)/V(y)
where:
V(y) is the restricted predicted variance,
V(y ) is the predicted vadance, and
112 is the Pearson relation coefficient.
(8)
.
(
as a measure to address the uncertainty of predictive models. Wallach and Genard (998)
used mean squared error of prediction (MSEP) with Taylor series expansion La study the
effect of uncertainty in input and parameter values on model prediction error. The authors
used the following equation:
(
where:
MSEP =E{ {yO f(u,p,q(p»)}2}
 13 
(9)
y. is the value of the output of interest for an individual chosen at random from the
population of interest.
f(u,p,q(p») is tbe corresponding model prediction.
In this paper, Taylor series expansion was used to adjust the effect of parameters. The
error propagation was estimated based on the error term of the Taylor series expansion.
The Taylor series expansion is shown in the following equation (Chapra and Canale,
2002):
(10)
where:
rn+l)(r) represents the (n+ 1)lh derivative of the function f(x) when x=r,
h = xr, r is in very close proximity to x.
Their work showed that the uncertainty in model inputs always increases the model
variance contribution to MSEP. However, one of the properties of these above methods is
that while adequate data are required for the uncertainty analysis, sufficient data are
usually not present. In addition, these equations (7 through 10) assume that the variables
or parameters have normal distribution. This decreases accuracy if the variable does not
have a normal distribution. The importance of the input distribution has been proven to
be very critical in uncertainty analysis by Frey et a1. (1995; 1998; 1999). Roe and
Reisman (1998) recently summarized the methods of uncertainty analysis. In their
summary, all of the methods typically can be categorized as qualitative, semiquantitative
and quantitative. For the qualitative method, sources of uncertainty are listed and
discussed, emission factors are listed and subjectively ranked, and then the uncertainty is
estimated in this way. For the semiquantitative method, some statistical properties of the
• 14 
"
I
1
'I
1
'I
(,
(
data are studied, such as tbe probability distribution, mean, and standard deviation, and
then the error propagation is estimated with standard statistical techniques, i.e. Taylor's
series expansion. For the quantitative method, simulation techniques, i.e" Monte Carlo,
Latin Hypercube, or Bootstrapping. are used to estimate the confidence intervals of
factors of interest. Also, modeling methods based on adequate field measurements are
used to achjeve high accuracy. A summary table is produced in Table 1 (Roe and
Reisman, 1998).
In a general sense, all of the methods used for uncertainty analysis can fit one of the
categories discussed above. But to some extent, it is not easy to follow. Hence, if the
analysis methods were based on the estimation of emissions from operational units, the
uncertainty and variability analysis would be more straightforward. The uncertainly
analysis would be conducted based on the variation of the estimation of emissions. The
analysis methods then include direct emission analysis, emission factors analysis with
simple equations, inputoutput table analysis, and GFMs (general fate models)
uncertainty analysis. Of these methods, the most accurate one is direct emission analysis.
However, it requires adequate field emission elata, which most facilities do not usually
have. Though it is technically feasible, it is not always costeffective.
Some authors have used the inputoutput table method to estimate the emissions, i.e.
CO2, natural gas and other air po1]utants, from power plants, manufacture factories and
wastewater plants (Jin, 1986; Caloghirou, et al., ]996; Ni, et al., 2001). The idea behind
this method is application of the principles of mass balance and energy movement.
However, based on this method, the uncertainty analysis can be done only when the
emissions have a linear and simplified relationship with the input and output. Under such
 15
(
Table 1. Summary of methods used to estimate emissions uncertainty
Methods Description Relative Level of
Effott
Qualitative Methods
Qualitative Source of uncertainty are listed and discussed;
Discussion general direction of bias and relative magnitude of Low
imprecision are given, if known.
Subjective Data Subjective rankings based on professional Low
Quality Rating judgment are assigned to each emission factor or
Darameter.
SemiQuantitative
Methods
Data Attribute Numerical values representing relative uncertainty Moderate
Rating System are assigned through objective methods.
Emission distribution parameters (e.g., mean,
standard deviation, and distribution type) are
Expert Estimation estimated by experts. Simple analytical and
Method graphical techniques can then be used to estimate Moderate
confidence limits from the assumed distributional
data.
Emission parameter means and standard
deviations are estimated using expert judgment,
Propagation of measurements, or other methods. Standard
Errors Method statistical techniques of error propagation Moderate
typically based on Taylor's series expansions are
then used to estimate the composite uncertainty.
Quantitative
Methods
Monte Carlo, latin hypercube, bootstrap, and other
Direct Simulation numerical methods are used to estimate directly
Method the central value and confidence intervals of Moderate to High
individual emission estimates.
Direct or Indirect Direct or indirect field measurements of emissions
Measurement are used to compute emissions and uncertainty High
(Validation Method) directly.
Receptor Modeling Can be used to provide a measure of the relative
Method contribution of each source type but noL absolute High
emission eSlimates.
Air quality simulation models are used In an
[nverse Air Quality inverse, iterative approach to estimate the High
Modeling Method emissions required to produce observed
,, concentrations.
(Roe and Reisman, 1998)
 16
"
~,:
a scenario, the uncertainty of the emission can be expressed with variance. For example,
the propagation of errors of input and outputvariables contributes to the uncertainty and
variability of the emissions.
Suppose:
(11)
where F(x) stands for the emission and is a function of input and output.
It can easily be expressed as:
Emission =Input  Output
then,
where:
Var(x) represents the variance of x,
Cov(Xj,Xj) represents the covariance of Xi and Xj.
If the variance of input and output are independent, then
Var(Emi) =Var (Input) + VOIr (Output)
(12)
(13)
(14)
....I
I
'1
I
.\
(
If the appropriate measurement of the statistics of the population is chosen, the
uncertainty can then be easily estimated by statistical calculation.
The theoretical treatment of inference is handled in most introductory statistics
books such as Ross' (1998) book "A First Course in Probability (5th ed.)". Models are
usually simplified when using emission factors and engineering equations to perform
uncertainty analysis of an emission inventory, and the models are called engineering
equations and comprised of most subjective factors. Roe (1998) provided an example of
this method, while he used a Monte Carlo simulation to estimate the emission inventory
 17
uncertainty. As addressed above, GFMs are widely used for estimation of emission
inventories. For its uncertainty, simulation techniques, i.e. Monte Carlo simulation, Lati.n
Hypercube simulation and Bootstrapping simulations, etc., are widely used. Frey (1998)
presented a paper on uncertainty analysis of air pollutants using the Bootstrapping
simulation method, where he also used a twodimensional approach to probabilistic
simulation. When using simulation techniques to estimate the uncertainty, knowledge of
the random variables' distributions is very critical. Wrong assumptions concerning the
distribution of variables, especially the input distribution, will result in unrealistic
outcomes. For the purpose of filling in the gaps due to lacking knowledge of input
distributions, Frey and Cullen (1995) published a methodological handbook as practical
guidance for uncertainty analysis. Since most of the methods mentioned above are simply
based on a single distribution to represent variability and uncertainty in the input of the
model, Zheng and Frey (2001) presented another method. They applied mixed
distributions to represent the variability in the input of the model. Frey, et at. (1999)
summarized the general steps of unceltainty analysis:
1. Assemble and evaluate a database,
2. Visualize data by developing empirical cumulative distribution functions for
individual variables and preparing scatter plots to evaluate dependencies among pairs
of variables,
3. Select, fit, and critique alternative parametric probability distribution models for
representing variability in activity and emissions factors,
4. Characterize uncertainty in the distributions for variability,
 I R
,
,
1
I
I"
I.,
5. Evaluate tbe effect of averaging, over both time and space, on variability and
uncertainty, and
6. Propagate uncertainty and variability in activity and emissions factors to estimate
uncertainty in emissions.
2.5 Distributions in Uncertainty Analysis
Usually, inputs for models are randomly sampled from specific distributions during
a simulation. Normal distribution is the most prevalent, the best known and most widely
used distribution in the world (Hahn and Shapiro, 1967). Normal distributions are
symmetric with scores more concentrated in the middle than in the tails. They are known
to have a bell shape. They are defined by two parameters: the mean (m) and the standard
deviation (s). Many kinds of behavioral data are approximated well by the normal
distributi.on. Many statistical tests assume a normal distribution. Most of these tests work
weB even if the distribution is only approximately normal and in many cases as long as it
does not deviate greatly from normality. The density function of a normal distribution can
be expressed as:
).J
.)
:I
.. \..~
. .:J
;~
.~i3
I"~ ..J)
<l:I Eo
.J.::
"'~
.~
I'
f(x) = .~ 2ncr 2
where:
Il =the mean of samples
(15)
cr = the standard deviation of samples
 19 

When 11=0, (J' = I, the normal distribution is called a standard normal distribution as
shown in Figure 1.

Normal (0, I)
X <= 1.6448
5.0%
X <= 1.6448
95.0%

0.4 .
0.3
0.2
0.1
0
3 2 I 0 2 3
Figure 1. Standard normal distribution
When the scores are not concentrated in the middle, but skewed in the left tail, the
distribution will not be symmetric or bell shaped. The distribution will not be a normal
distribution, but might be an exponential distribution, Pearson distribution. elC. The
following figures are some examples of these skewed distributions.
 20
.\
).j
)
J
X <=0.051
5.0%
0.8 
0.6 .
0.4 .
0.2 
.'
1.5 2 2.5 3 3.5 4 4.5 5
o +JL,r.::.....,~~_.=~~~~
0.5 0 0.5
Figure 2. Standard exponential distribution
PearsonS (3, 3)
)..
s..
x <= 0.476
5.0%
x <= 3.669
95.0%
2 } 4 567 8
0.9 ,.,..r,
0.8
0.7
0.6 0.5
0.4
0.3
0.2 0.1
o \tJ.~,__r_,....::::J~~~_t_'__t_1
\ 0
Hgure 3. PearsonS (3,3) distribution
The general t~onnula for the probability density function of the exponential distribution is
.
....~
f(x) = ~.e(\J.lJI /3, X ~ ~; ~ > 0
 2\ .
(16)
where Il is the location parameter and ~ is the scale parameter (the scale parameter is
often referred to as A. which equals l/~ when the scale parameter is expressed as A). The
case where 11 =0 and ~ =1 is called the standard exponentiat distribution. The equation
for the standard exponential distribution is
f(x) =eX, for x ~ 0 (17)
Pearson5 distribution is the transition form of the Pearson distributions (Jeffreys, L961).
1..\
Each family in the Pearson system can be generated as a soLution to the differential
equation
df (x) (x  <1>',)1 (x)
=
dx <1>0 +<1>IX+<1'>2X2
(18)
.'. .;
For the random variable x, a probability density function f(x) is determined by proper
choice of the four parameters <1>0. <1>1. <1'>2, and 4>3. Different choices of the four parameters
win lead to different distributions. The solution of this equation leads to a large number
of distribution families, including normal, Pearson I, and Pearson3 (Pearson distribution
families) distributions (Hahn and Shapiro, 1967). Equation 18 can be transformed to
..
. ~
(
" ..J
(19)
When the roots of the denominator in Equation 19 are equal, real and finite, then
Equation 19 can be written in the form
1df(x) =a+ 'f3::
f(x) dx xc (XC)2
 22
(20)

whence
where
f(x) = A(xct'exp[~/(xc)] (21)
A will be fixed by the condition that the integral of f(x) is 1,
C is the zero of the denominator in equation (18)
( \ a, pare shape parameter and scale parameter, respectively.
For detailed discussions of these distributions the reader is referred to Jeffreys (1961) and
Hahn and Shapiro (1967).
Once the distribution is determined, random numbers are sampled from the specific
distribution. Models are then run hundreds of times based on the required accuracy of the
results. But there is a prerequisite to achieving this goal. The model must be easily
mathematically expressed as an equation or the applicable models have the integrated
simulation function. All of the methods for uncertainty analysis mentioned above have
the same prerequisite, which is the models can be easily entered into Excel or other
spreadsheet, and then used with some statistics software, i.e., @RISK (Palisade, Inc.),
Crystal Ball (Decisioneering, Inc.), etc. However, estimating software for emission
inventories from wastewater treatment facilities such as WATER9, TOXCHEM+V3
focus mostly on air emissions based on mass balances, so that there is no such uncertainty
function built into the software. Under such a scenario, there are two alternatives. One is
coding the uncertainty function into the software, which is a complicated task. The other
is to assume the model is certain based on the expert advice and validation experience,
·23 
).S
..
)
j.
)
J

hence assuming the uncertainty and variability of the input parameters contribute all to
the uncertainty of the emissions. Since WATER9 (a newly version of WATER8), and
TOXCHEM were selected by USEPA as appropriate models for estimating air emissions
from water and wastewater resources (Card, 1995), TOXCHEM+V3 (a newly version of
TOXCHEM) and WATER9 were selected for this thesis and assumed to be certain for
the process of uncertainty analysis. The operation units were exactly modeled as they are
surposed to be in these two programs.
<
The uncertainty analysis of emissions from wastewater treatment plants has not been
adequately studied and discussed in literature. Instead, the uncertainty analysis of air
emissions was seen to date mostly to focus on coalfired power plants. Since GFMs were
developed as a deterministic tool to estimate the emissions from wastewater treatment
systems, it makes the uncertainty analysis of emissions from wastewater plants, or the
probabilistic analysis of emissions from wastewater plants, difficult. This thesis strives to
present a method to conduct uncertainty analysis of air emissions from wastewater
treatment plants.
 24 
s.,
)
j
....
Chapter III
Methodology
3.1 Emission estimating models
Emission estimating models mainly refer to general fate models (GFMs). Models used
to estimate emissions of volatile organic compounds from wastewater plant include
<';OXCHEM+, WATER8, BASTE, CORAL (Corsi and Olson, 1998). Since
TOXCHEM+ and WATER8 are considered as appropriate models for estimating air
emissions from water and wastewater resources (Card, 1995), and TOXCHEM+V3 and
WATER9 are the newly versions and similar in estimation of emission from process
drains (Corsi and Olson, 1998), WATER9 and TOXCHEM+V3 were chosen to be used
in this thesis.
3.2 Monte Carlo Simulation
Monte Carlo techniques have been used since the J940's, when they were first
developed by physicists working on the Manhattan project (Hammersley and
Handscomb, 1964). Only recently, however, have personal computers become
sufficiently powerful and widespread for Monte Carlo techniques to be widely applied for
health risk assessments.
Modem spreadsheet programs, such as @RISK (Palisade, Inc.), Crystall Ball
(Decisioneering, Inc.), and Xlsim (AnalyCorp, Inc) now provide a range of critical
factors to illustrate and order a model, including advanced statistical functions, charting,
 25 
.'
).~
)
~
etc. The ongm of the name "Monte Carlo" relates to the famous gambling city in
Monaco, but the relation to gambling applies only to the probability of a given event
occurring over the long term. Although one cannot know precisely which number will
appear on the next roll of a craps die or the spin of a roulette wheel, one can predict over
the long term (and as precisely as desired) the frequencies associated with each outcome
(Vose, 1996). Monte Carlo simulation techniques similarly cannot predict exactly which
exposures will occur on any given condition to any specific individual, but can predict
the range of potential exposures in a large population and each exposure's associated
probability. Monte Carlo simulation is conducted by randomly sampling from input
probability distributions for sufficient times to produce an output distribution, which
reflects the expected range and frequency of exposure. Once the model and distribution
are determined, the sampling method becomes the most critical issue because it will
affect the efficiency and accuracy of the simulations. Monte Carlo sampling and Latin
Hypercube sampling are two techniques widely used for Monte Carlo simulations.
3.2.1 Monte Carlo Sampling
Monte Carlo sampling is entirely random (with replacement). Sampling with
replacement means that a value might be sampled twice or more since sampling is totally
random. Over the whole range of the input distribution, samples are drawn randoml y
within the range. Therefore, some areas would have higher probabilities of occurrence,
which inevitably causes uneven sampling. Only after sufficient sampling can it be
assured that the sampled input stands for the input distribution. Figure 4, shows this
sampling method (Palisade Corporation, 2001).
 26
.)
1
:...:;
.....
.... '
....
Five Iterations of Monte Carlo Sampling With Clustering
MaXilrum
DJstributlon
VlIlJUO
Values
sampled
Random
number
generate
r
,~5               
.51               _
.52
.46
....1 _
.2
.3
o
Mininurn
Distribution
Value
1.0
.9 t B
.7
CumUl!ltlve
.6
Probability
.5
Figure 4. Monte Carlo Sampling. (Palisade Corp., 2001)
3.2.2 Latin Hypercube Sampling
Compared to Monte Carlo sampling, Latin Hypercube sampling is pseudo random
..'k
(without replacement). Sampling without replacement means that a value might not be
sampled twice. Depending on the numbers of samples, the whole range of the " ~
distribution is divided into even intervals. Then the values are sampled from these
intervals without replacement.. Obviously, it reduces time for sampling to achieve a
,
simulated distribution that stands for the input distribution. Figure 2, Latin Hypercube "
sampling, shows this sampling method (Palisade Corp., 2001).
 27
Five Iterations of Latin Hypercube Sampling
Figure 5. Latin Hypercube Sampling. (Palisade Corp., 2001)
l,
I,
.1
)
.'j.
Meximum
Dislribution
Value
IIIIIII
I I I :
Values Swttplea
 IP;
.1
o
MinimuM
Distrl bUtton
Value
1.0
.9
.B
.7
Cumulative
Probability .6
Distribu tion .5
Slcalif,ed
Into .4
Five
.3
Intervals
.2
Since Monte Carlo sampling is entirely random, unless sampling is sufficient, it is not
even over the distribution. Since Latin Hypercube sampling is pseudorandom, samples
"
are drawn equivalently without replacement. Hence, even a few samples can represent the
distribution. Therefore, Latin Hypercube sampling is more effective than Monte Carlo
sampling (Palisade, 2001). The common nature of these two sampling techniques is that
sampling is based on the known or assumed distributiono If the distribution is not known,
acceptable assumptions are then required.
OJ
.~
3.3 Bootstrapping
The Bootstrapping technique was introduced by Efron in 1979 for the purpose of
estimating confidence intervals for a statistic using numerical methods (Frey and
Burmaster, 1999). The key advantage of this technique is that it can provide estimates of
confidence intervals in situations for which analytical mathematical solutions may not
exist. Therefore, it is used for confirmation of a distribution by fitting the assumed
 28 
statistics (i.e. mean, standard deviation) to the confidence intervals that the bootstrapping
simulation produces. It is very helpful when there are only limited data and a specific
distribution is assumed for a variable. As defined by Efron and Tibshirani (1993) and
summarized by Frey and Burmaster (1999), bootstrap simulationts based upon drawing
multiple random samples, each of size fl, with replacement, from an empirical
distribution F. This approach is referred to as resampling. Each random sample of size n
is referred to as a bootstrap sample. The empirical distribution is described by an actual
dataset. If the original dataset is:
x =(Xl, X2, .•• , XII) (22)
~.S
)
.~
The asteri.sks indicate that x* is not the actual dataset x, but rather a randomized or
The bootstrapping sample is drawn from the sample set equi vaJently with the same
probability of lin. The bootstrapping sample of size n is denoted by
X* = (X*l ,X*2, •.. ,X*II)
resampled version of it. The resampled data describe an empirical distribution,
(23)
.,
)
.~
....
........
any given bootstrap sample.
Calculation of a statistical value (8) for each bootstrap sample, i.e. mean, standard
Since the sampling is done with replacement, it is possible to have repeated values within
F(:"'"\'"* I, x*2,···, x*II ) (24)
deviation, 95th percentile is done such that
8 =s(X*) (25)
where s(X*) is a statistical estimator applied to a bootstrap replication (bootstrapping
sample) of the original dataset. To estimate the uncertainty in the statistic, n bootstrap
samples may be simulated to yield n estimates (replicates) of the statistic.
 29 
8N=S(X*N), where n = 1,2, ... , n (26)
The n estimates of the statistic may be used to construct a sampling distribution for the
statistic (Frey and Burmaster, 1999).
Bootstrapping simulations are usually applied to provide the confidence intervals
and prove the assumed distribution of unknown variables, which have limited field data
and unknown distributions. As a rule of thumb, the assumption would be acceptable if the
expected value lies in the range of statistics of the bootstrapped samples.
3.4 Maximum Likelihood Estimators (MLEs)
The maximum likelihood estimators (MLEs) are the parameters of the function that
maximize the likelihood of the distribution given a set of observations. The MLEs are
derived from the input data set and are different for each distribution function. For any
density distribution f(x) with one parameter, a., and a corresponding set of observational
'.
data, Xi, define an expression called the l.ikelihood: (Palisade, ]997)
L =ITf(Xj , a)
To find the MLEs, simply maximize L with respect to a:
dLldo. =0
and solve for a.
(27)
(28)
'0,..
" ..
Based on this method, all of the parameters that fit the distribution would be estimated.
The distribution wm be fitted using goodness of fiL Usually, there are three methods for
achieving goodness of fit. (Palisade, 1997)
1. The ChiSquare Test
2. The KolmogorovSmimov (KS) Test
·30·
3. The AndersonDading (AD) Test
The Chisquare statistic is defined as:
where
Pi = the observed probability value for a given histogram bar,
(29)
Pi =the theoretical probability that a va]ue will fall with x range of the histogram
bar.
The KolmogorovSmirnov statistics is defined as:
Dn=Sup[lFn(x)  F(x)1l
where:
Sup(x) refers to the maximum value of the function, it equals to Max(x)
n =total number of data,
F(x) = the hypothesized distribution,
N
Fn(x) = _x ,
n
Nx =the number of x/s less than x.
The AndersonDaring statistic is defined as:
+«>
A/ =n J[Fn (x)  F(X)]2l1'(x)f(x)dx
where:
? 1
''1'= ,
F(x)[l F(x)]
f(x) =the hypothesized density function,
F(x) =the hypothesized distribution,
 3 ~ 
(30)
(3l)
s
.)
~
.'.i•
:j
:li
j
.~
'".,
.~
F  N,
n  ,
n
Nx =the number of Xi'S less then x.
Each of these methods has its advantages and weakness. The ChrSquare is the most
common goodness of fit test (Palisade, ]997). The weakness is that there are no clear
guidelines for selecting intervals. A different conclusion might be reached depending on
how the intervals are specified. The KolmogrovSmirnov test does not depend on the
number of intervals, which makes it more powerful than the ChiSquare test (Palisade,
1997). However, it does not detect tail discrepancies very well. While the Anderson
Darling test is similar to the KolmogrovSmirnov test it places more emphasis on tail
values, and can only be used with actual sampling data (Law and Kelton. 1993; Walpole
and Myers, 1993). In selecting the fitted distribution, the AD method was used in this
work.
3.5 Analysis Tools
Although a proper model for a specific prohlem is the first step of running a Monte
Carlo simulation, the variable distribution used in conjunction wilh the selected model is
the most critical item during a simulation event because all the simulating values will be
sampled from the distribution. Analysis software @RISK® (Palisade, Inc.) was used for
this thesis. @RISK® is a professional analysis software for risk assessment that uses the
Excel® spreadsheet. @RISK® uses Monte Carlo simulation techniques to combine all the
1.mcertainties identified in the modeling situation to predict the likelihood of occurrence
for each possible value. Minitab® (Minitab Inc.) is another professional software used for
 32
.)
"..i
.J,
:~
', .\
~,'"
,I'
statistics. In this thesis, Minitab® was used for bootstrapping sampling to produce
confidence intervals for chemicals in the influent and the ambient conditions based on the
actual data measured.
3.6 Model Setup for IWTP with TOXCHEM+V3 and WATER9
The treatment facility was mainly modeled as units in series and parallel. These units
included a primary paint chip clarifier, two covered blending basins, two oilwater
separators, two storage tanks, two aerated equalization basins, three mixing basins and a
solid contact clarifier, a sand filter and a chlorine contact chamber. Among these units,
the blending basins, oilwater separators, equalization basins and aerator basins were in
parallel while the others were in series (Veenstra, 200 I). The layout of the plant
configuration is presented in Appendix Ai.
For the two different estimating softwares, the plant configurations were slightly
different. Since there are several parallel units in the plant, the configuratlon should
assure the flow into each of parallel units is exactly half of the lotal i.nfluent. WATER9
and TOXCHEM+V3 have different configurations to set up parallel units. Figures 6 and
7 below are the modeled plant configurations by WATER9 and TOXCHEM+V3,
respectively.
·33 
.)
,:1
:";.
".
Figure 6. IWTP Configurations with WATER9
Notes:
1: Primary municipal clarifier modeling stripping waste clariifier
2: Storage tank modeling DIblending tank
3: Storage tank modeling D2blending tank
4: Open sump modeling venturi pipe from Building 3001
5: Weir/waterfall modeling the mixing unit of two influents
6: Covered separator modeling oilwater separatorl
7: Covered separator modeling oilwater separator2
8: Storage tank modeling storage tankl
9: Storage tank modeling storage tank2
10: Equalization basin modeling equalization basinl
11: Equalization basin modeling equalization basin 2
12: Mix tank modeling mixing basinl
13: Mix tank modeling mixing basin2
14: Weir/waterfall modeling diversionary structure (DS) between mixing basin 2 and 1
15: Mix tank modeling mixing basin3
16: Weir/waterfall modeling diversionary structure (OS) between mixing basin 2 and
SCC
17: Circular clarifier modeling solid contact clarifier (SCC)
18: Weir/waterfall modeling diversionary structure between SCC and wet well
19: Open sump modeling wet wen
20: Hard piped, no space unit modeling sand filter
21: Open sump modeling chlorine contact chamber
22: Circular darifier modeling thickener
23: Porous solids unit modeling filter press
24: Oil film unit modeling the final unit of sludge treatment, plate and frame filter
 34
~...,I
.~ '.
..o•J
..,
.,
y
;~:
"..,
'..
::, ",
f!!!
:..~
..
)
.i
Figure 7. IWTP Configuration with TOXCHEM +V3
Notes:
B5: Primary clarifier modeling stripping waste clarifier
C5: Blending tank modeling DIblending tank
D5: Blending tank modeling D2blending tank
B9: Channel modeling venturi pipe from building 3001
C9: Wastewater mixer modeling the mixing unmt of two influents
D9: Drop structureopen modeling diversionary structure (DS)
DlO: Drop structureopen modeling diversionary structure (DS)
F9: Dissolved airfloatation modeling oilwater separatorl
G9: Dissolved airfloatation modeling oilwater separator2
H9: Equalization basin modeling storage basin
19: Equalization basin modeling storage basin
19: Equalizationmixed/aerated basin modeling aerated equalization basinl
K9: Equalizationmixed/aerated basin modeling aerated equalization basin2
L9: Equalizationmixed/aerated basi n modeling mix ing basinl
M9: Equalizationmixed/aerated basin modeling mixing basin2
09: Equalizationmixed/aerated basin modeling mixing basin3
P9: Drop structureopen modeling weir (diversionary structure) between mixing basin
and SCC
R9: Secondary clarifier/sludge thickener modeling solid contact clarifier (SCC)
S9: Channel modeling wet well
T9: Force main modeling sand filter
U9: Drop structureopen mode.ling pump station between sand filter and chlorine
contact chamber
V9: Channel modeling chlorine contact chamber
 35 
.,
'.j
'" ,I..
~:
r'
' ..
:~
.~t '.
3.7 General procedures
The work steps in this thesis can be summarized as follow:
1. Field data acquisition
There are two potential SOUl'ces of field data for this thesis. One is from the record
logs (consumption of chemicals) over the past several years, from 1999 through 2001.
The other source is field sampling and measurement. A flux chamber was used for
field sampling of emissions. The flux chamber is one of the most promising
technologies for direct measurement of VOCs emissions devdoped by USEPA (Shen,
et al., 1993). SUMMA canisters, where the contents were analyzed by a gas
chromatography and mass spectrometer (GCMS), and an online FTIR (Fourier
Transform Infrared Spectmscopy) were used to measure the concentration of
chemicals in gas phase samples taken from the flux chamber. The liquid phase
constituents were measured by GCMS. The data used in this thesis is historic data. In
addition, methyl.ene chloride, acetone and 2butanone must be backsolved to the
concentration at the fmnt end of the plant since they were sampled at the transfer pit
which is in the middle of the plant.
2. Sensitivity study
TOXCHEM+V3 has an integrated subroutine for conducting a sensitivity study. The
subroutine was used to detennine which variables in the model, such as site
parameters, influent characteristics, or chemicals' physiochemical properties, would
have significant effect on the air emissions. The results were then applied to
WATER9 since WATER9 does not contain such a subroutine.
 36
)
.j
.
...
"
3. Fit distribution
One of the subroutines ill @RISK®, called "Best Fit", was used to fit the distribution
for the input data. The fitted distributions were ranked using ChiSquare test,
KolmogorovSmirnov test and AndersonDarling test based on the least square
method. The AndersonDarling test was used to select the distribution that had best
fit. The lower the value of the AD statistic of a distribution that the test ends up with,
the better the distribution fits the actual data. For each distribution fitting, the fitting
results were presented by a histogram, ProbabilityProbability (PP) plot and
QuantileQuantile (QQ) plots. For the best fitted distribution, these two plots would
be nearly linear.
4, Distribution confirmation
Since the distributions of most variables used in the model were not known,
assumptions concerning these variables had to be made or they must be fitted from
actual data. Bootstrapping simulation was used to make the assumed or fitted
distribution more appropriate and acceptable. Bootstrapping constructed a confidence
interval for each variable of interest based on actual data. If the statistical values of
the assumed or fitted distributions fell in the range that the Bootstrapping simulation
constructed, the assumption or fitted distribution was then acceptable.
5. Generation of random sample data sets for each model, WATER9 and
TOXCHEM+V3
Generating random numbers was an important step in simulation. Since the
WATER9 and TOXCHEM+V3 cannot generate random numbers themselves,
random numbers were generated externally. At this step, @RISK@ generated
 37
.'
.' ,.•
.1~
:;
'.'.
•1 '.
",
',
"',
hundreds of random number sets with the Latin Hypercube sampling method based
on the assumed or fitted distributions.
6. Running of models
All of the generated random number sets were substituted into the models manually.
The models were run with these data sets for two hundred times. The number of
simulations was determined based on the study of adequacy of simulations.
7. Results analysis and uncertainty estimate
Data analysis had as a focus the construction of a cumulative distribution for each l~
variable, on which uncertainty estimates were based. Furthermore, comparisons of ..
'.
detenninistic and probabilistic results helped evaluate the model itself, i.e.
applicability. The general procedure is illustrated in Figure 8.
"
,)
Rel;ultl;
Analysis
Probability Sampling Model
distribution
Field
Dala
Figure 8, Diagram of simulating procedures
 38 
Chapter IV
Results
4.1 Sensitivity study
Sensitivity studies of all of the chemicals in the influent by TOXCHEM+V3 were "
performed by ranging the values of parameters from 0.2 to 5 times their values. These
studies showed that VOCs emissions are mainly sensitive to influent flow rate,
concentration and slightly sensitive to influent water temperature, but not sensitive to air
temperature, pH, plant elevation, oil/grease concentration in the influent, volatile
I
\.
.I
distributions, of the other parameters in the models can be substituted.
of the flow rate and concentration, the averages over years, also called the single value
(Figure 9, Figure 10 and Figure 11). Thus, when the models were run, with the exception
i
 ._.
.. acetone
* benzyl alcohol
 butanone
+ phenol
2  .==11=. • • • • • ==
8.,
7 1,1
~ 6 @ 5 11
.~ 4 t===:::~;~~~~~~::==J
. 3  A''''
~
suspended solid ratio and wind speed. The results are presented in the following figures
methylene
chloride
5 JO 15 20 25
]11
o . .... ~=*=*~=_~~t:;=t;;;;;;;;;t:;~__l
o
Wind Speed (mph)
Figure 9. Sensitivity of wind speed on emission rate.
 39
35 r,
30 +.
::0 25 +1 ::c
:::: 20 )1 'c": o
o~ 15 _1 _
·s
t:.J;;l 101
5 + w.'!'!?"'''j
+ phenol
butanone
methylene chloride
~ benzyl alcohol
 Acetone
20 40 60 80 100 120
o~....,....~~~~~~__,1
a
Concentration (mgll)
Figure 10. Sensitivity of chemical concentration on emission rate. .I
+ phenol
....
.~.....
f)
d
methylene
chloride
~ benzyl alcohol
butanone
I_ODEOS 2_00E05 3.00E05 4.00E05 5.00E05
I.OOE02 r,
Ir, 9.00E03 +
,.. 8.00E03 '11
 
~6 7.00E03
~ 6.00E03 i ..£.11
:I:
~ S.OOE03 \
.~ 4.00E03 ) ;;.A'_.._""',__""""ll
o§ 3.00E03 1 ~~~::o:..==t_='_=II
t:.J;;l 2.00E03 t""V""':::AF==~~II
I.OOE03 t:A~"""'11
O.OOE+OO fJ,=,.~.__1
O.OOE+OO
3 Flow rate (Mm Is)
Figure 11. Sensitivity of flow rate on emission rate.
 40
4.2 Models inputs setups
Both models were run two hundred times. Based on the results of the sensitivity
studies, the yearly averages were substituted into the models for those insensitive
parameters, i.e., total solids, suspended solids, oil/grease concentration, dissolved solids,
volatile S5 ratio, radius of pipe, air temperature, wind speed, pH, wastewater temperature
and plant elevation. The influent wastewater temperature had an average of 68.7 OF with
a standard deviation of 8.5 OF. As such, it appears to be relatively stable and uniform.
These parameter values mentioned above remained unchanged during the simulations.
The inputs of these parameters are listed as follows:
Table 2. lnputs of variables for WATER9 and TOXCHEM+V3 simulation
I Parameter 1Units IWATER 91 TOXCHEM 31 Note (TOXCHEM) I
Total solids pom 827 724 Suspended solids
Oil/qrease mall 14.04 14.04
Dissolved solids ppm 588 24% Volatile SS ratio
Radius of pipe em 15.24 15.24
Air temperature F 61.02 61.02
Wind speed MIs 3.25 3.25
pH SU 7.79 7.79
Wastewater temp OF 68.7 68.7
Elevation ft 1250 1250
Note: SUstandard unit
Since benzyl alcohol was nol actually sampled and measured, its concentration was
estimated from the concentration of phenol according to a proportionality between these
two chemicals. The benzyl alcohol and phenol were assumed to follow the same
distribution. The equation for calculation of benzyl alcohol by this method could be
expressed as foHows (Veenstra, et al., 2001):
 4) 
where
Cba is the concentration of benzyl alcohol,
Cp is the concentration of phenol,
Xba, Xp is the composition of benzyl alcohol and phenol in water, respectively,
Uoa, Up is the amount of benzyl alcohol and phenol used, respectively,
Hba,Hp is the Henry's constant of benzyl alcohol and phenol, respectively,
VPba. Vpp is the vapor pressure of benzyl alcohol and phenol, respectively, and
Dba. Dp is the diffusivity coefficient of benzyl alcohol and phenol, respectively.
Whether the diffusivity of benzyl alcohol and phenol is considered or not, it depends on
value of the last item on the right side of equation 32. If the diffusivity is considered, the
value of the last item on the right side of equation 32 is 0.29289. If the diffusivity is not
considered, the value equals 1 (Veenstra, et aI., 2001). In equation 32, the estimated
concentration of benzyl alcohol in the influent would be higher if diffusivity of benzyl
alcohol and phenol was not considered. Then benzyl alcohol air emission rate would be
higher as a result. This scenario can be considered a worse case and was the main
approach used in this work. When diffusivity was not considered, the estimated influent
concentration of benzyl alcohol is about 3.64 times that of phenol (Veenstra, et aI, 2001).
4.3 Fitting Distributions
Distributions were developed for the individual compounds based on the analysis of
wastewater constituents in the transfer pit at Tinker's IWTP in 1999 and 2000. Since
 42·
"
.;,./ ..
~ ...
:':s i:~ .~:
I"
.~
these chemicals' probability density curves have a strong tail shape, the Anderson
Darling method, which places more emphases on the tail, was used for data distribution
fitting and the results are presented as follows:
I. Benzyl alcohol and phenol
Benzyl alcohol and phenol were assumed to fonow the same distribution since
the concentration of benzyl alcohol was estimated based on the concentration of
phenol. Phenol data was best fitted to be a Pearson5 distribution. Pearson5 (a,~)
distribution is defined by the shape factor a and the scale factor ~. Its statistics are
calculated as follow:
Mean =~/(aI)
Variance = ~2/(a1)\a2) if a> 2
(33)
(34)
i
"Ill
j~
These calculations are the theoretical values of mean and variance when samples
are continuous. The fitting figures are shown as follows:
..1
,,to,
rl
j
"I,
" ....
Pearson5 (1.3910, 30.945) Shift =3.4339
x<= 4.8922 X <= 212.46
3 .__~~.02!,o/o' ~~°I!L ,
o 100 200 300 400 500 600
01+
100
2
2.5 "
0.5
~
<o,....
)(
~ 1.5
:a
~o
Ii:
Concentration (mgll)
Figure 12. Fitted distribution of phenol
The PP and QQ plots of phenol are shown as Figures 13 and 14.
 43 
Pearson5 (1.3910, 30.945)
Shift =3.4339
v 0.8 ) 106 /
]u:: 0.4/ 0.2
0 ,
o 0.2 0.4 0.6 0.8
Input pvalue
Figure 13. PP plot of fitted distribution for phenol
Pearson5 0.3910,30.945)
Shift =3.4339
I
~ l ••
.ii.
"
. I
200
800
.~ 600
:'"l
CT
0 400
~
u::
a 200 400 600 800
Input quantile
"."...
Figure 14. QQ plot of fitted distribution for phenol
ProbabilityProbability (PP) graphs plot the distribution of the input data vs.
the distribution of the result. If the fit is "good", the plot will be nearly linear.
Similarly, QuantileQuantile (QQ) graphs plot the plot percentile values of the
input distribution vs. percentile values of the result. If the fit is "good", the plot
will be nearly linear. The linearity of all the PP and QQ plots in this thesis is the
 44
best among all fitted distributions for each compound. The mean calculated with
equation (33) is 79.14 mgll. A summary of the statistics is presented as follows:
Table 3. Summary of statistics for phenol
Fit Input
Distribution Pearson5 (1.3910, 30.945\ N/A
Shift 3.434 N/A
a 1.391 N/A
b 30.945 N/A
Left X 4.892 4.892
Left P 5.00"10 6.67%
Riqht X 212.457 212.457
Riqht P 95.00% 95.00%
Diff. X 207.565 207.565
Diff. P 90.00% 88.33%
Minimum 3.434 2.5
Maximum Infinity 540
Mean 75.711 55.408
Mode 9.508 15.000 rest
Median 25.339 25.5
Std. Deviation Infinity 87.265
Variance Infinity 7488.25
Skewness 6.240 [est1 3.608
Kurtosis 48.650 [est1 18.239
Note: [est] is the abbreviation of estimate. This note is applicable to all of
the tables in this thesis.
The mean of 79.14 mg/l calculated from equation (33) stands for the theoretical
value of the Pearson5 distribution with certain a and ~ values when samples were
continuous, while the mean of 75.711 mg/l from Table 3 was obtained from actual
data being fitted. Using different test methods, i.e., ChiSquare, AndersonDarling
(AD) and KolmogorovSmirnov (KS), the fitting will present different statistical
results. Table 4 presents different statistics of different te t method, i.e. test value,
P value, critical value at 0.05 significance level, number of bins and data points.
For each distribution fitted, a similar table will be produced.
 45
Table 4. Comparisons of test methods for phenol
ChiSO AD KS
Test Value 15.6 0.4153 0.1099
P Value 0.0485 N/A N/A
Critical Value 15.5073 N/A N/A
# Bins 9 N/A N/A
Data points 60 60 60
2. Acetone
Acetone data was best fitted to a Pearson5 distribution. The filling figures
are presented as follows:
Pearson5 (4.1981, 1200.6) Shift =+66.950
J
:»
l
I.
".
i
.I •.
0.9 1.1
X <= 876.97
0.2 0.3 0.4 0.5 06 0.7 0.8
0.5
3
3.5 r2.;~~~_____,
M
<I 2.5
o....
~ 2
~
:c 1.5
(lI
.0o 1
~
Concentration (rogJl.
Figure 15. Fitted distribution of acetone "
~..:.."
The PP and QQ plots of acetone are shown in Figure 16 and Figure 17.
46
0.2 0.4 0.6 0.8
Input quantile
Values in Thousands
Pearson5 (4.1981,1200.6)
Shift =+66.950
1
~ 0.8 // ~6.. 0.6 (
] 0.4  /~j
[; 0.2 r
o \'1+++1
o 0.2 0.4 0.6 0.8
Input pvalue
Figure 16. PP plot of fitted distribution for acetone
Pearson5 (4.1981, 1200.6)
Shift =+66,950
*l 1
<l) §
3 :g 0.8
I:l 0
~ ESO.6
'2 .5004
....... en i£ .2 0,2 
ro > 0 +l+++1
o
Figure 17. QQ plot of fitted distribution for acetone
According to equations (33) and (34), the mean and variance are 375.41 f..lgll and
64118.8, respectively. A summary of the statistics is presented as follows:
·47 
I"
i
~,
:i
I.
"'.
"
.~.: :,
Table 5. Summary of fitting statistics of acetone
Fit Input
Distribution Pearson5 (4.1981, 1200.6\ N/A
Shih 66.950 N/A
A 4.198 N/A
B 1200.594 N/A
Left X 216.374 216.374
Lett P 5.00% 4.17%
Riaht X 876.970 876.970
Riaht P 95.00% 95.83%
Diff. X 660.596 660.596
DiU. P 90.00% 91.67%
Minimum 66.95 189
Maximum Infinitv 988
Mean 442.36 437.29
Mode 297.92 267.00 restl
Median 377.19 395
Std. Deviation 253.22 208.11
Variance 64117.9 41504.29
Skewness 2.371 restl 1.205
Kurtosis 10.962 rest] 3.690
Table 6 lists the different statistics of different test methods.
Table 6. Comparisons of test methods for acetone
ChiSq AD KS
Test Value 2.6667 0.2485 0.1173
P Value 0.6151 N/A N/A
C.Val @ 0.05 9.4877 N/A N/A
# Bins 5 N/A N/A
Data points 23 23 23
3. Methylene chloride
Methylene chloride was best fitted to be a Pearson5 distribution. The fitting
figures are shown as follows:
 48
.,
••
1
.\
:~
:3,
.)
',j
..,.
,.
"
Pearson5 (4.4409,836.91) Shift =+13.988
X <= 113.90 X <= 528.39
o
o
50 100 150 200 250 300 350 400 450 500 550
6
5
<;l 4  <0
...
)( 3 'C/l
Q)
::J
>'iii 2
Conncentration (ugll)
Figure 18. Fitted distribution of methylene chloride
The PP plot and QQ plot are shown in Figures 19 and 20.
i
•i
Pearson5 (4.4409, 836.91)
Sh ift = +13.988
~.
I
v 0.8 
::::l
';;" 0.6
I
0..
v 0.4  ~
u::
0.2
0
0 0.2 0.4 0;6 0.8
Input pval,ue
Figure 19. PP plot of fitted distribution for methylene chloride
 49
Pearson5 (4.4409, 836.91)
Shift =+13.988
700
~ .!! 600
C 500
«l =' 400  0
"0 300 
~ 200 /.J
~ J
100
0
0 10 20 30 40 50 60 70
0 0 0 0 0 0 0
Jnput quantile
Figure 20. QQ plot of fitted distribution for methylene chloride
According to equations (33) and (34), the mean and variance are 243.224 J.lgl1 and
24236.14, respectively. A summary of the statistics is presented as follows:
Table 7. Summary of fitting statistics for methylene chloride
Fit Input
Distribution Pearson5{4.4409, 836.91) N/A
Shift 13.988 N/A
a 4.441 N/A
b 836.906 N/A
Left X 113.896 113.896._
Left P 5.00% 4.17%
Riaht X 528.394 528.394 
Riqht P 95.00% 100.00%
Ditt. X 414.498 414.4..9_8 
Diff. P 90.00% 95.83%
Minimum 13.988 94.2
Maximum Infinity 523
Mean 257.21 253.84
Mode 167.8 143.00 [est]
Median 217.49 223.5
Std. Deviation 155.68 126.13
Variance 24234.87 15245.39
Skewness 2.270 [estl 0.891
Kurtosis 10.318restl 2.592
 50·
Table 8 lists the different statistics of different test methods.
Table 8. Comparisons of test methods for methylene chloride
ChiSq AD KS
Test Value 2.25 0.3379 0.1065
P Value 0.6899 NJA N/A
C.Val @ 0.06 9.4877 N/A N/A
# Bins 5 N/A N/A
Data points 24 24 24
4. 2Butanone
2butanone, however, was best fitted by an Exponential distribution. For the
Exponential distribution [Expon (~)], ~ is the only parameter and refers to mean.
Then
Mean = ~
Variance = ~2
The fitting plots are presented in the following figures.
Expon (139'.03) ShUt = +3.5972
X <= 10.728 X <= 420.08
(35)
(36)
8
7 .
6 
"?
< 5 . 0...
)0( 4' II!
Gl
~ jij 3
>
2
1
o·
100 0 100 200 300 400 500 600
Concentration (ugll)
Figure 21. Fitted distribution of 2butanone
 51 
Expon (139.03) Shit! =+3.5972
I .r;:,
0.6
0.4
/
0.2 ~.~
Q) 0.8:
J c;l
;>
I
0
"'0
...Q......).
u::
o +++I'+~
o 0.2 0.4 0.6 0.8
Input pvalue
Figure 22. PP plot of fitted distribution for 2butanone
.,
Expon (139.03) Shift =+3.5972
,I
600
.~
500 / .c..:. 400
c;l
;:I
0" 300
0
<U ........ 200 
LI:
100  )
/
_.
a 
a 100 200 300 400 500 600
Input quantile
:',
11
II
..
.1
'j
".,
Figure 23. QQ plot of fitted distribution for 2butanone
In term of equations 35 and 36, the fitted mean and variance for 2butanone
are 139.03 Ilg/I and J9329.34, respectively. A summary of the statistics is
presemed as follows:
 52 
Table 9. Summary of fitting st.atistics for 2but.anone
Fit Input
Distribution Expon(139.03) N/A
Shift 3.597 NIA
B 139.026 N/A
Left X 10.728 10.728
Left P 5.00% 4.1,7%
RiQht X 420.083 420.083
RiQht P 95.00% 91.67%
Ditf. X 409.3541 409.354
Ditt. P 90.00% 87.50%
Minimum 3.597 9.39
Maximum Infinity 567
Mean 142.162 148.42
Mode 3.5972 114.00 1est
Median 99.963 104.45
Std. Deviation 139.03 147.94
Variance 19328.3 20974.11
Skewness 2 1.378
Kurtosis 9 4.215
Table1 0 lists the different statistics of different test methods.
Table 10. Comparisons of test methods for 2butanone
ChiSq AD KS
Test Value 0.5833 0.3135 0.1243
P Value 0.9649 >0.25 > 0.25
C.Val @ 0.05 9.4877 0.8572 0.2128
# Bins 5 N/A N/A
data points 23 23 23
5.. Flow rate
Flow rate was best fitted to a Nonnal distribution. A normal distribution has
two parameters; one is the mean (11), the other is standard deviation Co). The fitting
plots are shown in the foHowing figures.
 53 
•
x <= 0.50048
Normal (0.74348, 0.14773.)
x <= 0.98648
0.5 0.6 0.7 0.8 0.9 11
4 ..?JI.Po.~()'J<>~__,
0.5
0+0.4
35
3
g 2.5
E2 2"
f?
a. 1.5
Flowrate (MGD)
Figure 24. Fitted distribution of flow rate
(Note: MGDmillion gallon per day)
The PP and QQ plots of flow rate are shown as Figure 25 and Figure 26.
Normal (0.74348, 0.14773)
1....,
<l.)
::l 0.8"@
f 0.6
0
~........ 0.4
ii: 0.2
Ollt++1
o 0.2 0.4 0.6 0.8
Input pvalue
Figure 25. pop plot of fitted distribution for flow rate
 54
Normal (0.7434R. 0.14773)
1.2
<l) 1 
:.:.:..:. c 0.8 _.
('J
~
0' 0.6
'0
..<...l.). 0.4  .......
LL:
0.2
0
0 0.2 0.4 0.6 0.8 1 1.2
Input quantile
Figure 26. QQ plot of fitted distribution for flow rate
The fitted mean and variance are 0.743 million gallons per day (MOD) and 0.148,
respectlvely. A summary of statistics is presented in Table 11.
Table 11. Summary of fitting statistics for flow rate
Fit Input
Distribution Normal (0.74348, 0.14773 N/A
m 0.743 N/A 
s 0.148 N/A
Left X 0.500 0.500 .
Left P 5.00% 9.78%
Riqht X 0.986 0.986
Riqht P 95.00% 93.48%
Diff. X 0.486 0.486
Diff. P 90.00% 83.70%
Minimum Infinity 0.5
Maximum +Infinity 1
Mean 0.743 0.743
Mode 0.743 0.900 [est]
Median 0.743 0.7
Std. Deviation 0.148 0.148
Variance 0.022 0.022
Skewness 0 0.0158
Kurtosis 3 1.865
 55 
Table 12 lists the different statistics of different test methods.
Table 12. Comparisons of test methods for flow rate
ChiSq AD KS
Test Value 101.696 2.76 0.160
P Value 0 < 0.005 < 0.01
C.Val @ 0.05 18.307 0.746 0.093
# Bins 11 N/A N/A
data points 92 92 92
As an alternative to using concentration, the mass of each contaminant 10 the
influent was used to generate the bestfitted distribution. The results showed that the
same distribution was arrived at for each contaminant using either mass or
concentration. This result is detailed in the discussion section that follows.
4.4 Bootstrapping simulation
To confirm these fitting distributions, the Bootstrapping simulation was also
applied. Samples were bootstrapped from the actual field data (historic concentration
data from Tinker's IWTP) with replacement, and the mean and standard deviation of
each bootstrapped sample were calculated and then the cumulative distribution
function was plotted for all bootstrapped samples. A bootstrapping sample consisted
of 20 samples taken from field data. In total, 100 bootstrapping samples were
developed. The results are presented in the following section.
The bootstrapping results for flow rate are presented in Figure 27. The mean of
each bootstrapping sample was calculated. Figure 27 shows the cumulative
distribution of the means of lOa bootstrapping samples for flow rate.
 56
 ...•
1.6
.. I...... I
1.0
0.9
0.8
~ 0.7
:: 0.6
.c
~ 0.5
o Q: 0.4
0.3
0.2
0.1
0.0
o 0.2 0.4 0.6 0.8
Flow rate (mgd)
Figure 27. Cumulative distribution of flow rate with bootstrapping
The standard deviation of each bootstrapping sample was calculated and the
cumulative distribution of the standard deviation is presented in Figure 28.
,.,
 • ~
.. ..•  II'
r 
• •
 .. ~
.....1..1 l..wl
1.0
0.9
0.8
0.7
~= 0.6 :c
(II 0.5 .c
~ 0.4
Q.
0.3
0.2
0.1
0.0
a 0.05 0.1
Standard deviation (mgd)
0.15 0..2
Figure 28. Cumulatlve distribution of the standard deviation of tlow rate with
bootstrapping
The bootstrapping results for acetone are presented in Figure 29. The mean of
each bootstrapping sample was calculated. Figure 29 shows the cumulative
distribution of means of all bootstrapping samples for acetone. [n Figure 30, the
·57·
standard deviation of each bootstrapping sample was calculated and the cumulative
distribution of the standard deviation is presented.
,.1.1'
II
~ •
J•,
~
I"
I'
J
,..J
I~
1
0.9
0.8
0.7
~ 0.6
:cIV 0.5
.0e 0.4
Q.
0.3
0..2
0.1
o
o 100 200 300 400 500 600
Concentration (ug/I)
Figure 29. Cumulative distribution of acetone with bootstrapping
r IF
Ilil lII _. I 
~~
tt1
~,... I
'll t
Ij 
J • r f
~
.  I'" ..•  ._ ..,. t  ~ t 
4
....
1.0
0.9
0.8
0.7
>~
0.6
:cIV 0.5
.0e 0.4
a.
0.3
0.2
0.1
0.0
a 50 100 150 200 250 300 350
Standard deviation (ugll)
Figure 30. Cumulative distribution of standard deviation of acetone with
bootstrapping
 58·
The bootstrapping resutts for methylene chloride are presented in Figure 31. The
mean of each bootstrapping sample was calculated. Figure 31 shows the cumulative
distribution of means of all bootstrapping samples of methylene chloride.
,IW
I'"
· IF ....
~
· ..
" · r
I.
,'IT
po'
· ~""
1.0
0.9
0.8
0.7
~. 0.6
:0
('Q 0.5
.t:J e 0.4
D.
0.3
0.2
0.1
0.0
o 50 100 150 200 250 300 350
Concentration (ugll)
Figure 31. Cumulative distribution of methylene chloride with bootstrapping
The standard deviation of each bootstrapping sample was calculated and the
cumulati ve distribution of the standard deviation is presented in Figure 32.
 59
.. ..., ...
~
w
• .. I"
~ • ..•~•II' r
~ ..... ...........
1.0
0.9
0.8
0.7
>
== 0.6
.0,
~ 0.5 e 0.4
D..
0.3
0.2
0.1
0.0
o 50 100
Standard deviation (ugll)
150 200
Figure 32. Cumulative distribution of standard deviation of methylene chloride with
bootstrapping
The Bootstrapping results for 2butanone are presented in Figure 33. The mean of
each bootstrapping sample was calculated. Figure 33 showed the cumulative distri.bution
of means of all bootstrapping samples of 2butanone.
I.
...
... ..... II
I    I .. • _.  •  I

II I   
r .,.. I
•" .~
J    1 " ~
I
,JI
......~
I.... ~
1.0
0.9
0.8
0.7
:>=0.6
:0ca 0.5
.0e 0.4
ll.
0.3
0.2
0.1
0.0
o 50 100 150 200 250
Concentration (ug/I)
Figure 33. Cumulative distribution of 2butanone with bootstrapping
The standard deviation of each bootstrapping sample was calculated and the
cumulative distribution of the standard deviation is plotted in Figure 34.
 60
,
., ;I
..I' • I
! II •.
~
I ,
~
I.,..4l ...,
~
1.0
0.9
0.8
0.7
~ 0.6
:0
cu 0.5
.0e 0.4
0.
0.3
0.2
0.1
0.0
o 50 100 150 200 250
Standard deviation (ugll)
Figure 34. Cumulative distribution of standard deviation of 2butanone with
bootstrapping
The bootstrapping results for phenol are presented in Figure 35. The mean of
each bootstrapping sample was determined.
20 40 60 80 100 120 140
1.0 .rr'''''rT,rr.,...,.rTrr.rrrIT,,,,IAI::~''''.rr,,.,
O. 9 +!++11+++++HlH++t+J<HMIIl""+1++I+tt+I+J
IHH+t++t+I+tltlItII
..:..J.o"....."..'I+I+++III  ~  
0,8 +t+HI+ttt+t+l+tt+t+,l,rlJf+tlHHj  ~;
1++++1++++1+41'I\f,""IItl  ,Ir
0,7 ~r 1
> I\+I++I+tIt+I\+++.I.~..Itt  ItI+I >     
~ O. 6 +t++li/H+tl+HHI:;.++t+++II
:0 r cu 0.5 H+HI+IlbI++t++...I.+. 1\II 1  I  rr roO
e 0.4 1+l++I+t++1+h.I.~P"'It+l++Itl++H+t+I++H
a. 1+It+l++I+tlI..+. "IIfHlII++Itt++I+1 IH++I
0,3 ++++t++t++t~'1rl++++t++t+H+HIH+t+t
0.2 +t+HI+t++h••r+I+tt++ItltlltHt++ttttH
O. 1 +t++1I++++..JI"I"4+Htt+IlrltltIIttl+tl++flO.
0 1,'.L.....I..t.'.I...ltl~....L...J+'L.J....'fl''I.l...L....Jl.'fL'L.J....t''''!
o
Concentration (mgll)
Figure 35. Cumulative distribution of phenol with bootstrapping.
 6) 
...
•
.r. •
I,. '
.. ...I • " r .. ..... 
...".
 IolII'
1.0
0.9
0.8
0.7
~:.:: 0.6
:0
C'Cl 0.5
.0e 0.4
a.
0.3
0.2
0.1
0.0
o 50 100
Standard deviation (mglJ)
150 200
Figure 36. Cumulative distribution of standard deviation of phenol with
bootstrapping.
The standard deviation of each bootstrapping sample was calculated and the
cumulative distribution of the standard deviation is presented in Figure 36.
The results of fitting distributions and bootstrapping simulation are presented in the
tables below (Tables 13 and 14):
Table 13. Comparisons of fitted and bootstrapped means
Parameters Fitted distribution Bootstrapping Mean
Distributions Mean 5% 50% 95%
Phenol/Benzyl Alcohol PEARSON5 (a, ~) 79.14 34.02 60.58 97.8
Acetone .PEAIRSON5 (a, ~) 375.41 376.04 438.75 521.35
Methylene Chloride PEARSON5 (a. ~) 243.22 208.97 250.82 297.35
2Butanone EXPON (~) 139.03 93.14 146.53 189.50
Flow rate NORMAL (11, 0:) 0.74 0.69 0.75 0.8
 62 
Table 14. Comparisons of fitted and bootstrapped standard deviation
I Parameters I Fitted distribution ' Bootstrappinq standard deviation
Distributions SO 5% 50% 95%
Phenol/Benzyl Alcohol !PEARSON5 (lX. ~) N/A 27.7 82.7 157.9
Acetone PEARSON5 (lX, ~) 253.2 147.2 203.4 274.6
Methylene Chloride PEARSON5 (lX, ~) 155.7 84.5 122.4 152.8 I
I
2Butanone EXPON (~) 139.03 81.6 147.2 189.7
Flow rate NORMAL (~, cr) 0.15 0.12 0.15 0.17
Note: the standard deViatIOn of Pearson5 dlstnbutIOn for phenol IS not available because
ex is less than 2 in equation 34.
In summary, the distributions of all the chemicals and the flow rate are listed in Table 15.
Table 15. Summary of fitted distributions for compounds and flow rate
Parameters Fitted distribution
Phenol PEARSON5 (1.3910.30.945)
BenzyI Alcohol PEARSON5 (1.3910,30.945)
Methylene Chloride PEARSON5 (4.4409, 836.91)
Acetone PEARSON5 (4.1981,1200.6)
2Butanone EXPON (139.03)
Flow rate NORMAL (0.74348,0.14773)
Note: Benzyl alcohol follows the same distribution as phenol, but the values are
3.64 times more.
4.5 Backsolver study
Since acetone, methylene chloride, and 2butanone were sampled in the transfer pit,
which is almost in the middle of the treatment plant, these values are not representative of
the influent concentrations to the pl,ant. Therefore, the backsolver routine inside
TOXCHEM+V3 was used to project the influent concentrations of methylene chloride, 2
butanone, and acetone. This method was applied and backsolving was done during phase
1 of the Tinker project (Veenstra, et aI, 2001). The influent concentrations of acetone,
methylene chloride and 2butanone were backsolved again during this work. Sampling of
 63 
 64
37. These random number sets were used for the simulat.ions.
~Acelone
 phenol
Methylene
chloride
')( 2butanone
+ Benzyl alcohol
50 JOO 150
cone. al transfer pit (ug/l)
I:: 1500 +1~/';
U
I:: 8 LOOO
500 /..1r1
O+.,....rj
o
4000
3500 1....__1
~ 3000 ~11''1
~E 2500
1lJ'
~
'+: 2000 1171/7''''1
I::
Figure 37. Linearity of backsol ver for chemicals of concern
result the compoundsI dislribUlions in both locations are the same because the
transfer pit. But what was needed is the influent concentrations of compounds. Therefore,
from the fitted distribution stood for the random concentration of compound in the
multiplying generated random numbers by the slope of each compound's curve in Figure
random number sets for the influent concentrations of compounds were produced by
basins following the primary clarifier. These phenol concentrations were regarded as the
distributions were determined based on the field data in the transfer pit. The influent
proportional to those in the transfer pit. The results are shown in Figure 37. The fitted
concentrations in both locations are simply related linearly. Random numbers generated
concentration was obtained based on samples taken from DI1D2, which are mixing
concentrations have a linear relationship with the concentrations in the transfer pit. As a
influent concentration. The backsolver showed the concentrations in the influent were
the transfer pit was performed during October and November in 1999. For phenol, the
Table 16 presents the details of the backsolving studies and the slope of each
compound. The slope of each compound will be used for estimating the influent
concentrations.
Table 16. Backsolving studies of all the compounds
Chemicals Concentration (uo/l) SloDe
10 20 35 50 100
Benzyl alcohol 346.1 692.13 1211.31 1730.32 3460.96 34.61
Acetone 272.03 544.11 952.19 I 1360.14 2720.32 27.20
Methylene chloride 376.46 752.98 1317.67 1882.35 3764.28 37.64
2butanone 308.63 617.23 1080.19 1543.12 3086.28 30.86
Note: In the tabular area of concentration, the top numerical row stands for the
concentrations of compounds in the transfer pit. The rest are the backsolved influent
concentrations of compounds.
Then, the influent concentration for each compound was calculated using the following
equation:
Cin =Ctp x Slope
where:
Cin =concentration in the influent
C.p = concentration in the transfer pit.
4.6 Simu~ation adequacy study
(37)
The adequacy of the number of simulations was studied to determine whether more
simulations are needed. The results are presented in the following figures (Figures 3R, 39,
40, 41 and 42). Due to the different scales of the yaxis, the adequacy studies are
presented in Figures 38 and 39 for TOXCHEM+V3; the adequacy studies are presented
in Figures 40, 41 and 42 for WATER9. These results showed that running 200
65 
simulations was adequate, because the results began to converge at a constant value for
each compound.
2..
+Benzyl
alcohol
Phenol
:A Acetone
* 2Butanone
0+++++1
a 50 100 150 200 250
Number of simulations
Figure 38. Adequacy of simulation using TOXCHEM+V3 for benzyl alcohol.
phenol, acetone and 2butanone.
Methy'''' I chloride
'
50 100 150 200 250
~ 3.5 J
~ 3  j
E~
.2 2.5
g 2 +1
c
Il)
g 1.5 +1
o
u
~ I
~
~ 0.5 11
a .+tl++~
o
umber of simulations
Figure 39. Adequacy of simulations using TOXCHEM+V3 for methylene
chloride
 66
l~benZYI alcohol I
\
\.
~
I..ow
IIJ'II
,
0.014
0.012
0.01
0.008
0.006
0.004
0.002
o
o 50 100 150 200 250
Number of simulations
Figure 40. Adequacy of simulations using WATER9 for benzyl alcohol
1Phenol I
50 100 150 200 250
5...,...
0+++++.4
a
35 .~
30 +_1
25 1.......1
Number of simulations
Figure 41. Adequacy of simulations using WATER9 for phenol
The curve in Figure 41 starts to tail up at the end is probably due to the
variation of the average concentration of phenol substituted into the model.
 67 
However, it tends to flat again after 200 simulations. The tendency is clearer in
Figure 45 where the standard deviation is used.
Methylene chloride
* 2butanone
~ Acetone
3   01> § 2.5
c0
2 
~
bc 1.5 OJ
uc0
u 1
OJ co
'.".. 0.5  OJ
;>
<: 0
0 50 100 150 200 250
Number of simulations
Figure 42. Adequacy of simulations I.lsing WATER9 for methylene chloride,
2butanone and acetone
Simulation adequacy studies were also performed on the standard devjation of
aU the simulated results. These results are shown in Figures 43, 44 and 45. These
results show that 200 simulations is adequate for lhis study because the standard
deviation of these compounds start to converge on a constant value.
 68 
~ BeTl,Z~l alcohol
~meth ene cl1loride
aceta e
ohenol
~ z~Dutanone
5,.
4.5 +1
41aI
c0
.<.;.=;,> 3
" 2.5 "....
'" "0 2 cS
Ul 1.5
0.5
50 100 150 200 250
Number of simulali.on
Figure 43. Adequacy of simulations using TOXCHEM+V3 with standard deviation
~benzyl alcohol
~ 2butanone
methylene chloride
~acetone:: 1
+II~I___'c...___ . __ ~
1.4 1'~t
1.6 r.
0.2 Jti!1!:..1
c
o 1.2
.~ ......
ro:
>
Q.)
'"0
'"0 0.8
~
"0
C
C':l r;;5 0.4
50 100 150 200 250
Number of simulation
Figure 44. Adequacy of simulations using WATER9 with standard deviation
 69
~
'............... .. 
.......
~
r
120
c 100
.2 .....
.~ 80
;>
oU
"0 60 "0 ;a
"c0 40 .r..o..
en
20
o
o so
+phenol
100 ISO 200 250
Number of simulation
Figure 45. Adequacy of simulation using WATER9 with standard deviation for
phenol
4.7 Model simulation results
Two hundred random number sets were generated from the fitted distributions. Each
set included flow rate, as well as benzyl alcohol, phenol, methylene chloride, 2butanone
and acetone concentrations. The concentrations of these compounds were backsolved and
then each random number set was substituted into TOXCHEM+V3 and WATER9 to
serve as model input. The two models were then run. The pooled results from these
simulations are called probabilistic results. In addition, the determini.stic runnings of
these models were also performed. The actual field data for flow rate from logs of record
and the compounds concentrations from the transfer pit were also averaged and the
averaged concentrations were backsolved. The backsolved concentrations were then
substituted into the models and the models were run. Then, the deterministic results were
obtained. Since the Henry's constants in both models were different, simulations using
TOXCHEM+V3 were run four hundred times, two hundred times with the Henry's
70 
constants in its library, the other two hundred times with the same Henry's constants as
WATER9. TOXCHEM+V3 does not have a Henry's constant for benzyl alcohol in its
library, the Henry's constant of 1.60E5 for benzyl alcohol was derived from the
Michigan Environmental Response Division website. These results are presented below.
a. TOXCHEM+V3 simulation results with Henry's constants from its library (except
benzyl alcohol)
Simulation results are presented in Figures 46, 47,48, 49 and 50.
100.00% 'I"~===""""""'I
90.00% ,
,.... 80.00% ...,. 1
~'' 70.00% +aiI
:.~....'. 60.00%
c 50.00% .1
o 40.00% .....j
~
Q:) 30.00% ...1
0.. 20.00% __j
10.00% .....j
0.00% .,......,,.....r..,,.....,j
0.000 2.000 4.000 6.000 8.000 10.000 12.000 14.000 16.000
Emission rate (Ibid)
Figure 46. Simulation results of benzyl alcohol by using TOXCHEM+V3.
2.000 4.000 6.000 8.000 10.000 12.000 14.000 16.000
100.00% r=~~~~.....++___,
90.00% 1...I'~,,~~_i
~ 80.00% +..rEl
'' 70.00% ~rt
~ 60.00% +11'1
'.::J 50.00% ,...1
~ 40.00% .
~ 30.00% .1
~ 20.00% ...1
10.00% ..1
0.00% .~r_____,__,_____,_r_____,__,___~
0.000
Emission rate (Ibid)
Figure 47. Simulation results of phenol by using TOXCHEM+V3.
 71 

~
I
,I
J
I
I
I
.~
100.00%
90.00%
80.00%
70.00%
60.000/0
50.00%
40.00%
30..00%
20.00%
10.00%
0.00%
0.000 2.000 4.000 6.000 8.000 10.000
Emission rate (Ibid)
Figure 48. Simulation results of methylene chloride by using TOXCHEM+V3.
~
~
11"
, ,J"
I
I
I
 ,   I.
r
100.00%
90.00%
,.. 80.00%
~ 70.00%
'"
60.00%
50.00%
40.00%
30.00%
20.00%
10.00%
.00%
0.000 1.000 2.000 3.000 4.000 5.000 6.000 7.000
Emission rate (Ibid)
Figure 49. Simulation results of2butanone by using TOXCHEM+V3.
 72
~
~
,,JIIIII'
I
I
I
 J
I
'
100.00%
90.00%
80.00%
70.00%
60.00%
50.00%
40.00%
30.00%
20.00%
10.00%
0.00%
0.000 1.000 2.000 3.000 4.000 5.000 6.000 7.000 8.000
Emission rate (Ibid)
Figure 50. Simulation results of acetone by using TOXCHEM+V3.
b. WATER9 simulation results using Henry's constants contained in its database
Simulation results are shown in Figures 51, 52, 53,54 and 55.
~
,r~


100.00%
90.00%
80.00%
70.00%
60.00%
50.00%
40.00%
30.00%
20.00%
10.00%
0.00%
0.000 0.020 0.040 0.060 0.080 0.100 0.120 0.140
Emission rate (Ibid)
Figure 51. Simulation results of benzyl alcohol by using WATER9.
·73 
100.00%
90.00%
80.00% ~
~ 70.00%
"" 60.00%
50.00%
40.00%
30.00%
20.00%
10.00%
0.00%
~.....
Ir

, ,
0.000 1.000 2.000 3.000 4.000
Emission rate (Ibid) (X 101\2)
5.000 6.000
Figure 52. Simulation results of phenol by using WATER9.
.....
~ r
I
I
I
I
I
 I
J 
100.00%
90.00%
80.00%
70.00%
60.00%
50.00%
40.00%
30.00%
20.00%
10.00%
0.00%
0.000 1.000 2.000 3.000 4.000 5.000 6.000 7.000 3.000 9.000
Emission rate (Ibid)
Figure 53. Simulation results of methylene chloride by using WATER9.
 74
.,."~
I
 ,
I
I
I
Ir' , ,
100.00%
90.00%
80.00%
~ 70.00%
'' 60.00% ~'g 50.00%
~ 40.00%
~ 30.00%
20.00%
lO.OO%
0.00%
0.000 0.500 1.000 1.500 2.000 2.500 3.000 3.500 4.000 4.500
Emission rate (Ibid)
Figure 54. Simulation results of 2butanone by using WATER9.


~
I
I •r

 .J
100.00%
90.00%
80.00%
70.00%
60.00%
50.00%
40.00%
30.00%
20.00%
10.00%
0.00%
0.000 1.000 2.000 3.000 4.000 5.000 6.000 7.000 8.000
Emission rate (Ibid)
Figure 55. Simulation results of acetone by using WATER9.
C. TOXCHEM+V3 simulation results with Henry's constants from WATER9 library
Simulation results are presented in Figures 56, 57, 58, 59 and 60.
 75 
,. ~
r
,
100.00%
90.00%
80.00%
~ 70.00%
~ 60.00%
C 50~OO%
tl.l 2 40.00%
& 30.00%
20.00%
1.0.00%
0.00%
0.000 1.000 2.000 3.000 4.000 5.000
,~......
.
,
Emission rate (Ibid)
Figure 56. Simulation results of benzyl alcohol with WATER9 Henry's constant
using TOXCHEM+V3
100.00%
90.00%
80.00%
70.00%
60.00%
50.00%
40.00%
30.00%
20.00%
10.00%
0.00%
0.000 100.000 200.000 300.000 400.000 500.000 600.000
Emission rate (IbId)
Figure 57. Simulation results of phenol with WATER9 Henry's constant using
TOXCHEM+V3
 76

 ~
 f
I • I
I
 I
I
J
J ,
......'
 ,I'"
 r
I  I
 I
Ir
100.00%
90.00%
80.00%
~~
70.00%
~ 60.00% .r:: 50.00%
g.... 40.00%
~ 30.00%
20.00%
10.00%
0.00%
0.000 2.000 4.000 6.000 8.000 10.000 12.000 14.000 16.000
Emission rate (Ibid)
Figure 58. Simulation results of methylene chloride with WATER9 Henry's constant
using TOXCHEM+V3
100.00%
90.00%
80.00%
~
~ 70.00%
; 60.00%
r:: 50.00%
o 40.00% ~
~ 30.00%
20.00%
10.00%
0.00%
0.000 1.000 2.000 3.000 4.000 5.00n 6.000 7.000
Emission rate (Ibid)
Figure 59. Simulation resulls of2butanone with WATER9 Henry's constant using
TOXCHEM+V3
77 
~ ,
 ,I
j ,
j ,
100.00%
90.00%
80.00%
r
~ 70.00%
~ 60.00%
..c...:. 50.00%
o 40.00%
~
~ 30.00%
20.00%
10.00%
0..00%
0.000 2.000 4.000 6.000 8.000 10.000 12.000
Emission rate (1b/d)
Figure 60. Simulation results of acetone with WATER9 Henry's constant using
TOXCHEM+V3
With different Henry's constants, the simulations produced different results. The
detailed discussions are presented in the Discussion chapter. The Henry's constants
used during the simulations and a summary of these results is presented in Tables
17, 18, 19,20 and 21.
 78
.I
\D
Table 17. Summary of probabilistic and deterministic simulations of emission rates
(TOXCHEM+V3 and WATER9 using their own Henry's constants)
Emission rate (Ibid)
Simulation Models Benzyl alcohol Phenol Methylene chloride 2Butanone Acetone
5% 0.048 0.070 1.087 0.061 0.540
TOXCHEM+V3 Probabilistic 50% 0.196 0.236 2.445 0.685 1.108
95% 1.171 2.365 5.577 2.648 3.348
TOXCHEM+V3 Determi nistic 0.38 0.538 2.873 1.044 1.648
5% 0.001 1.904 0.995 0.035 0.317
WATER9 Probabilistic 50% 0.002 6.723 2.194 0.419 0.679
95% 0.009 66.199 4.906 1.695 2.122
WATER9 Deterministic 0.003 16.352 2.590 0.629 1.029
Table 18. Summary of probabilistic and deterministic simulations of emission rates
(TOXCHEM+V3 and WATER9 using same Henry's constants)
Emission rate (Ibid)
Simulation Models Benzyl alcohol Phenol Methylene chloride 2Butanone Acetone
5% 0.011 1.979 1.079 0.060 0.495
TOXCHEM+V3 Probabilistic 50% 0.045 6.904 2.429 0.682 0.682
95% 0.314 68.020 5.519 2.644 2.644
rOXCHEM+V3 Deterministic 0.074 16.362 2.822 0.999 1.570
5% 0.001 1.904 0.995 0.035 0.317
WATER9 Probabilistie 50% 0.002 6.723 2.194 0.419 0.679
95% 0.009 66.199 4.906 1.695 2.122
WATER9 Deterministic 0.003 16.352 2.590 0.629 1.029
00 o
Table 19. Cumulative frequency of deterministic results on the simulated cumulative distribution
(TOXCHEM+V3 and WATER9 using their own Henry's constants)
Benzyl alcohol Phenol Methylene chloride 2Butanone Acetone
TOXCHEM+V3 74.0% 76.5% 65.7% 65.0% 77.0%
WATER9 72.8% 72.1% 65.8% 65.0% 76.3%
Table 20. Henry's constants used in both models (dimensionless)
Benzyl alcohol Phenol Methylene chloride 2Butanone Acetone
TOXCHEM+V3 1.60E05 5.32E05 1.21 E01 5.32E03 1.50E03
WATER9 4.54E06 0.522 1.02E01 9.80E03 1.36E03
Note: the Henry's constant for benzyl alcohol in TOXCHEM+V3 is obtained from Michigan Environmental
Response Division.
Table 21. Cumulative frequency of deterministic results on the simulated cumulative distribution
(TOXCHEM+V3 and WATER9 using same Henry's constants)
Benzyl alcohol Phenol Methylene chloride 2Butanone Acetone
TOXCHEM+V3 66.6% 71.3% 64.1% 64.3% 77.8%
WATER9 72.8% 72.1% 65.8% 65.0% 76.3% .
Chapter V
Discussion
5.1 Distributions fitting
As discussed in the section on Monte Carlo simulation, a proper model for a specific
problem is the first step of running a Monte Carlo simulation. Variables distributions
used in conjunction with the selected model are the most critical elements during a
simulation event because all the simulating values win be sampled from the distribution.
Based on the principle of mass balance, the emission rate of VOCs is proportional to
the concentration in the influent. If the exact amount of these chemicals is known in the
influent, the emission rates can be determined with few uncertainties. However, the
chemicals in the influent fonow speci.fic distributions. It is critical for accurate
simulations to know the chemical's distribution. If there were adequate concentration
data available to fit a distribution for these chemicals, the fitted distribution would match
the actual distribution better. The fitted results would be more accurate. As the shapes of
the fitted distributions for the chemicals of interest were generally skewed, the AndersonDarling
test method, which places more emphases on the tail, was utilized to select the
proper fitted distribution. Distribution fitting of these chemicals was performed using the
sampled concentrations and calculated masses. The selection of the bestfit distributions
was based on the concentration data. The distribution developed using masses was used
for reference, or to confirm the distribution that was selected based on concentration, and
to try to help limit scatter in the data. The results were consistent except that there was a
slightly different ranking of the fitted distribu60n for 2butanone. An Exponential
(EXPON (~) distribution was ranked best and a LOGNORM2 (Il, d) distribution ranked
 81 
b
second using the AD test with concentration, while the ranking were reversed when they
were fitted with mass. The mass was calculated by multiplying the flow rate and the
concentration of each compound on the same day at the sample unit The sample unit is
the transfer pit for acetone, 2butanone, and methylene chloride. It is DlID2 for phenol
and benzyl alcohol. The flow rate used for fitting distribution was obtained from the
effluent record of the first quarter 2001 of the treatment plant. Since the flow
measurement is at the tail end of the plant and there is a lag time within the plant, the
EXPON (~) distribution was still considered to be the best fit for 2butanone in the
influent. For phenol, the fitted distribution using both concentration and mass were the
same. Both data sets were best fit by the PEARSON5 (a., ~) distribution. For acetone, the
top three fitted distributions were INVGAUSS (~,A), LOGNORM2 (Jl, 0'2) and
PEARSON5 (a., ~) based on mass and concentration. For methylene chloride, the top
three fitted distributions were INVGAUSS (Il,A), LOGLOISTlC (y, ~, a.) and
PEARSON5 (a., ~) based on mass and concentration. To pick a distribution thought to be
the best among these fitted ones, Bootstrapping simulation was used to provide a range of
the mean and standard deviation of the actual field data of chemicals. The theoretical
values of the mean and standard deviation of each distribution were calculated and
compared to the Bootstrapping results. The closer the mean and standard deviation of the
fitted distribution were to the mode of the bootstrapped values, the better the fitted
distribution would represent the true distribution. The theory behind this procedure was
discussed in detail in the Methodology section. On the other hand, the more data that is
used for fitting, the more accurate the result would be (Palisade, 1997). Since the quantity
of the fitting data is so limited, it is not possible to select an absolutely representative
 82 
distribution for each compound. Among the four chemicals for which data exist, phenol
had the most data so that the fitted distribution could be considered to be the most
repres.entati ve. Therefore, among the distributions thal fit the input data well, the one that
was similar to that of phenol was preferentially selected. Benzyl alcohol, phenol,
methylene chloride. and acetone aU followed the PEARSON5 (a, ~) distribution.
However, 2butanone followed the EXPON (~) distribution.
Tables 3, 5, 7, 9, and 11 present summaries of the statistics of the fitted distributions
for all compounds. It was found that the mean and standard deviation calculated from
Equations 33, 34, 35 and 36, are slightly different from those listed in these tables for
each compound. For example, the mean phenol concentration calculated from Equation
33 was 79.14 mgll while Table 3 lists a mean of 75.71 mg/I. The difference exists
because 79.14 mg!1 calculated from Equation 33 stands for the theoretical value of the
Pearson5 distribution with certain a, I~ values (n =1.3910; ~ =30.945), while 75.71 rngll
was obtained from samples being fitted. The theoretical equations assume that the
samples are continuous. However, as a matter of fact, samples being fitted are discrete.
When the bootstrapping results were used to confirm the fitted distribution, the
theoretical values of each statistics for each distribution were used.
5.2 Interpretations of probabilistic results
Simulations provided cumulative distributions for each compound. When USIng
TOXCHEM+V3, with 90% confidence interval (5% to 95%), the OCALC/IWTP was
estimated to release benzyl alcohol ranging from 0.048 IbId to 1.171 IbId (i .C., 5% and
95% values, respectively), phenol ranging from 0.070 IbId to 2.365 Jb/d, methylene
 83 .
chloride ranging from 1.087 IbId to 5.577 Ibid, 2butanone ranging from 0.061 IbId to
2.648 IbId, and acetone ranging from 0.54 IbId to 3.348 IbId. From WATER9, with a 90%
confidence rnnterval (from 5% to 95%), the OCALC/IWTP was estimated to release
benzyl alcohol ranging from o.om IbId to 0.009 IbId (i.e.. 5% and 95% values,
respectively), phenol ranging from 1.904 IbId to 66.199 IbId, methylene chloride ranging
from 0.995 IbId to 4.906 IbId, 2butanone ranging from 0.035 IbId to 1.695 IbId, and
acetone ranging from 0.317 IbId to 2.122 IbId.
Table 17 also presents the simulation results using average concentrations and flow
rate, which is the deterministic mode of the models. By using TOXCHEM+V3, the
estimated annual VOCs emission is 2366 lblyr (6.48 IbId), while WATER9 predicts an
annual VOCs (five compounds discussed in this thesis) emission of 7520 Ib/yT (20.60
IbId). [f the Henry's constants used in TOXCHEM+V3 were adjusted to the same as
those used in WATER9, the annual estimated emissions of VOCs would be 7966 lb/yr
(21.83 IbId). Card (1995) conducted a study on comparisons of mass transfer models
with direct measurement for free liquid surfaces at municipal wastewater plants and
found that WATER7 over estimated by a factor of fi ve and TOXCHEM over estimated
by a factor of two compared to the actually measured results. In Card's study it was
shown that the average measured emission rate for those IWTPs was about 70 pounds per
year per MOD of the influent flow. Compared to Card's (1995) results, TOXCHEM+V3
and WATER9 both predicted significantly larger emissions from the [WTP for this study,
which was 70008000 lb/yr with an average flow rate of IMGD. The huge difference
between these two estimations of VOC emissions might be due to the different industrial
processes in the plants of Card's study and this thesis. Different industrial processes
 84
would use different chemicals with different amount. Hence, more attention should be
paid to reduce the emissions from Tinker's IWTP and protect the residents from
exposure. However, these deterministic simulation results only represent a point on the
cumulative distribution curve of each compound. Information consisting of onepoint is
not enough to make a good decision relative to protective action. To make a better
decision, not only the emission rate is needed, but also the probability of the emission
rate is needed. The cumulative frequency of these values corresponding to the cumulative
distribution curve are shown in Tables 19 and 21. Comparing to probabilistic results, the
deterministic results could only tell there is about a 70% probability that the VOC
emissions rate from the IWTP would be or less as Tables 19 and 20 both indicate, or tell
that there is still 30% probability that the emission rate would be more. In this study for
the OCALC/IWTP, the total emissions simulated by WATER9 are significantly higher
than that by TOXCHEM+V3 if each model uses Henry's constants contained in its
library. Henry's constants for aU compounds in each model are different from each other.
These Henry's constants are listed in Table 20. If the Henry's constants in
TOXCHEM+V3 are adjusted to be the same as those used in WATER9, it was found that
both models have very similar estimations of emission rates. The simulation results are
presented in Table J8. The total estimated emission rate with TOXCHEM+V3 is 21.83
Ibid (3.98 tonlyr) while it is 20.60 Ibid (3.76 ton/yr) with WATER9. TOXCHEM+V3 has
slightly higher estimation of emission rate than WATER9.
Table 17 showed that the deterministic results predicted by WATER9 and
TOXCHEM+V3 are very different for benzyl alcohol and phenol. For benzyl alcohol, the
prediction derived from TOXCHEM+V3 is about 100 times greater than that predicted by
 85
WATER9. But for phenol, the simulation result from WATER9 is about 30 times greater
than that from TOXCHEM+V3. One of the reasons for these differences is the different
Henry's constant used in these two models, especially for phenol where there is four
orders of magnitude of difference in the dimensionless Henry's constant. These Henry's
constants are listed in Table 20. These Henry's constants are integrated in each model's
data library, with the exception of benzyl alcohol for TOXCHEM+V3. The data library in
TOXCHEM+V3 does not have Henry's constant for benzyl alcohol. The Henry's
constant was obtained from the Michigan Environmental Response Division website,
(website source, 2002). More work is needed to figure out why different Henry's
constants are used in WATER9 and TOXCHEM+V3. This area of research is partly
beyond the scope of this thesis.
When the Henry's constant for phenol in TOXCHEM+V3 was greatly increased after
it was adjusted to be the same as in WATER9, the emission rate increased. The emission
rate of phenol ranged from 1.079 ibId to 68.020 IbId with 90% confidence interval (from
5% to 95%) with TOXCHEM+V3. The prediction range of emission rates for henzyl
alcohol was changed to 0.012 Ibid to 0.314 Ibid with 90% confidence interval (from 5%
to 95%). While for the other compounds the emission rates were almost the same as those
using their own Henry's constants in the library of TOXCHEM+V3. For methylene
chloride, it ranged from 1.079 IbId to 5.519 Ibid, for 2butanone, it ranged from 0.060
IbId to 2.644 Ibid, and for acetone it ranged from 0.495 Ibid to 2.644 Ibid. The percentage
of the detenninistic estimation on the cumulative distribution curve is presented in Table
21. The simulated cumulative distribution curves are presented in Figures 53, 54, 55, 56
 86
and 57. The deterministic simulations were done by using the same Henry's constants in
both models. Furthermore, only the Henry's constants in the library of TOXCHEM+V3
were changed is because the database in WATER9 was not editable.
Due to different unit's definition in TOXCHEM+V3 and WATER9, the actual units
were modeled differently in each model. The emissions were different for each unit in
both models, though the total emissions for the whole plant were close. For example, the
diversionary structure ahead of the oil/water separator was modeled as an open dropstructure
in TOXCHEM+V3, while it was modeled as weir/waterfall in WATER9. Thus,
the emissions from both units were different. The emissions were 206 lb/yr from the open
dropstructure in TOXCHEM+V3 and 48.7 lb/yr from weir/waterfall in WATER9 with
the same model's setups for deterministic simulations. The detailed comparisons of these
units with different models are listed in Table 22, Figure 61, and Figure 62.
 87
co
00
Table 22. Comparisons of unit emissions in TOXCHEM+V3 and WATER9 (deterministic simulation with same Henry's constants)
Units TOXCHEM+V3 WATER 9
BenzYl alcohol Phenol Methvlene chloride 2Butanone Acetone Cpds Sum BenzYl alcohol Phenol Methylene chloride 2Butanone Acetone Cpds Sum
IbId Ibid IbId IbId IbId IbId IbId ibid IbId IbId IbId IbId
Primary Clarifier 0.01 2.68 0.47 0.16 0.24 3.560 0.003 2.000 0.400 0.152 0.362 2.917
lBlending Tank DI/D2 0.04 3.48 0.6 0.24 0.46 4.820 0.000 3.162 0.571 0.229 0.533 4.495
!oi versionarv Structure 1.23E05 0.26 0.02 7.19E04 5.1IE04 0.281 0.000 0.133 0.000 0.000 0.000 0.133
Pil Water Separators 8.64E04 0.26 0.04 1.41£02 0.02 0,335 0.000 9.905 0.914 0.076 0.038 10.933
Storage Tanks 2.50El1 1.08E06 4.55E08 1.39E09 9.28EIO 0.000 0.000 0.610 0.076 0.000 0.000 0.686
Equalization Basins 0.02 9.2 1.64 0.56 0.8 12.220 0.000 0.514 0.495 0.076 0.038 1.124
~ixing BasinI 2.97E04 0.Q3 4.95E03 3.34E03 6.23E03 0.045 0.000 0.000 0.038 0.019 0.000 0.057
~ixing Basin2 2.97£04 0.03 4.75E03 3.3IE03 6.2IE03 0.045 0.000 0.000 0.038 0.019 0.000 0.057
~ixing Basin3 5.26E04 0.06 0.01 8.48E03 0.01 0.089 0.000 0.019 0.038 0.019 0.000 0.076
Iweir fTom the Mixing Basin 6.07E05 0.06 5..53E03 1.29E03 1.84E03 0.069 0.000 0.000 0.000 0.000 0.000 0.000
~CCN 1.88E03 0.01 3.15E03 5.69E03 0.02 0.041 0.000 0.010 0.019 0.038 0.038 0.J05
Iwet Well 9.7 IEOS 2.43E04 5.35E05 1.76E04 9.26E04 0.001 0.000 0.000 0.000 0.000 0.000 0.000
~and Filter 0 0 0 0 0 0.000 0.000 0.000 0.000 0.000 0.000 0.000
rump staLion 3.36E05 0.03 2.86E03 7.08E04 1.02E03 0.035 0.000 0.000 0.000 0.000 0.000 0.000
/;::hlorine Contact Chamber 2.28E04 5.24E04 1.21E04 4. 12E04 2.17E·03 0.003 0.000 0.000 0.000 0.000 0.000 0.000
!oiversionarv Structure 1.23E·05 0.26 0.02 7. 19E04 5.IIE04 0.281 0.000 0.000 0.000 0.000 0.000 0.000
hickener 1.17E05 9.54E04 I 86E04 1.76E04 2.96E04 0.002 0.000 0.000 0.000 0.000 0.019 0.019

Belt filter 0 0 0 0 0 0.000 0.000 0.000 0.000 0.000 0.000 0.000
Sum 0.074 16.362 2.822 0.999 1.570 0.003 16.352 2.590 0.629 1.029
Total 21.827 20.603
Note: For blending tanks D l1D2, oil/water separators. storage tanks and equalization basins. the emissions are a combination of (WO parallel units.
When the Henry's constants were adju ted to be the same, TOXCHEM+V3 has higher
estimation of emission rates than that. of WATER9. Figure 61 shows the determini t.ic
results with same Henry's constants from WATER9.
iJ Benzyl alcohol
o 2Butanone
Phenol
Acetone
o Methylene chloride
25 ~,
~ 22.5 t==~~~~=====~;;~===1 "0 20 
:0
o 17.5 +j
~ ]5 +
I<
c;) 12.5 +t:
o 10 1
'U:)
c;) 7..5 1 ·s UJ 5 j
2.5 +0+
TOXCHEM+V3 WATER9
Figure 61. Prediction of VOC emissions with GFMs
From Table 22, it is obvious that the top four units that have the most emissions are
the equalization basins, blending tanks, primary clarifier, and oil water separators for
TOXCHEM+V3, while the sequence for WATER9 is the oilwater separators, blending
tanks, primary clarifier, and equalization basins. It is also clearly found that phenol is the
major VOC emitted from the Tinker'slWTP and mainly controls whi.ch unit is the major
emitter for WATER9 and TOXCHEM+V3. For TOXCHEM+V3, the equalization hasin
has the most emission of 9.2 Ibid while the oilwater separator has the most emissi,on of
9.9 Ibid for WATER9. More work is needed to figure out the reason for the big difference
in the future. A comparison of emission rates from major emitter for both models is
presented in Figure 62. The comparison is based on the deterministic simulation with the
same Henry's constants from WATER9.
 89
Primary clarifier
DOil water separators
Blending tank o Equalization basins
25 ,~
20 1~
=!2
.D
';; 15 11
c;:;
'c:
o
.~ 1011
'E
u.l
51
TOXCHEM+V3 WATER9
Figure 62. Estimated VOCs emissions by process.
During the Phase 2 of this project, field sampling (March 8, 20(2) was performed
by the Department of Civil and Environmental Engineering of Oklahoma State
University and Southwest Lab of Oklahoma. Three covered units (primary clarifier,
blending t.anks and oil water separator) and two uncovered units (storage tanks and
equalization basins) were sampled and the total emission rate was estimated. The results
of liquid phase samples analyzed were not received in time to be used in this thesis; only
the gas phase data was available. Moreover, benzyl alcohol and phenol was not
measured. The field measurements showed that the emissions of methylene chloride,
acetone and 2butanone were relatively low compared with the predicted emissions. The
predicted emission rate of all five compounds with TOXCHEM+V3 and WATER9 was
around 20 Ibid based on the historic liquid phase concentration data, while the actual
estimation of the emission rate was only 0.885 IbId. The field test results are listed in
Tabk 23. The compounds concentrations were analyzed and determined by GC/MS.
·90·
'0
Table 23. Emission rate of compounds based on field measurement
Sampling Compounds Emission
velocity Opening Flow concentration (Jlg/m
3
) rate (Ibid)
Covered Methylene Methylene
units at vents Area Rate chloride 2Butanone Acetone chloride 2Butanone Acetone
(fpm) (ft2) m3/d (llg/m1
) (llg/m3
) (1J.g/m3
) Ibid IbId Ibid
DI/02 50 7.8 15905 5630 47.85 594.3 0.197 0.002 0.021
Primary
clarifier 10 2.56 1044 7995 1104 54.1 0.002 0.00 0.00
Oil water
separator 15 2.56 1566 5590 24.5 217.2 0.019 0.00 0.001
SUM 0.219 0.002 0.021
Subtotal 0.242
Surface Flux Flow Compounds Emission
area of chamber concentration (ug/m3
) rate (Ibid)
Uncovered uncovered Xcross Rate Methylene Methylene
units unit area chloride 2Butanone Acetone chloride 2Butanone Acetone I
(£12) (£12) Umin (ug/m3
) (l!g/m3
) (ug/m3
) Ibid Ibid Ibid
Storage
tanks 4570 3.043 22 3,025.00 161.75 607.00 0.317 0.017 0.064
Equalization
basins 5005 3.043 22 1425 271.75 440 0.164 0.031 0.051
SUM 00481 0.048 0.114 I
Subtotal 0.643
Tala] 0.885
Note:
The calculation equations of emission rate are:
For Covered units, Emission rate (mass/d) :: Velocity (mid) * Opening Area (m2
) * Concentration (uglm3
)
For uncovered units, Emission rate (mass/d) =Flow rate (m3/d) * Concentration (uglm3
) * Surface area (fe) I Flux chamber area ere)
Chapter VI
Conclusion and suggestions
Uncertainty analysis of an emission inventory from an IWTP is a new issue, though
uncertainty analysis is widely applied in risk assessment. This thesis has presented a
method to perform this type of uncel1ainty analysis and its application in a case study.
Model uncertainty. parameter uncertainty and input uncertainty are three important
sources of uncertainty. Once the source of uncertainty is to be studied, the distribution of
the variables of interest becomes the most critical factor in the study. This thesis used
simulation techniques to fit the distribution for the variables of interest. Since the source
of uncertainty in this study was categorized as input uncertainty, all of the input variables
were discussed. The distributions for the influent flow rate and chemicals concentrations
were fitted and selected. The emissions rate of each compound was simulated using
TOXCHEM+V3 and WATER9. The simulation results presented a range of VOCs
emissions from the IWTP. When the same Henry's constants (from WATER9library) are
used, with WATER9, emission rate of benzyl alcohol ranged from 0.001 IbId to 0.009
IbId wilth 90% confidence interval (5%95%); phenol ranged from 1.904 IbId to 66.199
IbId; methylene chloride ranged from 0.995 IbId to 4.906 IbId; 2butanone ranged from
0.035 IbId to 1.695 IbId, and acetone ranged from 0.317 IbId to 2.122 IbId; with
TOXCHEM+V3, benzyl alcohol ranged from 0.012 IbId to 0.314 IbId, phenol ranged
from 1.979 IbId to 68.020; methylene chloride ranged from 1.079 IbId to 5.5 I9 IbId, 2butanone
ranged from 0.060 IbId to 2.644 IbId, and acetone ranged from 0.495 IbId to
2.644 IbId. The average emission rate of Tinker's IWTP was 21.83 IbId (3.98 tonlyr) with
 92 
TOXCHEM+V3 and 20.60 Ibid (3.76 tonlyr) with WATER9. However, there is only
about a 70% probability for both models to tell that the emission rate would be or less
based on these deterministic simulation results.
Uncertainty analysis of the emission inventory of the industrial wastewater treatment
plant can provide a range of VOC emission rates with confidence intervals. Monte Carlo
simulation was used to analyze the emission uncertainty. Based on the simulation results,
the emission rate and its probability can be easily read from the cumulative distribution.
This idea would be very useful for uncertainty analysis of environmental systems and
would be of great significance in environmental risk assessment.
The emission estimating models, TOXCHEM+V3 and WATER9, do not have a
simulation subroutine, i.e., Monte Carlo simulation; this makes it extremely difficult to
take other sources of unceltainty into consideration. More work is needed to solve this
problem, i.e. coding a simulation subroutine. In addition, more study is needed in finding
other possible reasons for the differences between these two estimating models, with the
exception of Henry's constants, such as the fundamentals of how the individual process
units are modeled in both models. FUfther exploration on these differences would be very
useful in the application of these two models.
Due to the limited time of study, the liquid phase concentrations of compounds of
interest sampling on March 8, 2002 are not analyzed yet, and thus those concentrations
are not substituted into models for deterministic simulation. The deterministic simulation
results with the field samples liquid phase cOll1centrations Inight provide useful
information for comparisons between field measurement and models estimates.
BlBLIOGRAPHY
1. Burrowes, P.A. et al. (1995). "Emission Assessment at a Large Canadian Wastewater
Treatment FacilityWhat are the Impacts?". Proceedings. Air and Waste
Management Association, San Antonio, Texas.
2. Caloghirou, Y.D. et 811. (1996). "Macroeconomic Impacts of Natural Gas Introduction
in Greece". Journal of Energy, 21(10), 899909.
3. Card, T.R. (1995). "Compm1son of Mass Transfer Models with Direct Measurement
for Free Liquid Surfaces at Wastewater Treatment Facilities". Proceedings.
Air and Waste Management Association, San Antonio, Texas.
4. Chapra, S. and Canale, R. (2002). Numerical Methods for Engineers, Fourth Edition.
McGrawHill.
5. Corsi, R.L. and Olson D.A. (1998). 'EmiSSIon Models' in Odor and VOC Control
Handbook. Ed. Rafson, Harold J. McGrawHill. New York, NY, ppS.ll5.25.
6. Curto, M.E and Daly, M.H. (1995)