United States
Environmental Protection
Agency
Environmental Sciences Research EPA-600/4-79-035
Laboratory May 1979
Research Triangle Park NC 27711
Research and Development
vvEPA
Systematic !;
Sensitivity
Analysis of Air -
Quality Simulation
Models
-------
RESEARCH REPORTING SERIES
Research reports of the Office of Research and Development, U S Environmental
Protection Agency, have been grouped into nine series These nine broad cate-
gories were established to facilitate further development and application of en-
vironmental technology Elimination of traditional grouping was consciously
planned to foster technology transfer and a maximum interface in related fields.
The nine series are'
1 Environmental Health Effects Research
2 Environmental Protection Technology
3 Ecological Research
4 Environmental Monitoring
5 Socioeconomic Environmental Studies
6 Scientific and Technical Assessment Reports (STAR)
7 Interagency Energy-Environment Research and Development
8. "Special" Reports
9. Miscellaneous Reports
This report has been assigned to the ENVIRONMENTAL MONITORING series.
This series describes research conducted to develop new or improved methods
and instrumentation for the identification and quantification of environmental
pollutants at the lowest conceivably significant concentrations. It also includes
studies to determine the ambient concentrations of pollutants m the environment
and/or the variance of pollutants as a function of time or meteorological factors.
This document is available to the public through the National Technical Informa-
tion Service, Springfield, Virginia 22161.
-------
EPA-6QO/4-79-Q35
Way 1979
SYSTEMATIC SENSITIVITY ANALYSIS
OF AIR QUALITY SIMULATION MODELS
Robtert J, Geltnas and J. Peter Vajk-
Science Applications, Inc,
1811 Santa Rita Road, Suite 104
Pleasanton, California 94566
Contract No. 68-02-2942
Project Officer
Robert E, Eskridge
Meteorology and Assessment Division
Environmental Sciences Research Laboratory
Research Triangle Park, North Carolina 27711
ENVIRONMENTAL SCIENCES RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
RESEARCH TRIANGLE PARK, NORTH CAROLINA 27711
-------
DISCLAIMER
This report has been reviewed by the Environmental Science Research
Laboratory, U.S. Environmental Protection Agency, and approved for
publication. Approval does not signify that the contents necessarily
reflect the views and policies of the U.S. Environmental Protection
Agency, nor does mention of trade names or commercial products constitute
endorsement or recommendation for use.
11
-------
ABSTRACT
This report reviews and assesses systematic sensitivity and
uncertainty analysis methods for applications to air quality
simulation models. The discussion of the candidate methods
presents their basic variables, mathematical foundations, user
motivations and preferences, computer implementation properties,
and costs (including both human and computer resources). Both
deterministic and sampling methods have been evaluated with
illustrative examples of "What you pay"and "What you get" with
present-generation systematic sensitivity analysis methods and
air quality models. Deterministic methods include the time- and
user-honored method of arbitrary parameter adjustments by trial
and error, variational methods, and a newly formulated Green's
function method (for kinetics systems only). Sampling methods
include Monte Carlo sampling of the outputs of air quality
models to compute variances and a Fourier analysis method of
sampling model outputs to compute expectation values of sensitiv-
ity coefficients.
Computational economics, inclusive of both programming
effort and computer execution costs, are the dominant governing
factors for the effective application of systematic sensitivity
and uncertainty analysis methods to arbitrary air quality
problem scenarios; several reasonable options and several un-
reasonable options emerge from the present evaluations. Recommen-
dations outline how the EPA should, today, most effectively
apply available multi-parameter systematic methods of sensitivity
analysis to arbitrary 3-D transient (PDE) air quality simulation
models. The report concludes with a discussion of key break-
throughs which would lead to the next major advances in system-
atic sensitivity analysis. In certain instances, these would
simultaneously provide major advances in air quality simulation
modeli ng , per se.
This report was submitted in fulfillment of Contract
No. 68-02-2942 by Science Applications, Inc., under the sponsor-
ship of the U.S. Environmental Protection Agency. This report
covers a period from May 12, 1978, to November 12, 1978, and
work was completed as of December 4, 1978.
-------
CONTENTS
Abstract iii
Figures vi
Tables vii
Acknowledgments viii
1. Introduction 1
2. The Problem 6
Air Quality Models and Their Sources of Uncertainty . . 6
Motivations for Sensitivity Analysis 12
Analytical Basis of Sensitivity Analysis Methods .... 14
3. Systematic Sensitivity Analysis Methods 30
Deterministic Methods 30
Sampling Methods 48
4. Evaluation of Candidate Sensitivity Methods for Air
Quality Model Applications 61
Basic Cases 63
Chemical Model Evaluation 66
1-D Time-Dependent Model Evaluation 76
3-D Air Quality Model Evaluation 82
5. Today's Approach 85
Conclusions 85
Today's Practical Application of Systematic
Sensitivity Analysis 87
Needed Breakthroughs in Systematic Sensitivity
Analysis Methods 92
References 95
-------
FIGURES
Number Page
1 Effects of parameter variations on the solutions 19
2 Schematic representation of solution surfaces 22
3 Interpretation of probability distributions of
uncertain parameters 24
4 Construction of first-order sensitivity measures and
coefficients 26
5 Schematic illustration of the Monte Carlo method 49
6a Schematic illustration of the Fourier Amplitude
Sensitivity Test (FAST) 52
6b Functional output of FAST 54
7 Arbitrary nature of functional forms used in FAST 56
vi
-------
TABLES
Number Page
1 Characteristics of candidate systematic sensitivity
and uncertainty analysis methods 60
2 Base cases for illustration of execution costs of
candidate sensitivity analysis methods 83
3 Illustrative computer execution costs for candidate
systematic sensitivity analysis methods 84
vii
-------
ACKNOWLEDGEMENTS
It is a pleasure to thank Dr. Richard Stolarski and
Prof. John Seinfeld for in-depth discussions of their
sensitivity work prior to publication. We are also
extremely grateful to our colleague, Dr. Said Doss, for
his exacting-review and suggestions regarding some of
the complex issues in numerical analysis which were
encountered in the course of this study.
VI 1 1
-------
IN
ECTION 1
RODUCTION
Uncertain parameters and sparse field data are
intrinsic elements of air quality modeling. This state
of affairs will persist indefinitely, despite the nation's
ever more urgent need for more detailed understanding of
cause and effect mechanisms in air quality. Thus a key
question is, "To what extent do uncertain input factors in
air quality models lead to impaired accuracy and precision
of the resultant outputs of these models?" This basic
question and many of its variants are examined in this final
report on "Systematic Sensitivity Analysis for Air Quality
Simulation Models ".
The principal objectives of this study were:
(1) to identify those systematic sensitivity
analysis methods which are applicable to
air quality simulation models;
(2) to evaluate critically the potential of
candidate methods of systematic sensitivity
analysis for effective application to key
air quality issues;
(3) to identify those breakthroughs most critically
needed for each candidate method; and
(4) to determine the most effective immediate
means of conducting sensitivity analyses on
existing air quality simulation models,
assuming no further developments in sensitivity
analysis methods.
-------
It is frequently unrecognized that sensitivity analysis
is a basic ingredient in any application of the classic
scientific method in which alternative hypotheses--including
models and input data--are tested and retested against
selected physical realizations. Historically, this process
has most frequently been conducted on a trial and error
basis. In the development of air quality models, sensitivity
testing is a particularly laborious process if all potentially
significant input parameters are to be varied, individually
or simultaneously, for every potentially significant spatial
1ocati on.
Somewhat more systematic approaches to sensitivity analysis,
including some semi-automated methods, have appeared recently
in the study of purely chemical kinetics systems. Enormous
benefits would undoubtedly be realized in air qual-ity analysis
and policy making if semi-automatic sensitivity methods could
be extended to economic applications in arbitrary space- and
time-dependent multivariate models. Such extensions of
systematic sensitivity analysis thus represent both a challenge
and a promise for the entire field of air quality work. The
principal focus of this study is the evaluation of this duality
--the challenge versus the pnomise--in the application of
candidate sensitivity analysis methods to arbitrary air quality
models.
Perhaps the key consideration for air quality applications
is the computational economics of candidate sensitivity
methods. On the one hand, we shall see that several dedicated
methods would require that a specific sensitivity program be
written for each specific application to an air quality model.
On the other hand, we shall see that several non-dedicated
methods can be applied, without additional programming effort
for either the air quality model or the general sensitivity
2
-------
method, to any air quality model. All sensitivity
analysis methods considered--in their current state
of development--wi11 be seen to have both strengths
and weaknesses in regard to potential applications
to 1-D, 2-D, and 3-D time-dependent air quality models.
The extent to which the weaknesses will be disabling
and/ or the extent to which the strengths can be ex-
ploited will necessarily depend upon the needs and
expectations of potential users of these systematic
methods. The motivations of air quality policy makers
and of air quality researchers, while closely coupled,
are sufficiently distinct that they may be served best
by different methods of viewing uncertainties and sen-
si ti vi ti es.
In the present study, it has been necessary to em-
phasize issues relating to the sensitivities of air quality
models With respect to those physical, chemical, or meteorological
factors which are most readily quantifiable. Thus we have
considered variations in model outputs such as the concen-
trations of Oj. NO, N02» and other air-borne species which
result from uncertain and variable input factors such as
photochemical rates, emission sources, pollutant sinks, and
winds. We have not considered sensitivities to computational
factors since no practical standards exist for comparative
analysis. It would be impossible, for example, to determine
the dispersion in oxidant concentrations attributable to
numerical diffusion and/or to operator-splitting in the
Systems Applications, Inc., regional model^ ' vis-a-vis
the LIRAQ model' ' vis a vis the MADCAP model^ for a specific
"typical day" in, say, St. Louis, even if the same data bases,
initial conditions, and boundary conditions were incorporated
-------
in each model. The multitude of other uncontrolled differ-
ences among air quality models precludes meaningful con-
sideration of model-to-model sensitivities at this time.
We have thus surveyed and evaluated those sensitivity
analysis methods which can, for any given air quality model,
determine systematically and quantitatively the dispersion
in the output variables--such as space- and time-dependent
concentrations of ozone, PAN's, NO , acids, aldehydes,
J\
peroxides, and SO -- which result from uncertain photo-
X
chemical rate data, winds, diffusivitives, initial and boundary
values, and emission sources and sinks, and from deviations
of actual conditions from "typical day" conditions.
Four types of systematic sensitivity methods are discussed
in this report:
(1) "time-honored methods," which examine variations
of selected parameters, either individually or
simultaneously, by trial and error;
(2) variational methods (including direct methods);
(3) Fourier amplitude sensitivity test (FAST) methods;
(4) Monte Carlo methods.
The report is presented in four sections. In Section
2, The Problem, we discuss the motivations, logical bases,
and key issues involved in the application of systematic
sensitivity analyses to air quality simulation models. In
Section 3, Systematic Sensitivity Analysis Methods, we de-
scribe the candidate methods in terms of their basic variables,
mathematical formulation, computer implementation, and execu-
tion properties.
-------
In Section 4, Evaluation of Candidate Sensitivity
Methods for Air Quality Model Applications, we discuss
"what you pay" , "what you get", areas of promise, and
problems. Section 5, Today's Approach, indicates the most
effective means, in our opinion, of applying systematic
sensitivity methods to air quality models in the face
of the current adverse constraints of computer economics.
Finally, we indicate some of the needed breakthroughs
which would lead to dramatic advances in the practical
utility of candidate sensitivity methods.
-------
SECTION 2
THE PROBLEM
AIR QUALITY MODELS AND THEIR SOURCES OF UNCERTAINTY
Air quality simulation models attempt to describe the
physical, chemical, and photochemical processes occurring
in specific portions of the lower atmosphere. Variables
of principal interest include the concentrations of sus-
pended particulate matter and of such gas phase chemical
species as ozone, nitrogen oxides, and carbon monoxide,
for which regulatory standards have been set by various
governmental jurisdictions. Although no such regulatory
standards have been set for species such as aldehydes,
peroxides, nitrates, and acids, these and other species
are also the subject of extensive current interest. The most
thoroughly developed and calibrated air quality models
available today describe, for the most part, the photo*
chemical kinetics and the transport processes for gas
phase species in certain limited atmospheric regions.
Several dissimilar sources of uncertainty in our know-
ledge of atmospheric processes affect the predictions of
air quality simulations. Such uncertainty factors include:
(i) Natural Variabilities
winds t temperature
0 solar flux atmospheric turbulence
e humidity atmospheric stability
natural sources and particulate matter
sinks of ai r-borne
chemical species
-------
It is important to note here that significant
natural variabilities occur in the real atmosphere
over scales which are often disparate with respect
to the scales used in practical computations. For
tropospheric models, typical computational scales
are tens of meters to tens or hundreds of kilometers.
Mixing of chemical species, on the other hand, occurs
on microscopic scales comparable to molecular mean
free path lengths. Small scale eddies and chemical
mixing are thus, in most instances, treated by sub-
grid scale characterizations rather than by detailed
calculations of reactive turbulence. Winds at specific
grid points in air quality models are generally inter-
polated from discrete field data from monitoring stations
separated far wider in space and in time than the horizontal
or vertical grids or integration time-steps of the models.
Variable attenuation and multiple scattering
effects associated with particulate matter, gaseous
components, cloud cover, and surface albedo all contribute
to significant local variations in solar flux. Coastal
areas, in particular, show extreme variations within a
single geographic area and within a single day.
Naturally occurring chemical sources include forest
fires, ocean spray, dust, local geogenic emissions from
oilfields and coalfields, eutrophication of inland waters
and oceanic processes releasing hydrocarbons, N^O, CO,
NH3, and H2S, as well as a nearly universal background
release of methane. To these are added widespread
sources attributable to agriculture and deforestation
around the globe. Natural chemical sinks include hetero-
-------
geneous reactions with suspended participate matter,
rainout, absorption at the air/ocean interface, and
other surface deposition processes.
(i i ) Photochemical Factors
0 upper and lower atmospheric compositions
(are determinants of solar flux intensity
and spectral distribution)
e diurnal variations in solar flux
o absorption and photodissociation cross-
section data
o hydrocarbon emission inventories
o photochemical characterizations of aggregate
hydrocarbon emissions*
Photodissociation rates are described by wave-
length integrations over products of local solar flux
and photodissociation cross-sections. These rates thus
depend strongly on the spectrum of the local solar flux
which is determined by atmospheric composition and time
(6)
of dayv . Because the concentrations of minor reactive
species--such as OH, 0, HUO^, NO-,, and N20g--have seldom
been measured simultaneously with solar flux and the con-
centrations of other important atmospheric species, much
of the validation of photochemical mechanisms remains to
be done.
Atmospheric hydrocarbon inventories, of course,
consist of huge numbers of individual organic compounds
which must be aggregated according to their reactive
properties in a photochemical mechanism. Additionally,
*See, for example, the articles of Hecht, Seinfeld, and
DodgeU), Hogo and Whitten (5), and Gelinas and Skewes-Cox(6)
for detailed characterizations and categorizations of reactive
hydrocarbons.
-------
detailed emission inventories and local hydrocarbon
distributions are difficult and costly to obtain for
regions as extended as those spanned by air quality
simulation models.
(i i i) Field Data
Measurements of such meteorological variables
as wind, temperature, stability, and transported
chemical species are relatively sparse in both time
and space. Because these quantities are associated
with exceedingly complex fluid dynamics processes
they have usually not been calculated from fundamental
equations but have instead been imposed upon air quality
simulation models using interpolations of the measured
data onto the spatial grids and time points of the
model. This imposition introduces uncertainties which
are associated with the specific interpolation processes
as well as uncertainties which are associated with non-
representative measuring stations and with the measure-
ment methods per se. An additional form of uncertainty
arises because meteorological conditions are necessarily
categorized to represent conditions which occur on
selected "typical days". Actual meteorological conditions
on any specific day will naturally differ from these
"typical day" conditions, thus constituting yet another
possible source of dispersion in the outputs of air
quality simulation models.
-------
(i v) Computational Factors
Air quality simulation models are formulated from
non-linearly coupled partial differential equations in
space and time describing physical transport and chemical
kinetics processes. Very simple models are sometimes
derived from systems of non-linearly coupled ordinary
differential equations in time only. In either case,
numerical integration methods are the most effective
means for obtaining solutions of the basic modeling
equations. Such numerical solutions are inherently
approximate, so that simulation model results inevitably
include dispersion due to these computational (numerical)
approximations.
Computational dispersion usually arises from
the well-known dilemmas of numerical modeling, most
notably numerical diffusion (in Eulerian formulations)
and aliasing of individual species under rezoning (in
Lagrangian formulations). Such computational features
as choice of differencing methods, representations of
sub-grid scale physics, simplifications in chemical
kinetics mechanisms, and imposed profiles for pollutant
sources result in additional computational dispersions
in the model outputs.
10
-------
Model descriptions of the processes mentioned above
are incorporated in systems of partial differential equations
(PDE's) which are generally solved by numerical integration
on large digital computers. In the present study, we con-
sidered air quality simulation models based upon PDE's of the
form
= R.(c1,c2,...,cli] + S.(x,t),
3t i = 1,...N, (2.1)
where c- = c-(x_,t) represents the mean concentration of the
i-th chemical species; _y_ is the mean fluid velocity; Kn is
a diffusion tensor; R. is the net sum of all chemical reaction
rates for production or destruction of the i-th species;
and S-(x_,t) represents the net sum of all sources and sinks
for the i-th species. This equation is a reduced version of
more fundamental (microscopically exact) fluid dynamics
equations, and thus incorporates well-known approximations
and assumptions which are peripheral to the present discussion
The solutions c-(x,{p},t) of Eq. (2.1) are functions of
I
space and time (_x, t), but they are also functions of the
parameters { p} = {y_, KD, ka, Jg, £ , c. ^t. g}. The parameters
k^ (T) represent temperature-dependent rate constants deter-
mining the thermal chemical reaction rates; J0 represents
P
coefficients for the photochemical reactions included in the
terms R-; and A represents parameters in the expressions for
source and sink rates included in the terms S-. Since the
solutions for the species concentrations are functionally
dependent on the initial and boundary values, c. and c.
'»l '»o ,
respectively, it is often useful to include these values
in the parameter array {p} as well.
11
-------
The values for this host of parameters needed to
describe important atmospheric processes are obtained from
laboratory experiments, field measurements, and theoretical
or empirical prescriptions. The source of the entire problem
under discussion in this study is the inevitable fact that
all of these parameter evaluations are uncertain to a degree.
In the face of such uncertainties, it is natural to ask,
"For which of the parameters do the uncertainties in the
parameters produce the greatest dispersion in the model out-
puts (especially for the most important atmospheric observables)?"
MOTIVATIONS FOR SENSITIVITY ANALYSIS
Air Q u a 1 i ty Policy
Air pollution control agencies and other governmental
bodies or officers are charged with the responsibility of
regulating pollutant emissions in order to meet or to main-
tain standards of air quality established for certain atmos-
pheric regions. To the extent that models are used responsibly
in decision making, it is both fair and necessary to ask,
"What are appropriate error bars for the outputs of the
model?" This question is appropriate whether the models used
are analytically based, detailed simulation models incorporating
as much of the current scientific knowledge as is practical
or the models are qualitative mental models which incorporate
intuition. In asking this question, we would like to obtain
a statement of the un certainty--or, equivalently, the credibi-
lityof the model results which are used in regulatory decisions.
The principal emphasis is then given to determining the expected
uncertainties in important atmospheric variables, given the
manifold of uncertainties present in the input parameters of
12
-------
the model. (Analytic definitions will appear in later
sections of this report.) As an example, one
might want to know the expected uncertainty in ozone
concentrations predicted by a given simulation model at
all significant geographical locations and times resulting
from the combined effects of rate data uncertainties,
meteorological uncertainties, emission uncertainties, etc.,
in order to assess whether proposed regulations are reason-
able and proper.
Air Quality Research
In contrast to regulatory policy applications of air
quality simulation models in which (expected) uncertainty
analysis or credibility analysis is emphasized, research
applications tend to emphasize sensitivity--cause-and-effeet-
analysis. The two concepts are, to a significant degree,
complementary. In research applications, the scientist seeks
to determine the magnitude of the effects which each and
every parameter input to the model produces in the outputs
of that model. This is a more explicitly detailed approach
than that of expected uncertainty analysis because, in any specific
simulation calculation, the model must use a single, specific
numerical value for each parameter of the model, not a dis-
tribution of possible values. Best case and worst case air
quality estimates are very conveniently made in terms of indi-
vidual, unaveraged parameter variations.
In studying the effects of varying an individual para-
meter (or a limited set of parameters) on the outputs of the
model, all other model parameters are generally held fixed,
13
-------
either at nominal "best" values or at "test values" selected
by the researcher. (It is also possible to consider simultane-
ous variations in all other parameters in order to obtain
best or worst case estimates.) Sensitivity analysis is thus
used by researchers as a means for identification of cause-
and-effect mechanisms and for ranking various mechanisms in
order of importance in the overall system of processes deter-
mining the state of a specific atmospheric region. The goals
of this approach are: (i) to increase basic scientific under-
standing of atmospheric processes; (ii) to focus scientific
research on the most critical issues, ultimately reducing un-
certainties in existing models; and (iii) to assist regulatory
applications, both directly and indirectly, by progress toward
the first two goals
ANALYTICAL BASIS OF SENSITIVITY ANALYSIS METHODS
The methods and language of statistics have been extensively
adopted in the development of systematic sensitivity methods,
undoubtedly because interest centers on deviations in model
solutions which result from deviations in governing input
parameters. In terms of air quality modeling equations which
are represented by Eq. (2.1), a major question is, "What
variations occur in the solutions, c-(x, (p},t), relative to
1
c.j (_x,{ p°), t), when model input parameters, {p}, are subjected
to variations relative to their nominal values {p°}?" The
input and internal parameters for Eq. (2.1) are written most
generally in terms of their nominal values and deviations as:
14
-------
r = r ° + >
ci,B = ci,B + «ci)B
v = v° + 6v
k = k + 6k
a a a
(2.2)
At this stage nominal values, {p°>, are not necessarily
mean values, nor are deviations, {6p}, necessarily standard
or random deviations. Considerable care must thus be exer-
cised to distinguish nominal values from mean,values or
central values, and also to distinguish arbitrary deviations
from standard deviations, variances, and random variations.
The logical question of, "How are nominal values to
be assigned in air quality models?", has several answers.
In principle, measurable parameters such as reaction rate
data, initial and boundary values, winds, diffusivities , etc
can be represented ultimately by mean values and associated
standard deviations; it is frequently true, in practice,
that available data bases do not allow reliable statistical
treatment. For example, in photochemical kinetics some
limited number of rate coefficients have been repeatedly
measured by independent experiments which have been well-
prepared (in the statistical mechanical sense) and thus
15
-------
have reliable statistical properties. But this tends to be
the exception rather than the rule for many reactions
which are objects of current research. Indeed, we can not
now profess to know all of the potentially important
atmospheric reaction processes, much less their mean rates
and statistical distributions. In such cases, nominal
values of known and hypothetical reaction rates are assigned
on the basis of each individual researcher's discretion.
Parametric deviations are also arbitrarily assigned for the
purpose of testing the sensitivity of photochemical mechan-
isms to specific reactions. The history of atmospheric
photochemistry developments emphatically illustrates that
meaningful probabi1i ty distributions of poorly known reac-
tion rates are simply not available at times of greatest
need. (Poorly known chemical rates have always been subject
to very large, frequently unpredictable changes which emerge
when markedly new or more carefully prepared experiments are
performed.) This discussion of meaningful probability distri-
butions for physical parameters is particularly pertinent to
our consideration of those sensitivity and uncertainty analysis
methods which use probability distributions.
Measures of Sensitivity and Uncertainty
Several measures are used in sensitivity and uncer-
tanity analysis which should be distinguished at the outset.
The most elementary measure is a deviation (variation)
in state variables such as species concentrations at specific
coordinates (_.x_K,t' ) ,
fic^xKpht1) = c.(x',{p0 + 6p},t') - c.(x;{p°}. ,t'), (2.3)
in the neighborhood of nominal parameter values, {p°}.
(Recall that the parameters can be either random or non-
random vari ables. )
16
-------
By Taylor's theorem
M
c(x',{p° + 6p},t) = c.(x',{p°}t) + Z
i k= 1
0((max6pk)2), -j = i,N; k = ijM, (2.4)
where the sensitivity coefficients are defined by
Zi'k = 8ci/9pk > 1 = l.N; k = 1,M .
Sensitivity methods which neglect the terms 0((max dp^) )
are referred to as first-order or linear methods. In space
and time dependent models, which are described by partial
differential equations (in contrast to ordinary differential
equations for purely kinetics models or for steady state
purely fluid dynamics models), species concentration varia-
tions associated with the k- parameter variations are
given by the functional derivatives,
p 9c, (x.1 ,t' ) .,
6^(21',t1) = I 3P /x tt \ <5pkU,t') djx (2.6)
and
6c,(x',t') =| ^ 6pk(x',t) dt. (2.7)
~ J 3pk(x',t) k
The sensitivity variable in Eq. (2.6) represents variations
in c- at (_x'»t') with respect to variations in p^ at locations
1 K
x and time t1. Likewise, the sensitivity variable in Eq. (2.7)
represents variations in c^ at (x',t') with respect to variations
in p. at times t and location x/. These functional derivatives
and the sensitivity coefficients, Z. ., form the basis of
I 5 K
17
-------
variational methods which will be described in Section 3.
These variables provide convenient measures of sensitivities
and uncertainties when probability distributions are poorly
defined and for "best case" and "worst case" air quality
considerations, which were discussed above in connection
with AIR QUALITY RESEARCH MOTIVATIONS.
The effects of parameter variations on the solution
c-(_x,t) are shown schematically in Figure 1. Two cases
must be distinguished here. In Figure l(a), the solution
for the "best value" of all parameters is shown by the solid
line. If some one of the parameters, say p., were changed
from its "best value" p°. to p*? + <$p-, the solution curve
J J J
would be different, although it would start from the same
initial value, c- T, for c.. Likewise, if the parameter
1,1 i
were changed to pQ -<$p., the solution curve would be different
J J
again. These altered solutions are shown by the dotted lines
in Figure l(a). If, on the other hand, we wanted to investi-
gate the effects of changing the initial value c. T itself,
i j i
the results might look like the dotted line curves shown in
Figure l(b).
It is important to realize that it is not enough to
consider the variations in the solutions ci at a single time
t1; we must consider the variations themselves as functions
of time if we are to properly evaluate the sensitivity of the
system to parameteric variations or uncertainties.
If the set of coupled differential equations describing
the system is non-linear, not only the value of c. at time
t may in general be altered, but the shape of the solution
curve itself may change. In certain extreme cases, changing
some parameters p (s f j) by small amounts may result in
18
-------
t
^.,..
x|
o
Figure 1 (a)
The effect on the solution for c^ as a function of time
when one of the parameters p- is altered by +_ 6p-, where
Pj is one of the rate coefficients or boundary conditions.
( For example, c^ may represent the concentration of ozone
and p. may represent a local photochemical rate coefficient.)
The parameter p. may also represent the initial value of a
different concentration c«,
x|
i,I=Pk
A\
'" "c1(x,pj+6pk,t)
Figure 1( b).
The effect on the solution for c- as a function of time
when the parameter p, representing the initial value
of c.j is altered by HH 6p. . Because of the non-linearity
of the system of coupled equations, the solution curves
are not the same shape as the solution for the "best value"
with a constant displacement. Changing the initial value
can alter the qualitative shape of the concentration-versus-
time curve significantly.
19
-------
changing the solution curve from a smooth, gradually
varying function of time into a very rapidly oscillating
curve of small amplitude which follows generally the
trajectory of the1 unperturbed solution, in such cases,
the magnitude of the variation 6c.].(_x,t) may always be
small, but we would certainly want to describe the system
as highly sensitive to variations in p .
Other measures are obtained by introducing statistical
averages taken over distributions of parameter values.
For example, mean species concentrations, » are
I
defined by
E /.../ ci(2clp1%,pM; t-)P (P1,...,PM^dpr ..dpM }
i = 1,N, (2.8)
where P (PJ >P2 > »PM) represents the probability distribution
for the parameters. Uncertainties in species concentrations
are often expressed in terms of variances, a^ ()(_', t'), defined by
a^x'.f) E< [cjCx'.t1) - ] 2> (2.9a)
or
Oi(x' ,t') = - 2 i = 1,N, (2.9b)
where
.(xI,t')> = /.../ c.(x';Plj.,.p ,t')
2
i = 1, N (2.10)
Variances provide convenient measures of expected uncertainties,
which were discussed above in connection with AIR QUALITY
POLICY motivations. As for expectation values of sensitivity
coefficients, sampling methods have been developed to evaluate
measures of , which will be described in Section 3.
I , K
20
-------
GRAPHIC INTERPRETATION
A graphic interpretation of the preceding discussion
may clarify the meanings of, and differences among, the
various measures of sensitivity and uncertainty. For
simplicity, suppose that the set of coupled differential
equations (2.1) depends on only two uncertain parameters,
pi and PO» whose values lie between p, and p, and
between Poln and p^*, respectively. For any given position
and time (x,t), the solution c- can be considered to be
I
a function of p, and p2 as shown in Figure 2. The "best
values" of p-^ and p~ are, respectively, p° and p2; in
general, these are not in the exact centers of the respective
domains of uncertainty.
The point Q on the solution surface represents the
magnitude of the solution c- at time (x.»t) computed with
both uncertain parameters assumed to have their "best values.'
If p, were then varied over its domain of uncertainty,
holding p2 fixed at its "best value," the solution point
would trace out the curve DQD1 ., Conversely, if P2 were
varied over its domain of uncertainty, holding p, fixed at
p°, the solution point would trace the curve AQA1. Varying
both of the parameters over the full domain of uncertainties
in both then generates a two-dimensional solution surface
for each variable c- for each position _x. This solution
surface changes as t itself is allowed to vary.
In the general case, far more than two uncertain para-
meters are involved in an air quality model. If M parameters
are uncertain, the solutions trace out M-dimensional hyper-
surfaces; the features of the analysis discussed below apply
equally well to the M-dimensional case as to the two-
dimensional case.
21
-------
Range of
uncertainty
Domain of uncertain
in the parameters
Figure 2. Schematic 'representation of the solution surface
c- over the domain of uncertainty in the parameters.
Tne resulting range of uncertainty in the value of
the solution c^ about the nominal value corresponding
to the "best values" of the parameters is also
indicated.
22
-------
The solution surfaces c-({p}) may have very complicated
topographies, with peaks, pits, undulations, sharp ridges
or gorges, saddles, and even such discontinuities as sheer
cliffs. The purpose of sensitivity analysis, from this
point of view, is to understand the shape of the solution
surfaces as we vary the uncertain parameters away from the
immediate vicinity of the current "best values." As new
field data or new experimental evaluations of parameters
become available, our estimation of the "best values" may
change, moving us across the complex topography of the
solution surface. Understanding the shapes of the solution
surfaces provides advance warning about what kinds of changes
in the parameters will most drastically change the predictions
of our simulation models.
As we mentioned earlier, for policy-making purposes
uncertainty analysis is generally of greater interest than
sensitivity analysis. The purpose of uncertainty analysis
is to assay the range of uncertainty in the solutions c-
resulting from the present degree of uncertainty in the
parameters, as indicated in Figure 2. A crude measure of
this uncertainty is the difference in elevations between the
nominal case (point Q) and the extremes (points B and C').
Although the extreme values of c^ are shown in the figure
to lie at corners of the domain of uncertainty in the para-
meters, the extreme values will more commonly lie somewhere
inside that domain.
If we have some idea about the probability distributions
for the uncertain parameters, we can in principle compute
the probability distribution for c. over the domain of uncer-
tainty in the parameters, as shown in Figure 3. From the
probability distribution for c^ , quantities such as mean
species concentrations and variances a.= can be computed
23
-------
Ci(x,{p},t)
D1 C1 max
-c-r- -
Figure 3. Given assumed probability distributions for
each uncertain parameter, the probability
distribution for the solution c- is obtained
by computing a probability-weighted distribution
of elevations of the solution surface. Note
that the "best value" of a parameter may differ,
in general, from both the most likely value
and the mean value of the parameter. Likewise,
> the mean value of the solution c^,
may differ from the most likely value of the
solution and from the nominal value c.(Q)
corresponding to the "best values" of the
parameters. We can never be absolutely certain
that the true value of a parameter lies within
a given domain of uncertainty. Thus P(p-i) and
P(p«), in general, have small but finite values
beyond the domain shown. In this illustration,
the uncertainties in p, and p,-, are assumed to
be independent and uncorrelated.
24
-------
by standard statistical methods as defined in Eq's (2.8),
2.9) , and (2.10).
Several important points should be noted in Figure 3.
First, the investigator's estimate of the "best
value" of a parameter may or may not coincide with either
the most likely value or the expected (mean) value of the
parameter as computed from the probability distribution.
(Note, in particular, p, in the figure.) Second, the shape
of the solution surface c- is independent of all assumptions
about the shapes of the probability distributions for the
parameters; the shape of the solution surface depends only
on the structure of the set of coupled differential equations
(2.1). Third, the probability distribution for the value
of the solution c^ (.x.t) depends on both the assumed probability
distributions for the parameters and the shape of the solution
surface. Thus the mean value may or may not coincide
with either the most likely value of c^ or the nominal value
of c- corresponding to the "best values" of the parameters.
Fourth, the variance o. of the probability distribution for
the solution c- will in most cases be significantly smaller
than the range of uncertainty in GI , that is, 2ai<_(c!jiax-crjn n).
The significance of Eqs. (2.3), (2.4), and (2.5) can be
seen in Figure 4. If we make small changes in the para-
meters, one at a time, and compute c- for the altered para-
meters, we find in general that c- has changed by a small
amount Sc-. Changing Pj by a small amount 6pj moves the
solution point from Q along the curve DQD1 toward D. (See
figure.) Changing p2 by a small amount 6p2 moves the solution
point from Q along the curve AQA1 toward A1. The differences
in elevation resulting from these small displacements within
the domain of uncertainty are given by Eq. (.2.3).
25
-------
rain
max
Figure 4. Construction of first-order sensitivity
measures and coefficients. First-order
sensitivity analysis approximates the
true solution surface bounded by the curves
BDB', BAG, B'A'C', and CD'C' by the tangent
plane at point Q which is bounded by the
straight line segments 363', Bay, B'a'y',
and Y^'Y'J respectively. Using this approx-
imation severely underestimates the range
of uncertainty in the solution c. and provides
no warning of such structural features as the
deep valley at V.
26
-------
In the limit of infinitesimal displacements Sp^,
the ratio of the deviation in c- to the displacement in p.
is the slope of the solution surface in the direction of the
p. axis at the point Q. The slopes thus defined are the
J
first-order sensitivity coefficients Z- . defined by Eq.
»J
(2.5). First-order sensitivity analysis, in effect,
approximates the shape of the solution surface bounded by
the curves BOB1, BAG, B'A'Ci and CD'C1 in Figures 2, 3,
and 4 by the plane bounded by the straight lines 36g',
BCCY'J B'a'y'' and Y^'Y' which is tangent to the solution
surface at the nominal solution Q. (See Figure 4.)
For very small displacements about the "best values,"
the true solution surface differs from the tangent plane
by only small differences in elevation. In this regime,
comparison of first-order sensitivity coefficients tells
us to which parameters the solutions are most sensitive.
But first-order sensitivity analysis provides no warning
whatever about such significant features as the deep valley
at V, only a modest distance in parameter space from the
"best values." Nor would this analysis indicate the full
range of uncertainty in the solutions c- over the domain
of uncertainty of the parameters p . For the case shown
in Figure 4, the range between extreme values on the tan-
gent plane-- [C-J(Y) - c-j (31)] -- is less than half the
range for the true solution surface-- fc.(C') - c.j(B)j .
Moreover, the values of the parameters corresponding to the
extreme values on the tangent plane are completely different
from the values of the parameters for the true extremes of c-
I
If we wish to learn something about the shape of the
true solution surface at finite distances from Q, we can
perform a Taylor series expansion about c^Q), evaluating
not only the first derivatives of ci with respect to all
the parameters p., but also higher-order derivatives
J
(higher-order sensitivity coefficients). Including second-
27
-------
order coefficients would be an improvement over first-
order sensitivity analysis, but it would be inadequate
for the case shown in Figure 4 since Q is an inflection
point of the curve DVQD1. Thus the second derivative
2 2
9 c-/9pi vanishes at Q, and we would still have no intima-
tion of the valley at V.
In general, a truncated Taylor series approximation to
the solution surface can be characterized by a radius of
approximate convergence , a small distance in parameter space
about the point Q within which the series approximation differs
from the true solution by less than some specified e. Regard-
less of how many terms are included in the expansion, however,
the difference between the approximation and the true function
depends on the still higher-order terms which have been omitted,
Unless we have sound theoretical reasons to show that higher-
order derivatives become vanishingly small (and this is almost
never the case for realistic air quality models), we have no
assurance that we will not overlook such features as the
valley at V a very short distance away from Q. Consider,
as an example, a function whose first and second derivatives
are of order unity, while all higher derivatives vanish,
except the 20tn, which has a value of 10 . The terms
involving this derivative are then of order (10 ) <5p ,
- 5
so that to achieve an e ~ 10 , the radius of convergence
would still be $10~ , even if we included all terms up through
the 19th order.
An alternative approach is to systematically examine the
solution surface by computing c. from the simulation model at
a significant number of locations in the domain of uncertainty
of the parameters {p}. Unless the functional forms involved
in the various coupling, source, and sink terms in the set of
differential equations (2.1) are simple expressions, it is
difficult to be sure we have not accidentally overlooked
28
-------
significant structures hidden between the locations actually
examined. Several of the sampling methods discussed in a
later section attempt to deal with this problem in one way
or another.
We noted earlier that the shapes of the solution surfaces
c.j ({p}) depend only on the structure of the basic set of
coupled differential equations (2.1). If the functional
forms for the sources, sinks, photochemical rates, or chemical
kinetics rates are changed, or if new reactions are added to
the system, or other reactions are removed from the system,
all of the solution surfaces are changed. Peaks, pits,
undulations, ridges, gorges, and cliffs in the surfaces
corresponding to any particular chemical species c., for
example, will all be altered; some may disappear, while new
topographic features may appear. The shapes of the distribu-
tion functions P(c.) deduced for the values of the solutions
c- (see Figure 3) and the range of uncertainty for c- (see
Figure 2) will also be changed. Thus, whatever method of
sensitivity analysis or uncertainty analysis is used, the
evaluation of sensitivity or uncertainty must be repeated
ab i n i t i o each time the system of equations (2.1) is altered.
29
-------
SECTION 3
SYSTEMATIC SENSITIVITY ANALYSIS METHODS
Concise descriptions are given in this section of the
candidate systematic sensitivity analysis methods, starting
with the logically simplest "time-honored" method and
proceeding next to variational methods and, finally, to
sampling methods. We subdivide this section into two
major subsections: Deterministic Methods and Sampling
Methods. The deterministic methods include both the "time-
honored" and the variational methods.
Arbitrary, unaveraged deviations, , are the
I IK
measures used in sampling analyses. Table 1 summarizes
the salient features which appear in this section.
DETERMINISTIC METHODS
The Time-Honored Method
Trial-and-error testing of model parameters has stood
the test of time in all scientific disciplines. This time-
honored method continues to be the most widely used sensitivity
testing method because it is flexibly applied to arbitrary
problem applications and, most significantly, because
most researchers understand the trial-and-error method and
are comfortable with it. Developments in our present study
suggest that the time-honored method will continue to have
an important role to play in conjunction with applications
30
-------
of more-highly automated sensitivity tests to complex air
quality simulation models.
The time-honored sensitivity method is applied by
ily altering the value of model par
(x. ,t) and then noting the deviations,
simply altering the value of model parameter, p., at any
J
6Ci(x ,{p},t) = c.j(x ,p9 + 6p.j,t)- CjCx ,Pj,t), (3.1)
which are obtained from the respective model solutions
for c.j(p° + SPJ) andc-j(pp. (see Figure 4.) The major
features of this method are:
Arbitrary deviations 6ci are economically obtained for
applications having poorly defined probability
distributions- of model parameters.
Variances a,(x't') and expectations of sensitivity
(2)
coefficients < Z .>, < Z. '. >, etc., can be evaluated
' » K 1 , J K
but may be extremely costly to compute. (Examples
will be given in Section 4.)
Sensitivities and uncertainties which are continuous
in time, t, are obtained for the entire problem time
domai n.
It is a non-dedicated method; that is, the time-
honored method applies to arbitrary air quality
simulation models without any need to alter the air
quality simulation model itself, or to develop special
programming for performing sensitivity analyses.
31
-------
Variations in photochemical kinetics and fluid
dynamics parameters are treated with equal facility
(whether the chemical kinetics and the fluid dynamics
are treated separately or in fully coupled form)
over the entire space and time domains of interest.
Applications to date include the full range of coupled
kinetics and hydrodynamics.
Clearly, the time-honored method has a tremendous
combination of simultaneous features which other sensitivity
methods do not yet have; and, depending upon the craft
and specific motivations of the practicing scientist, the
method is also both systematic and economical.
Variational Methods
As we discussed in Section 2 above, the basic task of
sensitivity analysis is to understand the shapes of the
solution surfaces defined by considering the solutions
c-(x, (p},t) at any given position and time to be functions
I
of the uncertain parameters of the system. If we consider
such a solution surface (as was illustrated in Figures 2,
3, and 4) at a sequence of times, the surface itself
changes shape with time. If we are interested in just a small
neighborhood of the "best values" of the parameters (p>,
we could in principal compute the solutions c^C^.t) for a
large number of sets of values of {p}in the vicinity of the
"best values" (p0) and then laboriously construct the solution
surfaces for each time of interest to obtain numerical
(2)
approximations to the sensitivity coefficients Zi k, Zj ^ ,
etc., by straightforward but tedious means.
32
-------
But a far simpler method is to treat the sensitivity
coefficients themselves as dynamic variables. In other
words, we want to obtain equations of motion describing
the change in time of the slopes and curvatures, etc.,
of the solution surfaces. These can be derived analyti-
cally and then solved once by numerical methods instead
of solving the original set of equations (2.1) many, many
times for different sets of values of {p} in the vicinity
of {p°}. Variational methods, in brief, are based on
this latter approach.
(i) Chemical Kinetics with Constant Rate Coefficients
Most automated sensitivity analyses implemented thus far
have been applied to the case of purely chemical kinetics in
which rate coefficients are held constant in time. Variational
methods are also sometimes referred to as direct methods
when problem parameters are held constant. By any terminology,
variational methods are computationally dedicated in the sense
that a set of equations for sensitivity coefficients must
be derived and solved in conjunction with the original model
equati ons.
Photochemical kinetics model equations are nonlinearly
coupled ordinary differential equations of the form
ti = dc./dt = RjU.-.fp}, t), i,j
33
-------
In the absence of spatial dependences, the parameter array
{p> contains all of the variables listed in Eq.(2.2),
except c.j B , _y_ , and KD-
The sensitivity coefficients are defined as' partial
derivatives of the solutions cn- with respect to the uncertain
parameters:
Zi ~- - ' i=1»N> k=l,M. (2.5)
Dynamic equations of motion for these coefficients can
be obtained by differentiating (2.5) with respect to the
time. Assuming that the order of differentiations may
be interchanged, we obtain
3R: N 3R. 7
' *"= L ^i k
J » K »
i,k dtu. r ,D v"i; ,D L
3pk 9pk j = l
i,j=l,N, k=l,M. (3.3)
Considering the solutions c- as elements of a column
vector of length N, we can think of the arrays Z. , r, Z-,u and
i 1 K. " 1 K
R. , E 3R./3p. either as matrices of N by M elements or as sets
1 5 K i K.
of M column vectors each of length N. The Jacobian matrix
J. . = 3R-/8C- is a square matrix of N by N elements. If
' 9 0 ' \J
we view the set of N by M equations as M vector equations,
we can write
+ J- (Z I . , k = 1,M, (3.4)
k
where it is evident that the columns of Z. . are completely
I , K
uncoupled from each other.
34
-------
Important features of this type of sensitivity analysis
application are:
Eqs. (3.3) can be automatically compiled and solved
for arbitrary chemical mechanisms in the same manner
as the chemical kinetics equations (3.2)^ '. This
eliminates tedious and error-prone hand-programming
and thus allows very rapid turnaround for solving !
arbitrary chemical mechanisms, in both the kinetics
and the sensitivity programs. ,
I
Although the M vector equations (3.4) are decoupled <
from one another, both R1 and J depend explicitly on
the current values of c^ which are determined by
solving (3.2). To obtain fully reliable solutions,
each of the M equations (3.4) should be solved
simultaneously with (3.2), resulting in a system of
2N simultaneous differential equations to be solved
to obtain each of the M column vectors making up Z. ^.
I , N
The cost of simultaneous solutions, using modern
/ o \
stiff ODE solvers,^ 'for the combined set of kinetics
and sensitivity equations greatly exceeds the cost
of solving the kinetics system separately from the
sensitivity equation system. It appears possible
to reliably accomplish this separation in many
photochemical problems, although some additional
research is still needed on this aspect of the
problem.
The method is first order. Although equations can
can be similarly derived for higher order sensitivity
coefficients, ZJ^, ZJjkr etc., no consensus has yet
emerged on the value which would be derived from
higher order evaluations. (For many research users,
first order terms provide sufficient information, whereas
for some policy applications the contributions of ;
higher-order terms may also be desirable. (See Eq.(2.4).)
35
-------
Sensitivity coefficients are evaluated continuously
over the entire time domain of interest.
t When statistically significant information about
the distribution of input parameters is available,
output variances, valid to first-order, can be
(9 )
evaluated.v '
(i i) Chemical Kinetics with Time-Varying Rate Coefficients
When rate coefficients are varying in time (e.g.,
diurnally varying photochemical rates) the sensitivity
coefficient is defined as
Z = O/3e) fc-Uj, ... ,CN; . PR( t )+egk (t) ;t] , (3.5)
e = 0
where the function gk(t) measures the uncertainty in pk(t).
The sensitivity equations are then given in analogy to Eq.
(3.4), by
Z = R(t) + J- Z. (3.6)
The functional derivatives Rk(t) can be evaluated, for the
kth parameter, in variational calculus, by
Rk(t) = (3/3e) [R(ca,..,cN; pk(t) + egk(t); t]
e=0
(3.7)
This method has been applied for diurnal variations in
a simplified atmospheric ozone example by Dickinson and
Gelinas^ , and specific choices for the function 9i<(t)
have been indicated. The important features of the
diurnally varying examples were qualitatively much the same
as in cons.tant rate coefficient case which was discussed
above, but with some operational distinctions.
The evaluation of functional derivatives requires
some dedicated programming of the functions gk(t)
in the sensitivity program.
36
-------
0 Kinetics and sensitivity solutions can be much more
costly to obtain when rate coeffievents are varying
with time than when alternative constant rate
coefficients are used. For example, experience in
diurnal phothchemistry has shown^ ' that diurnal
kinetics integrations can be ten to fifty times
more costly than corresponding constant (daily
averaged) rate photochemistry. This proves to be
a serious practical constraint for all dedicated and
non-dedicated sensitivity methods.
(i i i) Non-Reactive Hydrodynamics
Examples of variational sensitivity methods, as well
as the Fourier Amplitude Sensitivity Test, which is a
sampling method, have recently been presented and compared
by Koda et al^ , they present several computational examples
for the atmospheric diffusion of inert species while also
discussing reactive systems. The governing non-reactive
hydrodynamics equations are (in one spatial dimension).
3c.(x,t) 8 r 3c,(x,t) .
= [KD(x) !- 1 . i = i.N,
8x J
where KD(X)represents spatially varying diffusion coefficients,
which are assumed to be identical for all species. Boundary
conditions are represented by
- KD(0) - = S.(t) at x = 0, (3.8)
where S^ represents all sources and sinks of the i-th chemical
species at time t at the boundary x = 0, and
Hi.
3x = ° at x = H(t), (3.9)
where H(t) is a time varying boundary height. Initial
conditions, c- T(x) = c-(x, 0), are usually assumed to be
1,1- i
known with no error.
37
-------
Sensitivity equations are obtained by first discretizing
Eq. (3.7) in the spatial domain, such as by the method of
lines. ^ The partial differential equations (3.7) are there-
by reduced to a large, finite set of ordinary differential
equations in which c.(x,t) and KD(X) are evaluated'only at
discrete spatial locations x£, £ = 1 ,L. For a single
species the discretized form of Eq. (3.7) is (in vector
notation, and including boundary conditions)
dC(t)
= A(Kn)-C(t) + S(t)5 (3.10)
dt U
where C(t) = [c(Xpt), C(x2,t) C(xL,t) ] T, KQ = [ KD (xj),
KD(x2),. . . ,KD(xL) ] T, and S(t) = |s (Xj ,t), 0 ,. . ,oj T. Initial
conditions are denoted by CQ = C(o)» and the matrix A(KD) is
a finite difference representation of the diffusion terms in
Eq. (3-7). For constant KD, the analytic solution for Eq.
(3.10) is
t
C(t) = e A(KD}t [c + C e"A(KD)1:l S(t')dt'] . (3.11)
L ° 0
o
Sensitivity coefficients for this case are defined by
T
3C(x,,t ) 3C(x9,t) 3C(x, »
Z,(t) = l 2 L
O
9KDj
j=l,L. (3.12)
where KD- sKD(x-). A sensitivity equation (in vector form) is
obtained by a derivation analogous to that used in the develop-
ment of Eqs. (3.3) and (3.4). This development yields the
sensitivity equations
dZ-(t) 3A(K)
. C(t), j = 1,L,
(3.13)
= A(KJ ' Z,(t) + - ^- . C(t),
J
which can also be solved analytically.
38
-------
Important features of this method are similar to those
mentioned in (i) above, namely:
The solutions to Eq. (3.10) and (3.13) can be
evaluated either analytically or numerically
for arbitrary applications.
The diffusion coefficient at each spatial grid
point is treated as an individual, uncertain
parameter.
The method is first order. Higher order terms
have not yet been evaluated for this case, to
the best of our knowledge. (The utility of higher-
order terms is application-dependent.)
Sensitivity coefficients are evaluated continuously
in time over the entire space and time domain of
i nterest.
This method can be generalized to also include (constant)
wind advection terms. In such cases, the wind values
assigned to each grid point are also treated as individual,
uncertain parameters.
Output variances, valid to first order, can be
evaluated when statistically significant information
about the distribution of diffusion data is available.
(i v) Reactive Hydrodynamics with Constant Rate
Coefficients and/or Diffusion Coefficients
Sensitivity coefficients are defined for the i-th species
(i n vector form) by
Zi,j *k.C1(t> - IP,
J J
1=1,N, j=l,(M+L), (3-14)
and the spatially discretized sensitivity equations are given
by:
3Z- .(t)
1 »J _
3t
3A(KD)
. +
A(Kn) + J
1=1, L, j=l,(M+L).
Z, .(t),
(3.15)
39
-------
In this case, p. spans all chemical rates k,...kM and
diffusion coefficients KD ^ -,..., KD L at all spatial grid
points. The total number of parameters is (M + L).
To our knowledge, applications at this level of complexity
have not yet been done for air quality simulation models.
The important features follow the same general lines as
in the cases which have previously been discussed. Clearly,
the amount of dedicated programming and data management can
approach prodigious proportions for this level of applications
where there can easily be 30 chemical species, 75 photochemical
reactions, and hundreds of grid points.
(v) Reactive Hydrodynamics with Space and Time Varying
Coeffi cients
Cases in which Kn(x), S.(t), H(t), and k- can vary
" (n) J
have been considered by Koda et al^1. Equations have been
developed for the functional derivatives, 6C - (.x1 , t' )/6 K[)(x)
6ci (x1 ,tJ)/£S.j(t), &c.(x' ,t')/6k.., and &c- (x', t' )/
j A j j i j
(3.16)
40
-------
subject to initial conditions
6ci (x, 0) = 0, (3.17)
and boundary conditions
36Ci(x.,t) 3c?(x»t)
-K°(0) - - 6KD(0) = 6S1(t) at x = 0 ,
3x 3x (3.18)
and
36c,(x,t)
J
3x
= 0 at x = H, i=l,N. (3.19)
suggested that input parameters be expressed
as expansions of orthogonal basis functions, which leads
to the following development of sensitivity equations:
(a) Select a set of functions {<]; (x_,t)}.
(b) Multiply Eq. (3.16) by \l>- and sum over i = 1,K.
(c) Integrate the resultant equation in (b) over the
intervals ^x = [o,HJ and t = fo,fj .
(d) Apply the initial and boundary conditions, Eqs.
(3.17) - (3.19).
(e) Specify a defining equation for ij>.(x,t)
3^. ~ Sty. "n°
7T " ^ (K°D(-)"97')"
wi th
ax
i= o, x = [O,H]. (3.21)
41
-------
(f) Specify a terminal condition on the { ty.}, e.g.,
^(x.T) = 6ijjS( x - x') , 1.j = l,N,
(3.22)
where 6.. and 6(x - x1) are, respectively, the Kronecker
' J ~ ~
and Dirac delta functions.
The results of performing operations (b) and (c) and solving
Eq. (3.20) gives N H J q
Sc.(x',T) = - I S S*D (^ S at1 Iw^1 dt dX
1 n r\ 0 X OX
N T
+ E S 6S..(t) ^j^O,
.j=l 0
M N H T aDo
?s r,. o K ^
6k,
r- .»> r> on--i
J I -J i ^*i dt dx- <3-23)
=100 8kj
The functional derivatives which we seek are
N J °
= r^ \
6c-(x1,T)
' r^ \ J J
6Kn(x) n 9x
D ~ j=l U
6c,(x',T)
- - =
il'j^O.t) , (3.25)
N H T 3Ro
S ^o-i dt dx . (3.26)
6kj £=1 0 0 3kj
42
-------
A similar exercise can be carried out in order to
obtain Sc (x ' , t ' )/6H ( t) when H is a function of time.
1
Sensitivity functions (3.24) to (3.26) are obtained by
solving a decomposed equation for the nominal values
c?(_x,t), and solving Eq, (3.20) for ^£j(x_,t) in order
to carry out the operations in (3 . 24)- (3 . 26) f Koda
et al ^ ' have obtained solutions for
<5c(x't) 6c(x' ,T) , 6c(x',T)
_ _ , d 1 1 U _ ~ __
<5KD(x) ' 6S(t) 6H(t)
for a single non-reacting species in a 1-D (vertical)
boundary layer simulation and have related the resulting
functional derivatives at specific points to sensitivity
functions obtained in a discretized grid calculation which
was discussed in (iii) above. The relationship at spatial
location x.' between sensitivity coefficients which are
developed from functional derivatives and those which are
developed from discretized grid point formulations i's
x
'a £, ( x ' , t ) ^ 6 d . ( x ' t )
-.(x',t) EL-Z - = \ 1 ~
(3.27)
ac,(x',t) 6c,(x',t)
_J -- - = AX - ^ = -
(3.28)
where Ax represents an appropriate grid spacing and 6p,-(x)
J
is assumed to be zero at grid points other than x-
43
-------
The major features to note for the variational methods
involving functional derivatives are:
t The solution of the system equations for nominal
values, <:(>(,t), and of the adjoint equations for
tyo-(x.\t) are complex calculations, even with no
JC 1
reactions present; but the method produces highly
rigorous sensitivity measures.
This highly dedicated method will be difficult to
automate sufficiently for reliable, rapid turnaround
calculations in arbitrary model applications.
The method is first order and subject to the same
considerations which were previously discussed. No
higher order evaluations have yet been made.
Sensitivity coefficients are evaluated at specific
space and time points (xjt1) but rigorously contain
correlative effects which have occurred along
characteristic trajectories which lead to (_x',t').
Arbitrary deviations are the measure of sensitivity
in this method. It is not yet apparent that
expectation values of observables could be readily
obtained when meaningful parameter statistics are
available.
The Green's Function Method
Hwang, Dougherty, Rabitz, and Rabitz^ 'have developed
a variant of the direct variational method of sensitivity
analysis for chemical kinetics systems. The basic develop-
ment of the nethod is similar to that described in our
discussion of the variational method for chemical kinetics
with constant rate coefficients.
Starting from the kinetics model equations (3.2), the
dynamic equations for the evolution of the sensitivity
44
-------
coefficients Z^ k are derived as before, obtaining equations
(3.3) or (3.4). As we remarked in the earlier discussion,
Eq. (3.4) can be viewed as M uncoupled equations, each des-
cribing the evolution of a column vector of length N:
(Z)k =(R')k + J'(z)k ' k = ^M. (.3.4)
In the most general applications of the direct method,
each one of these equations is solved numerically together
with the basic kinetics equations (3.2) in order to allow
the stiff equations system solvers to adjust the time steps
optimally for the numerical integration of each of the vector
equations in (3.4), without introducing numerical errors due
to interpolation from an independent solution run for the
basic kinetics equations (3.2).
The Green's function method continues its development
by assuming that each of the vector equations in (3.4) can
be reliably and safely solved by itself, interpolating
values of c. obtained by a careful separate integration of
(3.2). Since Eq. (3.4) is linear in Z, the solution of
the i nhomogeneous equation can.be expressed in terms of a
Green's function K.j^(t,T) as
Tl R;R UI,...CN; {p};T )dt,
0 i=l
i,£ = l,N, k = l,M, C3.29a)
if Pk is one of the initial conditions, where R'. = 8R. ,3p ;
or as
N
E
R=l
ZiJ<(t) = K(t,r) R (c ..... C; {p};r) dt
(3.29b)
if pk is any other parameter of the system. In either case,
45
-------
the Green's function is the solution of the homogeneous
equation
N
Kik(t,T) - JT Ju (c15..., CN; (p};t) Ku(t,T) = 0,
1=1 i,*5k = 1,N, (3-30)
where J ^ E 3^.7 3cA is the Jacobian matrix as before, subject
to the causality condition
I (t,r) = 0 for all t N), which is usually the case for realistic simulation
model s .
Practical implementation of the Green's function method
proceeds as follows:
(a) The simulation model equations (3.2) are integrated
numerically using the nominal values 'j> } f or the uncertain
parameters for times t=0 to t=T, recording or storing values
for the solutions c^ at numerous intermediate points in time,
(b) The time interval T is subdivided by a coarse arid of
1 + 1 quadrature times xm, TQ= 0
-------
(c) Eq. (3.30) is solved numerically for each of the
subi nterval s x ,< t< T with initial conditions
ffl J. HI
K-uU ., Tm ,) = 6... The values of c-(t) needed to evaluate
i K m -r I in i IK i
J,o(t) are interpolated as needed from the stored solutions.
(d) The sensitivity coefficients themselves are computed
by using the appropriate form of Eq. (3.29), replacing the
integrations over time by numerical quadratures on the coarse
grid T . It can be shown that the Green's function
for any two points on the mesh may be obtained as a
product of the Green's functions obtained for the inter-
vening subi nterval s :
N N
m,nr, (3_33)
(e) If higher order sensitivity coefficients are
desired as well, they may also be computed by quadratures
with the same Green's function:
vt>T)
(3.34a)
where the differential operator D- is defined as
N ^
3 + T 7 ^
Dj = - L Zq,j (3.34b)
^p,- q=i acn
j q
The method is clearly a dedicated method, requiring
new programming for each problem, but it is highly amenable
to automatic programming methods just as the chemical
kinetics equations are.
47
-------
SAMPLING METHODS
Monte Carlo Method
Stolarski et al^ ' have recently presented an atmospheric
uncertainty analysis which is based upon Monte Carlo selection
of trial values for chemical rates in a 1-D stratospheric
model which gives steady state solutions of Eqs. (2.1) with
v^ = 0. Many of the photochemical reactions of stratospheric
interest have been measured often enough and well enough
that it is reasonable to assign probability distributions for
the parameters describing these reactions. Most rates which
are known to be significant in the stratosphere are uncertain
within factors ranging from 1.1 to 10.
Hypothetical reactions must simply be included and
assigned nominal values and uncertainties on a tentative basis,
The Monte Carlo method is executed as follows (see Figure 5):
(a) A log-normal distribution is assigned to the
probability that true reaction rates lie between
(k°j/r.) and rjkj°, where r->l is the factor by which
k-j isjuncertain. '
J
(b) The atmospheric model is run, using a chosen "best
set" of reaction rate values and center line values
for other parameters.
(c) A random number generator is used to select trial
values for the entire set of reaction rates, simul-
taneously.
(d) The atmospheric model is run again, using the set
of simultaneously perturbed trial reaction rates
from (c).
(e) The processes (c) and (d) are repeated until enough
output samples are accumulated to evaluate
statistically meaningful mean values and variances
of the output variables. Stolarski et aZ(l5) ran
sufficient cases so that the accumulated 10 im-
precisions on both sides of the output distributions
converged to within 1 or 2%.
48
-------
Figure 5. Schematic illustration of the Monte Carlo
method of sensitivity analysis. A random
number generator is used to select values
of all the uncertain parameters within the
domain of uncertainty. The simulation model
is then run using those values for the para-
meters. Values thus computed for c^(x,{p},t)
are then tallied and analyzed by standard
statistical methods. The distribution of
values obtained in this way is illustrated
above as a histogram, with mean value
and standard deviations a~*~ and a". Stolarski
et al assume log-normal distributions for
the uncertain parameters, but any distribution
function could be used, including a joint
probability distribution for two or more
correlated parameters.
49
-------
The important features of the Monte Carlo method
are:
The method is non-dedicated and thus readily applies
to arbitrary air quality models. To date only
chemical rate uncertainties and boundary condition
uncertainties have been considered.
The method has, so far, been applied as an
uncertainty analysis method which evaluates output
variances. Questions of how a shift in one parameter
value would change the model response to all other
parameters have not yet been indicated.
Dependences of the method upon selected probability
distribution functions have not yet been indicated.
Best and worst case estimates can be identified and
retained from the accumulated test cases.
Initial studies are underway forcorrelated parameters.
Computer economics presently require that the 1-D
stratospheric model be physically simplified and
solved for steady-state solutions rather than the
far more costly transient solutions. Stolarski
et al emphasize that the feasibility of using the
Monte Carlo method depends most heavily upon reducing
the computational cost of the atmospheric model to
an absolute minimum. In their case, complete repro-
gramming and physical simplifications were required
for their original, more general atmospheric model;
this went far beyond simply altering the original
model and obtaining nominally improved economies.
Sensitivities were determined separately by the
time-honored method applied to individual ratesr
The rates were varied, one at a time, by +_ la and by+2a.
The number of test cases which were required to obtain
statistical convergence in a-j(x^,t) was proportional,
in some sense, to the number of sensitive parameters
rather than to the total number ofuncertain
parameters. Convergence was more rapid for those species
and locations which were most affected by rates with
small imprecision; relatively poor convergence was
experienced for species such as the methane oxidation
chain, for which important reaction rates are quite
uncertai n.
50
-------
Stolarski et al emphasized that the uncertainty
analysis must be redetermined when:(i) a new reaction
is included in the model (ii) substantial changes
occur in parameter values or their imprecisions,
or (iii) substantial changes occur in the basic
understanding of atmospheric processes.
FOURIER AMPLITUDE SENSITIVITY TEST METHOD
The Fourier Amplitude Sensitivity Test was first developed
by Cukier and Shuler ejt a_]_ to provide a means for
determining the sensitivity of the solutions c^Cx., (p},t) of
a large set of coupled rate equations to uncertainties in
each parameter p when averaged (in a certain sense) over
uncertainties in all the other parameters.
Figure 6(a) i11 ustrates the basic idea of the FAST method.
For simplicity in illustration, consider a system of coupled
rate equations depending on only two parameters,pi and p^.
Although the investigator can select some "best value" for
each of these parameters, the true value lies within some
range of uncertainty about that "best value".
If only one parameter were uncertain, we could profitably
ask, "What is the slope of the solution surface as that single
parameter is varied over its range of uncertainty?" (See
line AA1 in Figure 6(a).) But other parameters are also un-:
certain; thus in the general case we do not know where on the
solution surface we are. It is thus of some value to ask,
"What is the average slope of the solution surface in the
direction of parameter p- when we average over the full range
\J
of uncertainties in all other parameters as well?" In other
words, what is the average slope of curves AA', BB', CC1,
etc.?*
*
In some cases, the average slope may be zero or very small for parameters
which in fact are sensitive. The average slope transverse to a symmetrical
trough or ridge, for example, vanishes, even if the trough is very deep
(or the ridge very high). To avoid being mislead by such cases, a
squared norm is sometimes used to indicate averaged sensitivities.
51
-------
Ci(x,{p},t)
>P2
Schematic illustration of the Fourier Amplitude
Sensitivity Test (FAST). The FAST algorithm
generates a closed search curve across the
domain of uncertainty of the parameters which
begins at the "best values" of all parameters.
Values of c.^ at each sampling point (marked by
large dots in the figure) are tabulated against
the search parameter s along the search curve
to produce a periodic function of s as shown
in Figure 6(b). Note that the sampled points
shown here do not include the extreme values of
c^ (points B and C').
In principle, one could determine the slopes of the
solution surface along each of the parameter axes at a large
number of points on the solution surface above the domain of
uncertainty of the parameters and compute suitable
52
-------
averages. Such brute force methods, indeed, have been used
on occasion for systems involving only a few parameters and
small numbers of coupled differential equations. The FAST
algorithm, on the other hand, provides a technique for
selecting much smaller numbers of sampling points. The
method is most likely to be successful for those systems in
which the solution surface is smooth and preferably mono-
tonic, with minimal structure in the way of deep valleys,
sharp ridges, or peaks close to potholes.
The method consists of four steps. First, a set of
M integer frequencies o)£ is selected subject to certain
constraints to be discussed below. Second, one frequency
is assigned (arbitrarily) to each of the parameters considered
to be uncertain, setting
P£ = p° exp (u^Csin w^s ) } (3.35)
where s is a parameter of variation; p° is the "best value"
of p?; and u? is a continuous function of its argument
selected so that its extreme values result in values for
p. at the extremes of the range of uncertainty for p^. Note that
p.£ is then periodic in the parameter s; as s varies, the point
£ determined by (3.35) moves about in the domain of uncer-
tainty of all the parameters. If the set of frequencies
co£ are mutually prime, the point £ traces out (exactly once)
a closed curve (a Lissajous figure) within the domain of
uncertainty as s ranges from 0 to 2ir.
Third, the values of p^ resulting from setting
s = 2-nq/H*, q = 1, 2, 3,...., N*, are used in running the air
quality simulation model from given initial conditions to
time t. The values of the solutions c . (_x,{ p (s ) } ,t) are
stored after each simulation run. II* is proportional to
the largest frequency to .
A/
53
-------
w
O
max
jnin
i
t/2
TT
S
3TT/2
2lT
Figure 6(b).
The solution c.(-{p(s)}) is a periodic function
of the search parameter s. The Fourier amplitudes
of this function at the frequencies assigned to
any parameter p. provide an estimate of the
average slope of the solution surface in the
direction of the p.-axis.The search curve does
not in general reach the extreme values of c.
Fourth, the solutions c. for a given position and time
(x.»t) are now periodic functions of s itself (Figure 6(t>).) Fourier
analysis of this function yields amplitudes corresponding to
each of the integer frequencies to^ :
2ir
Ad> = fc ^ ci({p(s)})sin (iHs) ds . (3.36)
£ 0
Only sine components occur in the expansion of c^(-p(s) )
because of the symmetries of Eq. (3.35). This integral can
be approximated by a simple summation over the accumulated
soluti ons:
A
2
^ £ ci({p(sq)))sin
q=i
sq)
(3.37)
The Fourier amplitude A measures the contribution made by
variations in p^ to the total variations in c.,- as all the
parameters p oscillate back and forth across the domain of
uncertainty.
54
-------
It is thus not surprising that the Fourier amplitudes
defined above can be shown to be proportional to a weighted
average of the partial derivatives of c^ with respect to each
parameter:
max max max .
Pi Po PM 3c.
A A " P(p) L dp, dp,... dpM
t
v)
p-n pjm pnnn
a < 7 > I 3 38)
J0' \ -J
-------
mm
max
nan
'(P,)
Figure 7. The functional forms used to generate the search path
in the FAST algorithm are highly arbitrary, resulting
in a great deal of ambiguity about the probability
distributions P(p^).
56
-------
The number N* of simulation runs required is determined
by the highest frequency in the set {w^K since each cycle
of a sinusoidal wave must be sampled at several points in
order to determine the amplitude of that Fourier component.
For economy, we would like N* to be as small as possible; for
accuracy in evaluating the amplitudes, we would like N* to be
as large as possible. If we make N* larger, then we can also
increase tom, . the largest frequency in the set, which has
in a x
the effect of folding the search curve back and forth in the
domain of uncertainty more times, increasing the likelihood
of detecting fine structures in the solution surface and of
the search curve coming very close to the highest and lowest
points (the extremes) of the solution surface.
Various algorithms have been developed for generating
sets of frequencies for use with FAST, The frequencies should
be chosen to minimize ambiguities in th.e determination of
Fourier amplitudes which arise from interferences between
frequencies or harmonics, As a minimum, the frequencies must
be linearly independent;
M
t a^ f 0 (3,39)
for integer coefficients an- subject to
M
£_ a. < Mj + 1, (3.40)
where Mj is an integer. Such a set of frequencies will then
be free of interferences to order Mj,
Because the Fourier amplitudes are computed from a finite
number N* of equally spaced sampling points, "aliasing" can
also confuse the picture. Aliasing occurs whenever a sum of
two or more frequencies (and/or their harmonics) is equal
to an integer multiple of the number of sampling points;
M
Y a', to. = JN* (3.41)
f=l n n
57
-------
where the coefficients a- and J are also integers. To
minimize such aliasing, we require
M
ai eoj j JN* (3.42)
wi th
M
L
< M+1. (3.43)
The frequency setdo^} will then be free of aliasing to order
My. For MT = 4, it has been found empirically that the maxi-
mum frequency and the number N* of sampling points necessary
? s
are proportional to M . No analytic expression has yet been
obtained for N* in terms of M and MJ-, such a formula would
have to be derived by the techniques of combinatorial algebra
which would have to provide first of all an explicit algorithm
for generating the set of frequencies {w,}.
The important features of the FAST method are:
It is a non-dedicated method which requires no
alterations of the air quality models.
It provides relative measures of expected sensitivities.
To date, only first-order measures, , have been
I , K
obtained in applications, although higher order informa-
tion is potentially available should there be a need
or interest in abstracting it.
Applications to date have emphasized pure kinetics or
pure non-reactive hydrodynamics.
t Consideration of correlated parameters is in the early
stages .
Computer economy of individual air quality simulation
models is a major factor in whether or not a sufficient
number of cases can be accumulated as adequate input
to the FAST method.
58
-------
Variances obtained by summing all specific frequency
components are generally inexact since only a finite
number of terms can ever be computed numerically.
The method can ultimately serve the interests of
both policy-makers and of scientific researchers.
Each application of the FAST method provides
sensitivities at a specific time point, in snapshot
fashion, rather than continuously in time as do the
deterministic methods. (The same is true for other
sampling methods, as we saw for the Monte Carlo
method. )
In difficult problem regimes, the results obtained
depend upon specific sampling distributions, as
is evident from considering Figure 7 and as
(19)
observed in practice by Penner;- ' If the solution
surfaces are very irregular, different search patterns
yield different output results.
Given a sufficiently dense search pattern, the method
can in principle resolve subtle system sensitivities
associated with complex wrinkles or pathological features
in the solution surface. (The extent to which this
possibility is ever realized in practice is problematical
due to considerations of computer economics as discussed
in Section 4.)
59
-------
TABLE 1
SUMMARY TABLE
CIIAKACTtmSTICS OF CAMDIDATt STSTtHMIC SCHSIlltllt AMD IIHCmAlll Tt ^HAHSIS METHODS
METHOD VARIABLES
OCUHHiHlSllC KTHQOS «C( (x.t).
(A.) 1!ME-HO»0«£D Z( j, z"{
(B ) VARMTIOHAL 7, ,, Ac, ( first
'J 'order)
Curstant parameters
(direct method)
(1) pure Kinetics
(11) hydrodynamics
(reactive and
non- reactive)
variable para-
cters 7 functional
derivatives]
(1) pure kinetics
(11) pure hydro-
4yna«ics
SAtf-HHG KTHOOS
(A) Monte Carlo o.(x.t)
1 ~
PROGRAMHIHC FEATURES
Nun-Dedfcattd-
Rpsults are obtained as
continuous functions of
time for all variables.
Dedicated- Results ar*
obtained as continuous
functions of tine for
all variables. Sensiti-
vity vodel should be de-
coupled fron air quality
model.
(1) Highly automated
(11) only 1-D transient
sensitivity models
have be"r<4fls are r.uw
reasonably near at
hand with Highly
automated features
due to recent ad-
vances In PDC
ethods.
(1) Seal-automated
(1i) Not-automated
Non- Dedicated.
Results ar« obtained at
single time points for
all variables
APPLICATIONS TO DATE
Models frcMi essentially all
scientific and engineering
disciplinespure kinetics,
! f>, 2-0, and 3-0 transient
rodel*.
Only puie kinetics and 1-0
models. Ho Myher dimensional
sy.Uas.
(1) atmospheric chemistry.
chemical laser*, and com-
hustlon rnenfstry
(11) 1-D atmospheric dynamics.
further developments can
be readily made In 10,
whereas scwcwhat nor*
extensive developmental
f f 01 ts art 3 Oil re-
quired For 2-D applica-
tions.
(1) diurnal photochemistry for
the simple Chapman
neclianisB.
(11) Test proble«s for a
dlffuslvt boundary layer.
1-D ste*4? state stratospheric
ozone Model
REQUIRED INPUT AND
EUCUTION fOH M PABA-
H£Tf US
Hlnlnuff of ( MM) itr
qualfty node! Integra-
tions-Includes both
nominal and perturbed
paianeter values. All
paidineters can be tested.
Che integration of th*
air quality rwxlel using
nominal parajtetcr values
plus M Integrations of
sensitivity equations.
The above Integrations
plus solutions of adjoint
sensitivity equation! are
are required, plus some
double Integrals *u$t be
computed. In addition to
being more costly, these
operations are difficult
to automate for arbitrary
problem scenarios.
On tbe order of (y) ^
nodel Integration*,
where Y Is the nunfcer of
saaf>Te points per para-
«eter anrt K. 1$ ch«
nukaber of sensitive para-
ncters. Generally Y'% 2
and K^ 10 for computer
conol.y. All »tode1 paca-
COHMEHTS
This fs the favored wthod of mat
piactfcal tisers. No special pro-
nr*iw; are nrqulnd. The trlal-
and-error concept Is Intuttlvcl/
appealing and logically trans-
parent to » jjotlty of users.
No statistical constdeiatlons art
required. All »od*l paramters.
whether constant or varying can be
tested. Results are a»*.iablc ta
best and worst case Identification.
Coming Into widespread use as at/to-
Mted proqiaw becune opplirahle
to arbitrary problem ictnarlos.
Requires no pre-screenfng of para-
neters for atmospheric kinetics.
No statistical ccms l.lei ationi arr
required, and 2. , and or, (f1n,t
order) are conceptually trans-
parent » Apptars to be mure ecorm-
nlcal than all other alternative
nethixJs (In both human effort >*4
computer costs). Results are
amenable to best and worst case
Identification.
Functional derivatives require
significant hands-on control by
the user, and are considerably
more costly to solve than cases
with constant parameters.
In practical applications It wovld
be Biost reasonable to approxlaut*
time variations by plecewlse coat-
slant variabilities to that a nra
efficient direct method can be «.-.*d
In place of functional derivatives.-
Results an* expectation values
whlcfl ate generally dependent u»«m
assuiwd paraaicter distributions.
Pre-screenlng Is necessary to
estimate In advance the nuaf>«r ef
sensitive parameters which In tun
govern the rate of convergence of
the method (istvjls). Results an
good to all orders and are concett-
can be tested, but unity transparent. R«uHs art mot
applications to datt
eftfriiaslze cht*«!cal rates.
amenabU to bsst and worst ca»
Identlftcatlon uiilcss raw samp It-
data are additionally (n-ocessed.
'Even with prc-screenlng estlMtes,
tlte otmfcer of satuples required U
give convergence Is not fnown
exactly for each specific ?|(x.O;
nonitorlng Is retjulrtd during
evaluation of o '«,
partial varla
and lt >
Non-Dedicated. Results
art obtained at slngli
time points for all
vat tables.
Hodali from several On tlte order of
-------
SECTION 4
EVALUATION OF CANDIDATE SENSITIVITY METHODS
FOR AIR QUALITY MODEL APPLICATIONS
This section illustrates "what you get" and "what
you pay" with systematic sensitivity and uncertainty
analysis methods in their present states of development.
Evaluations of the candidate systematic sensitivity
and uncertainty analysis methods are given in the
discussion below in terms of our detailed assessments of
the performance characteristics of each method for:
(i) chemical models, (ii) 1-D transient models and (iii)
3-D air quality models. The results are self-evident and
indicate a mixed picture--- one in which both reasonable
and unreasonable options are avaiable, Accordingly, several
criteria are used for our evaluations.;
Moti vat ions. Fortunately, "what you get" can be very
much a function of what you want or need in choosing
systematic sensitivity and uncertainty analysis
methods. Optional methods are available for those
with poli cy-maki ng motives and for those in atmos-
pheric research. These respective motivations will
largely determine the selection of candidate analysis
methods for air quality simulation models.
Cost. This can be the make-or-break factor, in-
volving not only computer execution costs but also
human programming costs. Human programming costs
include the development of new sensitivity programs
(for dedicated methods) as well as the essential
tasks of preparing input data and writing preprocessor
and graphics routines. For large, complex regional
air quality models the human costs alone can prove
to be prodigious- depending upon specific applications,
as can computer execution expenses. This is true
61
-------
both for dedicated sensitivity programs, per se,
and for the multitudes of air quality model runs which
may be needed as input to non-dedicated sensitivity
methods. For smaller, purely photochemical kinetics
box models or for purely hydrodynamics models,
computer costs are more reasonable as are human
programming costs when automated compilation
features are used.
Useful ness. We particularly focus on criteria which
underlie whether or not potential users may sooner
or later be motivated to use candidate systematic
sensitivity and uncertainty analysis methods. In
addition to labor and computer economy, it is essential
that a method be conceptual1y transparent and that its
results be readily interpretable. (This factor alone has
accounted for a majority of practical users continuing
to stand by the time-honored method. The overwhelming
appeal of the time-honored method derives from several
factors; it is highly intuitive in nature, it requires
no special preparations to operate initially, and no
special operations are needed to interpret output
results.) An additional criterion of utility is
the potency of candidate analysis methods. Here,
the major promise of some systematic sensitivity
methods is their ability to identify more effectively
and to evaluate quantitatively subtle sensitivities
in complex solution surfaces. Potent methods also
have the capacity to anticipate unexpected events
[new discoveries or other surprises.)
Whether or not potency of a method should be associated
with the availability of higher-order sensitivity measures
is a subject of continuing debate. On the one hand,
policy makers desire reliable variances (uncertainty
measures) for the results of air quality models. The
62
-------
most reliable means of obtaining such information
is by sampling results from repetitive runs of the
air quality model (basically, the time-honored or
Monte Carlo method). It is clearly out of the range
of practicality to construct variances by first
evaluating and then summing terms involving
Z- -, Z- .. , etc,,from either the deterministic or the
1 »J 1 » J K
FAST methods. On the other hand, users who require
sensitivity measures are generally quite satisfied
with first-order sensitivity coefficients (regardless
of the methodology which is used to generate the
coefficients) because first-order coefficients contain
all of the information which most users want and need
*
on a value-for-yalue basis \n practical applications,
With some reflection, all of these criteria are
recognized to be branches of the same tree, the tree of user
acceptance. This should not be surprising because these
issues of user acceptance have perpetually governed the rate
of development and the pattern of use of systematic sensitivity
and uncertainty analysis methods. Tables 2 and 3 provide a
summary of our evaluations of the candidate analysis methods.
Detailed discussions of these evaluations follows directly below,
BASIC CASES
Three basic examples are used for purposes of illustrating
evaluations in this section. In all cases, computer cost
* This view has been expressed repeatedly by the overwhelming
majority of scientific practitioners; it represents, in fact,
the hard reality of the situation. This does not say that
higher-order sensitivities are insignificantwe have clearly
shown the contrary in previous sections--but rather that, on
a value-for-value basis, first-order sensitivity coefficients
provide a major part of the information which practical users
seek from sensitivity analysis.
63
-------
estimates assume that a Control Data Corporation (CDC) 7600
computer is used for sensitivity and uncertainty computations
at a cost of $1200 per execution hour ($20/minute). This is
obviously an arbitrary basis, but it can be converted approxi-
mately to the appropriate, equivalent basis for other computer
systems. The main points of our evaluations should be quite
evident within this selected context despite the great vari-
ability of computer execution times for the various models
and computers which are used by different investigators.
A. Chemical Models (purely kinetics)
This base case assumes that:
Thermal and photochemical rates are constant.
(Time-varying rates are far more costly to integrate
using existing ODE system solvers)
External chemical sources are constant.
f The system has N chemical species and M reaction
rate parameters.
Numerical integration of the ODE system of kinetic
(species) equations has costs which are proportional
to Na . It is reasonable to assume that a^2-3, which
is representative of implicit solution methods for
stiff ODE systems, such as variants of the Gear method "'
Generally, a system of 20 species would require
approximately 15 seconds of CDC 7600 execution time
at a cost of about $5 per kinetics run; a system of
50 species requires about 4 minutes CDC 7600 execution
time at an approxinate cost of $80 per kinetics run.
Kinetics execution times are not very sensitive to the
number of reaction rate parameters, M.
Automatic kinetics programs are available which
allow users to select, compile, and solve arbitrary
reaction sets (with arbitrary rates), requiring
virtually no dedicated programming. "
64
-------
B. 1-D Transient Atmospheric Models
One-dimensional transient models have been applied |
primarily to reactive plumes and to global stratospheric
ozone modeling. In addition, 1-D models provide a convenient i
foundation for subsequent considerations of 2-D and 3-D |
regional air quality simulation models. This base case j
assumes that: i
thermal and photochemical rates are constant over
intervals of one hour or longer. (Instantaneous
variations greatly increase execution costs).
i
t External chemical sources are constant over
intervals of one hour or longer.
The system has N chemical species, M parameters,
and L spatial zones.
The cost of numerical integrations depends upon the
choice of solution methods. (For example, simultaneous
space-time solutions using conventional matrix
methods are far more costly than space-and-time
split solutions using Gauss elimination methods.)
Current models tend to require on the order of
several minutes of CDC 7600 execution time per
simulation when simultaneous space-time solution
methods are used.
Automatic compilations of arbitrary chemical
mechanisms, velocities, and diffusion profiles
are possible, but are not yet in widespread use.
N = 20, L = 50, M = 150. (We assume there are
50 chemical rates, 50 velocity values, and 50
diffusion coefficient values at the respective
grid points, resulting in M = 150.) Initial and
boundary conditions will be considered in later
examples.
65
-------
C_. 3-D Regional Air Quality Models
This base case assumes that:
Thermal and photochemical rates are constant over
intervals of one hour or longer
The spatial grid is 25 x 25 x 5, for a total of
3125 elements.
N = 10
t The system has 25 chemical rates, of which 5 are
photodi ssociation.
t 3 velocity components (measured at 3 hour intervals)
Diffusion is characterized by horizontal and vertical
components KH and KV,
e The model treats 5 emission factors (NO , CO, 3
J\
hydrocarbons), measured hourly.
Initial values are available for all species, the
temperature, and the relative humidity (R.H.).
Boundary values are available at all surfaces.
0 While numerical integration costs are highly
dependent upon choices of solution methods, current
models of this type (which split space and time
integrations) require approximately 30 minutes of
CDC 7600 execution time for a 24 hour simulation.
CHEMICAL MODEL EVALUATION
A. The Time-Honored Method
(i) Sensitivities
What can be simpler than merely altering the value of
any parameter p-, while holding the remaining parameters
J
at their nominal values, and directly noting the resulting
variations in all solutions, 6c-(t) for itl,..N, con-
tinuously in time? This is the essence of the time-
66
-------
honored method. The performance characteristics of
this method are:
Statistical properties of parameter distributions
are not required.
Values of 6c-(t) are accurate to all orders.
Values of 6c-(t) can be obtained with equal- facility
for time-varying rate coefficients as for constant
r^tes. (The cost of diurnal kinetics runs can exceed
constant parameter runs by factors of 5 to 50.)
A minimum of (M + 1) kinetics runs are required--
one run for the nominal parameter set and one run
for each rate which is varied. (Assume all M are
to be varied.) In practice, several trial values
may be required in order to locate the sensitive
range of some of the parameters, which frequently
leads to a cost of ^ 2M kinetics runs. For 20 species
and 50 rate parameters, the cost of 2M runs would
be approximately $5 x 2 x 50 = $500. For 50 species
and 100 rate parameters, the corresponding cost
would be $80 x 2 x 100 -$16,000 (for testing all
100 parameters).
Potentially new discoveries or unknown rates can
be tested by including hypothetical reactions in the
basic mechanism and then determining values of their
rates which would be required in order for the
hypothetical reaction to become significant in the
overall mechanism.
Whenever new nominal rate values are used or new
reactions are included in a mechanism, the sensitivity
analysis must be repeated in toto.
67
-------
Sensitivity coefficients can be obtained as
Z. . = Urn A-c-i ,
1,J Apj+0 Apj '
A2C.
l.(2} - lim 1- , etc.
1 >Jk Apj ,Apk->0 Ap,Apk
Recall, however, that it is very poor economy to
evaluate sensitivity coefficients as a means of
obtaining 6c.j (see Eq. 2.3) to high-order accuracy.
It is far more effective to get Ac-, directly by simply
introducting any desired perturbation Ap- and using
J
the time-honored method.
Expectations of sensitivity coefficients,
, can be obtained by accumulating a huge
1 , K
number of outputs from kinetics runs in which all
permutations of simultaneous rate perturbations
have been introduced. This requires on the order
r ]M
of [ Yj kinetics runs where y is the number of
trial perturbations to be applied to each para-
meter. If only two trial values are used for each
of 50 parameters, |2| kinetics runs would be
required, which is clearly out of the range of
p r a c t i c a 1 i ty .
(i i) Uncertainties
0 Variances a- can be evaluated from data accumulated
in the same fashion as was mentioned immediately above
In this case Stolarski et aP15' have found that the
" i M
number of samples is much less than JY| '' Tn1's wil1
be discussed further below under Monte Carlo Methods.
68
-------
Best case and worst case estimates for 6c1- can be
evaluated by directly applying selected combinations of para-
meter perturbations. The choices of parameter perturbations
are arbitrary, and the process of finding best and worst
case parameter combinations is not unique; exhaustive trials
appear to be required.
B. Variational (Direct) Methods
These methods are most appropriate for evaluating sen-
sitivity coefficients Z- *. The performance characteristics
i »J
for first-order sensitivity computations are:
The sensitivity equations (3.3) can be automatically
compiled and solved for arbitrary chemical mechanisms
The solutions can be obtained either simultaneously
with the chemical kinetics equations or decoupled
from them. Decoupling must be carefully performed
and thoroughly validated against the simultaneous
case; this is the subject of current research.
The costs for decoupled sensitivities are about
3
the same as a kinetics run, which goes as N . This
o
amounts to (MN ) for a system of M parameters
plus one kinetics run at nominal parameter values
3
which gives a total cost proportional to (M+1)N .
For 20 species and 50 parameters, this implies a
cost of approximately $250; for 50 species and 100
parameters, the approximate cost is $8000. The
costs of solving the sensitivities and kinetics
o
simultaneously is proportional to M(2N) for all
its parameters, or about 8 times as expensive as
the corresponding decoupled case. The approximate
costs of the simultaneous case are $2000 for 20
69
-------
species and 50 parameters, and $64,000 for 50 species
and 100 parameters. Clearly the redundant kinetics
solutions which characterize the simultaneous case
should be avoided.
Costs of evaluating systems with time-varying
coefficients are inordinately greater, both in
machine execution time and in dedicated programming
of functional derivatives.
Whenever new nominal values or new reaction sets
are used, the sensitivity analysis must be reperformed.
Expectations of sensitivity coefficients
i »j
are not practical to evaluate. Estimates of 6c- and
Oj based on first-order sensitivity coefficients have
on "y local validity (see Section 1.)
C. Green s Function Method
This r-ethod is similar in many respects to the variational
(direct) nr^thods just discussed, but has not been tested or
applied as extensively as yet. Key performance characteris-
tics of th is method are as follows;
The homogeneous differential equations (3.30) for the
Gre -en's function are easier to solve (in terms of
cor-,putation time) than are the sensitivity equations
(3. 3) of the direct method, and can be automatically
corn-foiled and solved for arbitrary chemical mechanisms.
The numerical solution of (3,30) must be reinitialized
rep-- eatedly, adding some expense to the computation.
Add itional expense is incurred in computing the sen-
sit- ivitiesZi k(t) from the Green's function by
qua_ drature (trapezoidal rule or' equivalent). These added
cos- ts, however, should be smaller than the savings
aqh- ieved by decoupling the sensitivity equations from
the kinetics equations.
70
-------
Whether or not the decoupling itself can be done
without losing accuracy remains to be proven for
reasonably complicated kinetics systems.
If desired, higher order sensitivity coefficients can
be obtained at far less cost than by the direct
method.
Additional programming efforts totalling
perhaps one person-year should allow the develop-
ment of a user-ready package which could resolve
these questions.
D. The Monte Carlo Method
This is the extension of the time-honored method
which computes variances of model output data. Performance
characteristics are as follows:
t This method is conceptually transparent; it is
non-dedicated; and it applies to one time point
at a time.
["Y] kinetics runs are required as input, where
Y is the number of trial perturbations applied to
each parameter. For j=2, the following dollar
costs would be incurred for (constant parameter)
kinetics runs alone:
20 Species
50 parameters - $ 5.6 x 10^ ^ 1nf1n1te
100 parameters - $ 6.3 x 10 ,
50 Species -,
50 parameters - $ 9 x 10" ; ^ 1nf1n1te
100 parameters - $ 1 x 10 |
For all practical purposes, these costs are infinite,
suggesting that more selective parameter sampling
71
-------
must be used. One alternative would be the FAST
sampling scheme. So far as we can determine, this
possibility has not yet been explored in a practical
example. We will see in (E) directly below that
even the FAST sampling would be inordinately
expensive in kinetics integration for our basic
(IS)
examples. Stolarski et alv ' estimate in advance
the likely order-of-magnitude cost of variance
evaluations. They apply the time-honored method,
using trial perturbations of +0 and +_ 2 a for all
rates, in order to make a preliminary identification
of the most sensitive parameters in the system.
(Recall that this is a local evaluation, varying
one rate at a time about nominal parameter values,
which is good to arbitrary order and is continuous
in time.) In the event that the kinetics mechanism
is sensitive, say, to only 10 rates out of the
original set of 50 or 100 rates, then only (2) or
1024 kinetics integrations are required as input
for Y=2. In this event the cost of the input for
computing variances would be $5120 for 20-species
systems and $81,920 for 50-species systems. This
approach is now a naturally-evolved hybrid of sen-
sitivity and uncertainty analysis. This practical
alternative will prove useful for all of the systematic
methods . Whereas in the Monte Carlo method, prescreening
only serves to estimate execution costs, in the FAST
method, prescreening followed by rejection of insensitive
parameters may result in considerable savings of input
samples without compromising the basic results. This
point should be further pursued andeither validated
or invalidated by FAST users.
In cases where prescreening is used for parameter
rejection,the hybrid approach exacts a cost elsewhere,
72
-------
namely, some of the promise which systematic methods
hold for a more rigorous search of model solution sur-
faces is compromised by computer economics,
t The entire analysis must be redone if any nominal
values change. In some cases, changes in the magnitude
of parameter uncertainties (probability distributions)
may also require that the entire analysis be redone.
Sensitivities tested with perturbations of +a do
not sample the "corners" or extreme regions
of the solution surfaces. Best case and worst
case examples are more likely to escape notice
in these circumstances.
0 Stolarski et al found that output variances were
essentially independent of the insensitive parameters.
But the rate of convergence to reliable variances
was dependent upon the degree of sensitivity of
particular variables in various time (and space)
domains. Hence the user of the Monte Carlo method can
not know a priori the number of input runs which are
required to obtain convergence of variance values,
E. Fourier Amplitude Sensitivity Test
The performance characteristics for first-order
sensitivity coefficients are:
The FAST is a non-dedicated method which can be
applied to arbitrary model outputs at one time
poi nt of the model.
The number of kinetics integrations which are
required as input to the FAST is approximately
(M)2-5. Therefore, approximately 17,700 kinetics
runs are required for 50 parameters; and 100,000
73
-------
Nineties runs are required for 100 parameters.
The following dollar costs would be incurred for
(constant parameter) kinetics runs alone:
20 Species
50 parameters - $ 88,500
100 parameters - $ 500,000
50 Species
50 parameters - $1.42 x 10
100 parameters - $8. x 10
t It is clearly necessary to minimize the number of
parameters under consideration as well as to minimize
the cost of reliably solving kinetics systems.
A conundrum arises. The promise of systematic
methods is to search the model solution surfaces
more thoroughly in order to detect subtle as well
as obvious reaction coupling influences. Increasingly
thorough searching, in turn, suggests that more
and more complete reaction mechanisms should be
put to the test. But the computer economies argue
strongly in the opposite direction, thereby thwarting
some of the promise of systematic methods which
happen to be expensive.
t Variances can be obtained by two means:
(i) The output of the numerous kinetics runs which
are required by the FAST can be processed directly.
It is known that the parameter space is sampled non-
uniformly, depending upon the specific choice of
FAST distribution functions. The effect of sampling
choices upon such output variances has not yet been
assessed, but is likely to be significant for sen-*
sitive parameters.
74
-------
(ii) The FAST obtains the sensitivity measures
from partial variances a (p -) =a . (w . t)/a. . The
»J ' J ' J » i
terms c-(u>. t) are defined in terms of the Fourier
1 j j
amplitudes at frequency w. by
J
In practice only a finite number of to. and associated
J
harmonics are actually used, so that the resulting
variance values may or may not be reliable.
Unaveraged sensitivity measures Z- (analogous
> J
to measures obtained by the deterministic methods)
can be obtained in principle by using distributions
in the FAST of zero width (zero uncertainty) for
all Pi^i* but this would be very inefficient relative
the other sensitivity methods for kinetics systems.
In principle, second-order and perhaps higher-order
sensitivity information is available from the FAST.
In practice, this capability does not yet appear to
have been implemented for computations by general users
As discussed above, this additional information
will not provide a more effective means of obtaining
improved variances vis a vis other alternatives;
and the application of higher-order, averaged
sensitivity information appears to be a very low
priority need of most practical users, due to
problems of conceptual transparency and ease of
application. (It is our impression that this is
the case for higher-order sensitivity information
from all existing methodsj
75
-------
The sensitivity measures depend in general
1 5 J
upon the choices of distribution functions, as was discussed
fl.8)
above and also found by Penned ; in applications. The
sampling methods thus tend to be best-suited to problems in
which probability distributions are well-known for key parameters
1-D TIME DEPENDENT MODEL EVALUATION
From the evaluation of chemical models just above,
it is clear that computer economics will dominate our
evaluations of sensitivity methods for 1-D and higher-
dimensional transient air quality models. The cost of
integrating multi-dimensional transient models can be
highly variable and is particularly dependent upon the
numerical methods and programming efficiency in any specific
air quality model. To provide a basis for subsequent
evaluations we will first discuss briefly a set of well-
known numerical solution options which are used in multi-
dimensional transient models and then indicate approximate
scale relationships of the costs of these solution options
to the costs of purely chemical models which were previously
stated.
Transient 1-D chemical-transport models are generally
solved simultaneously in (.x,t) by employing the method of
lines to discretize spatial gradient operators. An implicit
stiff ODE integration method is then applied to solve the
overall system of size NL x NL| . When no special numerical
advantages are exploited, the cost of solving such a 1-D
model is proportional to rNL|a, where a has a value between
2 and 3. When central differences are applied to gradient
operators, resulting in an overall matrix having a block-
tridiagonal structure; Gauss elimination solutions of the
matrix can be used which lead to reduced solution costs on
the order of LNa.
76
-------
JBy similar analysis, simultaneous solutions in (x..y_»t)
in 2-D transient models (e.g., LIRAQ) have solution costs
ranging from L3Na to [L2Nj a. (For the remainder of this
discussion we can assume a ^ 2-3, but the structural aspects
are clearer if we Scale costs in terms of a.) The use of
ADI (alternating direction implicit solution) methods for
2-D transient models reduces the numerical solution cost
to the range of 2L2Na to 2LJLN]a. In 3-D air quality
models (e.g., the Systems Applications, Inc., regional model), ADI
is absolutely essential. Solution costs for 3-D air quality
f\
models would then scale to span the range of 3L Na to
3L LNJa (for Gauss elimination versus conventional matrix
solution options).!
For estimation purposes, we will assume that 1-D models
are solved simultaneously in (_x,t) at a cost of 5 minutes
of CDC 7600 time, or $100 per simulation. Again, this
estimate is representative of existing transient 1-D models
with which we are familiar; the main thrust of our con-
clusions will not be greatly affected by this particular
choice of cost estimate vis a vis execution costs of other
such 1-D transient atmospheric models actually in operation
today.
A. The Time-Honored Method
Assume that, in addition to the single model solution
which is required using nominal parameter values, only
one additional run will be required for each parameter in
order to evaluate 6cn-(x,t). (This might be the case for a
I "~~~
very experienced practitioner.) A minimum number of para-
meters would be M = 150, which is composed of:
50 chemical rates which are constant in time and
have the same value in every zone. Assume that the
same rate perturbation applies simultaneously to
every zone. '
77
-------
t 50 diffusion coefficients which are uncorrelated
and take specific values at each zone. Assume
that a single fractional perturbation at each
zone is appropriate over the entire time history.
50 wind velocities which are uncorrelated and
take specific values at each zone. Assume that
a single fractional perturbation at each zone is
appropriate over the entire time history.
In practice, many other parametric variations may
also be considered, including
t pollutant emission sources at each zone as a
function of time;
t boundary values of species concentrations or
species fluxes as a function of time;
initial conditions at every zone; and
photodissociation rate variations from zone-to-zone
due to variable cloud cover (even for constant rate
coefficients).
On this minimum basis, the cost of 1-D air quality
model runs alone would be (151 x $100) or $15,100 of com-
puter time. (The reader can readily make similar estimates
which are germane to his own computing circumstances and
experiences.) Major improvements in computer economy could
be realized from possible advances in four areas. First,
improved numerical techniques such as matrix solutions by
Gaussian elimination would probably bring the computational
cost under $1000 for 150 model runs. Second, intuitive pre-
screening for obviously insignificant sensitivities can reduce
the number of parameters to be varied; the utility of this
option depends strongly upon the craft of the user. Third,
correlated variables can sometimes be identified and varied
78
-------
simultaneously rather than individually. Fourth, optimized
programming, using automatic compilation of arbitrary problem
scenarios, machine coding, parallel processing, and other
specialized program optimization routines, can result in
(15)
substantial savings. For example, Stolarski et aP ' found
it necessary to write a completely new stream-lined production
version of their original 1-D atmospheric model and further,
to reduce it to a steady-state model rather than a transient
model. However, this latter option would probably not be
appropriate for multi-day air quality simulations which are
intrinsically transient in nature.
B. Variational Methods
Using the same basic example as above, several operational
possibilities can be dismissed immediately. For example,
simultaneous solution of the sensitivity and 1-D model equations
is out of the question. (Recall the discussion for Chemical
Models.) Similarly, the sensitivity equations would be solved
for one parameter at a time rather than solving the larger
matrix for, say, all M parameters simultaneously. Further,
automatic compilation for arbitrary problem scenarios is
essential for both the sensitivity programs and 1-D air
quality models. (Can you imagine the difficulty of writing
a new dedicated sensitivity program which corresponds per-
fectly to each air quality model and/or to each specific
application which is of interest?)
Following these suggestions and assuming that solutions
of the sensitivity equations would cost the same to obtain
as solutions of the 1-D air quality model, one would be solving
a system of size [ML x Nil (M+l) times. A practical approach
to calculating sensitivities for this problem (with constant .
chemical rates) might be to pursue the following lines.'
79
-------
For a given set of nominal v° and KR at each grid
point, solve the sensitivity equations MR times for
the reaction rates of interest. This gives Z- -L.
at all grid points continuously in time.
For a given set of nominal chemical rates, solve
the sensitivity equations Mv and M,, times for the
- =r>
velocities and diffusion coefficients at selected
grid points of interest. (v^ and KD may be constant
over, say, hourly intervals; this would pose no
difficulty to the direct methods.) This procedure gives
Zi,J Ivel. and Zi,J 'oiff. at a11 System 9rid
points, continuously in time.
o For a given set of nominal values for chemical rates,
winds, and diffusion coefficients, the sensitivities
Z- . can also be determined at all grid points, con-
's J
tinuously in time, for selected initial, boundary, and
source conditions.
So far as we can determine, no sensitivity analysis of
this level of generality has yet been performed with direct
methods. (Koda et al^ have considered much more limited
examples, such as non-reactive atmospheric diffusion, which
could be solved analytically.) In the event that M=150 in
the general case (due to any combination of values for Mn,
M , M., , and other parameters), the computer cost would be
D
identical ($15,100) to the minimum case discussed for the
time-honored method immediately above. Notice that the
possibility of using intuitive or other judgements to eliminate
obviously insensitive parameters from the systematic analysis
is tacitly left open, recognizing:
(i) the element of risk in such eliminations, and
(ii) the fact that some of the promise of systematic
methods is compromised to computational economy by such
eliminations.
80
-------
The functional derivative approaches for either time-
varying rates or non-discretized spatial variables are severely
limited by the following factors:
0 The methods are too highly dedicated, at present,
for multi-dimensional application.
1-D model execution costs increase dramatically
for diurnally varying rates and continuously varying
emissions sources.
0 We estimate execution costs for solving the air
quality and the adjoint sensitivity equations to
be on the order of solving an [ML x Nil system
of ODE's for the air quality equations plus an
[N2L2 x N2L2] system of ODE's for the adjoint
equations. Thus the execution cost would scale
r o o^
approximately as [N Lja, with a relatively large
proportionality factor applying in the case of
variable rates, sources, etc. In addition, N
double integrals at,say, L values of x^ would have
to be evaluated at each time point for each chemical
rate. Likewise, N double integrals for each
pair (_x' , xj and each time point would have to be
evaluated for sensitivities to winds and diffusion
coefficients. We conclude that, vis a vis other
methods, this approach is far too costly and difficult
to apply in air quality modeling for the foreseeable
future.
C. The Monte Carlo Method
For 150 parameters, (2) runs of the 1-D model is,
practically, an infinite input cost. To date this method has
been used primarily to obtain steady state atmospheric
model variances associated with uncertain chemical rates.
Even in the extremely unlikely event that only 10 parameters
were found (by some other means) to be sensitive, 1024 runs of
a 1-D transient model would be on the order of $100,000 in
computer execution cost. This is the reason that
f 15)
Stolarski et a 1 emphasize the need to take every
81
-------
reasonable step to reduce the execution cost of the atmos-
pheric model which generates the basic data for the Monte
Carlo method of evaluating variances. In the event that
20 sensitive parameters are present in an air quality
model, over a million input samples would be required to
evaluate variances, which is out of the question if air
quality model execution costs greatly exceed 1/100 cents
per run.
D. The FAST Method
The input cost of testing 150 parameters would be
% (150)2'5 % 275,500 runs of the 1-D air quality model,
costing on the order of $27.5 million for computer execution.
The non-dedicated features of this method remain very
attractive, but the number of parameter trials must be
markedly reduced either by pre-screening eliminations or
by considering only segments of the overall problem (such
as testing the chemical mechanisms independent of the
fluid dynamics or vice versa).
3-D AIR QUALITY MODEL EVALUATION
The general picture is abundantly clear by now. If
we consider only a 24-hour simulation with the 3-D model at
a cost of approximately $600 per run, the computer costs
relative to the 1-D examples present a grim picture, indeed.
To carry through this discussion, Tables 2 and 3 have
presented what we havetermed an "Oversimplified Case" and a
"Realistic Case" for the number of parameters which must be
considered in our basic 3-D air quality model.
The "Oversimplified Case" is simplified beyond
scientific usefulness. The "Realistic Case" is scientifically
acceptable but is well out of the range of practicality.
The results in Tables 2 and 3 are self-explanatory.
82
-------
TABLE 2
SUMMARY TABLE
BASE CASES FOR ILLUSTRATION OF EXECUTION COSTS
OF CANDIDATE SYSTEMATIC SENSITIVITY ANALYSIS METHODS
BASE CASE
DESCRIPTION
ASSUMED EXECUTION
COSTS (on basis of
"typical" CDC 7600)
COMMENTS
PURE CHEMICAL
KINETICS MODELS
(A.) 20 species, SO
parameters
(B.) 50 species, 100
parameters
$5/klnetlcs
simulation
$80/klnetlcs
simulation
The paramatera
are assumed to be
constant. Diurnal
and other time
dependent effects
would be treated by
using piecewise
constant parameters
over intervals of an
hour or longer.
1-D TRANSIENT MODELS
20 species, 50 grid
points.(Hence there
are 50 diffusion
coefficient values
and 50 wind values on
the basis of one value
per grid point and
uncorrelated y_(x) and
KD(x).)
$100/1-0 simulation
Parameters are
assumed at least
piecewise constant
(as above).
3-D REGIONAL AIR
QUALITY MODELS
10 species,
3125 grid points (based
on zoning which Is
dimensioned as 25 x 25
x 5),
25 chemical rates,
3 velocity components,
"Wir. and "Srert.!
5 emission factors (NO ,
CO, three reactive
hydrocarbons, measured
hourly),
Initial values for
species, temperature,
and relative humidity.
Boundary values at all
surfaces.
$600/3-D simulation
A. Over-simplified Case. 137 parameters, assuming
that:
v_ has identical values at all grid points.
* ^H* ^v nave identical values at all grid points.
No cloud corrections; thus J's do not vary from one
grid point to another.
k's do not vary from one grid point to another
jCconstant T, R.H. at all grid points)
Initial conditions for 10 species, T, and R.H. are
identical at all grid points.
Boundary conditions of v_, KU»KV. and 10 species are
identical over each surface. (15 components x 6
surfaces - 90 values)
B. Realistic Case. *\» 234,200 parameters, assuming that:
_v la variable at each grid point.
KH and Ky are variable at each grid point.
5 Jfs, which are variable at each grid point.
20 k's, which are variable at each grid point.
Initial conditions are variable at each grid point.
Boundary conditions are variable at each grid point.
83
-------
TABLE 3
SUMMARY TABLE
ILLUSTRATIVE COMPUTER EXECUTION COSTS FOR CANDIDATE SYSTEMATIC
SENSITIVITY ANALYSIS METHODS (REFER TO BASE CASES IN TABLE 3.1(a))
METHOD
TIME-HONORED
(6Ci(t))
DIRECT VARIATIONAL
(Z, ,(t))
i , J
FAST
(^ ,-(!')>)
1 »J
MONTE CARLO
(O.j(t'))
PURE CHEMICAL
KINETICS MODELS
(A. ) $500 (Assumes
two trials
per para-
meter)
(B. ) $16,000
(A.) $250 (Assumes
sensiti-
(B.) $8000 argyd^S'
coupled
from kine-
tics egs,)
(A.) $2000 (As^H6?
(B ) $64, 000*1 "-V"*
eqs. are
coupled)
(A.) $88,500 (Assumes
, no pre-
(B.) $8.10b screen-
ing)
(A.) $5.6 x 1015
(B.) $l.xl032 $82-200(assume
B One
(B.) $1.4xlOB trial /
para-
meter)
(13 yrs CDC 7600
-------
SECTION 5
TODAY'S APPROACH
The central objective in this section is to suggest
the means by which the EPA, today, could most effectively
apply systematic sensitivity analysis methods to arbitrary
air quality models, which we will assume to be transient
2-D or 3-D regional or reactive plume models. Our suggested
choices of method will be guided to a significant extent
by the conclusions which have evolved from our evaluations in
theprevioussections.
CONCLUSIONS
A. On balance, the full promise of systematic sen-
sitivity analysis methods cannot yet be realized in applications
to entire 2-D and 3-D transient air quality models. However, by
segmenting the models or by prescreening for sensitive para-
meters, very useful sensitivity information can be gained with-
out si gni fi cant alterations being required in the air quality
models.
B. A combination of time-honored and systematic sen-
sitivity analysis methods can be effectively used on individual
air quality model segments; for example just the chemistry can
be analyzed in arbitrary zones with the fluid dynamics remaining
decoupled, and conversely the chemistry can be frozen while
analyzing the fluid dynamics, the boundary conditions, the
initial conditions, sources, etc. (Without question pre-.
screening and segmenting compromise some of the promise of
systematic methods in exchange for computational economy;
this compromise is necessary until further advances occur in
systematic sensitivity analysis and/or in the efficiency
of execution of air quality models.)
85
-------
C. Reliable variances appear to be out of practical
range for multidimensional transient models, with a multitude
of sensitive parameters.
D. Dedicated sensitivity methods must be programmed
so that arbitrary problem scenarios can be automatically
compiled in a hands-off manner; otherwise, non-dedicated
methods are a far more practical alternative for analysis
of complex air quality models.
E. For better or for worse, practical users favor con-
ceptually transparent sensitivity methods over those which
require more subtle interpretation without providing com-
mensurately more useable information. (For example, higher-
order sensitivity coefficients hold important topological
information about complex solution surfaces; but the volume
of information which is needed to accurately map complex
surfaces grows quickly to unmanageable and conceptually
obscure proportions.)
F. Numerical efficiency of air quality models and
sensitivity methods assumes ever-increasing importance when
extensive sensitivity studies on the models are desired.
(Most air quality models would require considerable re-
programming to reach their maximum attainable efficiencies
for production purposes.)
6. Depending on specific circumstances, field or
laboratory measurements maybe a preferred alternative to
sensitivity analyses. At one extreme are the chemical models
where sensitivity analysis is very cost-effective for
indicating the most critical reactions which need to be
measured in the laboratory in order to markedly improve a
chemical mechanism. At the other extreme are 3-D regional
models where the costs for even a si ngle sensitivity analysis
86
-------
may easily exceed $100,000 for computer expense, alone.
But we have seen that, in addition, many such analyses
must be done to account for alternative parameter choices
which are furthermore subject to many permutations. In
this case, comprehensive uncertainty analysis of the 3-D
model is unreasonable and should be supplanted by more
meaningful experimental work. Between these extremes,
however, more subtle tradeoffs must be very carefully
considered on an individual case basis.
TODAY'S PRACTICAL APPLICATION OF SYSTEMATIC SENSITIVITY ANALYSIS
Variances of air quality model outputs can be economically
and reliably obtained for only very special circumstances
(such as models with less than about 10 sensitive variables
and for highly specialized models which require only seconds or
fractions of seconds of CDC 7600 execution time
(or equivalent) per simulation). This conclusion is
equally true for air quality model segments: For example,
variances from a purely chemical model (or submodel), in-
dependent of other segments of an air quality model, can
be reliably obtained only for cases in which about 10 or
fewer chemical rates are sensitive parameters. We thus
believe that policy-making applications will be best served
by shifting the emphasis away from studying variances of
limited scope and validity and, instead, to increased
understanding of the physical and chemical cause-and-effect
mechanisms with the aid of hybrid methods of sensi ti vi ty analysis,
Today's best approach to sensitivity analysis of complex
air quality models would first decompose the subject model
into a chemical component and a hydrodynamics component,
noting several additional factors which apply to the model
segmentation process:
87
-------
(i) Segmenting is most effective when it is applied
at the points of weakest mutual coupling in the air
quality model. The chemical-hydrodynamics decoupling
is an obvious first choice, based on present experience;
improved decoupling points should always be sought.
(ii) Carefully considered application of model segmenting
requires--and returnsgreater understanding of the chemistry
and physics in a given model.
(iii) If incorrect assumptions are made in model seg-
menting, quite often the result is that more--not less--
is learned about both the air quality model and the
simulated air quality problem.
Analysis of Sensitivities in the Chemical Component
The chemical mechanism of any large air quality model
can be readily simulated by stand-alone automated chemical
kinetics solvers such as those described in reference 17,
Arbitrary conditions (e.g.,initial conditions, species
fluxes, dilution factors, and pollutant sources and sinks)
can be abstracted from any spatial zone and time domain of
the subject air quality model and provided as input for
solution by the more efficient purely chemical kinetics
solvers. This can be done repeatedly for representative
conditions and circumstances in the original air quality
*
model.
Photochemical sensitivities can be determined most
economically if rate parameters are at least piecewise
constant in time (over one-hour intervals, for example)
rather than instantaneously varying on a diurnal basis.
* This process of abstracting such conditions from a more general air
quality model to a more specialized (photochemical) model is known as
"linking". Linking has repeatedly been demonstrated to be effective
in validating complex computational models in numerous scientific
disciplines.
-------
(Instantaneous diurnal variations exact such an extreme
cost in kinetics integrations that all sensitivity methods
can become unduly expensive.) Recognizing that the
application of a kinetics sensitivity method is very much a
matter of user preference, we will nevertheless indicate
as a guiding suggestion that our present choice would be
to use a direct variational method on the chemical com-
ponent for several reasons:
9 Economy of effort and computer cost. If decoupled
from the kinetics solutions, a suitably automated
direct sensitivity method is more economical than
the other alternatives in terms of human effort and
computer cost, particularly for new applications
to large chemical mechanisms. Data handling is
minimized as are arbitrary operational decisions
such as trial and error of selection of rates, rate
distributions, etc.
Conceptual transparency. The measures 6c-(_x,t)
of the time-honored method and Z- .(t) of the direct
' »J
method are both obtained continuously in time and
are independent of statistical assumptions, unlike
the measure of the FAST. The time-honored
> J
variable 6c-(x,t) has an advantage of being valid
to all orders. From experience with all of the
alternative sensitivity methods (except for the
89
-------
Green's function method which is not yet imple-
mented for general use) we have usually found that
the automated direct method gives the basi c sen-
sitivity information which is needed in a more
organized and efficient, hands-off manner than the
alternative methods. (Again, individual user pre-
ferences may very well differ from ours.)
Potency. Direct variational methods come closest
to fulfilling the full promise of systematic sen-
sitivity analysis methods in chemical applications.
No pre-screening (usually done by trial and error)
is required in atmospheric applications of direct
methods, whereas pre-screening is advisable for the
FAST when more than about 20 parameters are to be
tested. Pre-screening should always be minimized
to the greatest possible extent because any over-
sights in this stage will compromise all subsequent
evaluations. The measures Z- . and 6c-(x,t) of the
1 > J i
direct and time-honored methods can also be used for
best and worst-case identification, whereas
' »J
from the FAST suppresses such information by the
averaging process. By reverse reasoning, the FAST
is the most effective method for providing averaged
sensitivity coefficients in those special instances
where the averages have meaning. It may also occur
to some that the FAST can effectively compute Z-
> J »
as well as , by simply assigning zero uncer-
1 > J f h
tainties to all except the j~Ln parameter. But we
have seen that, in the absence of new advances, this
apparent option of the FAST is grossly inefficient
relative to the direct variational and time-honored
alternatives, particularly for moderate to large-
size chemical mechanisms.
90
-------
Ease of Application. A majority of practical users
prefer the trial and error approach and tend to resist
more formal sensitivity analysis methods which require
their own numerical programs. The great intuitive appeal
and immediacy of the time-honored approach ensure that
it will continue as a favored sensitivity method. On
the other hand, direct methods have now been developed
and automated to a sufficient degree^' ^' * ' for all
types of chemical model use that widespread application
should be near-at-hand. Koda et al^' are similarly
reducing the FAST method to semi-automatic general
applicability in atmospheric chemistry problems.
The sensitivity testing of the purely chemical component
is likely to be most meaningful if the systematic analysis (by
the direct method) of rate and initial condition uncertainties
is augmented by independent trial-and-error assessments of
factors such as instantaneously varying photodissociation rates,
more detailed characterizations of time dependent sources and
sinks, and any other factors which were either by-passed or
treated approximately by the direct method.
Analysis of Sensitivities in the Hydrodynamics Component
Whereas automated chemical kinetics solvers can be used
efficiently to assess essentially unlimited chemical scenarios
without running the full air quality model, the same is not
true for the hydrodynamics component. The problem is that,
even for non-reacting fluids, independent hydrodynamics solvers
are not feasible for higher than one-dimensional transient
models with arbitrary zoning, sources, sinks, and wind and
diffusivity resolutions. Similarly, direct sensitivity
methods have not yet been effectively automated for arbi-
trary scenario evaluations in higher than one-dimensional,
transient hydrodynamics. It is therefore necessary for the
91
-------
original air quality model, or some trivial derivative of it,
to be used to generate necessary output data for sensitivity
analyses of hydrodynamics factors. For example, the chemistry
can be turned off by assigning zero rates to all chemical
reactions in order to study the purely hydrodynamic (meteoro-
logical) characteristics of the air quality model.
With extensive basic understanding of the physical
problem, the time-honored method would then be used almost
exclusively to identify the key variables and their basic
sensitivities to space and time dependent meteorological
parameters, sources, sinks, initial conditions, and boundary
conditions. The time-honored method would serve also as a
means of pre-screening insensitive variables before applying
the FAST to small groups of carefully selected hydrodynamics
parameters. Correlated parameters (i.e. those which are re-
presented by common functional relationships at all spatial
zones) may also be treated effectively by the FAST in particular
circumstances. Clearly, basic understanding of the physical
system, coupled with sound intuition, is the most important
element for incisive sampling and successful analysis of
hydrodynamic sensitivities; for it is simply impossible to
test every potentially sensitive parameter at every spatial
grid point with any of the presently developed sensitivity
methods .
NEEDED BREAKTHROUGHS IN SYSTEMATIC SENSITIVITY ANALYSIS METHODS
Breakthroughs which would provide truly major advances
in sensitivity analysis applications fall into two major
categories; first, there are general advances which would
benefit al1 specific sensitivity and uncertainty analysis
methods, and second, there are specialized advances which
will benefit, exclusively, a particular sensitivity method.
92
-------
General Advances
Model execution speed. We have seen that several
sensitivity methods would become more attractive if the cost
of air quality simulations were on the order of 1/100 cents
per run. However, this possibility appears to be out of
reach for even perfectly optimized 2-D and 3-D transient
air quality models with envisioned computers of the.next
generation or two.
Advances in partial differential equations.
Significant advances in PDE numerical analysis are likely
over the next 5-10 years which may dramatically improve air
quality modeling. First, numerical integration of non-
linear systems of PDE's is.in the midst of major advances.
Finite element methods with adaptive gridding features '
have the capacity to use far fewer grid points, located at
optimal positions. These features tend to eliminate numerical
diffusion and the overshooting errors which have always
plagued the finite difference methods which are used in
current-generation models. These new approaches are, additionally,
very amenable to automatic compilation of arbitrary problem
scenarios in higher dimensions, in a manner which is very
much akin to the automated kinetics (ODE) solvers which were
mentioned in previous sections. These general advances un-
doubtedly will be many years in coining to fruition, unless
special initiatives are undertaken.
Correlated variables allow significant reductions in
execution costs for all systematic sensitivity analysis methods.
The identification of correlations in winds and diffusivities ,
for example, will most likely come from very fundamental considera-
tions of fluid dynamics. Such studies are far more basic and
rigorous than most air quality models, per se, and will thus
require many more years of highly detailed theoretical and
experimental analysis.
93
-------
Specialized Advances
The FAST method would be markedly improved by more
efficient parameter sampling which relates to the selection
of incommensurate Fourier frequencies.
The FAST method would run markedly faster if
correlated parameters can be identified and represented by
a common parametric variation, which in turn reduces the
required number of independent sample points.
0 The direct variational methods are markedly improved
when the sensitivity computations are decoupled from the
kinetics or hydrodynamics computations. A very modest amount
of research (about l person-year) will now accomplish the
development of reliable decoupling methods.
The direct variational methods will advance markedly
when 2-D hydrodynamics models can be automatically compiled
for arbitrary problem scenarios and solved by generalized
hands-off PDE solvers (such as is presently the case with
automatic kinetics solvers which use the Gear ODE solution
method^1* ^8'). Additionally, the sensitivity equations
must be amenable to the same measures. At Science Applications,
Inc., in Pleasanton, California, in-house efforts are, in fact,
remarkably near (within 3 person-years) of accomplishing this
advance using very potent finite element methods. (We are
aware of only one other such effort, in a different field,.
which is near a comparable point of advancement, so it is
difficult to give a more universal status estimate than the
one above which is based on our own work.)
Finally the Green's function method may ultimately
offer some advantage over alternative sensitivity methods
for purely chemical systems. It is still too early to tell.
(We estimate that 1-2 person-years may yet be required to
satisfactorily implement this method in a suitable' automatic
format for arbitrary mechanism evaluations.)
94
-------
REFERENCES
1. Reynolds, S.D., P.M Roth and J.H. Seinfeld, "Mathematical
Modeling of Photochemical Air Pollution -I. Formulation
of the Model.", Atmos. Environ., 7, 1033-1061. (1973).
2. Mac Cracken, M.C., D.J. Wuebbles, J.J. Walton, W.H. Duewer,
and K.E. Grant, "The Livermore Regional Air Quality Model:
I and II". J.A.M.. 17, No. 3, 254 (March 1978).
3. Sklarew, R.C., M.A. Joncich, K.T. Iran, "An Operational Air
Quality Modeling System for Photochemical Smog in San Diego
Air Basin", Fourth Symposium on Turbulence,Pi ffusion,and Ai r
Pollution, 301, (AMS Jan 15-18, 1979).
4. Hecht, T.A., J.H. Seinfeld, and M.C. Dodge,"Further
Development of a Generalized Kinetic Mechanism for Photo-
chemical Smog", Environ. Sci. Techno!., Vol. 8, No. 4,
pp. 327-339 (197TT
5. Whitten, G.I. and H. Hogo, "Mathematical Modeling of
Simulated Photochemical Smog", EPA Report EPA-600/3-77-011
(Jan. 1977).
6. Gelinas, R.J. and Skewes-Cox, P.O., "Tropospheric Photo-
chemical Mechanisms", J. Phys. Chem. . 81, 2468 (1977).
7. Dickinson, R.P. and Gelinas, R.J., "Sensitivity Analysis
of Ordinary Differential Equation Systems", Journal of
Computational Physics, Vol. 21, No. 2, 123-143, (June 1976).
8. Hindmarsh, A.C. and G.D. Byrne,"EPISODE: An Experimental
Package for the Integration of Systems of Ordinary Differential
Equation Systems",, UCID-30112 , Lawrence Livermore Laboratory
(May, 1975).
9. Atherton, R.W., R.B. Schainker, and E.R. Ducot, "On the
Statistical Sensitivity Analysis.of Models for Chemical
Kinetics", A.J_._C h_.J. J_. , 21., 441 (1975).
10. Gelinas, R.J., "Diurnal Kinetic Modeling", Proceedings of
the International Conference on Structure, Composition and
General Circulation of the Upper and Lower Atmosphere and
Possible Anthropogenic Perturbations. Vol. II, Jan. 14-25,
1974, IAMAP/IAPSO, Melbourne, Australia, (1974).
11. Koda, M., AH. Dogru, and J.H. Seinfeld, "Sensitivity
Analysis of Partial Differential Equations with Application
to Reaction and Diffusion Processes", preprint (1978).
95
-------
12. Walter, W., Differential and Integral Inequalities,
(Springer-Verlag, New York, 1970).
13. Porter, W.A., Int. J. Control. _5, 393 (1967).
14. Hwang, J.T., E.P. Dougherty, S. Rabitz, and H. Rabitz,
"The Green's Function Method of Sensitivity Analysis in
Chemical Kinetics", preprint (1978).
15. Stolarski, R.S., D.M. Butler, and R.D. Rundel , "Uncertainty
Propagation in a Stratospheric Model, 2. Monte Carlo
Analysis of Imprecisions due to Reaction Rates", J6R,
in press (1978).
16. Cukier, R.I., C.M. Fortuin, K.K. Schuler, A.G. Petschek,
and J.H. Schaibly, J. Chem. Phys. 59, 3873, (.1973).
17. Dickinson, R.P. and Gelinas, R.J., "SETKIN", in L. Lapidus
and Schiesser (eds.), Numerical Methods for Differential
Systems, (Academic Press, 1976) , pp 167-180.
18. Penner, R.C., private communication, also see A.A. Boni
and R.C. Penner, Combust. Sci. Technol. 15, 99, (1976).
19. Doss, S., private communication (_1978),
96
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)
REPORT NO.
EPA-600/4-79-035
3. RECIPIENT'S ACCESSIOf*NO.
4 TITLE AND SUBTITLE
SYSTEMATIC SENSITIVITY ANALYSIS OF AIR QUALITY
SIMULATION MODELS
5. REPORT DATE
May 1979
6. PERFORMING ORGANIZATION CODE
7 AUTHOR(S)
R. J. Gelinas and J. P. Vajk
8. PERFORMING ORGANIZATION REPORT NO.
9 PERFORMING ORGANIZATION NAME AND ADDRESS
Science Applications, Inc.
80 Mission Drive
Pleasanton, CA 94566
10. PROGRAM ELEMENT NO.
1AA603A AB-042 (FY-79)
11. CONTRACT/GRANT NO.
68-02-2942
12. SPONSORING AGENCY NAME AND ADDRESS
13. TYPE OF REPORT AND PERIOD COVERED
Environmental Sciences Research Laboratories-RTP, NC
Office of Research and Development
U.S. Environmental Protection Agency
Research Triangle Park, N.C. 27711
Final 1/78 - 12/78
14. SPONSORING AGENCY CODE
EPA/600/09
15. SUPPLEMENTARY NOTES
16. ABSTRACT
Systematic sensitivity and uncertainty analysis methods are reviewed and
assessed for applications to air quality simulation models. The candidate methods
in terms of their basic variables, mathematical foundations, user, motivations and
preferences, computer implementation properties, and costs (including both human and
computer resources) are discussed. Both deterministic and sampling methods have been
evaluated with illustrative examples of "What you pay" and "What you get" with present
generation systematic sensitivity analysis methods and air quality models. Determi-
nistic methods include the time- and user-honored methods of arbitrary parameter
adjustments by trial and error, variational methods, and a newly formulated Green's
function method (for kinetics systems only). Sampling methods include Monte Carlo
sampling of the outputs of air quality models to compute variances and a Fourier
analysis method of sampling model outputs to compute expectation values of sensiti-
vity coefficients.
Computational economics, inclusive of both programming effort and computer exe-
cution costs, were found to be dominant governing factors for the effective applica-
tion of systematic sensitivity and uncertainty analysis methods to arbitrary air
quality problem scenarios; several reasonable options and several unreasonable option'
emerge from the present evaluations. Recommendations are made outlining how EPA
should, today, most effectively apply available multi-parameter systematic sensitivity
j17 analysis metnoas TO arDitrar^^-^^n^^t^u^^d^^udj ,ty ^.nuiau.u,, ,,,uucia.
'a. DESCRIPTORS jb. IDENTIFIERS/OPEN ENDED TERMS
* Air pollution
* Mathematical models
* Analyzing
* Probability theory
;13 DISTRIBUTION STATEMENT
i
. RELEASE TO PUBLIC
19. SECURITY CLASS (This Report)
UNCLASSIFIED
20. SECURITY CLASS (This page)
UNCLASSIFIED
c. COS AT I Field/Group
13 B
12 A
14 B
21. NO. OF PAGES
105
22. PRICE
EPA Foi-m 2220-1 (9-73)
97
------- |