Improved Feasibility of Joint Uncertainty and Sensitivity Analyses Performance for Complex Scenario of Accidental Radioactivity Release into the Calm Atmosphere

. During several hours of the calm meteorological situation, a relatively significant level of radioactivity can be accumulated around the source. At the second stage, the calm situation is assumed to terminate and convective movement of the air induced by wind immediately starts. Random realisations of the input atmospheric dispersion model parameters for the CALM scenario are generated using LHS (Latin Hypercube Sampling) scheme. The resultant complex random radiological trajectories passing through both calm and convective stages of the release scenario represent inevitable prerequisite for prospective uncertainty analysis (UA) and sensitivity analysis (SA). Novel concept of Approximate Based (AB) solution approximates non-Gaussian sum of individual puffs at the end of the calm period by only one Gaussian “superpuff” distribution. Substantial acceleration of generation of sufficiently large number of random realisation makes further UA and SA feasible. Both procedures come from common mapping of the pairs of matrix random dependant output fields and vector of random input parameters realisations. Examples of 2-D random trajectories of deposited 137 Cs are presented. Global sensitivity analysis utilising random sampling methods is outlined.


Definition of the task
Deterministic analysis of the accidental release of radioactive substances into the motionless (calm) atmosphere followed by radioactivity dissemination by wind into the living environment is introduced in (Pecha and Kárný, 2020). This follow-up study provides basic form of probabilistic analysis of the CALM scenario. On basis on uncertainty and sensitivity studies is shown that thanks to novel "Superpuff " concept, the sampling-based Monte Carlo procedures are fairly feasible even for high numbers of random realisations.
Before anything else we should stress the difference between variability and uncertainty of a certain variable. Variability reflects changes of a certain quantity over time, over space or across individuals in population. Variability represents diversity or heterogeneity in a well characterized population. The term "uncertainty" covers stochastic uncertainties, structural uncertainties representing partial ignorance or incomplete knowledge associated with the lack of perfect information about poorly-characterized phenomena or models and input model uncertainties.
Within the frame of selected dispersion model, many uncertainties related to both conceptual model (parametrization errors, uncertain submodel parameters, stochastic nature of some measured input data,…) and computational scheme (step of computation net, averaging land-use characteristics, averaging time for dispersion parameters etc.) are involved. For particular definition of the limited group of input random parameters and their ranking we have accepted the former results of extensive literature review of the analogous international environmental codes, e.g. UFOMOD, COSYMA, MARC-2A, OSCAAR, NPK-PUFF, summarized in (Pecha and Pechova, 2005). The random characteristics should be regularly refreshed according to recent recommendations from the elicitation procedures of experts (the weak point of this study).
Let X  {X1, X2, …, XN } denotes a vector of N random input model parameters Xi with the corresponding sequence of random distributions D1, D2, …, DN which are usually selected (range, type of distribution, potential mutual dependencies) on the basis of commonly accepted agreement of experts (elicitation procedures). Output random fields from the CALM dispersion model can be schematically differentiated as alternatives: spatial fields Y k k=1, … , K generated simultaneously on the polar calculation grid  The number K from Eq. (1.c) is rather large and represents the possible (radiological) output results generated on the polar computational grid consisting of 42 concentric circles up to 100 kilometres from the source of pollution and 80 angular beams (clockwise, from North). These output matrices have dimensions 80  42 = 3 360 elements (see Figure 1). The examples of various fields are:  Radioactivity concentration of radionuclides in the air and deposited radioactivity on the ground (analysed in this article), their time integrals;  Doses of external irradiation;  Internal irradiation from inhalation of contaminated air;  Internal irradiation from ingestion of contaminated food;  …. and many others.
1 The index k stands here for 2-D radiological field of the deposited radioactivity of 137 Cs on the ground (in the discrete computational nodes according to Figure 1)

Generation of random trajectories of the output fieldssampling-based methods
Considerations in this article are limited to one output random field (matrix or "rose" on the polar computational grid ) of the radionuclide 137 Cs deposited on the ground. Detailed description of the radionuclide propagation within the accidental radioactivity release of the CALM scenario (dissemination into the motionless atmosphere suddenly ensuing by wind) is given in the above mentioned study (Pecha and Kárný, 2020). The radioactivity [Bq m -2 ] deposited on the ground is examined during the complex CALM scenario consisting of the 5 hours discharges into the motionless atmosphere immediately succeeded by 4 hours of ensuing windy transport.
Due to the excessive number of the input model parameters only those with significant effects of their fluctuations on the outputs are assumed and they enter the analysis as random. Necessary reduction of the input random vector dimension to N is demonstrated on the scheme: where x b stands for best estimated (nominal) values. Y is assumed as 137 Cs deposited on the ground in points of the computational grid , marked as (Y 137 ) .
Sampling-based method for UA and SA includes: (1) Definition of probability distributions to characterize uncertainty in analysis inputs, (2) Generation of samples from uncertain analysis inputs, (3) Propagation of sampled inputs through an analysis, (4) Presentation of uncertainty analysis results, and (5) Determination of sensitivity analysis results. In practice, it consists in multiple repetitions of calculations of outputs successively for each specific sample of random input vector, specifically:  Generation of a particular j-th sample of the input vector x j  x1 j , x2 j , … , xN j  where xi j are realisations of the input random parameter X j (successively for j=1,2,…,NREAL).  Propagation of the sample j through the model, it means calculation of the corresponding resulting realisation y j  of random output value Y (output vector on the polar computational grid see Figure 1) when running the model y j  =  CALM (x1 j , x2 j , … , xN j ).
It results in the mapping of pairs: The expression (3) represents the key scheme for uncertainty analysis and sensitivity studies. Statistical processing of the pairs (3) can determine the extent of the uncertainty on predicted consequences and yields various statistics such sample mean and variance, percentiles of the uncertainty distribution on the quantity given, uncertainty factors, reference uncertainty coefficients etc. Simulation of uncertainty propagation through the model has a cardinal importance for introduction of advanced methods in modelling, because:  offers essential data for transition from deterministic procedure of consequence assessment to probabilistic approach which enables to generate more informative probabilistic answers on assessment questions,  provides detailed analysis of the model error covariance structure thus making possible to improve the reliability of model predictions on basis of application of advanced statistical techniques of assimilation of mathematical prognoses with real measurements incoming from terrain.

Uncertainty and sensitivity analysis procedures
Resultant mapping of pairs given schematically by (3) provides mutual bases for:  Uncertainty analysis (UA)statistical processing of the pairs provides quantification of uncertainties and determines an extent of the uncertainty on predicted consequences. It yields various statistics such sample mean and variance, percentiles of the uncertainty distribution on the quantity given, uncertainty factors, reference uncertainty coefficients etc.  Sensitivity analysis (SA)determines how different values of an independent variable affect a particular dependent variable. Its strategies are applied depending on the settings. An embracive outstanding study of the SA basis is in (Saltelli A., K. Chan and E. M. Scott, 2001). Many authors use the term "sensitive" when referring to the degree to which an input parameter affects the model output. The basic methods include e.g. the differential sensitivity analysis, one-at-a time sensitivity measures, factorial design, sensitivity and importance index, importance factor or subjective sensitivity analysis. Significant role play methods of parameter sensitivity analysis utilising random sampling methods providing scatter plots, importance indices, relative deviation ratio, correlation coefficients, rank transformations or various regression techniques often used to replace a highly complex model with a simplified "respond surface".
Roughly speaking, uncertainty analysis quantifies the variation in the model output, while sensitivity analysis identifies the input sources potentially responsible for this uncertainty. UA is not same as SA, but, ideally, both analyses should be run in tandem (CC-MOD, 2020, Eur. Commission). A survey of sampling-based methods for calculations of the joined UA-SA analyses is reviewed in (Helton et al., 2006). UA and global SA often use similar mathematical techniques. Close relation of SA with UA of numerical models is emphasized in (Pianosi et al., 2016). Propagation of uncertainty by Monte Carlo simulation is also used to perform the initial steps of sampling-based sensitivity analysis. SA can address miscellaneous questions, e.g. what input factors cause the largest variation in the output, what factors are negligible or if some factors can amplify or dampen output variability (Pianosi et al., 2016). UA is recommended to precede a SA.
There are mainly two approaches to analysing sensitivity:  Local SA considers the output variability against variations of an input factor around a specific best estimated value x best (or x ). The derivatives are taken at a single point and an impact of one parameter on the cost function at a time is analysed keeping the other parameters fixed. Local SA investigates the effect of variation of uncertain inputs from a baseline point. The sensitivity measure S for the i-th input factor takes the form The partial derivatives are usually approximated by finite differences. Function g stands for the CALM model  CALM . The test of local sensitivity examines the change in output as each parameter is individually increased by a factor of its standard deviation (±SD).
 Global sensitivity analysis (GSA) is based on M-C sampling and considers variations within the entire space of variability of the input factors. Profound SA knowledge is given in the principal publication (Saltelli et al.,2001). GSA investigates the effects of variation of certain inputs across their entire variability space and naturally overcomes nonlinear relationships of the CALM scenario. Three GSA methods are discussed in (Sarrazin et al., 2016). Variance-based SA method comes out from variance decomposition, Elementary Effect Test (or method of Morris) improves performance when the computing time of a single model run is high, or when the number of factors is very large. Regional Sensitivity Analysis (also called Monte Carlo filtering) is mentioned in (Pianosi et al., 2016) as a family of methods mainly aimed at identifying regions in the inputs space corresponding to particular values (e.g. high or low) of the output, and that can be used for mapping and for dominant controls analysis.
Sensitivity analysis of environmental models is treated in (Sarrazin et al., 2016) from the points of sample size and threshold for the identification of insensitive input factors. In (Hamby, 1994), the fundamental SA techniques are overviewed for the various modelling situations. A method of full decomposition of variance of dependant output is introduced in , where the computational cost is relatively low.

Constraints in choice of method of sensitivity analysis.
The most common constraints are: the computational expense, completeness of the input random parameters space, correlated inputs, nonlinearity (inaccurately applied sensitivity measures of linear regression), model interactions during simultaneous input parameters effect or complexity connected with multiple outputs. A number of sensitivity analyses techniques are reviewed in the preceding publication (Hamby, 1994). The simplest one-at-atime method varies each parameter while all others are held constant. Further sensitivity analyses techniques examine parameter influence based on output variation while simultaneously changing input parameters. It enables to build respond surface that approximates complex model. A comprehensive review of the fundamental research related to introduction and utilisation of the random sampling methods is presented onward. (Razavi and Gupta, 2015) revisit the theoretical basis for SA, critically evaluate existing approaches in the literature and demonstrate their flaws and shortcomings through the published examples. A systematic review of sensitivity analysis practices is guided in (Saltelli et al., 2019). A question "Why so many published sensitivity analyses are false?" is opened. Many published SA fail the elementary requirement to properly explore the space of the input random factors. Warning before perfunctory settings of the sensitivity scenario is emphasized in (Saletelli, 2010). One should be aware of pitfalls and avoid the false procedures. (Pianosi et al., 2016) warns that conclusions drawn from GSA should be taken with care when the input variability space is poorly known.
When sensitivity indices cannot be computed analytically, sampling-based sensitivity analysis must be used. Typically, the number of model evaluations NREAL increases with the number of input factors N subject to SA. However, the ratio between NREAL and N significantly varies from one method to another and often also from one application to another. The combinations (NREAL,N) are usually taken as an compromise between the environmental model computational complexity and completeness of the input parameter uncertainty group.

What originality lies in our solution?
In this article we have proposed common basis for joint UA and SA analyses for the complex radioactivity release scenario CALM given by scheme (3). Following (Pecha and Kárný, 2020), at the calm period end CALM END T the radioactivity was gradually heaped up in the motionless atmosphere in the form of superposition of M Gaussian pulses in different stages of their dispersion. The superposition is definitely non-Gaussian which brings troubles to the further solution. Two alternative procedures of handling the considered ensuing convective transport are used: are projected on a single representative Gaussian super-puff distribution. Unlike BF, it requires only one-shot run modelling the convective transport. The benefit of the reduced computational load is evident, especially for a large number of puffs M and prospective demands on necessary application of more powerful but laborious dispersion codes. The reduction is paid by the introduced approximation error. Its extent, experimentally studied in (Pecha and Kárný, 2020), seems to be acceptable. The most important result lies in the conclusion: Thanks to AB approach, the UA and SA analyses for the CALM scenario are fairly feasible even for the relatively high numbers of random realisations NREAL, number of input random parameters N and number of puffs M.

Setting up the CALM scenario of radioactive discharges.
We are focusing on the near-field analysis below the mixing layer in a small spatial domain up to 100 km from the source of pollution. The considered area is demonstrated on Figure 1.

Definition of the specific uncertainty group of the atmospheric dispersion and deposition model of the CALM scenario
The computations mostly require to express the input random parameters X to be split in the form of two alternatives: for additive character: X = x best + f Random fractions f resp. f is random number with a suitable frequency distribution, x best is the best estimate (reference, nominal) value of the input (see figures for LHS sampling). The release scenario is composed of two parts: the accidental release of radioactive substances into the motionless (calm) atmosphere followed by radioactivity dissemination by wind into the living environment.

Input parameters selected as random for atmospheric dispersion model in the CALM region:
ADMXX(1) …. random fraction for horizontal dispersion parameter CALM y ADMXX(2) …. random fraction for vertical dispersion parameter CALM z ADMXX(3) …. absolute random value for gravitational settling velocity in the calm regionin LHS sampled from absolute intervalsee Panel_LHS_Part I. ADMXX(4) …. random fraction for scavenging coefficient for aerosol After 5 hours of the calm, the motionless atmosphere is lastly disseminated by wind to the surrounding environment. The radioactivity package accumulated locally during the last 5 hours is immediately drifted onward. We assume successive 4 hours of the convective transport. In the 4 th hour, atmospheric precipitation with intensity 1 mmh -1 occurres.

Latin hypercuble sampling (LHS) generation subsystemfor relevant input
parameters their random characteristics are specified filling up the PANELs LHS (see below) with additional more detailed options shown on the example of detailed random characteristics entering for ADMXX(7) -see example of detailed selection for random parameter ADMXX(7) from Panel LHS_Part I.
2. Special graphical subsystem for visualization -resulting output dependent fields y j  (from Eq. 2) are drawn as 2-D picturessee. Figures 2 and 3.
Adopted scheme of Monte Carlo modelling integrated into the CALM scenario uses stratified sampling procedure LHS ((Latin Hypercube Sampling) -see LHS Panel below). Multiple generation of NREAL samples from the various types of random distributions Dn of input random vector components Xn (n=1, … , N). A certain technique for correlation control between input components Xn is included. The LHS samples tend to produce more stable results than the common random samples.
Detailed options of the input parameters of the atmospheric and deposition model and its random characteristics are given in the following PANNEL. As evident, 5000 LHS samples of the sixteen values of the ADMXX(1:16) are generated.
The main contribution of this article relates to introduction of the "Superpuff" concept, leading to the substantial reduction of computational load of UA and SA of the CALM scenario. More attention should be devoted to other critical areas. The space of input random factors space should be properly explored with respect to their completeness 2 . Sometimes inputs can be strongly correlated or model responses can be nonlinear with respect to its inputs. Setting of the group of random input factors given on the Panels pictures above are coming out from the former international codes COSYMA, UFOMOD, MARC-2 and others analyses. The data were formerly collected and used in (Pecha and Pechova, 2005). Quality of the input parameter uncertainty group should be actualised of basis of the new knowledge and elicitation procedures.

Results for UA and SA of the accidental scenario CALM
Near-field analysis of the dispersion and deposition of radioactive admixtures below the mixing layer in a small spatial domain up to 100 km from the source of pollution (Figure 1) is examined. Two consecutive stages of numerical experiment are designed for a hypothetical release from the category of the worst-case analysis WVA (Weather Variability Assessment). In the first five hours, the release dissipates into the motionless ambient atmosphere just until to the calm-period termination. The total inventory of radionuclide 137 Cs is set as n TOT Q = 6.0 E+07 Bq. It is discharged into the motionless ambient air during the calm conditions lasting , when the calm condition ends, immediately the second convective phase of the transport induced by wind is assumed to start. Meteorological data for each hourly transport (wind speed and velocity, rain episodes, atmospheric stability class) are provided by the regional meteorological service.

UA results
The results are illustrated in the form of complex 2-D radiological trajectories (more detailed in  ). The values calculated for the best estimate of input parameters are demonstrated on Figure 2. UA analysis comes from pairs  y j ;x j  given by pattern according to Eq (3) . Influence of the input random parameters realisations is demonstrated on Figures 3.1  3.5. The ordinal number of realisation of the vector x j is marked as NRE. The pictures of y j  stands for several selected realisations NRE from the whole range of NRE <1;NREAL>). The complex random 2-D radiological trajectories are valid for the CALM scenario with 5 hours of radioactivity release into the motionless atmosphere immediately succeeded by 4 hours of convective windy transport (with rain [1 mm.h -1 ] in the 4 th (last) hour of the convective transport). In other words, the pictures can be taken as "snapshots" of the 137 Cs deposition on the ground related precisely to the beginning of the 10 th hour just after the release start.

Comments on SA results: Parameter sensitivity analysis utilising random sampling method
The input parameters whose uncertainty makes a major contribution to overall uncertainty can be identified using the correlation coefficients PRCCs, PCCs and SRCs. The identification of important contributors to variations of consequences is done using a sensitivity measure, so called partial (rank) correlation coefficient. It is possible to calculate the percentage contribution of each uncertain model parameter to uncertainty of consequences by use of so-called coefficients of determination (R 2 ). The simplest type of SA perturbs the input factors entering the model g from their nominal values one at a time. The measure of output sensitivity to the i-th input factor is based on the partial derivative of the dependant value (g) /xi at the nominal values of other factors k i . Many other tests with different degrees of complexity are used.
The most common sensitivity analysis is a sampling-based method. The model is running repeatedly for a combination of factor values from the respective distributions. Samplingbased methods do not require access to model equations. Based on this technique, we shall show the benefit of our research in the field of feasibility of the global SA. The main problem of SA concerns the computational cost of the analysis defined in terms of the number of evaluations of the environmental model of the CALM scenario needed to complete the analysis. The convergence analysis assesses whether sensitivity estimates are independent of the size of the in the input-output sample. The problem could appear even if the obtained samples are of medium size. Because of substantial acceleration due to implementation of the AB procedure (see Chapter 5), the problem of robustness and convergence of the CALM scenario connected with the sample size is not so critical.
The dependent output values are arranged into the matrix consisted from N  = 3360 elements (42 concentric circles  80 angular beams). N  can be understood as dimension of the output vector field on computational grid . We shall somewhat simplify our consideration on Eq. (1a) for a single point node on  in the form: (y Cs137 (node,x k )) node = yk node (6) yk node represents the deposition of 137 Cs on the ground exactly in the discrete spatial point node of the calculation polar grid  (Fig. 1) given for the k-th realisation x k  x1 k , x2 k , … , xN k  of the input random vector x.
As stated above, the sampling-based methods for UA consist in NREAL multiple repetitions of calculations of outputs, successively for each specific sample of random input vector. The sensitivity evaluation itself can be done successively at any point of the spatial computational grid. Various techniques can be used providing different measures of sensitivity (scatterplots, regression and correlation analysis, rank transformations etc.).
As a demonstration example we shall outline determination of the Pearson correlation coefficient c(xj , y, node) which expresses a measure of strength of the linear relationship between xj and y at the local computational point node. We use the form: The larger the absolute value of c( xj, y, node) the stronger the degree of linear relationship between the input and output values. The points labelled as node =1 ; node=17 ; node=20; node=25 are demonstrated on Fig.4. All are laying on the angular beam b=69 (around the maximum nominal values). Even in the simplified case of this point-based SA, a more detailed consideration should be mentioned. It turned out to be necessary to discriminate the nature of the action of the random input parameters xj as local or global. For example, from this point of view, the input random factors ADMXX(1), ADMXX(2), ADMXX(3) and ADMXX(4) are assumed act independently inside the calm region. On the other hand, the random inputs factors in the convective region can impact the convective region either globally (same in all phases -) or locally (independently in each convective phase -precipitation and speed and direction of wind are available for each hour of propagation). The sample-based technique leading to Eq.
(3) can naturally comply with the local or global character of input random parameters impact. Finally, some input factors are correlated mutually. Specifically, the horizontal dispersion in the convective region is evaluated according to the expression (see Eq. (4.2) in (Pecha and Kárný, 2021) ) from two parts: The first part stands for dispersion of the puff m from its birth at tm until t= CALM END T . Randomisation is expressed consistently according to the random scale factor ADMXX(1). The second part is the contribution to the dispersion along the total convective transport length Lp+l. The responsible random scale factor ADMXX(5) is assumed to act globally (in all four hourly phases equally). We cannot simply anticipate the correlations of such composite variables.

Several selected numerical results
We have focused on the input factors where is a straightforward relationship with the dependant variable of radioactivity deposition on the ground. x7 and x8 have global character (same in all convective phases p).

Conclusion
Uncertainty analysis of the CALM release scenario is based on scheme (3). Sensitivity analysis comes out thereof and its quality directly depends on degree of implied details. The UA can inherently include global or local character of the random input parameters (same or different fluctuations in separate spatial regions) or account for their correlations. The more detailed UA the more sophisticated SA can be accomplished. Even point-based SA can benefit from the sampling based scheme (3). It provides fast and effective estimation of the sensitivity measures in the various computational nodes.
The above facts are supported by numerical results. Positive correlations between random scale factors of the input parameters of deposition processes (gravitational settling, dry velocity deposition, washout from the plume by rain) and dependant random variable y of the 137 Cs radioactivity deposition on the ground have been confirmed. At the same time, the negative correlations for input parameters ADMXX(1) ( CALM yhorizontal dispersion) and ADMXX(2) ( CALM zvertical dispersion) comply with the local character of their effect limited to the calm region only.
The main benefit of this study can be declared as: UA and SA analyses for the CALM scenario are fairly feasible even for the sufficiently high numbers of random realisations NREAL and the number of puffs M. Findings from calculations on common notebook shown, that for high NREAL>10 3 and number of puffs M  100  (AB ) solution in convective area is roughly about M-times faster (instead about tens of hours for the BF, the AB solution with one Superpuff takes about tens of minutes).