WATERRESOURCES

RESEARCH,
VOL.30,NO. 11,PAGES
3159-3169,
NOVEMBER
1994

Parameter
uncertainty
andinteraction
in complex
environmental models
Robert C. Spear
Environmental
Engineering
andHealthSciences
Laboratory,
Richmond
Research
Center
Universityof California, Berkeley

Thomas M. Grieb
Tetra
Tech,Inc.,Lafayette,
California,
andEnvironmental
Engineering
andHealth
Sciences
Laboratory,
Richmond
ResearchCenter,Universityof California,Berkeley

Nong Shang
Environmental
Engineering
andHealthSciences
Laboratory,Richmond
ResearchCenter
Universityof California, Berkeley

Abstract.Recentlydeveloped modelsfor the estimationof risksarisingfromthe

release
of toxicchemicals fromhazardous wastesitesareinherently complexboth
structurally
andparametrically.To betterunderstand theimpactof uncertaintyand
interaction
in the high-dimensional
parameterspacesof thesemodels,the set of
procedures termedregionalsensitivityanalysishasbeenextendedandappliedto the
groundwaterpathway of the MMSOILS model. The extensionconsistsof a tree-
structured
densityestimationtechniquewhichallowsthe characterization
of complex
interactionin that portionof the parameterspacewhichgivesrise to successful
simulation. Resultsshowthat the parameterspacecanbe partitionedinto small,
denselypopulatedregionsand relativelylarge, sparselypopulatedregions.From the
high-density regionsone can identifythe importantor controllingparametersas well as
theinteractionbetweenparametersin differentlocalareasof the space.This new tool
canprovideguidancein the analysisandinterpretationof site-specific applicationof
thesecomplex models.

Introduction In this paper we extend the approach to this type of

problem developed by R. C. Spear and G. M. Hornberger
Complexnumericaland analyticalmodelshave beenused [Hornbergerand Spear, 1980;Spear and Hornberger, 1980],
tosimulatethe release,transport,andfate of toxicchemicals who in 1978 carried out an analysis of a multiparameter
fromhazardouswaste sites [Mills, 1993; Huyakorn et al., eutrophicationmodel whose parameters were describedby
!987].Recently, a new set of these modelshas been devel- statistical distribution functions. The Monte Carlo method-
opedto simulate both the release of toxic chemicalsfrom ology that they developed to evaluate the relative impor-
severaltypes of waste management units (e.g., landfills, tance of individual parametersin determiningmodel perfor-
surfaceimpoundments,waste piles, undergroundinjection mance has now come to be called regional sensitivity
wells,andunderground storagetanks)andthe movementof analysis (RSA).
chemicals within differentenvironmentalmedia(groundwa- In subsequentyears the RSA procedurehas been applied
ter, surfacewater, soil, air, and biota). Examples of these to other water quality investigations [Van $traten, 198!;
multimedia exposuremodelsincludeMMSOILS [U.S. En- Whiteheadand Hornberger, !984; Lence and Takyi, 1992;
vironmentalProtectionAgency (EPA), 1992; Staskoand Jakernan eta!., 1990] and to quite different problems,
Fthenakis,1993],MULTIMED [EPA, 1990],and MEPAS includingthe dynamics of moth [Auslander, 1982] and mos-
[Whelanet al., 1992].Multimediamodelshavebeendevel- quito [Eisenberget al., 1994]populations,toxicology[Spear
opedas screeningtoolsto guidethe assessment of potential et al., 199!], nuclear safety [Cook and Gimblett, 1991], and
humanand ecologicalrisks and to formulatealternative control engineering[Auslander et al., 1982; Tsai and Aus-
actionoptionsfor regulatoryagencies.To ensuretheir lander, 1988]. However, Beck [1987, p. 1401], in his ency-
effectiveuse, however, there is a need to establishthe
clopedicreview of uncertainty in water quality modeling,
magnitude andsourcesof uncertaintyassociated with model
pointedout that a disadvantageof the approachwas that
predictionsin order to achieve a better understandingof
"the interpretation of the derived a posteriori parameter
simulated
systems,
to increase
thereliabilityof modelpre- distributions becomes more difficult as the dimension of the
dictions,andto definerealisticvaluesthat shouldbe usedin
subsequent risk assessments. parametervectorincreases,andfor all practicalpurposes,it
seems probable that any conclusionswill have to be re-
Copyright
1994bytheAmerican
Geophysical
Union. stricted to the properties of, at most, the univariate and
Papernumber94WR01732. bivariatemarginaldistributions."This has been the casein
0043-1397/94/94WR-01732505.00 all subsequentstudies carried out by our group using this
3159
3160 SPEAR ET AL.: UNCERTAINTY AND INTERACTION IN COMPLEX MODELS

methodology.In our analysesof complexnumericalmodels ping.If suchdifferencecanbe discernedby an appropriate

we havealsofoundthat there are extensivesetsof parameter statistical
testappliedto the sampledistribution
functions,
valuesthat result in plausiblemodelresultsthat are consis- then evidenceexiststhat xi is an "important" or a "sensi-
tent with experimental observations.Moreover, we have tive" parameter.However, it was realized from the outset
found linear statisticalmodelsto be ineffectivein identifying thatthisindexof sensitivity
isa sufficient,
butnotnecessary,
the underlyingstructureof parameterinteractions. conditionfor sensitivity.One canenvisionvarioustypesof
These resultshave led us to explorea particularclassof parametricinteractionsthat would not be observablefrom
computer intensive algorithms for multivariable data analy- the univariate marginal distributions, a fact that motivated
sis. In this paper, we describea new approachwhich has multivariable analysiseven in the earlieststudyusingthis
arisenout of theseexplorationsandgive preliminaryresults approach.In our experience,however, conventionalmulti-
arisingfrom an applicationto the MMSOILS model. variable analysis of the parameter vectors has not been
particularlyrevealing.

The RSA Concept

The most attractive feature of the RSA idea is its concep- ExperienceWith the RSA Approach
tual simplicity. A model G contains a set of constant An importantdifferencebetweenthe Peel Inlet studyand
parametersx and a set of inputs z which togetherproducean all of the other applicationsof RSA that our grouphas
outputy. Uncertainty and/or variabi!ityin the parametersis carried out concerns the fraction of the total number of
describedby assigningto each element of x a statistical simulations that resultin goodoutcomesor "passes." In the
distribution function which, taken together, constitute a PeelInlet casethisfractionwasabout45%. In all subsequent
multivariatedistributionfunctionf. This assignmentresults studies,it has been a struggleto achieve valuesas highas
in an ensemble of models, each of structure G and with a 5%. The worstcasethat we have encounteredwas 20 passes
parameter vector which is a random member of the multi- in 2.6 million simulations.In all of these subsequentcases
variate distributionf. The goal is often to ascertainwhich therewassomeminorpruningof the rangeusedin thevery
elementsof the parametervector are importantin producing first runs, but for the most part these ad hoc adjustments,
simulationswhich mimicthe essentialfeaturesof the system which were always based on the univariate marginaldistri-
being studiedas reflectedby the outputvector y. To this butions,were very modest, e.g., 10% to 20% reductionsin
end, a criterionfunctionC is defined,which classifiesanyy the range of one or two distributions out of 20 to 25. The
as representing a "good" simulation or a "bad" one. Good interesting
observationhas beenthat even beforethe prun-
simulationsare often referred to as "passes" or "behav- ing, the univariatemarginaldistributionsof passingparam-
iors." The good/badclassificationschememakespossible etersextendedover the entire allowable range for almostall
the app!icationof the largeand powerfularmoryof methods parametersand that after the prune it extended acrossthe
of multivariatestatisticalanalysisin exploringdifferences
in entirerangeof all parameters.We have never seenevidence
the posterior distributionsassociatedwith good and bad that the passingparametersoccupieda single well-defined
simulationsor in exploringthe structureinducedinto the region internal to the hypercube, which might arise, for
joint distributionof parametersassociatedwith goodsimu- example, from a multidimensional normal distribution cen-
lations. tered at some interior point. However, we have made
The power of the method arises from the classification explorationsof the connectednessof the pass regionusing
notion. The situationcan be visualizedgeometricallyby nearest neighbor metrics of several sorts [Li et al., 1994].
consideringeach element of the p-dimensionalvector of These results suggestthat the pass region is generallya
allowableparametervalues to be independentlyand uni- singleconnectedregion.
formly distributedover the interval [0, 1]. When the initial The generalpictureof the passregionthat we haveformed
parameterdistributionsare not sensiblymodeledas uniform is oneof very complexcorrelationstructurebetweenparam-
or there is prior knowledge of a covariancestructure,this eterswhichvarieslocallyover the space.That is, we arenot
informationis best dealt with in a secondary phaseof the dealingwith hyperellipsoids or othermacrogeometries of the
analysis.Hence for presentpurposes,the prior parameter shortwhichunderlieprincipalcomponentsanalysisor other
spaceis theunit hypercube.Any pointwithinthehypercube standardmultivariatestatisticalprocedures.This observa-
canbe identifiedas leadingto a goodor badsimulation by tion suggeststhat there is much information of interest
runningthe modelwith the corresponding parametervector contained in the correlation structure, but it has been a
and applyingthe classification criteriato the outputso challengeto describeand analyze these interactionsin the
produced. Hence the model and the classification criterion 20- and 30-dimensional spacestypical of the appliedprob-
providea meansof dividingthe hypercube intotwo regions, lemsthat we have addressed.In our work on toxicological
one associatedwith good simulationsand the other with bad. modelswe had occasionto explorethe classification and
The informationavailablefrom the approachis containedin regression trees(CART) procedure[Breimanet al., 1984]to
the problem-specificinterpretationof the observablefea- definethe regionsin the spacethat are most usefulin
tures of these two regions. classifying a parametervectoras one leadingto a goodora
As wasmentionedat the outset,mostapplications of the badsimulation. CART is a nonanalytic,computerintensive
RSA concepthavebeenconfinedto viewingthe good/bad procedure which leads to classification rules based on in-
resultsfrom the perspectiveof the univariatemarginal equalityconstraints appliedto individualparametervalues
distributions.For example,if F{xi) is thepriordistribution or to linearcombinationsof parameters.While CART ledto
of theith elementof the parameter vectorxi, thenoneasks usefulinsights,it was not the ideal too1 since it usedboth
whetherF(x•lG) = F(x•IB), thatis, whethertheconditional goodandbadparametervectors,andthe badweregenerally
distributionsshowanydifference underthegood/bad map- in muchgreaternumbersthan the good. However,the
 SPEARET AL.: UNCERTAINTYAND INTERACTIONIN COMPLEXMODELS 3161

potential
oftheCARTapproach
ledustoinvestin exploring m

variationson the tree-structured,computerintensivestrat- (4)

egywhichis its foundation. i=1

For any partitioningof the space, the value of the loss

Tree-Structured
DensityEstimationApproach functioncanbe calculated.As the partitioninggetsfiner,the
The challengeof describingthe passregionbeyondwhat valuegetssmaller.Hence the challengeis to devisea search
canbe observedfrom the'univariatemarginaldistributions, procedurethat findsa seriesof splitsthat resultin a smallL
wecontend,is best approachedas a problemin multivariate while retainingrelatively simple structure. The series of
densityestimation[Silv•,rrnan,1986]. Here we treat the splits is referred to as a tree for reasons that will become
processof gettingpasspointsasa randomprocess.Usingan apparent.The procedurebeginsby splittingS, the parameter
appropriatescalingfactor, the densityof passesat a given hypercube, intotwosmaller hypervolumes S½andSr. The
parameter pointis definedastheproportion of passpointsin rangeof eachparameteris equally divided into k bins, and
a sufficientlysmall neighborhood of the point or, equiva- the grid pointsare usedto splitS. Hence thereare a total of
lently,the chanceof gettingpassesper unit volume of the k x p ways to make a split. If no informationabout the
parameterspaceif samplesare taken randomlyfrom the parametersand their relationshipsexisted, i.e., if the data
neighborhood. Hence the densityvalue is 1 or 0, depending pointswere uniformlydistributedin the originalspace,the
onwhetherthe point is containedwithin the passregionor estimated
density,
f0, wouldbe I/S andthecorresponding
thefail region. For thoseboundarypoints, the densitieswill L, L0, is equal to -l/S, as can be seen from (1) and (4).
bebetween0 and 1. The underlyingconceptof a new method Otherwise,the unevendistributionof the densitysuggests
that more informationabout parameterscan be obtainedby
developed by oneof usinvolvesthe extensionof the concept
of the variable length histogramto multiple dimensions splittingthe regioninto two subregions.Actually,for a given
[Shang,1993].That is, the complexdensityof passesin the split,supposeS isdivided
intotwoareas,Sq and$r, withnq
parameterspace is approximated by uniform densitiesin andrtr pointsin eachandwherertq+ nr = n. Thenthenew
density function is
localregions,just as the histogramwith variablelengthdoes
in a singledimension.The idea is based on the following rtq
philosophy:if the density is uniformly distributedin a local fl(x) = dq: • (Sq)-I xGSq
region(includingthe extreme cases where the density is (5)
equal to 0 or 1), then no further information about the
n I.
parameterand parameter interactions can be extractedfrom
fl(x ) = dr =_n (Sr)- I x {• Sr
thelocal region. The first challenge, however, is to find the
localregions, determine their extent, and arrange them in
somesystematic way. /lq F/r
To gain an overview of the approach, we begin with a L 1= dq- • dr (6)
samplespaceS and a set of observations,x 1, x2, ß' ß, Xn,
whichare independentand identically distributedsamples It is easyto showthat L 1 -< L0, and L0 - L• measuresthe
froma population with unknowndensityf(x). Eachof thexi increasein accuracy when a split is made. A search is then
is a vectorof p parameters.The object is to constructan conductedto findthe bestsplitby maximizingL0 - L1. The
estimate,
f(x), so thatf is a goodapproximation
of f. splittingprocesscontinuesfor the newly generatedregions
SupposeS is partitionedinto m subspaces, S = U/n=•Si, in until the data in each region are judged to be uniformly
whichthe n observations are splitas n = •1 hi. Thenthe distributed or the number of observations in the region is
mostnatural estimate of density in subspaceSi is smaller than some critical value. The splitting process and
the resultingregionswhich can be describedby a binary tree

fi(x)
ni(•i)
= -- (1) structureare presentedin the following example.
Figure I showsa test casein whichf comprisesa weighted
sum of four bivariate normal distributions of different means
The accuracyof estimationis defined by an integrated and variances. A random sample of size 200 was drawn from
squareerror (ISE): this complexbivariatedistribution,and thesepairsof values
were input to the densityestimationprocedure.The task is
to divide the spaceinto regionsof relatively uniform density
ISE
=f.r{f(x)
-j'(x)}
Es 2dx (2)
within the constraintsof the precisionallowed by the sample
size. The scatterplotfor the 200 samplesis shownin Figure
wherethe smallerthe ISE, the betterthe fit. Note that 2. The irregularrectangularbox structure superimposedon it
representsthe boundariesof 13 subregionsdeterminedby
the estimationprocedure.Figure 3 showsthe corresponding
tree. The density values in the nodes are relative densities,
•S •S rl i= 1 which are calculatedby assumingthat the parameter space
(3) has a unit volume. From (!), it is easy to show that the real
densityisjust the relativedensitydividedby the real areaor
It canbeshownthatforanyf withstructure
(1),the!SEcan volume of the parameterspace. The first split is at x2 -<
'• estimated
to withina constant by a lossfunctionL, -9.7, As is notedon the tree, only 5 of the 200 points lie in
definedas that region, with 195beinggreater than that value, as can be
3162 SPEARET AL.: UNCERTAINTY AND INTERACTION IN COMPLEXMODELS

comparing
trees,adjusting
thesplittingcriterionto avoidend
cuts, and developingpruning strategiesfor generatingthe
best subtrees.To give some senseof the concept, we define
the roughnessparameter as

1 m I

Rs
=--
rn i•
1--
Si (7)
The definitionhastwo importantproperties'R -> rn for any
tree with rn terminalnodesand R = rn if and only if Si =
S/rn. Hence equal splittinggives the least roughness.Con-
versely,the more uneventhe splitting,the larger the rough.
nessparameter.In general,R is an increasingfunctionof m.
That means that the roughness of a tree comes from two
aspectsof a partition:finenessand evenness.Many criteria
in the tree-generating
processneedto be critically examined.
For example,the modifiedcriterionfor makinga split is to
Figure 1. Perspectiveplot of the real density. maximize

iX(s)= (L o - L1)/(R• - Ro) (8)

seenin Figure 2. The correspondingdensitiesare 0.07 and
1.48.The proceduredeterminesthat it hasno furtherinterest thereforethe criterion is trying to find the split with a small
in the low density region and concentrateson the higher estimation error and, simultaneously, a small roughness
densityportion of the space.The next split is at x 1 > -7.8. parameter.
Five of the remaining195 points go to the right into a region The roughnessparameter plays a central role in "tree
of density0.16, and the other 190 lie in a regionof density pruning,"which is basicallythe processof arrivingat the
1.89. The processcontinuesuntil the proceduredetermines besttree amongall subtreeswith the sameroughnessparam-
that it does not have the sample size to proceedfurther. With eter. The secondconceptimportant to arriving at a goodtree
200 samples in two dimensions, however, the method has is the notion of cross validation, which is used to arrive at a
discernedthe four separateclustersof pointscenteredon the defensibleestimate of the ISE. Cross validation is a repeti-
subpopulationmeans, which is consistentwith the true
density shown in Figure 1.

Tree Variability and Complexity

The foregoing discussion and example suggesta decep- 0.69
tively simple picture of the art and scienceof successful 0,09

tree-structured density estimation. In practice, there are ,

many difficulties.Among them are the following. o _

.

1. For a given sample, how doesone decidewhen a node •,16'

is a terminal node? Any partitioning of S correspondsto a ß

specialsmoothingof the data. The finer the partitioning,the 0.16 1.06 0.17
rougher the estimation. Both oversmoothing and under- ß

smoothingresult in bad estimations. .- ,

ß , , .....

o -
2. How doesone verify the accuracyof a particulartree?
If the samesampleis used to constructa tree and to estimate '.. ':'
the ISE, the error is called a resubstitution error, and it is ,,

,
.

not, in general, a reliable estimate of ISE. 0.52.

3. How doesone avoid bad splits?The simplifiedpicture -
.... ß •.:•...-•
..
.........

givenabove for selectingsplitsfavors "end cuts," i.e., splits

closeto the boundariesof regions. This can be seenfrom (6),
since smaller volumes result in higher densities and smaller
estimation errors.
0.07
4. Becausea tree is built from a random sample,the tree
itseft is subject to random variability. How does one deter-
mine thosefeatureswhich persistfrom sampleto sample?
5. How can one find hidden interactions?The simplified ,,,

version looks for the largest decreasein the loss function at

any stageof the splittingprocess,and as a result, it can lead
-15 -10 -5 0 5 10 15
one to fall into a local minimum.
The key notion in dealing with many of these problemsis Variable 1

to introduce a roughnessparameter, which measuresthe Figure2. Partitioning

of parameter
spaceanddensity
esti-
complexity of a tree. This measure provides a basis for mates based on 200 observations.
 SPEAR ET AL.' UNCERTAINTY AND INTERACTION IN COMPLEX MODELS 3163

6(•Intermediate
Node 6.49 = Density
5

X2 < -9.7
195

0[• Terminal
Node 0.07 = Density
0.34 o.3.._•4
= Area X1 > -7.8

183
Xl • 5.4

137

X2• 1.1

X2 > -4.3 Xl •;o

51 11 34 41 11 i S0 4

Xl •-1.2 Xl < 3.0 x2 • lO.1 X2 • 12.8

/ \ / \

0.026 0.026 3 0.0096 0,.,,0.

52 0.0.$7 6 0.029
Xl •-1.8 X2 • 6.5

/ \

0.024 0.019 0.022 0.025

Figure 3. Tree diagramfor exampleproblem.

tiveprocedurein which the complete data set is randomly includinglandfills, surfaceimpoundments,waste piles, un-
dividedinto v subsets(v = 10 was used). In each of the v dergroundinjection wells, and undergroundstoragetanks.
repetitions,v - 1 subsets are combined into the learning The physicalcharacteristicsof each type of waste manage-
samplethat is used to constructa tree, and the remaining ment unit are different, but in each case source characteris-
subsetis usedas a test sampleto estimatethe ISE. For each tics are defined and mass balance calculations are made on
valueof the roughnessparameter,the estimatedISEs are an annual basis to ensure that mass is conserved in waste
averaged over the v repetitions.The resultingquantityis a managementunits that have multiple release pathways.
reasonableestimation of the real ISE and is referred to as the Several mathematical models are used to represent fate
crossvalidationerror. The best roughnessparameteris that and transport processesin five pathways: atmospheric,
surfacewater, groundwater,soil erosion, and food chain
onewhich has the smallest cross validation error, or approx-
imatelythe smallestISE. The final tree is constructedby bioaccumulation.The uncertaintyanalysesconductedin this
usingthe best roughnessparameterand the completedata studyfocusedonthe simulationof contaminantmovementin
set. the groundwaterpathway.The simulatedgroundwaterpath-
way includesthefollowingprocesses:net recharge,leaching
of contaminantsfrom the soil, transportthroughthe partially
ModelDescription saturatedzone, and contaminantdispersionwithin an under-
MMSOILS was used to explorethe applicationof the lying aquifer.A seriesof equationsand simplemodelsare
tree-structured
densityestimationmethodology to investi- used to link the movement of water into the soil, fluid flow
gatingmodeluncertaintyandparameterinteraction.It is a and contaminanttransport through the partially saturated
multimediacontaminant fate, transport, and exposure zone, and mixingat the boundaryof the saturatedzone. The
model.It wasdevelopedas partof an assessmentmethod- fate and transport of a contaminantin the aquifer are
ology[EPA, 1992]to estimateexposures andhealthrisks calculatedby using a quasi-analyticalsolution to a three-
with the release and subsequentfate and trans- dimensional
associated transientequation.It incorporatesthe effectsof
longitudinal
port of chemicals from contaminated soils and hazardous advectionandthree-dimensional
dispersion,
re-
chemical tardation, and first-order decay.
wastesites.The modelis capableof simulating
•leasesfrom severaltypesof wastemanagement units, MMSOILS was designedto 'run on a personal computer
3164 SPEAR ET AL.: UNCERTAINTY AND INTERACTION IN COMPLEX MODELS

Table !. Model Parameters

Range of
Parameter Definition Values

Saturated Zone Compartment

GRAD* aquifer
gradient,
m/m 1 x 10-5-0.11
COND* aquiferconductivity,
m/yr 0.31-6300
AQUIDP* aquiferdepth,m 5.0-13.0
Unsaturated Zone Compartment
DEPLYR2* soil layer thickness,cm 83-1090
FIELDC* field capacity 0.01-1.0
ROOTZN rootzonedepth, cm 1 x 10-7-100
KS1 layer1soilsaturated
hydraulic
conductivity,
cm/h 1 X 10-4-25
KS2 layer2 soilsaturated
hydraulic
conductivity,
cm/h 1 x 10-6-35
POERGMI* layer 1% organicmatter 2.5-10.0
POERGM2 layer 2 % organicmatter 1.0-4.0
CN runoff curve number 72-92
RAIN ! * monthly rain, cm/yr 78--117
Landfill
AREAl* layerarea,m2 58-232
CONDWB! hydraulic
conductivity,
cm/h 1 x 10-4-35.1
CONDWB3 hydraulic
conductivity,
cm/h 1 x 10-4-0.78
THICK3* waste layer thickness,m 0.38-1.5
ICLOSE* waste site closure date 25-75
AREAINC incremental
wastearea,m2/yr 1.2-3.4
SOLB* benzenesolubilitylimit, mg/L 315-1260
CSL benzenesoil concentration,mg/kg 3.0-30.0

*Parametersthatexhibitedstatistically
significant
differences
in theirdistributionsbetweenpassand
nonpass groups.

system, and the model output consistsof the predicted The initial analysisinvolvedthe applicationof regional
concentrationin eachenvironmentalpathwayanda selected sensitivityanalysisto investigatewhich modelparameters
exposure point. For a previous study [Tetra Tech, Inc., make the largest contribution to the ability to meet the
1992],the model has also been incorporatedinto a Monte performance criterion.The methodinvolvesthe application
Carlosimulation packageto randomlyselectvaluesof input of univariatestatisticaltechniquesto examine differences
in
parameters from specified distributions. The simulation the marginaldistributionsof each parameterbetweenthe
packageallows up to !80 of the 300 MMSOILS input passand nonpassgroups.In the primary set of analyses the
parameters to be varied in each model run. tree-structureddensity estimation technique was usedto.
examinethe interactionsbetween parametersin the pass
Analyses region of the parameterspace.
The analysesdescribedin thispaperexaminethe relative
importanceof model parametersin predictingchemical Results
concentrations
in the groundwaterat a solidwastemanage- UnivariateAnalysis
ment site regulated under the Resource Conservationand
RecoveryAct of 1976.Modelparameter
valuesrequired
for The model simulationsproduced 1090 passesand 3910
the model simulationshad previouslybeen obtainedfrom nonpasses.The hydraulicconductivityand hydraulicgra•-
numerousdocumentspreparedduringthe courseof remedial ent showedthe largestdifference in the cumulative distribu-
investigationand feasibility studiesconductedat the site tionfunctionbetweenthe two groups.Nine additionalm•el
[Tetra Tech, Inc., 1992].Twentymodelparameters, de- parametersexhibitedstatisticallysignificantdifferences
in
scribedin Table 1, were randomlyvaried in 5000 model their distributions betweenthe passand nonpassgroups.
simulations.These20 parameterswere selectedon the basis Differenceswere exhibitedbetweenthe two groupsin t•½
of theirrelativeimportance
in preliminary sensitiv- typesof distributionbut not in the range of valuesthatmet
regional
ity analyses[TetraTech,Inc., 1992]using54parameters. the performance criteria. Correlation coefficientswere a!s0
Eachmodelrun wasseparated intotwo groups on the calculatedto evaluatethe relationshipbetweenmodelpa-
basisofthemodel performancecriterion: ground- rameters. The results, summarized in Figure 4, indicateaa
predicted
water concentration of benzeneat a specified receptor absenceof importantinterparameterrelationships.Thelarg-
well located505m fromthe center estcorrelationcoefficient(-0.41) occurredbetweenconduc,
location(a drinking-water
of the site)andtime(77 yearsafteroperation of the solid tivity and gradient.The absolutevalue of the majorityof
waste management site began). Model runs in which the correlation coefficients was less than 0.10.
predictedbenzeneconcentration
exceededthe U.S. Envi-
ronmental Protection
Agency's applicablePartitioningof the PassPoints
chemical-specific
or relevantandappropriate requirement for benzene(1 Thediagram producedin the evaluation
of thedensity
:of
/zg/L)wereclassifiedaspasses.
All othermodelrunswere thepasspointsis shownin Figure5. The relativedensity
classifiedas nonpasses. valueis shownfor eachnode,and the proportion of •
 SPEARET AL.: UNCERTAINTYAND INTERACTIONIN COMPLEXMODELS 3165

n= !65
d)/N. This informationon the density in each node is also
summarized in Table 2. In terminal nodes 4, 10, 1, 3, and 14,
the relative density values were greater than 2.0, which
meansthat we expect that at least 2.0 x 1090/5000= 44% of
the pointsin the local regionsare passpoints. The expected
number of pass points in terminal node 4, for example, is
n =23 95%. Overall, 858 of the 1090passes,or 79% were found in
six of the terminal nodes,which representedapproximately
28% of the volumeof the parameterspace.Examinationof
the other terminal nodes shows that 51% of the parameter
space has a relative density of less than 0.5. The expected
n=l n=l proportion of passes in these regions is less than 10%.
m • m
0.1 0.2 0.3 0.4 0.5
Generally, these results show that the parameter spacecan
be partitioned into relatively small, densely populated re-
gionsand large, sparselypopulatedregions.
The model parametersthat are responsiblefor the parti-
Figure4. Correlationbetween model parametersin the tioning of the parameterspaceand the valueson which the
passregion(n = 1090). parameter space is partitioned are also shown in the tree
diagram(Figure 5). The five parametersthat definedthe high
density terminal nodes are COND (aquifer conductivity),
volume of the parameter space occupied by the terminal GRAD (aquifer hydraulic gradient), SOLB (chemical solu-
nodeis also indicated. The relative density value indicates bility), ICLOSE (waste site closure date), and DPLYR2
the density of pass points in each node. In general, if we (unsaturatedzone thickness). These parameters are a subset
haven passesout of N simulations,then a relative density, of those identified in the regional sensitivity analysis that
d, in a local region can be explained as the expectationof accounted for the partitioning of the pass and nonpass
gettingn x d passesif we took N samplesfrom the local regions of the total parameter space (Table 2).
region,or the proportion of passes in the region is (n x The parameterrangeswithin the high density nodes of the

776 3!4

CONDUCT > 630.3

501 275 72 242

GRAD> 0.01 TN1 TN: SOLB.• 12.45

472
SOLB > 9.75 0.0.._•.9 0.,035 88 154
I,CLOSE =• 5O

TN4

102 370
26 3
CONDUCT > 3150.2
CONDUCT • 3465.1
0.033 o 033

74 28 86 284

GRAD 5 0.033 GRAD > 0.061

O.lO! O.lOl

1
0.07_•5 2
I I 0.1,35, 11 3
iii i iiiiiiiim iii iiiiiiiiiiiiiii iiii ii iiiiii i iii i i ii
AREA1 > 101.5 ICLOSE < 55

0•o•,• 1 0.054
Q Intermediate
Node
TN3
1.26 = Density

Terminal Node 3
AQUIDP • 9 DPLYR2 • 586.5
TN13/ • 14 2,48 = Density
0.033 = VolUme
0'033(Highdensity
termtna!
nodes
are shaded)
0,041 o,o•_!
o_..•.• o,o9•

Figure5. Partitioning
ofpass
points
using
density
estimation
approach
(number
of cases
= 1090).
3166 SPEAR ET AL.: UNCERTAINTY AND INTERACTION IN COMPLEX MODELS

Table 2. Tree DiagramStructureand ImportantParameters

Terminal Relative Number of
Node Density Cases Volume Parameters
4 4.35 154 0.033 COND, SOLB, ICLOSE
10 2.94 173 0.054 COND, GRAD, SOLB, ICLOSE
1 2.80 275 0.09 COND, GRAD
3 2.48 88 0.033 COND, SOLB, ICLOSE
14 2.17 96 0.041 COND, GRAD, SOLB, ICLOSE, DPLYR 2
2 ! .89 72 0.035 COND, SOLB
7 0.91 74 0.075 COND, GRAD, SOLB
8 0.58 86 0.135 COND, GRAD, SOLB
13 0.34 15 0.041 COND, GRAD, SOLB, ICLOSE, DPLYR 2
5 0.24 26 0.101 COND, SOLB, GRAD
11 0.22 24 0.098 COND, GRAD, SOLB
12 0.037 4 0.098 COND, GRAD, SOLB, AREAl, AQUIDP
6 0.027 3 0.101 COND, SOLB, GRAD
9 0.00 0 0.066 COND, GRAD, SOLB, AREAl

pass region provide an effective method for explainingthe The tree structure and these parameter ranges within the
tree diagram. Figure 6 depicts the numericalrangesfor each terminal nodes also provide information on the relative
of these five parameters within the six terminal nodesthat importanceof individualparameters. The initial split of the
exhibit high densities(i.e., relative densitiesgreaterthan 1). tree is made on the basis of values of COND, and the
This figure shows that the high-densityregions take on importanceof COND is indicated by the fact that several
valuesover the full range of each of thesefive parameters. subsequentsplitsof the tree refine the acceptablerangesof
The passregion is not a singlehigh-densityarea describedby COND within additional high-density terminal nodes. In
restrictedrangesof values for the five parameters.Instead, contrast,DPLYR2 appearsto be less important in defining
there are several distinct high-density areas defined by the high density nodes. Full ranges of values were observed
differentcombinationsof parameterranges. in all but one of the high density terminal nodes (terminal
Examination of these combinationsprovides a means of node 14).
comparingand contrastingthe effectof model parametersin Another view of the importance of these parametersin
differentportionsof the passregionaswell as informationon differentportionsof the pass region is provided in Figure7.
the interactionbetweenparametersthat was obscuredin the Terminal node 4 has the highest relative density, contains
initial correlationanalysis(Figure4). For example,very low 14% of the observedpasseswithin 3.5% of the parameter
values of COND and unrestricted values for GRAD charac-
space, and is defined by very low values of COND •d
terize terminal nodes 2, 3, and 4. Terminal node 1, on the relatively high values for SOLB and ICLOSE. Terminal
other hand, is characterizedby very low values of GRAD node 1 encompassesa relatively large portion of the pass
and all values for COND except thoseincludedin terminal region (approximately 9%) and is defined on the basisof
nodes2, 3, and 4. The importanceof these combinationsin restrictedvaluesofjust two parameters,GRAD and COND.
definingthesefour high-densitynodeshelpsto explainthe Terminal nodes 10, 13, and 14 have the same values for
relatively strong linear relationshipthat was observedbe- COND, GRAD, and SOLB but different values for ICLOSE
tween these parameters in the initial correlationanalysis and DPLYR2 at separatenodes, which results in different
(r = -0.41). The tree-structured density estimation tech- densities. Terminal nodes 10 and 14 are defined on the basis
nique, however, also explainswhy a strongerrelationship of differences in the values of ICLOSE and DEPLYR2. Both
does not emerge from the linear correlation analysis;while are parametersthat affect contaminantavailability, andit
terminal nodes 1, 2, 3, and 4 exhibit restrictedrangesof seemslikely that theseparametersfunction as surrogates.
eitherCOND or GRAD, the remaininghigh-density terminal Differencesin the densitybetweenterminal nodes13 and14
nodes (10 and 14) are characterized by wider ranges of showthe effectsof a singleparameterin defininga high-
moderate values for both COND and GRAD. densityarea within the passregion.

Waste Site Unsaturated Zone

Conductivity Gradient ChemicalSolubility ClosureDate Thickness
TN 2
(COND) (GRAD) (SOLB) (ICLOSE) (DPLYR2)
TN 3

TN4

TN1

TN 10

TN 14

Parameter Parameter Parameter Parameter Parameter

Ranges Ranges Ranges Ranges Ranges
Figure 6. Numerical
ranges
of fiveparameters
withinhighdensity
terminalnodes(TN).
 SPEAR ET AL.' UNCERTAINTY AND INTERACTION IN COMPLEX MODELS 3167

Terminal
Terminal Node Density
Node Densi_t• I 10 2.94
T" 4.35 '..•
::. !4 2.14
TN1 2.80 • 13 0.34
Parameter Parameter
Conductivity

Conductivity
(COND)
Gradient
(COND)

Gradient
(GRAD) (GRAD)

Solubility Solubility
(SOLB) (SOLe)

Waste Site Waste Site

Closure Date Closure Date
(ICLOSE) (ICLOSE)

Unsaturated Zone Unsaturated Zone

Thickness Thickness
(DPLYR2) (DPLYR2)
•• .• •,•.__.. •_..•
• 'v"'
ParameterRanges Parameter Ranges

Figure 7. Parameterrangeswithinselectedterminalnodes.

Direet• Search terminalnodes.Noneof theseparameters wereinvolvedin the

The terminal nodesin the tree structurediagram(Figure 5) originaldelineation
of terminal node 5 (Figure5), andtwo of
canbe individuallysampledby rerunningthe Monte Carlo the parameters(monthly rainfall,RAIN1, and wastelayer
simulationsusingthe restrictedrangesof the parameters that
thickness,THICK3) played no role in the original
parta•oning
definethe boundaries of the terminal node of interest. These of the parameterspace (Figure5). The remaining threeparam-
directedsearchesof local regions of the parameter space etersappearto haveincreased
importance
in thislocalized
wereconductedwith the goalof extractingadditionalinfor- region.Landfillarea (AREAl) was only importantin the
mationon the importanceof individualparameters. partitioning
of lowdensityareasin the initialanalysis.
Unsat-
uratedzone thickness(DPLYR2) was importantin defining
The terminalnodesalready describedin Figure 5 can be
onlyoneof thesixhighdensityterminalnodesin the original
separated
intotwo categories basedon informationcontent.
Thoseterminal nodesthat exhibit either low or high relative analysis(Figure6). Whilethe wasteclosuredate(ICLOSE),
densitiesd contain little additional information. The low-
which defines the numberof yearsthat wasteswere addedto
the landfill,was importantin definingthree of the original
densityportionsof the parameterspacecontainlittle infor-
mationbecausethey are sparselypopulated.The high- high-density terminalnodes,it is the singlemostimportant
densitynodes(d _>2.0) containlittle additional information parameter in thepartitioning
of thelocalizedregion.
While new parameterswere identifiedin the directed
because almostall pointsarepasspointsandthe important
search,a patternin parameterrelationships similarto that
parameters as well as their rangesare alreadywell defined.
observedin the originalanalysisexists: localizedinterac-
However,thoseterminalnodeswith a relativedensityof
tionsareimportantin describing thesehighdensityregions.
about1.0 representan unclearregionat the boundary
Forexample, high-densityregionsincludebothlowandhigh
betweenthe low- and high-densityregions.They are of values of the landfill closure date (ICLOSE), but distinct
specialinterest becausefurther partitioningof this space
interactionsare observablebetweenother parametersin these
mayidentify other high-densitysubspaces.
differenthigh-density
regions.In TN5-7, low valuesin the
Terminalnode5 from the originaltree structurediagram
(Figure by parameter
5), witha relativedensityof 0.91,wasexamined range
ofICLOSEareassociated withhighvaluesin
theparameter rangeof wastelayerthickness (THICK3)and
usingthe directedsearchstrategy.A new set of 5000
thefull rangeof valuesfor monthlyrainfall(RAIN1).High
simulations
wererunusingrestrictedrangesforthreeparam- valuesof ICLOSE in TNS-10andTNS-11, on the otherhand,
eters:conductivity
(COND), hydraulicgradient(GRAD), areassociated withthefullrangeof valuesin THICK3 butwith
•d chemicalsolubility(SOLB).Theresultsof thisanalysis a relativelylimitedrangein the valuesof RAIN 1.
aresummarized
in Figures8 and9. Figure8 showsthenew
treestructurein which the portion of the parameterspace
definedby terminalnode5 has beenpartitionedinto 11 Discussion
terminal
nodes.The relativedensityin threeof thesenodes
(TNS-7,TNS-10,andTN5-11)wasgreaterthantwicethe Monte Carlo simulationis a basic tool that is now used
relative
density
intheparentterminal node(TN5,Figure5). routinelyto quantifyuncertaintiesassociated with model
Whilethese
nodesrepresented ofthe predictions.
only14%ofthevolume It hasalsobeenusedin examiningtherelative
l•a! parameter
space, 42%,or448of 1067, importance
theycontained of modelparametersin affecting
modelperfor-
ofthenew modelrunsdesignated
as passes. mance.In the mostelementarysensitivityanalyses,single
Figure
9 shows
theranges
offiveadditional that modelp•ametersarevariedusingMonteCarlosimulation
paramete•
to the delineationof the three new high-density whileall otherparameters
contribute
ß are heldconstant.Regionalsen-
3168 SPEAR ET AL.: UNCERTAINTY AND INTERACTION IN COMPLEX MODELS

3OO 767

!CLOSE < 52.5

252 48 66 701

DPLYR2 > TNS.1 AREA1 ,• 92.8

! 93 59
2 46 0.0.._.•9 56 645
THICK3 > 0.938 DPLYR2• 385.1 TN5-S DPLYR2< 284.4

3 56
I 0.111 0.165
0i-;
0.165
353
RAIN1 • 1.08
292

AREA1 _<136.3

=•5-6
0.0495
TN5-7
,,,
0.0605
261

THICK3
92

> 0.881
117

DPLYR2 < 737.55

o
175

TNS-11

1(•)Intermediate
Node
ms.?
1,64 = Density
0.111 0,091
0.0,,9 0.03

Terminal Node 5-7

2.42 = Density
0.03 = Volume
0.03 (High
density
terminal
nodes
are shaded)

Figure8. Results
of directed
search
ofterminal
node5 (number
of cases= 1067).

sitivityanalysis(RSA) techniquesextendthe basicsensitiv- observed

relationships
are discontinuous
andarenotnecessar-
ity analysis
by identifying
therelativeimportance ofindivid- ily monotonic.Theseresultsalsoshowthatsincetheimpor-
uaIparameters throughout the multidimensional parametertanceofindividual parameters varieswithindifferentregions
of
space.However,sinceRSA techniques havebeenlimitedto theparameter space,notonlyaretheresults of theelementary
theexamination of the multidimensional parameter space sensitivity
analyses limitedto particularregions,butthegen-
oneparameterat a time, they haveprovidedlittle informa- eralizationof theseresultsmayalsobemisleading.
tionaboutthe interactionbetweenparameters. Oneof theprincipalbenefitsof thetree-structured density
The tree-structureddensityestimation technique de- estimationtechniqueis the informationthat is obtainedto
scribedinthispaperextends theabilityoftheRSAapproachguidedirected searchesof theparameter space.Theprelim-
toelucidate multivariable
interaction intheparameter space. inaryanalysespresented inthispapershowtheabilitytouse
In the preliminaryanalyticalresultspresented,we have the parameterboundariesof individual terminal nodesm
shown theabilityto identifyrelativelysmall,densely popu- guidetheselection of theinputparameter rangesin orderto
latedregionsand large,sparselypopulated, unstructureddramaticallyincreasethe fraction of model simulationsthat
regions of the parameter space.Fromthe high-density meetthespecified modelperformance criteria.Thisdirected
regions we are ablenot onlyto identifythe important or searchwas used to identify a secondaryset of model
controllingparameters,but also to demonstratethe interac- parameters with well-defined,localizedimportance,butit
tionbetweenparameters in particularlocalizedareasof the couldalsobeusedto improvetheefficiency of MonteCarlo
parameterspace. The high-densityregionsidentifiedare techniques in sensitivity
analyses. In theexample analysis,
characterizedbydistinctlydifferent
combinations ofparam- thenumberof passesout of 5000simulationswasincreased
eter ranges.High values of one parameterare associated from74atterminal node5 intheinitialanalysis
(Figure
5)to
withlowvaluesof another parameterin onehigh-density 1067in the directed search.
region,andin a differenthigh-density
regionthecorrelation The preliminary results from the examination of the
betweenthe sameparameterscan be reversed.In both MMSOILS model are intended to show the use of the
cases, otherparameters
govern theobserved interparameter tree-structured
densityestimationtechniquein uncertain•
relationships.
Theseresultsexplainthe low correlationthat analyses.Only a singlesubmodelof the MMSOILS model
was observedin the originalanalysisof linear correlation (groundwater fate andtransport)was evaluated in •ese
betweenmodelparametersinthepassregion
(Figure
4).The preliminaryanalyses.Although theanalyses described
rep-
parameter
spaceof thepassregioncannotbedescribedby resentonly a cursoryexaminationof the selectedsubmodel,
simple
linearrelationships
between
parameters.
Instead,
the theresultsshowhowthistechnique canbeusedto ident'ffy
 SPEAR ET AL.: UNCERTAINTY AND iNTERACTION IN COMPLEX MODELS 3!69

Terminal Beck,M. B., Waterqualitymodeling:A reviewof the analysisof

Node Density uncertainty,Water Resour. Res., 23(8), 1393-1442, 1987.
Breiman, L., J. Friedman, R. Olshen, and C. Stone, Classification
I 5-11 4.34 and RegressionTrees,WadsworthInternationalGroup, Pacific
:".,'• 5-10 2.26 Grove, Calif., !984.
• 5-7 2.42 Cook, I., and C. G. Gimblett, A risk perspectiveon fusion safety
phenomena,FusionEng. Des., 17, 301-306, 1991.
Eisenberg,J. N., W. K. Reisen, and R. C. Spear, Sensitivity
Parameter analysisof a dynamic model comparingthe bionomicsof two
Waste Site
isolatedCu!extarsalis(Diptera: Culicidae)populations,J. Med.
Closure Date Entomol., in press, 1994.
(ICLOSE) Hornberger,G. M., andR. C. Spear,Eutrophicationin PeelInlet, I,
The problem:Definingbehaviorand a mathematicalmodelfor the
Landfill Area phosphorusscenario,Water Res., 14, 29-42, 1980.
(AREA1) Huyakorn,P.S., M. J. Ungs,L. A. Mulkey, and E. A. Sudicky,A
three-dimensionalanalytical method for predictingleachate mi-
gration, Ground Water, 25(5), 588-598, !987.
Unsaturated Zone
Thickness
Jakeman, A. J., F. Ghassemi, and C. R. Dietrich, Calibration and
(DPLYR2) reliabilityof an aquifersystemmodelusinggeneralizedsensitivity
analysis. IAHS Pub!., 195, 45-51, 1990.
Lence, B. J., and A. K. Takyi, Data requirements for seasonal
Monthly Rain
(;AIN•)
dischargeprograms:An applicationof a regionalizedsensitivity
analysis, Water Resour. Res., 28(7), 1781-1789, 1992.
Li, H., K. Watanabe, D. Auslander, and R. C. Spear, Model
Waste Layer parameter estimation and analysis: Understanding parametric
Thickness
(THICK3) structure,Ann. Biotaed. Eng., 22, 97-!11, 1994.
Mills, W. B., Development of soil cleanup levels at a former
Parameter Ranges
manufacturedgas plant site in Sacramento, California, in Pro-
ceedingsof DevelopingCleanup Standardsfor Contaminated
Soil, Sediment and Groundwater--How Clean Is Clean?, pp.
Figure9. Parameterrangeswithinselectedterminalnodes 277-288, Water Environment Federation, Alexandria, Va., 1993.
defined in directed search.
Shang, N., New developmentsin tree-structuredmethodology,
Ph.D. dissertation, Univ. of Calif., Berkeley, 1993.
thoseparametersthat controlmodel outcomesand to exam- Silverman, B. W., Density Estimation for Statistics and Data
incthe complexinteractionsthat exist betweenparameters. Analysis, Chapman and Hall, London, 1986.
This techniquealso has far greater applicabilityin model Spear,R. C., andG. M. Hornberger,Eutrophicationin Peel Inlet,
evaluationthan just the definition of parameterinteractions. II, Identificationof critical uncertaintiesvia generalizedsensitiv-
ity analysis,Water Res., 14, 43-49, 1980.
Fromthe perspectiveof the groundwaterfate and transport Spear,R. C., F. Bois,T. Woodruff,D. Auslander,J. Parker,andS.
submode!, this techniquecan be usedto comparealternative Selvin,Modelingbenzenepharmaco-kinetics acrossthree setsof
submodels or the relative importance of the different com- animal data: Parametricsensitivityand risk implications,Risk
partments of thesemodels(i.e., saturatedzone,unsaturated Anal., 1I, 641-654, 1991.
Stasko, S., and V. M. Fthenakis, Software review: MMSOILS,
zone,andlandfillparameters).From the perspectiveof the version 2.2, Risk Anal., 13(5), 575-579, !993.
overall risk assessmentframework which is provided by Tetra Tech, Inc., Uncertainty analysisof the Poppy and Petunia
MMSOILS and similar models, this techniqueprovidesa RCRA facilities for RCRA corrective action regulatory impact
toolwith which to evaluate the relative importanceof other analysis,U.S. Environ.Prot. AgencyEnviron.Res. Lab., Ath-
exposuresubmodels(e.g., sourcecharacterization,surface ens, Ga., 1992.
Tsai, K. C., andD. M. Auslander,A statisticalmethodologyfor the
waterpathways,atmosphericpathways,and food chain designof robustprocesscontrollers,J. Dyn. Syst. Meas. Con-
transport)
andto comparealternativemodelconfigurations. trol., I10, 126-133, 1988.
Risk assessment frameworks such as MMSOILS are ulti- U.S. EnvironmentalProtectionAgency(EPA), A subtitleD landfill
matelyintendedto guide the clean-upor remediationof applicationmanualfor the multimediaexposureassessment
model(MULTIMED), preparedby Aqua Term Consultantsand
hazardouswaste managementunits. The use of the tree- ComputerScienceCorporation,U.S. Environ. Prot. Agency
structured
densityestimationtechniquein conjunction
with Office of Res. and Der., Athens, Ga., 1990.
therisk assessmentmodel framework providesnot only a U.S. EnvironmentalProtectionAgency(EPA), MMSOILS: Multi-
toolforidentifying
thosesitesthatrepresentthegreatest
risk media contaminantfate, and exposure model, Environ. Prot.
andthereforeshouldbe cleanedup first, but alsoa method AgencyOfficeof Res. andDev., Athens,Ga., 1992.
Van Straten,G., Analysisof modeland parameteruncertaintyin
for identifyingwhich sites are most likely to respondto simplephytoplanktonmodelsfor Lake Balaton,in Progressin
availableclean-upstrategies. EcologicalEngineering
andManagementby MathematicalMod-
elling, edited by D. M. DuBois, pp. 107-134, Editions
Acknowledgments. Partialfinancial
support
forthisresearchwas CEBEDOC, Liege, Belgium, 1981.
provided
by a ChevronPostdoctoral Fellowshipfor RiskAnalysis Whelan,G., J. W. Buck,D. L. Strenge,J. G. DroppoJr., andB. L.
(NS)andby the U.S. Environmental ProtectionAgency,Officeof Hoopes,Overviewof the MultimediaEnvironmental Pollutant
Researchand Development (contract68-CO-0019). The authors Assessment System(MEPAS), Hazard. WasteHazard. Mater.,
gratefully
acknowledge
threeanonymous
reviewers
for theiruseful 9(2), 191-208, 1992.
commentson the manuscript. Whitehead,P. G., andG. M. Hornberger,Modellingalgalbehaviour
in the River Thames, Water Resour. Bull., 18(8), 945-953, 1984.
References T. M. Grieb, N. Shang,and R. C. Spear, EnvironmentalEngi-
Aus!ander,
D. M., Spatialeffects ofa food-limitedneering
onthestability andHealthSciences
Laboratory,
Richmond
Research
Cen-
mothpopulation,
J. FranklinInst., $14, 347-365,1982. ter, Universityof California,Berkeley,CA 94720.
Auslander,
D. M., R. C. Spear,andG. E. Young,A simulation
basedapproach to thedesign
of controlsystems withuncertain (ReceivedJanuary24, 1994;mvis.ed
June6, 1994;
parameters,
J. Dyn. Syst.Meas.Control,!04, 20-26,1982. acceptedJune27, 1994.)