Second-order Monte Carlo uncertainty/variability analysis
using correlated model parameters: application to
salmonid embryo survival risk assessment
Fu-Chun Wu
a,b,∗, Yin-Phan Tsang
aaDepartment of Bioenvironmental Systems Engineering, National Taiwan University, Taipei 106, Taiwan, ROC
bHydrotech Research Institute, National Taiwan University, Taipei 106, Taiwan, ROC Received 12 September 2003; received in revised form 17 February 2004; accepted 26 February 2004
Abstract
This paper presents a Monte Carlo uncertainty/variability analysis applied to assessment of salmonid embryo survival affected by the deposition of fine sand in the spawning gravel riverbed. For the system being modeled, the uncertainty originates from the lack of perfect knowledge about the spawning habitat, while the variability arises from the natural diversity in the amount of sand deposit. A second-order Monte Carlo simulation (MCS) procedure consisting of multiple realizations of model parameters and full-range iterations of input variable is used to separately propagate the knowledge uncertainty and natural variability of the system. Given that the three parameters of the present model are highly correlated (correlation coefficients >0.8), we use four multivariate methods (Iman–Conover method, standard normal transformation method, normal copula method, and maximum-entropy copula method) to generate correlated model parameters. The scatter plots of the correlated parameters gen-erated with various methods demonstrate different correlation patterns. An important finding of the present uncertainty/variability analyses and risk assessments is that the outcomes resulting from the four correlated scenarios are very similar, regardless of the diverse correlation patterns. The output of the second-order MCS is a collection of cause–effect relations (spaghetti plot) that can facilitate the analyses of survival-rate uncertainty/variability. The favorable habitat conditions are also defined on the basis of the spaghetti plot. It is revealed that significant values of information are associated with the parameter correlations, especially when higher levels of embryo survival are targeted. The results further indicate that the outcomes of the independent (less realistic) scenario substantially differ from those of the correlated ones, indicating that the effect of parameter correlations should be incorporated into the MCS whenever possible. This work provides a useful paradigm for the systematic uncertainty/variability analyses in the context of ecological risk assessments.
© 2004 Elsevier B.V. All rights reserved.
Keywords: Monte Carlo simulation; Uncertainty/variability analysis; Ecological risk assessment; Model parameter; Correlation; Salmonid embryo survival; Spawning habitat
∗Corresponding author. Tel.:+886-2-2369-2466; fax:+886-2-2369-2466.
E-mail address: [email protected] (F.-C. Wu).
1. Introduction
The total uncertainty (unpredictability) of a system is the combination of two components, namely, uncer-tainty and variability (Hattis and Burmaster, 1994; Rai 0304-3800/$ – see front matter © 2004 Elsevier B.V. All rights reserved.
et al., 1996; Vose, 2000). Uncertainty is the analyst’s lack of knowledge (or level of ignorance) about the system, which may be reduced by further measure-ment or study. Variability represents the random (or stochastic) heterogeneity in a population (people or events) that is irreducible by additional measurements or studies but may be reduced by changing the sys-tem (Anderson and Hattis, 1999; Vose, 2000). Both uncertainty and variability can be described by prob-ability density functions (PDFs), one might therefore mix them in the Monte Carlo risk assessment (in which some distributions reflect the uncertainty about model parameters, while the others reflect the inher-ent stochastic nature of the system) and claim that the output would take account of all uncertainty and vari-ability. However, the resultant single distribution is technically difficult to interpret and information is lost regarding the individual contributions of uncertainty and variability to the total uncertainty. Moreover, a much larger problem than mixing uncertainty and variability can occur when a variability distribution is used as if it were an uncertainty distribution, in this case a variable could be mistakenly treated as an uncertain parameter such that the output that should be a single figure is replaced by a meaningless distri-bution (Vose, 2000). Therefore, it has been increas-ingly emphasized the need to separate uncertainty and variability in risk assessments (Hoffman and Hammonds, 1994; Frey and Rhodes, 1996; Rai and Krewski, 1998; Nauta, 2000; Mitchell and Csillag, 2001). The separation of uncertainty and variability allows us to identify the steps that can be taken to reduce the total uncertainty, to gauge the value of in-formation or of potential system changes, and to avoid the serious errors induced by mixing variables with parameters.
In recent years, the second-order Monte Carlo simulation (MCS), also known as two-phase or two-dimensional MCS, has been extensively used to separately propagate uncertainty and variability in risk assessments (e.g., Bogen and Spear, 1987; Hoffman and Hammonds, 1994; Hession et al., 1996; Burmaster and Wilson, 1996; Nauta, 2000; Vose, 2000; Sanga et al., 2001; Moschandreas and Karuchit, 2002; Pouillot et al., 2003). The second-order MCS involves a double looping (or nesting) procedure consisting of multiple realizations of model param-eters and iterations of input variables. The outcome
of the second-order MCS is a collection of cumu-lative distribution functions (CDFs) or cause–effect relations that simultaneously display the uncer-tainty and variability in the results. In this study, we apply the second-order MCS in the context of salmonid embryo survival risk assessment. Based on the results of such analyses, the implications to the management of spawning habitat are also ex-plored.
Parameter distributions in a model often have to be correlated in order to ensure that only meaningful sce-narios are simulated (Vose, 2000; Todd and Ng, 2001). It has been recognized that when the parameters of a model are correlated, incorporation of such corre-lations in the MCS is important, especially when the correlation coefficients exceed 0.7 (Smith et al., 1992; Bukowski et al., 1995; Haas, 1999; Pohlmann et al., 2002). The correlation (or dependence) between the random numbers is generally indicated by such mea-sures as the product-moment correlation coefficient (Pearson’s r) or rank–order correlation coefficient (Spearman’sρ or Kendall’s τ) (see, e.g.,Nelsen, 1999 andVose, 2000). Parameter correlations are incorpo-rated into the MCS by generating many sets of corre-lated model parameters and using different set in each realization. In this study, we employ four multivari-ate methods to genermultivari-ate correlmultivari-ated model parameters. The parameter vectors produced with these meth-ods are compared. Further, the effect of parameter correlations and the associated value of information are evaluated. These results provide valuable insights into the problem and useful guidelines for habitat improvement.
Monte Carlo methods have been widely used in the context of ecological modeling (Neave et al., 1997; Harmon and Challenor, 1997; Annan, 2001; Qian et al., 2003; Waller et al., 2003). However, Mitchell and Csillag (2001)pointed out that ecosystem mod-els that include both variability of driving variables and uncertainty in their predictions are rare. The second-order MCS is employed herein to address this issue. Moreover, various model approaches to survival risk assessment have been presented, such as the stochastic model (Todd and Ng, 2001) and vitality-based model (Salinger et al., 2003). The present survival model is based on deterministic cause–effect relations validated by field data, thus eliminates the process-based stochasticity.
2. Salmonid embryo survival model
Field studies have indicated that salmonids (salmon and trout) explore suitable locations in the gravel riverbed for spawning, where the female excavates a pit and releases fertilized eggs into the bottom. After spawning, the female resumes digging up-stream to bury the area of egg deposition (so called “egg pocket”). The embryos are therein protected against bedload motion and scour during high flows. The typical morphology of the resultant nest (redd), with a pit upstream and a tailspill downstream, is shown in Fig. 1. The hydrodynamics induced by the bed topography promotes a downwelling flow into redd. The pressure gradient exerted between the up-stream and downup-stream faces of the tailspill forces streamwater to flow through the gravels. This flow benefits the exchange of dissolved oxygen (DO) and removal of metabolic wastes (MW) such that a suit-able intragravel environment for embryo development is maintained. Among the factors adversely affect-ing embryo survival, accumulation of fine sand in the gravel substrate is one of the most detrimental. Fishery researchers generally agree that intrusion of sand into spawning gravels substantially reduces substrate permeability and intragravel water veloci-ties, thereby restricting supply of DO and removal of MW. Modeling of embryo survival is based on the above-stated mechanisms (for details see Wu, 2000), which were validated by field data (Tappel and Bjornn, 1983). A brief overview of this model is as follows.
Fig. 1. Typical morphology of salmonid spawning nest (redd), with a pit upstream and a tailspill downstream. The bed topography promotes a downwelling flow into redd.L1 is the length of flow path through surface sand seal (layer 1),L2is the length of flow path through surrounding gravels (layer 2) (modified from Wu, 2000).
2.1. Model components
2.1.1. Relation between sand deposit and substrate permeability
A nonlinear relation between the reduced substrate permeability and sand content is given by the follow-ing expression derived from the hydraulic radius the-ory (Wu, 2000): K K0 = 4.54 × (0.42 − 1.54σ)3 (0.58 + 1.54σ)2 + 3.66 ds Dg 2 σ (1)
whereK0is the permeability of the clean gravels;K is the reduced permeability of the polluted gravels;σ is the specific sand deposit, defined as the volume ra-tio of sand deposit to bulk gravels;ds/Dgis the ratio of sand to gravel grain sizes. For a clean gravel bed (σ = 0), the second term on the RHS ofEq. (1)is in-effective, whereas the first term vanishes as the intra-gravel pores are saturated with sand (corresponding to
σ = 0.273).
2.1.2. Relation between substrate permeability and seepage velocity
Streamwater enters the redd gravels in the high-pressure region (upstream face of tailspill) and leaves the gravel bed in the low-pressure region (downstream face of tailspill). The typical flow path is illustrated in Fig. 1. The substrate gravels tend to filter out the sand particles carried in the water, thus forming a sand seal in the surface layer of the gravel bed that inhibits fur-ther sand infiltration. The permeability of the sand seal is relatively low compared to the surrounding gravels. On the basis of Darcy’s law, the apparent seepage ve-locity through the two-layer redd gravels can be de-termined by:
V = (h/L1)K2
L2/L1+ K2/K1 (2)
where L1 is the length of flow path through layer 1 (sand seal); L2 is the length of flow path through layer 2 (surrounding gravels);K1andK2are the per-meability of layers 1 and 2, respectively; h is the total pressure-head drop between the upstream and downstream faces of tailspill. The ratio K1/K2 used in Eq. (2) is simply the K/K0 value obtained from Eq. (1), for whichσ is taken to be the specific deposit in layer 1.
2.1.3. Relation between seepage velocity and embryo survival
The apparent seepage velocity can serve as an in-dicator for embryo survival assessment. The experi-mental data ofCooper (1965)were used to develop an empirical relation between seepage velocity and em-bryo survival rate (R2= 0.99, n = 8):
S = −17.6(log V)2− 39.6(log V) + 68.7 (3)
whereS is percentage survival; V is in cm/s. 2.2. Model parameters
2.2.1. PDFs of model parameters
The embryo survival model described above in-volves three physical parameters, namely, grain-size ratiods/Dg, pressure-head ratioh/L1, and flow-path ratio L2/L1. To determine these values, Chinook salmon (Oncorhynchus tshawytscha) is selected as the target species because more data are available for this species. The ranges of these parameters and their most probable values, listed inTable 1, were obtained through a comprehensive survey of literature (Wu, 2000). However, the currently available data are in-sufficient to determine the PDFs of these parameters. Normal and log normal distributions are among the most widely adopted PDFs for modeling uncertain pa-rameters, but they are unbounded at two and one side, respectively, thus are inappropriate for the bounded parameters. Nevertheless, since we have the possible ranges and the most probable values of parameters, it seems reasonable to use the bounded distributions, such as triangular distributions, as a first
approxima-Fig. 2. Probability density functions (PDFs) of three model parameters. The minimum and maximum observed values are used to infer the lower and upper bounds, and the most probable values are taken to be the modes of triangular PDFs.
Table 1
Ranges and most probable values of habitat characteristics and model parameters (modified fromWu, 2000)
Item Range Most probable value
Habitat characteristic
Observed data range Most frequently observed
Dg/ds 7–30 14
h 1–10 cm 4 cm
L1 5–10 cm 8 cm
L2 100–350 cm 280 cm
Model parameter Inferred parameter range
Inferred ds/Dg 0.03–0.14 0.07 (=1/14)
h/L1 0.1–2 0.5 (=4 cm/8 cm)
L2/L1 10–70 35 (=280 cm/8 cm)
Target species: Chinook salmon (Oncorhynchus tshawytscha).
tion to the PDFs of model parameters. The minimum and maximum observed values are used to infer the lower and upper bounds of the triangular PDFs, and the most probable values are taken to be the modes of the PDFs (Fig. 2). Although the observed ranges of parameter values are used herein, there is no support for setting such limits on the ranges of uncertainty, thus one could allow slightly wider ranges if desired. 2.2.2. Correlations among model parameters
The experimental data of sand infiltration into a gravel bed (Wu and Shen, 1999) are used to evalu-ate the correlations among the parameters. A total of 13 data sets are analyzed, resulting in a correlation matrix as given in Table 2. Because the calculated product-moment and rank–order correlation coeffi-cients are very close, the values listed inTable 2can
Table 2
Correlation matrix of model parameters
ds/Dg h/L1 L2/L1
ds/Dg 1 0.8 0.9
h/L1 1 0.8
L2/L1 1
be used in either context. It turns out that the cor-relations among the parameters are all positive and greater than 0.8, with the most significant correlation existing between ds/Dg and L2/L1. These values also suggest that we need to incorporate parameter correlations in the present study.
3. Generation of correlated random numbers
There are numerous ways of obtaining correlated random numbers (Burmaster and Anderson, 1994). Herein, we employ four methods appropriate for multivariate correlations to generate a large number of parameter sets using the marginal PDFs and cor-relations specified in the preceding sections. These methods include a distribution-free approach that in-volves the orthogonal transformation of a van der Waerden random score matrix (i.e. Iman–Conover or IC method); an inverse scheme based on the or-thogonal transformation in the standard normal (SN) space (i.e. standard normal transformation or SNT method); and two multivariate copula methods re-spectively employing the continuous normal and piecewise-uniform copula densities (i.e. normal cop-ula and maximum-entropy copcop-ula or NC and MEC methods). The theoretical background and procedure of these methods are included inAppendix A.
For illustration,Fig. 3 presents the scatter plots of
ds/Dg versus L2/L1 generated with various meth-ods (with rank or moment correlation= 0.9), where different sample sizes are the result of convergence, goodness-of-fit, and correlation tests. The convergence test is to check if the means and standard deviations of the generated random numbers converge to the pre-scribed marginal values. If the test criterion is not met we generate 100 more observations, primarily due to the concern that the generated random numbers do not sufficiently span the parameter space. As such, among the four methods, only IC method can generate a fixed
sample size of parameter vectors meeting the speci-fied criteria. The independent scatter plot and the ex-perimental data (Wu and Shen, 1999) are also shown for comparison (Fig. 3), where the independent ran-dom numbers were generated using Latin hypercube sampling (LHS) technique. In those correlated scatter plots, an obvious positive-correlation scatter pattern is demonstrated, although their shapes are different. The IC, SNT, and NC scatter plots are elliptical in shape, whereas the MEC plot is piecewise “pinched” in four blocks. The elliptical-shaped scatter plot of IC method is the result of using van der Waerden scores in the generation of correlated random numbers, which was the original intention ofIman and Conover (1982)to use such scores for producing “natural-looking” cor-relations. The NC scatter plot appears to be in good agreement with the IC plot (Haas, 1999). However, due to the acceptance/rejection procedure employed, the NC method needed more iterations and computa-tion time (in the order of several hours) to meet the specified criteria for the three tests. The sample size of SNT method was the greatest, making its scatter plot spindlier in shape. Nonetheless, owing to the straight-forward inverse transformation scheme employed, the computation time taken by the SNT method was in the order of one minute.
The piecewise pinched scatter plot of MEC method is unique among others, which is a combined result of partitioning each uniform marginal on the unit hyper-cube into four equal-size pieces for constructing the piecewise-uniform copula densities and the high de-gree of correlation adopted. For example, the random numbers in the second interval of ds/Dg are associ-ated with the majority of random numbers in the sec-ond interval ofL2/L1, with only rare in the first and third intervals but none in the fourth. Note that the unequal-size intervals on both marginals of the MEC scatter plot are attributed to the transformation be-tween uniform and triangular distributions. Although the result of MEC method does not look very realistic, it does provide proper correlation, convergence, and goodness of fit. While these matched moments are not sufficient enough to encourage use of this approach, as will be shown later on, all these four methods yield convergent results in the uncertainty/variability anal-yses and risk assessments.
The results shown in Fig. 3 indicate that the cor-relation pattern between the random variates can be
Fig. 3. Scatter plots showing the correlation between ds/DgandL2/L1. The random parameters are generated with four correlation methods (with rank or moment correlation= 0.9) and independent Latin hypercube sampling. Experimental data (Wu and Shen, 1999) are also shown.
rather different even specified with the same marginals and correlations. The correlation pattern observed in real data may not match with that produced (Fig. 3). Previous researchers have pointed out that to specify
only the marginal distributions and correlation coef-ficient might not be sufcoef-ficient when it is desired to generate random numbers of a certain correlation pat-tern (e.g.,Haas, 1999; Vose, 2000). The correlation
patterns among the present model parameters are be-yond our current knowledge. Nevertheless, the effects of using different correlation methods on the outputs of various convolutions and, of more direct interest, on the outcomes of Monte Carlo uncertainty/variability analyses and risk assessments are investigated in this study.
4. Second-order Monte Carlo simulation
The salmonid embryo survival model presented in Section 2 and the correlated model parameters ob-tained inSection 3are integrated into the framework of second-order MCS to propagate the knowledge un-certainty and system variability. The components of uncertainty and variability involved and the procedure of second-order MCS are described below.
For the salmonid embryo survival model, the pa-rameter uncertainty originates from the lack of perfect knowledge about the physical conditions of the spawn-ing habitat, which include the gravel/sand grain sizes, flow hydraulics, and redd dimensions. Such knowl-edge uncertainty may be reduced through further
Fig. 4. Schematic diagram showing the procedure of second-order MCS and the outcomes. The inner loop consists of m iterations propagating the system variability; the outer loop consists of n realizations propagating the parameter uncertainty.
investigations. As mentioned earlier, with the current knowledge these parameter uncertainties are approxi-mated by the triangular PDFs. Accordingly, sufficient sets of correlated model parameters have been gen-erated; each set represents a meaningful and possible combination of habitat conditions. Propagation of such uncertainty is implemented by using a different set of parameters in each realization (Fig. 4). On the other hand, the amount of sand accumulated in the bed surface represents a natural diversity of the system, which varies as a complex function of many temporal and spatial factors (Wu and Shen, 1999; Wu, 2000). Such diversity is irreducible through further studies, and is the cause of survival-rate variability under a given combination of habitat conditions. Propagation of such variability is implemented by running multiple iterations with the inputσ-value incrementally varying from 0 (clean gravels) to 0.273 (saturated sand seal) to cover the full range of specific sand deposit (Fig. 4).
The double-looping procedure of the second-order MCS is demonstrated inFig. 4, where the outer loop consists of n realizations of model parameters to prop-agate the knowledge uncertainty; the inner loop con-sists of m iterations of input variable to propagate
the system variability. The ultimate outcome of the second-order MCS is a collection ofS versus σ rela-tions simultaneously displaying the survival-rate un-certainty and variability. The results concerning the uncertainty and variability of embryo survival, the risk arising from knowledge uncertainty and system variability, and the value of information associated with parameter correlations are discussed in subse-quent sections.
5. Uncertainty/variability analyses and risk assessments
The outcomes of the second-order MCS using the parameters generated with different correlation meth-ods are shown inFig. 5, where the “spaghetti-looking” graph demonstrates a number of possibleS versus σ relations, each one of which corresponds to a differ-ent habitat scenario. On the basis of these results, the uncertainty/variability analyses and quantitative risk assessments can be conducted, and the effect of pa-rameter correlations on the model outcomes can be evaluated. Currently a standard procedure for uncer-tainty/variability analyses is not yet available, although the second-order MCS is increasingly employed as a practical means to separately propagate uncertainty and variability. In this paper we systematically address these issues in an attempt to provide a paradigm for such analyses.
5.1. Comparison of the four methods for generating correlated random numbers
In this section, the uncertainty and variability of salmonid embryo survival are investigated, and the outcomes resulting from different correlated scenarios are compared. The four methods for generating cor-related random numbers are evaluated on the basis of these results.
5.1.1. Uncertainty of embryo survival rate
For each spaghetti plot inFig. 5, the uncertainty of embryo survival rate at a givenσ, denoted as Su(σ), is evaluated bySmax(σ) − Smin(σ), where Smax(σ) and
Smin(σ) are the maximum and minimum values of S
corresponding to the givenσ, respectively (as shown inFig. 4). The survival-rate uncertainty Su(σ)
gen-erally increases with σ. However, diverse results of
Su(σ) are associated with different correlation
meth-ods. For lower values ofσ (roughly for σ < 0.15), the SNT uncertainty is the smallest; the NC and MEC un-certainties are similar in magnitude and greater than other correlated ones. The survival-rate uncertainty associated with the independent parameters is con-sistently the largest among all results. These appear to indicate that, for those lowerσ-values, incorporat-ing parameter correlations into MCS reduces the un-certainty of the outcome due to limiting attention to only realistic combinations of parameter values. How-ever, at higher values ofσ (roughly for σ > 0.15), no consistent trend can be found. The spaghetti plot of the independent parameters appears less orderly than the correlated ones, implying that some unrealistic or meaningless scenarios could be generated. Although the spaghetti plots of the correlated parameters appear more orderly, some extreme conditions might have arisen from the random sampling procedure. For ex-ample, in the NC spaghetti plot, the lowestS versus
σ curve isolated from others is a good illustration of
such situation. Herein, the concept of uncertainty in-terval (Sanga et al., 2001) is adopted to address this issue.
To eliminate the extreme effect, the 5th and 95th percentile values ofS, along with the median S, are plotted againstσ inFig. 6afor the outcomes of differ-ent correlation methods. The difference between the 95th and 5th percentile values of S is the length of the 90% uncertainty interval (not to be confused with confidence interval), which is plotted against σ in Fig. 6bfor different methods. It is revealed inFig. 6a that the 5th, 95th percentile and medianS versus σ curves of all correlation methods are in good agree-ment. The independent 95th percentile and medianS versus σ curves coincide well with those correlated ones. However, the independent 5th percentileS ver-susσ curve is quite different from the correlated ones. The independent 5th percentile values of S are con-sistently smaller than the correlated ones forσ < 0.2, but then exceed the correlated ones forσ > 0.2. If the correlated scenarios are considered as more realistic, the implication of the results shown inFig. 6bwould be that by ignoring the parameter correlations the survival-rate uncertainty is overestimated for those less harmful conditions of sand deposit but is under-estimated for those worst sand deposit conditions.
Fig. 5. The outcomes of the second-order MCS using the parameter sets generated with different correlation methods. The spaghetti-looking graph presentsS vs. σ relations corresponding to a number of possible habitat scenarios.
Fig. 6. (a) The 5th and 95th percentiles and median ofS vs. σ relations derived from the spaghetti plots for different correlation methods. (b) The width of the 90% uncertainty intervals (= difference between the 95th and 5th percentile values) for different correlation methods.
5.1.2. Variability of embryo survival rate
For anyS versus σ curve in the spaghetti plot (i.e. for any set of model parameters), the variability of embryo survival rate induced by the system variabil-ity, denoted as Sv, is defined by S(0) − S(σsat), whereS(0) and S(σsat) are the S-values corresponding toσ = 0 and σsat(= saturated value), respectively (as illustrated inFig. 4). Because eachS versus σ curve is associated with aSv-value, theSv-values of all curves in a spaghetti plot are calculated and the his-togram is shown inFig. 7for the outcomes of different
correlation methods. The survival-rate variabilitySv ranges between 0 and 90%, mostly between 0 and 20% but rarely beyond 70%. To further examine the distribution of the survival-rate variability, the statisti-cal moments (mean, standard deviation, skewness and kurtosis coefficients) of theSv-values are calculated. FromTable 3we see that the mean values, skewness and kurtosis coefficients are similar for the indepen-dent and correlated scenarios, with the indepenindepen-dent meanSv-value being 13.7% and the correlated ones ranging from 10.2 to 16.1%. The skewness coefficients
Table 3
Statistical moments of survival-rate variability for different meth-ods
Method Mean (%) S.D. (%) Skewness Kurtosis
Independent 13.7 9.4 1.5 6.6
IC 15.6 13.8 1.6 5.7
SNT 16.1 14.2 1.7 7.5
NC 10.2 12.9 1.8 7.6
MEC 15.0 12.7 1.2 3.9
are all positive (ranging from 1.2 to 1.8) and the kur-tosis coefficients vary between 3.9 and 7.6. However, a significant discrepancy exists between the standard deviations for the independent and correlated scenar-ios. The correlated standard deviations are within a narrow range between 12.7 and 14.2%, whereas the independent standard deviation (=9.4%) is much less than the correlated ones. Given that the magnitude of standard deviation represents the uncertainty in a dis-tribution, the implication of such discrepancy would be that by ignoring the parameter correlations the un-certainty of survival-rate variability is underestimated. Based on the results presented above, it can be con-cluded that the outcomes corresponding to the param-eters generated with different correlation methods are very similar. However, the IC method can be selected as the simplest and most satisfactory one in view of its efficiency and overall performance. In the subse-quent sections, although four correlation methods are used in the analyses, only the outcomes of IC method are used to demonstrate the results associated with the correlated scenario.
5.2. Definition of favorable spawning habitat 5.2.1. Favorable habitat conditions
It would be natural for the salmon conservationists to ask: what constitute the physical conditions of a favorable habitat? To answer this, herein we define a favorable habitat as the habitat with S(σsat) > 70%. To examine the physical conditions that facilitate the favorable habitats, the model parameters correspond-ing to the favorableS versus σ curves are identified. Two findings are gained from these results. First, the distributions of the favorable conditions are ob-tained, as shown inFig. 8. For the three parameters, the favorable ranges for the independent scenario are consistently larger than those for the correlated one.
Fig. 8. Distributions of favorable habitat conditions for indepen-dent and correlated scenarios. (The original ranges of the three parameters are shown for comparison.)
The independent favorable parameters cover 87, 84, and 94% of the original ranges for ds/Dg, h/L1, and L2/L1, respectively, while the corresponding ranges of the correlated favorable parameters only cover 50, 44, and 62% of the original ranges, pri-marily within the middle-lower portions. These much smaller ranges of the correlated favorable parameters imply that the independent scenario might be overop-timistic, especially for the parameterh/L1. The fa-vorable range of h/L1 for the correlated scenario is
Fig. 9. Dependences among favorable habitat conditions for independent and correlated scenarios.
only half of the corresponding range for the indepen-dent one. Much of the favorable range overestimated with the independent scenario is in the uppermost portion.
Second, the dependences among the favorable pa-rameters are observed (Fig. 9). For the correlated sce-nario, apparent negative and positive dependences are demonstrated in the scatter plots ofh/L1versusds/Dg andL2/L1versusds/Dg, respectively, whereas no ob-vious dependences are shown in the scatter plots of
L2/L1versush/L1. For the independent scenario, an apparent negative dependence is revealed only in the scatter plot ofh/L1versusds/Dg. To make quantita-tive comparisons, the correlation coefficients between the favorable parameters are shown inFig. 9. The cor-relation coefficient betweenh/L1andds/Dgis−0.48 for the correlated scenario, while the corresponding value for the independent scenario is−0.82. The phys-ical interpretation of this negative correlation is that the finer sand deposit needs a greater pressure head to facilitate favorable seepage velocity and thus the sur-vival rate. From a management perspective, the much higher negative correlation associated with the
inde-pendent scenario would potentially lead to an over-estimated cost for habitat improvement. It is also in-teresting to note that the dependence between h/L1 andds/Dgis demonstrated by a band-shaped scatter plot rather than a natural-looking pattern as shown in the scatter plot ofL2/L1versusds/Dg, implying that the pressure head is more dominant than the length of flow path in affecting the seepage velocity and thus the survival rate for a givends/Dg. On the other hand, the correlation betweenL2/L1andds/Dgis 0.51 for the correlated scenario, whereas no obvious correla-tion is observed for the independent one. The physi-cal interpretation of this positive correlation is that the finer sand deposit needs a shorter flow path to facil-itate the favorable seepage velocity and thus the sur-vival rate. The much lower correlation associated with the independent scenario would probably lead to an ineffective design of habitat improvement measures. In general, the correlations betweenL2/L1andh/L1 are low. Such negligible correlations imply that the pressure head and the length of flow path are unlikely to constrain each other in the formation of favorable habitats.
Fig. 10. (a) Maximum tolerable sand deposit vs. target survival rate for independent and correlated scenarios. (b) Value of information vs. target survival rate.
5.2.2. Maximum tolerable sand deposit for target survival rate
To implement a project aimed to improve the spawning habitats for embryo development, it would be necessary to set a target survival rate. Once the target is set, the mitigating measures are taken to control the deposition of fine sand into the spawning gravels. What must be known at the planning stage is the amount of sand deposit tolerable for achieving the target survival rate. The spaghetti plots given inFig. 5 may be used to determine the maximum tolerable sand deposit corresponding to the specified target. For a specified target survival rate (represented by a horizon-tal line in the spaghetti plot), the maximum tolerable
specific deposit, denoted asσt,max, is determined by the first intersection (from the left side) withS ver-susσ curves. As each S versus σ curve corresponds to a possible habitat scenario, this intersection repre-sents a most conservative estimate of σ (and thus a maximum tolerableσ-value) for a least favorable com-bination of habitat conditions. Theσt,max-values corre-sponding to various target survival rates are shown in Fig. 10a. Theσt,max-values consistently decrease with target survival rate for both scenarios. The indepen-dentσt,max-values become smaller than the correlated ones as the specified target survival rates are higher than 40%. Since the habitat improvement projects are most likely aimed at those higher target values, the
underestimated σt,max-values associated with the in-dependent scenario would probably lead to overesti-mated costs.
To further demonstrate the value of information as-sociated with the parameter correlations, the abso-lute difference between the correlated and independent
σt,max-values, denoted as|σt,max|, is plotted against the target survival rate. It is shown in Fig. 10 bthat the value of information increases with target survival rate (primarily at the higher target values), indicating that incorporating the parameter correlations in risk
Fig. 11. (a) Variations of risk for independent and correlated scenarios. The risk for not achieving the target survival rate varies with the specific deposit. (b) Value of information vs. specific deposit for different target survival rates.
assessments is crucial for the planning of habitat im-provement, especially when the target survival rate is set to a higher standard.
5.2.3. Risk for not achieving target survival rate Given a target survival rate, the corresponding risk for not achieving the specified target is of important concern to risk assessors and decision makers. Such a risk, arising from the knowledge uncertainty in spawn-ing habitats and natural variability in sand deposits, may be evaluated with the spaghetti plots given in
Fig. 5. For a target survival rate, the risk at some
σ-value is defined as the percentage of the samples for
which the correspondingS-values are below the spec-ified target. The resulting risk diagram, expressed as a function of target survival rate and specific deposit, is shown in Fig. 11a. The risks associated with dif-ferent correlation scenarios consistently increase with specific deposit. Significant rises of risk occur atσ > 0.16, 0.18, and 0.2 for the target S = 80, 70, and 60%, respectively. The risks corresponding to the higher tar-get survival rate are consistently greater than those corresponding to the lower one. No consistent trend between the independent and correlated risks is ob-served.
To assess the value of information associated with the parameter correlations, the absolute differ-ence between the correlated and independent risks, denoted as |Risk|, is calculated and its cumula-tive value is plotted against σ in Fig. 11b. For the lower two targets (60 and 70%), significant values of information are associated with those worst spe-cific sand deposits (σ > 0.2). However, the value of information corresponding to the highest target (80%) is much greater than those corresponding to the lower ones, consistent with the result obtained previously.
6. Summary and conclusions
This paper addresses two important issues involved in the survival risk assessment. First, a second-order MCS is employed to separately propagate the knowl-edge uncertainty and natural variability of the system. Second, the effect of parameter correlations and the associated value of information are investigated. Four multivariate methods are used to generate correlated model parameters. The outcome of the second-order MCS, called “spaghetti plot,” is used to assess the un-certainty, variability, and risk of salmonid embryo sur-vival. The following conclusions can be drawn from this study.
(1) Despite the extreme sampling effect demonstrated in the spaghetti plot, the 90% uncertainty inter-vals are consistent for all correlated scenarios and thus provide a commendable means to analyze the survival-rate uncertainty. The survival-rate
uncer-tainty increases with the amount of sand deposit. The survival-rate uncertainty associated with the independent scenario is overestimated for the lower values of specific sand deposit but underes-timated for those worst sand-deposit conditions. (2) The distribution of the survival-rate variability is
positively skewed, mostly ranging between 0 and 20%. The uncertainty of the survival-rate variabil-ity associated with the independent scenario is un-derestimated.
(3) The ranges of favorable habitat conditions for the correlated scenario are significantly smaller than the corresponding ranges for the independent one. The favorable pressure-head and grain-size ratios are negatively correlated; the favorable grain-size and flow-path ratios are positively correlated, whereas no apparent correlation is observed be-tween the favorable pressure-head and flow-path ratios. The independent scenario overestimates the negative correlation between the favorable pressure-head and grain-size ratios, but under-estimates the positive correlation between the favorable grain-size and flow-path ratios, both of which would lead to ineffective planning of habitat improvement measures.
(4) The maximum tolerable sand deposit for achiev-ing the target survival rate decreases with the increasing target value. The results for the inde-pendent scenario are more conservative at those higher target values. The risk for not achieving the target survival rate increases with the speci-fied target value and the amount of sand deposit. Significant values of information are associated with the parameter correlations, especially when the higher levels of embryo survival are targeted. (5) The scatter plots of the random numbers gener-ated with the four methods demonstrate different correlation patterns. An important finding of this study is that the outputs of various convolutions involved in the uncertainty/variability analyses and risk assessments are consistent for all corre-lated scenarios, regardless of the diverse correla-tion patterns. The outcomes associated with the independent scenario, however, substantially dif-fer from those of the correlated ones, indicating that parameter correlations should be incorporated into the Monte Carlo risk assessments whenever possible.
In this paper we establish a paradigm for uncer-tainty/variability analyses in the ecological context. Presently a systematic framework for such analyses is not available, thus the procedure proposed in this work may provide useful techniques to ecological risk assessors. More field data should be collected in the future, and alternative PDFs of model parameters may be used to investigate the distribution effect.
Acknowledgements
This research was supported by the National Sci-ence Council, Republic of China through Ecological Engineering Grant. We appreciate comments from two anonymous reviewers and the editor, which helped im-prove the clarity of the paper.
Appendix A.
A.1. Iman–Conover method (IC method)
The procedure presented by Iman and Conover (1982)is the most well known method for generating correlated random numbers. Existing risk assessment software, such as Crystal Ball or @Risk, uses IC method to produce correlated random numbers with any set of marginal distributions and specified values of rank–order correlation. The theoretical basis of IC method is that the independent random numbers can be transformed into the correlated ones using the orthogonal transformation. Given N sets of k inde-pendent random numbers (X1, . . . , Xk), expressed by aN × k matrix X, one can transform them into the correlated random numbers corresponding to a given rank–order correlation matrix R∗ by X∗ = X · LT, where X∗ is a N × k matrix of correlated random numbersX∗1, . . . , X∗k, L is the lower triangular matrix of R∗, obtained by Cholesky decomposition R∗ =
L· LT. Generation of correlated random numbers by IC method also involves the use of van der Waerden scores a(i), defined as a(i) = Φ−1(i/(N + 1)) for
i = 1, . . . , N, where Φ−1 is the inverse SN CDF.
The stepwise procedure of IC method is as follows:
• Step 1: Generate a N × k random score matrix A,
each row of which consists ofk independent
com-ponents randomly sampled from van der Waerden scores.
• Step 2: Transform the independent score matrix A
into the correlated score matrix A∗by A∗= A · LT. Then the row components of A∗have a rank–order correlation matrix M close to the given R∗.
• Step 3: Generate a N × k matrix X of
indepen-dent random numbersX1, . . . , Xk, each column of which contains N random numbers sampled from the corresponding marginal distribution.
• Step 4: Rearrange the components in each column
of X such that they have the same rank order as the corresponding column of A∗. The resulting matrix
X∗containsN sets of random numbers X∗1, . . . , X∗k, which also have a rank–order correlation matrix M close to R∗.
• Step 5: Check if the means and standard deviations
of the generated random numbers X∗1, . . . , X∗k converge to the prescribed values. Throughout this study, the convergence criterion is set as
|ε| ≤ 0.015, where ε is the difference between the
specified and sample values. If the convergence test is passed, proceed to step 6. Otherwise, reject the random numbers, and return to step 1 with
N = N + 100.
• Step 6: Use Kolmogorov–Smirnov (K–S) test to
ex-amine if the generated random numbers preserve the given marginal distributions. Throughout this study, a significance level of 0.05 is used in the K–S test. If the goodness-of-fit test is passed, proceed to step 7. Otherwise, reject the generated random numbers, and return to step 1 withN = N + 100.
• Step 7: Check the correlations among the generated
random numbers. Throughout this study, the accep-tance criterion is set as|εc| ≤ 0.05, where εcis the difference between the specified and observed val-ues of correlation. If the correlation test is passed, then stop. Otherwise, reject the generated random numbers, and return to step 1 withN = N + 100. A.2. Standard normal transformation method (SNT method)
The SNT method is based on the transformation of random variates between the SN and marginal prob-ability spaces (Der Kiureghian and Liu, 1986; Todd and Ng, 2001). For two correlated random numbersXi andXj in the marginal space, their product–moment
correlation coefficient,rij, is related to the equivalent correlation coefficient in the SN space, rN,ij, as the following: rij= ∞ −∞ ∞ −∞ x i− µi σi x j− µj σj
× φij(yi, yj, rN,ij) dyidyj (A.1) where µi and σi are the mean and standard devia-tion of Xi, respectively; φij is the bivariate SN joint PDF;Yi andYj are two correlated SN random vari-ates, obtained by the SN inverse transformationYi =
Φ−1(Fi(Xi)), in which Fiis the marginal CDF ofXi.
The joint PDFφij can be expressed as:
φij(yi, yj, rN,ij) = 1 2π 1− rN,ij2 exp −y 2 i − 2rN,ijyiyj+ y2j 2(1 − rN,ij2 ) (A.2) for−1 < rN,ij < 1. Given the correlation coefficient
rij and marginal distributions, the equivalent correla-tion coefficient rN,ij can be solved using Eqs. (A.1) and (A.2), which can be expressed as a simplified semi-empirical form (Der Kiureghian and Liu, 1986):
rN,ij= Tij· rij (A.3)
whereTij is the transformation factor. For the corre-lations and marginals given in Section 2.2, the val-ues ofT12,T13, andT23 so obtained are 1.002, 1.004, and 1.002, respectively, in which the subscripts 1, 2, and 3 correspond to the parametersds/Dg,h/L1, and
L2/L1, respectively. The Tij-values given above are close to unity, indicating that the differences between
rij andrN,ij are negligible in our setting. Generation of correlated random numbers with SNT method is proceeded as follows:
• Step 1: Transform the given product-moment
corre-lation matrix R between marginals into the correla-tion matrix RN between SN usingEq. (A.3).
• Step 2: Generate N sets of independent SN random
numbers(z1, . . . , zk).
• Step 3: Transform the independent SN random
vari-ates Zi into the correlated SN random variates Yi by orthogonal transformation Y= L · Z, where L is the lower triangular matrix of RN, obtained by Cholesky decomposition RN= L · LT.
• Step 4: Transform the correlated SN random variates Yi into the correlated marginal random variatesXi, using the inverse transformationXi= Fi−1[Φ(Yi)].
• Step 5: Check the convergence of the means and
standard deviations of the generated random num-bersXi. If the test is passed, proceed to step 6. Oth-erwise, reject the generated random numbers and return to step 2 withN = N + 100.
• Step 6: Use K–S test to examine the goodness of
fit. If passed, proceed to step 7. Otherwise, reject the generated random numbers and return to step 2 withN = N + 100.
• Step 7: Check the correlations rij among the gen-erated random numbers. If the correlation test is passed, then stop. Otherwise, reject the generated random numbers and return to step 2 with N =
N + 100.
A.3. Normal copula method (NC method)
The multivariate normal copula has been used in the decision and risk analysis (Clemen and Reilly, 1999). Like other copula families, the multivariate normal copula allows any marginal distributions. In addition, the normal copula permits the use of multivariate pairwise correlations. The flexibility and analytical tractability of the multivariate normal copula make it an attractive method for generating correlated random numbers. The basic idea of a copula is that the joint distribution of random variates can be expressed as a function of the marginal distributions (Sklar, 1959), i.e.
F(x1, . . . , xk) = C[F1(x1), . . . , Fk(xk)] (A.4) where F(x1, . . . , xk) is the joint CDF for ran-dom variates X1, . . . , Xk with marginal CDFs
F1(x1), . . . , Fk(xk); and C, called a copula, is a joint
CDF with uniform marginals on [0,1]. Given that
Fi(xi) and C are both differentiable, the joint PDF
f(x1, . . . , xk) can be expressed as:
f(x1, . . . , xk) = f1(x1) × · · · × fk(xk)
× c[F1(x1), . . . , Fk(xk)] (A.5) where fi(xi) is the marginal PDF corresponding to
Fi(xi), and c = ∂kC/(∂F1· · · ∂Fk) is the copula
density. Eq. (A.5) reveals that the joint PDF is a product of the marginal PDFs and the copula density.
For independent Xi’s, c = 1 and f(x1, . . . , xk) =
f1(x1) × · · · × fk(xk), which implies that c encodes
information about the dependence (or correlation) among theXi’s. Hence,c is also called a dependence function.
The multivariate normal distribution is typically parameterized with a matrix RN of product-moment correlation coefficients rN,ij, while the multivariate normal copula is defined by a rank–order correla-tion matrix R∗ (Clemen and Reilly, 1999). Solving Eq. (A.5)for the normal copula density cN in terms of the multivariate SN joint PDF φ(y1, . . . , yk|RN) yields:
cN[Φ(y1), . . . , Φ(yk)|R∗]= φ(y1, . . . , yk|RN) [φ(y1) × · · · × φ(yk)]
(A.6) where Φ(Yi) and φ(Yi) are the univariate SN CDF and PDF, respectively, in which Yi is obtained by
Yi = Φ−1[Fi(Xi)]. Substituting the SN PDFs into Eq. (A.6) and using the linear algebra, Clemen and Reilly (1999)obtained an expression forcN:
cN[Φ(y1), . . . , Φ(yk)|R∗]
= exp( − yT(R−1N − I)y/2)
|RN|1/2 (A.7)
where y = (y1, . . . , yk)T, I is the k × k identity matrix, and |RN| is the determinant of RN. The nor-mal copula density given in Eq. (A.7) can be used to generate correlated random numbers by an accep-tance/rejection procedure (Gentle, 1998; Haas, 1999), as described below:
• Step 1: Transform the given product–moment
cor-relation matrix R between marginals into the corre-lation matrix RNbetween SN usingEq. (A.3).
• Step 2: Determine the maximum cN-value on the k-dimensional unit hypercube, and call thiscN,max, wherecN is calculated withEq. (A.7).
• Step 3: Draw N sets of uniform random numbers (u1, . . . , uk, ξ) from [0,1].
• Step 4: For each set, if ξ ≤ cN(u1, . . . , uk)/cN,max, then accept(u1, . . . , uk); otherwise reject this set, wherecN(u1, . . . , uk) is evaluated withEq. (A.7). The accepted uniform random numbers are trans-formed to the marginal usingXi= Fi−1(Ui).
• Step 5: Check if the means and standard deviations
of the generated random numbers Xi converge to the prescribed values. If the test is passed, proceed to step 6. Otherwise, reject the generated random numbers and return to step 3 withN = N + 100.
• Step 6: Use K–S test to examine the goodness of
fit. If passed, proceed to step 7. Otherwise, reject the generated random numbers and return to step 3 withN = N + 100.
• Step 7: Check the correlations rij among the gen-erated random numbers. If the correlation test is passed, then stop. Otherwise, reject the generated random numbers and return to step 3 with N =
N + 100.
A.4. Maximum-entropy copula method (MEC method)
The piecewise-uniform copula presented by MacKenzie (1994) is applicable to generating cor-related random numbers. Determination of the piecewise-uniform copula density is based on the principle of entropy maximization. Shannon (1949) proposed a measure of the expected amount of in-formation contained in an observed probability dis-tribution, which is called entropy and defined as the following for a discrete case:
entropy= −
k
i=1
piln(pi) (A.8)
wherep1, . . . , pkare the probabilities corresponding to a set of events {E1, . . . , Ek} such thatpi = 1. For a set of random numbers{X1, . . . , Xk} having a joint PDFf(x1, . . . , xk), the entropy is defined as: entropy(f) = − · · · ln(f(x1, . . . , xk)) · f(x1, . . . , xk) dx1· · · dxk (A.9) Maximizing entropy subject to the given constraints is considered by many to be an intuitively compelling paradigm. By doing so one would be maximizing the amount of information to be learned from ob-serving a manifestation fromf . Based on Eq. (A.5) and the definition of entropy given by Eq. (A.9), the following relation is obtained (MacKenzie, 1994), which expresses the entropy of the joint PDF with the entropies of the marginals and of
the copula, i.e. entropy(f) =
k
i=1
entropy(fi) + entropy(c) (A.10) where entropy(c) = − · · · ln(c(u1, . . . , uk))
× c(u1, . . . , uk) du1· · · duk (A.11) in which c(u1, . . . , uk) is the joint copula den-sity. As revealed by Eq. (A.10), if the marginal distributions are known, then the problem of max-imizing entropy reduces to the problem of finding the maximum-entropy copula subject to the given constraints on dependence. The piecewise-uniform copula provides a useful tool for solving such a problem.
The k-dimensional piecewise-uniform copula den-sity ck,m(u1, . . . , uk) is constructed by partitioning each marginal on the unit hypercube intom equal-size pieces. The resultingmkcombinations of intervals are denoted as Ii1,... ,ik, where the subscript ij indicates that the uniform random variableUjbelongs to the in-terval [(ij− 1)/m, ij/m], for ij= 1, . . . , m and j = 1, . . . , k. Each Ii1,... ,ik is associated with a constant density, as expressed by:
ck,m(u1, . . . , uk) = pi1,... ,ik (A.12)
The piecewise-uniform copula densitypi1,... ,ik is sub-ject to the following constraints:
pi1,... ,ik ≥ 0 (A.13) and pi,... mk = p... ,i,... mk = p... ,i mk = 1 m fori = 1, . . . , m (A.14) wherep... denotes the summation of piecewise cop-ula densities over the remaining mk−1 combinations of subscript.Eqs. (A.13) and (A.14)are to assure that
pi1,... ,ikis a copula density and has uniform marginals.
An additional constraint is provided to model the de-pendence among thek(k − 1)/2 pairs of random vari-ates, which is given by the following formula for Spearman’sρ in terms of the copula density (Nelsen, 1999). For the multivariate case, the correlation coef-ficientρjlbetween two random variatesUjandUlcan
be expressed as: ρjl= 12 1 0 · · · 1 0 uj· ul × c(u1, . . . , uj, . . . , ul, . . . , uk) du1· · · duk−3 (A.15) Using the piecewise-uniform copula density, we turn Eq. (A.15)into the following form:
ρjl= 12 1 0 1 0 uj· ul m i1=1 · · ·m ik=1 pi1,... ,ij,... ,il,... ,ik × i1/m (i1−1)/m · · · ik/m (ik−1)/m du1· · · duk ×dujdul−3 = 3 mk+2 m ij=1 m il=1 m i1=1 · · ·m ik=1 pi1,... ,ij,... ,il,... ,ik × (2ij− 1)(2il− 1) −3 (A.16)
Moreover, with the piecewise-uniform copula density, Eq. (A.11)can be rewritten as:
entropy(c) = − 1 mk m i1=1 · · · m ik=1 ln(pi1,... ,ik) · pi1,... ,ik (A.17) MacKenzie (1994)has shown that the optimization problem defined by maximizing Eq. (A.17) sub-ject to the constraints given by Eqs. (A.13), (A.14) and (A.16) has a unique global solution. The mk piecewise-uniform copula densitiespi1,... ,ik so deter-mined form a basis for generating correlated random numbers. Before we proceed to the generation proce-dure, it would be appropriate to give a criterion here for determining the minimum value ofm. It is shown that±(1 − 1/m2) are the limiting values of pairwise correlation captured by a piecewise-uniform copula (MacKenzie, 1994), thusm must be large enough to assure that all|ρjl| ≤ 1 − 1/m2. For the correlations given in Table 2,m = 4 may be taken. The accep-tance/rejection procedure used in the NC method is also adopted here for generating correlated random numbers.
• Step 1: Given the rank–order correlation matrix R∗
among the random variates, solve the optimization problem defined by maximizing the entropy given in Eq. (A.17) subject to the constraints given in Eqs. (A.13), (A.14) and (A.16).
• Step 2: Find the maximum pi1,··· ,ik-value obtained
from step 1, and call thisck,m,max.
• Step 3: Draw N sets of uniform random numbers (u1, . . . , uk, ξ) from [0,1].
• Step 4: For each set, if ξ ≤ ck,m(u1, . . . , uk)/
ck,m,max, then accept (u1, . . . , uk); or else reject this set, where ck,m(u1, . . . , uk) is the pi1,··· ,ik obtained from step 1.
• Step 5: Check the correlations ρjl among the ac-cepted uniform random numbers. If the correlation test is passed, proceed to step 6. Otherwise, reject the uniform random numbers and return to step 3 withN = N + 100.
• Step 6: The accepted uniform random numbers are
then transformed into the marginal random numbers usingXj= Fj−1(Uj).
• Step 7: Check the convergence of the means and
standard deviations of the generated random num-bersXj. If the test is passed, proceed to step 8. Oth-erwise, reject these random numbers and return to step 3 withN = N + 100.
• Step 8: Use K–S test to examine the goodness of fit.
If passed, then stop. Otherwise, reject the generated random numbers and return to step 3 with N =
N + 100.
References
Anderson, E.L., Hattis, D., 1999. Uncertainty and variability. Risk Anal. 19, 47–49.
Annan, J.D., 2001. Modelling under uncertainty: Monte Carlo methods for temporally varying parameters. Ecol. Model. 136, 297–302.
Bogen, K.T., Spear, R.C., 1987. Integrating uncertainty and interindividual variability in environmental risk assessment. Risk Anal. 7, 427–436.
Bukowski, J., Korn, L., Wartenberg, D., 1995. Correlated inputs in quantitative risk assessment: the effects of distributional shape. Risk Anal. 15, 215–219.
Burmaster, D.E., Anderson, P.D., 1994. Principles of good practice for the use of Monte Carlo techniques in human health and ecological risk assessment. Risk Anal. 14, 477–481. Burmaster, D.E., Wilson, A.M., 1996. An introduction to
second-order random variables in human health risk assessments. Hum. Ecol. Risk Assess. 2, 892–919.
Clemen, R.T., Reilly, T., 1999. Correlations and copulas for decision and risk analysis. Manag. Sci. 45, 208–224. Cooper, A.C., 1965. The effect of transported stream sediments on
survival of sockeye and pink salmon eggs and alevins. In: Int. Pac. Salmon Commun. Bull., vol. 18. New Westminister, BC. Der Kiureghian, A., Liu, P.L., 1986. Structural reliability under incomplete probability information. J. Eng. Mech. 112, 85–104. Frey, H.C., Rhodes, D.S., 1996. Characterization, simulation and analyzing variability and uncertainty: an illustration of methods using an air toxics emission example. Hum. Ecol. Risk Assess. 2, 762–797.
Gentle, J.E., 1998. Random Number Generation and Monte Carlo Methods. Springer, New York.
Haas, C.N., 1999. On modeling correlated random variables in risk assessment. Risk Anal. 19, 1205–1214.
Harmon, R., Challenor, P., 1997. A Markov chain Monte Carlo method for estimation and assimilation into models. Ecol. Model. 101, 41–59.
Hattis, D., Burmaster, D.E., 1994. Assessment of variability and uncertainty distributions for practical risk analyses. Risk Anal. 14, 713–730.
Hession, W.C., Storm, D.E., Haan, C.T., Burks, S.L., Matlock, M.D., 1996. A watershed-level ecological risk assessment methodology. Water Resour. Bull. 32, 1039–1054.
Hoffman, F.O., Hammonds, J.S., 1994. Propagation of uncertainty in risk assessments: the need to distinguish between uncertainty due to lack of knowledge and uncertainty due to variability. Risk Anal. 14, 707–712.
Iman, R.L., Conover, W.J., 1982. A distribution-free approach to inducing rank order correlation among input variables. Commun. Stat.—Simul. Comput. 11, 311–334.
MacKenzie, G.R., 1994. Approximately maximum-entropy multivariate distributions with specified marginals and pairwise correlations. Ph.D. Dissertation, University of Oregon. Mitchell, S.W., Csillag, F., 2001. Assessing the stability and
uncertainty of predicted vegetation growth under climatic variability: northern mixed grass prairie. Ecol. Model. 139, 101–121.
Moschandreas, D.J., Karuchit, S., 2002. Scenario-model-para-meter: a new method of cumulative risk uncertainty analysis. Environ. Int. 28, 247–261.
Nauta, M.J., 2000. Separation of uncertainty and variability in quantitative microbial risk assessment models. Int. J. Food Microbiol. 57, 9–18.
Neave, H.M., Cunningham, R.B., Norton, T.W., Nix, H.A., 1997. Preliminary evaluation of sampling strategies to estimate the species richness of diurnal, terrestrial birds using Monte Carlo simulation. Ecol. Model. 95, 17–27.
Nelsen, R.B., 1999. An Introduction to Copulas. Springer, New York.
Pohlmann K.F., Hassan A.E., Chapman J.B., 2002. Modeling density-driven flow and radionuclide transport at an underground nuclear test: uncertainty analysis and effect of parameter correlation. Water Resour. Res. 38(5), doi:10.1029/ 2001WR001047.
Pouillot, R., Albert, I., Cornu, M., Denis, J.B., 2003. Estimation of uncertainty and variability in bacterial growth using Bayesian inference: application to Listeria monocytogenes. Int. J. Food Microbiol. 81, 87–104.
Qian, S.S., Stow, C.A., Borsuk, M.E., 2003. On Monte Carlo methods for Bayesian inference. Ecol. Model. 159, 269–277. Rai, S.N., Krewski, D., 1998. Uncertainty and variability analysis
in multiplicative risk models. Risk Anal. 18, 37–45. Rai, S.N., Krewski, D., Bartlett, S., 1996. A general framework for
the analysis of uncertainty and variability in risk assessment. Hum. Ecol. Risk Assess. 2, 972–989.
Salinger, D.H., Anderson, J.J., Hamel, O.S., 2003. A parameter estimation routine for the vitality-based survival model. Ecol. Model. 166, 287–294.
Sanga, R.N., Bartell, S.M., Ponce, R.A., Boischio, A.A.P., Joiris, C.R., Pierce, C.H., Faustman, E.M., 2001. Effects of uncertainties on exposure estimates to methylmercury: a Monte Carlo analysis of exposure biomarkers versus dietary recall estimation. Risk Anal. 21, 859–868.
Shannon, C.E., 1949. The Mathematical Theory of Communi-cation. University of Illinois Press, Urbana.
Sklar, A., 1959. Fonctions de répartition à n dimensions et leurs marges. Publications de l’Institut Statistique de l’Université de Paris 8, 229–231.
Smith, A.E., Ryan, P.B., Evans, J.S., 1992. The effect of neglecting correlations when propagating uncertainty and estimating the population distribution of risk. Risk Anal. 12, 467–474. Tappel, P.D., Bjornn, T.C., 1983. A new method of relating size
of spawning gravel to salmonid embryo survival. North Am. J. Fisheries Manag. 3, 123–135.
Todd, C.R., Ng, M.P., 2001. Generating unbiased correlated random survival rates for stochastic population models. Ecol. Model. 144, 1–11.
Vose, D., 2000. Risk Analysis. John Wiley & Sons, Chichester. Waller, L.A., Smith, D., Childs, J.E., Real, L.A., 2003. Monte
Carlo assessments of goodness-of-fit for ecological simulation models. Ecol. Model. 164, 49–63.
Wu, F.C., 2000. Modeling embryo survival affected by sediment deposition into salmonid spawning gravels: application to flushing flow prescriptions. Water Resour. Res. 36, 1595– 1603.
Wu, F.C., Shen, H.W., 1999. First-order estimation of stochastic parameters of a sediment transport model. J. Hydraulic Res. 37, 213–227.
