F Cadini

J De Sanctis

I Bertoli

E Zio

UPSCALING OF A DUAL-PERMEABILITY MONTE CARLO SIMULATION MODEL FOR CONTAMINANT TRANSPORT IN FRACTURED NETWORKS BY GENETIC ALGORITHM PARAMETER IDENTIFICATION

Keywords: Radioactive waste repository, performance assessment, fracture networks, Upscaling, Monte Carlo simulation, genetic algorithms

The transport of radionuclides in fractured media plays a fundamental role in determining the level of risk offered by a radioactive waste repository in terms of expected doses. Discrete Fracture Networks (DFN) methods can provide detailed solutions to the problem of modeling the contaminant transport in fractured media.

However, within the framework of the performance assessment (PA) of radioactive waste repositories, the computational efforts required are not compatible with the repeated calculations that need to be performed for the probabilistic uncertainty and sensitivity analyses of PA. In this paper, we present a novel upscaling approach, which consists in computing the detailed numerical fractured flow and transport solutions on a small scale and use the results to derive the equivalent continuum parameters of a lean, one-dimensional Dual-Permeability, Monte Carlo Simulation (DPMCS) model by means of a Genetic Algorithm search. The proposed upscaling procedure is illustrated with reference to a realistic case study of Pu 239 migration taken from literature.

INTRODUCTION

The modeling of solute transport through networks of fractures is of extreme importance for the risk assessment of radioactive waste repositories [START_REF] Carr | Fractal characterization of joint surface roughness in welded tuff at Yucca Mountain, Nevada[END_REF][2][3][4] [START_REF] Reeves | Incorporation of conceptual and parametric uncertainty into radionuclide flux estimates from a fractured granite rock mass[END_REF].

Multiple-continuum models, based on advection-dispersion equations, allow to approximately describe the transport phenomena on large scales but may fail to capture certain non-Fickian features on small scales. [START_REF] Benke | Modelling conservative tracer transport in fracture networks with a hybrid approach based on the Boltzmann transport equation[END_REF] [START_REF] Painter | Upscaling discrete fracture network simulations: An alternative to continuum models[END_REF]. These models assume uniform dispersion, which is not always satisfied in fractured networks where flow and transport can take place in highly localized channels [START_REF] Benke | Modelling conservative tracer transport in fracture networks with a hybrid approach based on the Boltzmann transport equation[END_REF] [START_REF] Painter | Upscaling discrete fracture network simulations: An alternative to continuum models[END_REF] and may even give rise to apparent scale dependence of the dispersivity [START_REF] Benke | Modelling conservative tracer transport in fracture networks with a hybrid approach based on the Boltzmann transport equation[END_REF] [START_REF] Painter | Upscaling discrete fracture network simulations: An alternative to continuum models[END_REF][8] [START_REF] Vanderborght | Stochastic Continuum Transport Equations for Field-Scale Solute Transport: Overview of Theoretical and Experimental Results[END_REF].

Single-continuum models based on equivalent parameters may in general be exploited in case of very regular and well connected fracture networks, which can be approximated by a porous medium; however, significant biases may be introduced when more realistic fractures geometries are taken into account [START_REF] Chevalier | Stochastic delineation of wellhead protection area in fractured aquifers and parametric sensitivity study[END_REF].

A detailed characterization and mapping of the fractured fields is usually not feasible due to the high degree of complexity and heterogeneity of the network, and information on geometric and conducting properties can be obtained only in terms of probability distributions and correlation functions [START_REF] Benke | Modelling conservative tracer transport in fracture networks with a hybrid approach based on the Boltzmann transport equation[END_REF] [START_REF] Painter | Upscaling discrete fracture network simulations: An alternative to continuum models[END_REF][10] [START_REF] Long | From field data to fracture network modeling: an example incorporating spatial structure[END_REF] [START_REF] Reeves | Incorporation of conceptual and parametric uncertainty into radionuclide flux estimates from a fractured granite rock mass[END_REF].

Discrete Fracture Network (DFN) methods allow handling the high degree of complexity and heterogeneity of the networks of fractures [START_REF] Smith | An analysis of the influence of fracture geometry on mass transport in fractured media[END_REF] [START_REF] Long | Porous Media Equivalents for Networks of Discontinuous Fractures[END_REF]. These methods do not rely on volume averaging techniques, but solve the flow and transport problem in individual fractures by accounting for both anisotropy and stochasticity [START_REF] Long | Porous Media Equivalents for Networks of Discontinuous Fractures[END_REF] [START_REF] Cacas | Modeling fracture flow with a stochastic discrete fracture network: calibration and validation: 1. The flow model[END_REF]. Nevertheless, an important drawback of these methods is related to the high degree of complexity of the modeling which requires an intensive computational effort and thus limits their range of applicability to near-field scales [START_REF] Benke | Modelling conservative tracer transport in fracture networks with a hybrid approach based on the Boltzmann transport equation[END_REF] [START_REF] Painter | Upscaling discrete fracture network simulations: An alternative to continuum models[END_REF][14] [START_REF] Reeves | Incorporation of conceptual and parametric uncertainty into radionuclide flux estimates from a fractured granite rock mass[END_REF], and actually prevents the use in the probabilistic uncertainty and sensitivity analyses required in the PA of radioactive waste repositories.

Another approach proposes a stochastic continuum model which deterministically incorporates the major fracture planes and resorts to multiple statistical populations to model the statistical distribution of the hydraulic conductivity [START_REF] Hendricks Franssen | 3D inverse modelling of groundwater flow at a fractured site using a stochastic continuum model with multiple statistical populations[END_REF]; the model has been successfully applied to a threedimensional flow and transport problem in a fractured domain of the order of tens of meters, at the expense, again, of very large computational times and numerical problems due to the solution of the groundwater flow over a large domain in presence of very steep conductivity gradients.

An approach recently proposed to overcome these problems [START_REF]Rock Fractures and Fluid Flow, Contemporary Understanding and Applications[END_REF] is the "hybrid approach" [START_REF] Painter | Upscaling discrete fracture network simulations: An alternative to continuum models[END_REF]. It consists in carrying out detailed flow and transport simulations over small domains, thus limiting the numerical flow solution problems, and using the results to determine equivalent properties for continuum models to be used on larger scales [START_REF] Benke | Modelling conservative tracer transport in fracture networks with a hybrid approach based on the Boltzmann transport equation[END_REF] [START_REF] Painter | Upscaling discrete fracture network simulations: An alternative to continuum models[END_REF][15] [START_REF] Blum | Stochastic simulations of regional scale advective transport in fractured rock masses using block upscaled hydro-mechanical rock property data[END_REF]. When adopting this "upscaling" approach, the underlying assumption implicitly made is that the small domain is representative of the properties of the network and the results can be employed to upscale the conducting properties, i.e. to extend them to larger domains [START_REF] Benke | Modelling conservative tracer transport in fracture networks with a hybrid approach based on the Boltzmann transport equation[END_REF] [START_REF] Painter | Upscaling discrete fracture network simulations: An alternative to continuum models[END_REF].

In this paper, an upscaling procedure is developed to determine the equivalent parameters of a lean, one-dimensional Dual-Permeability, Monte Carlo Simulation (DPMCS) model [START_REF] Marseguerra | Modelling the Transport of Contaminants in Groundwater as a Branching Stochastic Process[END_REF] [START_REF] Ferrara | A comparison between the advection-dispersion and the Komogorov-Dmitriev model for groundwater contaminant transport[END_REF] [START_REF] Cadini | Monte Carlo simulation of contaminant migration in a dual-permeability fractured medium[END_REF] of radionuclide migration at the field scale. Monte Carlo simulation is adopted because of its flexibility, necessary to account for the complex phenomena and geometries typically encountered in realistic Performance Assessments (PAs) [START_REF] Marseguerra | Groundwater contaminant transport in presence of colloids: a stochastic nonlinear model and parameter identification[END_REF] [START_REF] Cadini | Challenges in Monte Carlo Simulation of Non-Linear Radioactive Contaminant Release Processes[END_REF]. On the other hand, the model onedimensionality is chosen in view of it utilization for fast dose predictions, within computationally demanding sensitivity and uncertainty analyses. Specified probability distributions of fracture characteristic parameters are taken as input and numerical experiments are performed to collect reference data for the parameter calibration procedure. Different realizations of fracture fields are generated on small sub-domains by random sampling from the assigned distributions and by mapping into continuum grids with constant cell size, according to a recently suggested methodology [START_REF] Reeves | Incorporation of conceptual and parametric uncertainty into radionuclide flux estimates from a fractured granite rock mass[END_REF] [START_REF] Reeves | Transport of conservative solutes in simulated fracture networks: 1. Synthetic data generation[END_REF]. Then, the MODFLOW and MT3DMS solvers [23][24] are employed to get detailed flow and transport solutions over these sub-domains [START_REF] Reeves | Incorporation of conceptual and parametric uncertainty into radionuclide flux estimates from a fractured granite rock mass[END_REF]. These are used to compute a mean breakthrough curve which serves as a reference for the identification of the equivalent parameters of the DPMCS model. The parameters identification is performed by a Genetic Algorithm search aimed at minimizing the mean squared error between the DPMCS estimated breakthrough curve and the mean one of reference. Once the DPMCS model parameters are identified, they can be implemented in a new DPMCS model upscaled to a larger domain, as demonstrated in the application concerning Pu 239 migration in a fractured medium.

The paper is organized as follows. In Section 2 the parameter identification strategy and the resulting upscaling procedure are illustrated with reference to the DPMCS model, whose principal features are recalled. The results of the application of the procedure to the PA problem under consideration are reported in Section 3. Conclusions on the capabilities of the proposed method are drawn in Section 4.

GENETIC ALGORITHM-BASED UPSCALING PROCEDURE OF A DPMCS MODEL FOR TRANSPORT IN FRACTURE NETWORKS

Fracture networks do not in general behave like a continuum [START_REF] Long | Porous Media Equivalents for Networks of Discontinuous Fractures[END_REF]: channeling of flow and other phenomena may occur, which cannot be accounted for in continuum models. Only when the density and the connectivity of the network are large enough, equivalent descriptions can be found. In this case each fracture, if connected to the others, contributes to the overall permeability of the network; however, this is not simply the sum of the permeability of each individual fracture [START_REF] Neuman | Universal Scaling of Hydraulic Conductivities and Dispersivities in Geological Media[END_REF][13] [START_REF] Hestir | Analytical Expressions for the Permeability of Random Two-Dimensional Poisson Fracture Networks Based on Regular Lattice Percolation and Equivalents Media Theories[END_REF].

Moreover, the high degree of spatial variation of the conductivity of the fractures influences the velocity field, which in turn causes an increment in the spreading of the dissolved contaminants. This phenomenon causes an apparent increase of the dispersivity [START_REF] Neuman | Universal Scaling of Hydraulic Conductivities and Dispersivities in Geological Media[END_REF], which can even be found to vary with the observed scale [8][9]. It is therefore clear that even if an equivalent description is possible, the evaluation of the equivalent parameters of the network to solve flow and transport problems is far from being trivial. In certain cases, analytic expressions for equivalent properties can be obtained by using equivalent media and percolation theories [START_REF] Hestir | Analytical Expressions for the Permeability of Random Two-Dimensional Poisson Fracture Networks Based on Regular Lattice Percolation and Equivalents Media Theories[END_REF][26] [START_REF] Sahimi | Flow and transport in porous media and fractured rock: from classical methods to modern approaches[END_REF], but the extension of these methods to realistic networks, whose parameters are described by realistic probability distributions, is burdensome. In such cases, if the network is not poorly connected, it is possible to resort to upscaling procedures which consist in carrying out detailed simulation of contaminants migration over small fractured sub-domains, assumed to be representative of the global properties of the network, and from these results determine the parameters for the equivalent continuum problem [6][7].

In this work, an upscaling procedure is developed in three main steps: i) a set of fracture networks is generated in a small sub-domain, by random sampling the characteristic parameters of each individual fracture, i.e. position, orientation, length and conductivity, from the corresponding distributions; these parameters are shown to be sufficient to determine a fracture network whose geometry is geometrically similar to that observed [START_REF] Chevalier | Stochastic delineation of wellhead protection area in fractured aquifers and parametric sensitivity study[END_REF] [START_REF] Long | From field data to fracture network modeling: an example incorporating spatial structure[END_REF]; ii) each fracture network is mapped into continuum grids with constant cell size [START_REF] Reeves | Incorporation of conceptual and parametric uncertainty into radionuclide flux estimates from a fractured granite rock mass[END_REF] and the corresponding flow and transport problems are solved by a detailed numerical model; and iii) the parameters of the DPMCS model are identified on the basis of the results of the detailed numerical model. The DPMCS model with the identified parameter values previously identified can then be applied to a larger, fractured domain.

In what follows, the details of the upscaling procedure are provided with reference to the migration of a single species of non-reactive contaminant, for simplicity.

Fracture network generation

Assuming that the individual properties are not correlated [START_REF] Reeves | Transport of conservative solutes in simulated fracture networks: 1. Synthetic data generation[END_REF], fractures are generated in two groups, of fixed orientation ( 1  , 2  ), with position uniformly distributed in a two-dimensional domain, and length ( f L ) and conductivity ( f K ) distributed according to power laws, such as: In order to limit the fracture spatial densities in the domain, the maximum density is controlled to decide whether the realization can be used or not. By choosing a small power law exponent for the fracture length distribution, the following criterion can be adopted to control maximum spatial density [START_REF] Reeves | Transport of conservative solutes in simulated fracture networks: 1. Synthetic data generation[END_REF]:

f L f a f L f f l w l L P    ) ( (1.a) f K f a f K f f k w k K P    ) ( (1.b)
5 . 2 1 1 ,    f N i i f l A (2)
where

A is the area of the domain, i f l , is the length of the th i  fracture and N is the total number of fractures in the domain.

The lower limit is not computed by evaluating the percolation threshold; rather, here, in analogy with [START_REF] Reeves | Transport of conservative solutes in simulated fracture networks: 1. Synthetic data generation[END_REF], it is instead determined indirectly by carrying out trial numerical simulations on the flow field in order to determine a value which guarantees a certain degree of connectivity and allows for a flow in the fractures from one side of the domain to the other.

According to the distributions and the constraint described above, a total number of fr M twodimensional fracture network realizations is generated by random sampling.

Flow and transport solution in a fractured sub-domain by MODFLOW and MT3DMS

The non-reactive single species flow and transport problems are solved for each fr M fractured network realization by means of the finite element solvers MODFLOW and MT3DMS [START_REF] Mcdonald | A modular three-dimensional finite-difference groundwater flow model[END_REF] [START_REF] Zheng | MT3DMS: a modular threedimensional multispecies transport model for simulation of advection, dispersion, and chemical reactions of contaminants in groundwater systems; documentation and user's guide[END_REF].

The corresponding rates of release outside the domain (i.e. the breakthrough curves) are subsequently averaged to obtain a mean breakthrough curve, whose the shape ultimately depends on the statistical properties assumed for the fracture field. An analysis is performed to assess whether the number of field realizations fr M guarantees the convergence of the mean of the breakthrough curve.

Parameter identification and upscaling

The mean breakthrough curve is used as reference for the identification of the equivalent flow and transport parameters of a continuum model. In this case, a DPMCS model is employed [START_REF] Cadini | Monte Carlo simulation of contaminant migration in a dual-permeability fractured medium[END_REF], whose main features are summarized as follows. In general the model can be one-, two-or threedimensional but here, for simplicity of illustration and without loss of generality, a one-dimensional spatial representation of the migration domain along the horizontal direction is considered.

The domain along which the transport occurs is subdivided in

x N discrete zones, x N x ..., , 2 , 1  . The th N x   1
compartment is an absorbing state representing the "environment". In each zone the amount of contaminant in the porous matrix and in the fracture system is to be determined. To this aim two categories of particles are introduced: the solutons, i.e. the particles of contaminant in the porous matrix, and the fracturons, i.e. the particles of contaminant travelling in the fracture system.

In an elementary time interval

dt , each one of these two particles can either travel to one of the neighboring zones, (for a fracturon becoming a soluton). Finally, each contaminant particle might undergo a radioactive decay or a chemical reaction with rate  . For simplicity of illustration, we assume that the daughter of the decay or of the chemical reaction is neither radioactive nor contaminant, so that it is of no further interest in the analysis.

By assuming that the process is i) linear (i.e. each particle is independent from the others), ii) markovian (i.e. the future state depends only on the present one) and by exploiting the properties of the probability generating functions (pgfs), a system of

) 1 ( 2  x N coupled Kolmogorov's equations
for the expected values of the numbers of solutons (S) and fracturons (F) at time t in each zone

1 ..., , 2 , 1   x N x can be obtained [17]: ) 0 , 1 | , ' ( ) 0 , 1 | 1 ' ( ) 0 , 1 | , 1 ' ( ) 0 , 1 | , ' ( ] [ ) 0 , 1 | , ' ( , ' , , , , ' , , x s x s f x s s x s s x s x f s s s x s t z F r x S f t x S b t x S l f b t t x S                 (3.a) ) 0 , 1 | , ' ( ) 0 , 1 | 1 ' ( ) 0 , 1 | , 1 ' ( ) 0 , 1 | , ' ( ] [ ) 0 , 1 | , ' ( , ' , , , , , , x s x f s x s f x s f x s x s f f f x s t x S r x F f t x F b t x F l f b t t x F                 (3.b) with x N x ,..., 2 , 1 ' 
and where we assume to have one soluton in zone x at time

0  t ( x s , 1 ).
In analogy to [START_REF] Ferrara | A comparison between the advection-dispersion and the Komogorov-Dmitriev model for groundwater contaminant transport[END_REF], the transition rates of the Kolmogorov-Dmitriev model may be related to the advection-dispersion model parameters. In this case, however, two sets of properties are needed, one describing the matrix and the other the fractures system to alternatively determine the solutons and the fracturons transition rates. For this reason, an analogy can be drawn with the dual permeability advection-dispersion model [START_REF] Barenblatt | Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks[END_REF][29][30] [START_REF] Gerke | Evaluation of a First-Order Water Transfer Term for Variably Saturated Dual-Porosity Flow Models[END_REF] (see Appendix A for more details).

The stochastic process of radionuclide migration across domains modeled by the system of partial differential equations (3.a) and (3.b) represents a continuous-time Markov process [START_REF] Papoulis | Probability, random variables, and stochastic processes[END_REF], whose solution is the probability of finding the m-th contaminant particle in compartment ' z at time t . To account for realistic conditions, e.g. to account for non-homogeneities in time and space [START_REF] Dentz | Time behavior of solute transport in heterogeneous media: Transition from anomalous to normal transport[END_REF] or non-linear phenomena [START_REF] Marseguerra | Groundwater contaminant transport in presence of colloids: a stochastic nonlinear model and parameter identification[END_REF][21], a solution by Monte Carlo simulation may be quite adequate, whereby the stochastic migration of a large number M of contaminant particles in the domain are simulated by repeatedly sampling their births, from the release sources and their transitions across the medium compartments, from the proper probability density functions.

Let us denote by ) (t X m the state variable representing the position at time t of the m-th

contaminant particle ( f s m , 
for a soluton or a fracturon, respectively) in the discretized domain, i.e.,

n t X m  ) ( implies that at time t the m-th particle is in cell 1 ...., , 2 , 1   x N n
. The random walk of the individual radionuclide is simulated either until

1 ) (   x m N t X
, i.e. until it exits the domain to the "environment", that is an absorbing state from which the particles are not allowed to come back to the migration domain, or until T t  , i.e. its lifetime crosses the time horizon T of the analysis. For computational purposes, the time horizon T

    M k n Count k P s s n , ,  (4.a)     M k n Count k P f f n , ,  (4.b)
The total probability that a contaminant particle occupies cell n at time k is:

        M k n Count M k n Count k n Count k P f s n , , ,    (5)
Similarly, the probability density function of the release into the environment (breakthrough curve), ) (k pdf env , can be estimated as:

  t M k N Count k pdf x env     , 1 ) ( (6)
As anticipated above, the aim is that of identifying the parameters of the DPMCS model so that it is capable of approximating the mean breakthrough curve obtained on the small sub-domain by the MODFLOW and MT3DMS suite. In particular, the DPMCS model parameters to be determined are: i) the equivalent conductivity of the fractures network, eq K , which allows determining the velocity in the equivalent-continuum fracture system, eq v , from Darcy's equation; ii) the equivalent dispersion coefficient in the fractures, eq D ; iii) and the mass exchange parameters ( ) describing the process of retardation due to solute interaction with the matrix. The values of these parameters allow computing the transition rates of the DPMCS model (see Appendix A for details). In general, the method allows considering advective transport also in the porous matrix, but, in this modeling, for simplicity and without loss of generality, this is neglected and only the retardation due to diffusion into the matrix pores is considered.

A Genetic Algorithm [START_REF] Goldberg | Genetic algorithms in search, optimization, and machine learning[END_REF] search is performed to find the equivalent DPMCS model parameters. The parameter values found are implemented in the upscaled DPMCS model, which can be applied to larger, fractured domains. An exemplary application of this procedure is described in the next Section.

UPSCALING APPLICATION ON A REALISTIC CASE STUDY

The proposed procedure is applied to a case study from literature [START_REF] Cadini | Monte Carlo estimation of radionuclide release at a repository scale[END_REF] concerning the release of Pu 239 from a near surface repository, whose design concept was studied by ENEA [START_REF] Cadini | Monte Carlo estimation of radionuclide release at a repository scale[END_REF][36] and has similarities with the currently operative disposal facility of El Cabril, Spain [START_REF] Zuloaga | Long term water scenario in low-level waste disposal vaults, with particular regard to concrete structures in el cabril, cordoba, spain[END_REF]. As the radionuclides exit the repository floor, they migrate through the unsaturated zone, whose most important effect is that of retarding the infiltration towards the saturated zone (refer to [START_REF] Cadini | Monte Carlo-based assessment of the safety performance of a radioactive waste repository[END_REF] for further details). Then, as in [START_REF] Chapman | Preliminary Feasibility Assessment for Near-Surface Engineered LLW Repositories at Two Sites[END_REF], radionuclide migration in the groundwater domain is described to occur within a highly fractured limestone layer exhibiting extensive karstic phenomena. In this work, the equivalent parameters of the fracture network are estimated by adopting the upscaling procedure described in Section 3. For the properties of the fractures, the distributions of equations (1.a) and (1.b) are assumed with the parameters reported in Table 1. The criterion of equation ( 2) is considered for the maximum spatial density. The corresponding networks are sparsely fractured domains dominated by a few long fractures connecting their boundaries [START_REF] Reeves | Transport of conservative solutes in simulated fracture networks: 1. Synthetic data generation[END_REF]. Moreover, the following correction is applied to the sampled value of the conductivity in order to account for the "stair step" pattern of the velocity field, which causes an apparent decrease of the hydraulic gradient from cell to cell in MODFLOW [START_REF] Reeves | Transport of conservative solutes in simulated fracture networks: 1. Synthetic data generation[END_REF]: in the vertical direction. These dimensions are comparable to those adopted in [START_REF] Benke | Modelling conservative tracer transport in fracture networks with a hybrid approach based on the Boltzmann transport equation[END_REF]. For simplicity's sake, in the determination of the equivalent parameters, the migration of a single species of a non-reactive tracer is considered, characterized by a retardation factor 1  R (no adsorptions). The transport is simulated over a period of h 300 and the injection is assumed to be impulsive, i.e. all the particles are released at time

|) | cos | | (sin      K K MODFLOW (7)
0  t . A total of 100  fr M
realizations of the fracture field are generated and, for each realization, the flow and transport are solved by MODFLOW and MT3DMS numerical methods. The hydrological parameters used for the simulations are summarized in Table 2. The mean breakthrough curve over the different realizations is computed, together with its 1 standard deviation bands (Figure 1, solid line). This curve is considered as a reference for the calibration of the parameters of a small scale DPMCS model, which will be later exploited for developing an upscaled DPMCS model.

Table 2: Hydrological parameters

Matrix conductivity

To promote model simplicity for repeated use in PA probabilistic analyses, a one-dimensional model is considered. The parameters to be determined are, then, the equivalent conductivity of the network eq K , the equivalent dispersion coefficient ) and dx dh is the average hydraulic gradient over the domain. As mentioned in Section 3, the transport in the matrix subsystem is neglected because of the low permeability value (see Table 2).

The range of variability of each gene for the search is determined by computing the minimum and maximum values of the corresponding physical quantities over all the generated fracture field realizations. As these values vary more than one order of magnitude, a logarithmic search is performed for all parameters apart from the reciprocal of the Peclet's number, whose variation is kept linear. The search ranges are set equal to . At the end of every simulation, the computed breakthrough curve, i.e.

) (k pdf

(see equation ( 6)) is compared with that obtained from the MODFLOW and MT3DMS simulations and the fitness parameter, i.e. the mean squared error, is computed. To enhance its variability, the fitness parameter is set in logarithmic scale.

The optimal curve found with the GA procedure, i.e. the one with the minimum mean squared error, is reported in Figure 1 (dashed line). For the small-scale DPMCS model, the optimal values of the parameters found by the GA search are: . The general behavior follows that of the mean curve obtained by MODFLOW and MT3DMS, lying well within the 1 bands. It shows good predictions for the peak's value and time, whereas it is less satisfactory for the tail of the curve. This is expected, because of the one-dimensionality of the DPMCS model. In fact, the shape of the peak predicted by the detailed MODFLOW and MT3DMS is mainly due to the fastest tracer particles which are transported by advection along the main fractures and rarely interact with the matrix nor change their fracture direction. To intuitively describe this process, we can assume that, in a sparsely fractured domain dominated by a few long fractures connecting its boundaries, as the one under consideration [START_REF] Reeves | Transport of conservative solutes in simulated fracture networks: 1. Synthetic data generation[END_REF], the advective transport is driven by a single fracture, e.g. with orientation

  21 8 1 
and average hydraulic conductivity

s m K f / 10 1 . 1 3   
; then, we can write the average time  one of these particles takes to exit the domain boundary as:

h dx dh K l q l v l f f f 24 ) tan( ) tan( ) tan( 1 1 1 1                        (9)
which is very close to the actual peak time obtained with MODFLOW and MT3DMS. Since the advection of the tracer particles along a main, straight fracture is well represented by a onedimensional process, then the DPMCS model, which is intrinsically one-dimensional, is capable of satisfactorily capturing this behavior. On the other hand, the shape of the tail of the curve is due to hydrodynamic dispersion, both at the microscale of the pores and macroscale of the fracture network, and to matrix interaction. The DPMCS model satisfactorily represents the matrix interaction phenomena, as demonstrated by the rather flat tail mainly due to those particles which are trapped in the matrix and then slowly released (see [START_REF] Cadini | Monte Carlo simulation of contaminant migration in a dual-permeability fractured medium[END_REF]); on the contrary, the differences in the tails predicted by the two models is mainly due to hydrodynamic dispersion, which is an intrinsically two-dimensional phenomenon (three-dimensional in more general applications). The larger spread of the particle release in the MODFLOW and MT3DMS due to dispersion is more evident at earlier times when the advection process is dominant, giving rise to the initial underestimation before 150 h; then, simple mass conservation considerations explain the successive overestimation. The equivalent dispersion coefficient eq D in the one-dimensional DPMCS model is capable of partially capturing this behavior, but cannot exactly reproduce However, DPMCS tests carried out with different values of eq D , not reported here for brevity' sake, show that the tail can be better reproduced at the expense of a decrease in the model peak accuracy and an increase in the total mean squared error [START_REF] Cadini | Monte Carlo simulation of contaminant migration in a dual-permeability fractured medium[END_REF] [START_REF] Cadini | Monte Carlo estimation of radionuclide release at a repository scale[END_REF]. . The Pu 239 particles are injected in the groundwater according to the release rate computed in [START_REF] Cadini | Monte Carlo estimation of radionuclide release at a repository scale[END_REF] and the migration in the domain is assumed to occur under linear isothermal conditions with retardation factor equal to 5667  R [START_REF] Chapman | Preliminary Feasibility Assessment for Near-Surface Engineered LLW Repositories at Two Sites[END_REF]. Radioactive decay of Pu ; all other hydrological parameters are left unchanged with respect to the previous case (Table 1). The spatial domain is discretized into 200 x 50 cells of m in the results of the calibration phase, since they are, again, mainly due to the one-dimensional average advective transport component. The shape of the curve is more symmetric than before, probably because the matrix interaction effects are "covered" by those due to dispersion, which is larger now due to the larger domain. The width of the curve is very similar to that of the reference one, which shows that eq D well describes the average longitudinal dispersion effects. On the other hand, the differences in the rising and descending parts of the curve are related to the asymmetry of the reference curve, which is the result of the intrinsic two-dimensional effects of the macroscopic dispersion at the full-scale domain due to the fractures. These effects cannot be captured unless a more complicated two dimensional DPMCS model is developed, but this is out of the scope of this work.

The apparent temporal upscaling is actually due to the fact that in the parameter calibration phase an instantaneous source of tracer particles with retardation factor , not shown here for brevity, does not show significant differences in the model behavior, although the time scales become comparable to that of the calibration simulations. The entire procedure has been carried out with a 6GB RAM computer with an INTEL® CORE™ i5 CPU M520 @ 2.40GHz processor. The calibration procedure requires a total of to calibrate the parameters with the Genetic Algorithm). This procedure, however, needs to be performed only once and the results can subsequently be employed to solve the upscaled migration problem with a much faster and "leaner" model, which lends itself to be used for the repeated calculations of probabilistic uncertainty and sensitivity analysis of PA. In the case of the by the MODFLOW and MT3DMS solver. Note that computational advantage increase as the domain becomes larger.

CONCLUSIONS

The assessment of the safety performance of a nuclear waste repository stands on the use of quantitative models for the prediction of the release and transport of radionuclides on extended time and spatial scales. In the transport process, fractures play a relevant role which must be properly evaluated and accounted for.

In this paper, we propose a novel approach for upscaling fracture networks properties into a continuum model. The approach consists in computing detailed solutions of the flow and transport problems in a two-dimensional sub-domain of the fractured medium with the finite difference solvers MODFLOW and MT3DMS, and using the results to determine the equivalent conducting properties to be employed in a one-dimensional Dual-Permeability, Monte Carlo Simulation (DPMCS) model. These properties are determined by means of a Genetic Algorithm search which minimizes the mean squared error between the breakthrough curve found with MODFLOW and MT3DMS and the one computed with the Monte Carlo scheme.

The approach is employed to carry out realistic simulations of Pu 239 release from a radioactive waste repository and migration within a fractured domain. By inserting the calibrated parameters in the Monte Carlo scheme, good agreement is found with the corresponding solution computed with a reference, detailed MODFLOW and MT3DMS simulation.

The proposed approach allows to drastically simplify the numerical model, by coarsening the spatial discretization and thus reducing the computational time. The procedure does require a certain effort but, once the equivalent parameters are found, they can be employed to solve transport problem on domains of any size, where the detailed solutions would become too expensive and even fail to reach convergence.

The choice of a one-dimensional model reduces the required computational times, so that the approach becomes attractive with respect to the fundamental PA tasks of uncertainty and sensitivity analyses on the physical parameters governing the radionuclides migration in groundwater. On the other hand, this choice affects the accuracy of the results, in particular due to the impossibility of fully capturing the two-dimensional (or, more generally, three-dimensional) dispersion patterns. A two-dimensional model would probably provide improved estimation accuracy. It is to be noted, however, that this would be at the cost of increased computational times, so that the model would not be effective from an uncertainty and sensitivity analyses viewpoint. Nonetheless, due to the flexibility offered by the Markovian modeling, investigations are currently being carried out in this regard.

APPENDIX A

In dual permeability and dual porosity models, the transport phenomenon is assumed to occur in two domains: the porous matrix and the fracture system. In particular, dual porosity models divide the medium in a mobile and an immobile region, assuming that flow occurs only in the larger permeability zone, whereas dual-permeability models assume the flow occurs in both subsystems.

The model describes the heterogeneous porous media as a system made up of two different coupled subsystems and therefore the properties of the bulk medium depend on those of the different porous subsystems. Thus, following the notation of [30][31], the transport of solute in this system can be described by two coupled advection-reaction-dispersion equations, where the subscript m refers to the matrix subsystem and f to the fractures:

f s f f f f f f f f f w c c c t c R                 ) ( ) ( ) ( f f q D (A1.a) m s m m m m m m m m m w c c c t c R                 ) ( ) ( ) ( m m q D (A1.b)
where D is the dispersion tensor,  is the decay or chemical reaction term, R is the retardation coefficient, w is the relative volume of the pore system (

f m w w   1
),  is the water content, q is the Darcy's velocity, c is the solute concentration and s  is the exchange term between the porous matrix and the fracture system. According to [START_REF] Gerke | A Dual-Porosity Model for Simulating the Preferential Movement of Water and Soluted in Structured Porous Media[END_REF] [START_REF] Gerke | Evaluation of a First-Order Water Transfer Term for Variably Saturated Dual-Porosity Flow Models[END_REF], assuming that only the diffusive contribution to the exchange process is taken into account and that the flow field is at steady-state, the exchange term s  can be expressed as:

) ( m f m m s s c c w      (A2)
The coefficient s  in equation (A2) is [START_REF] Gerke | A Dual-Porosity Model for Simulating the Preferential Movement of Water and Soluted in Structured Porous Media[END_REF][31]: ) and, subsequently, discretize them in space by a centered Euler method: The difference in the exchange rates depends mainly upon the different relative volume of each system: the larger the volume of the subsystem is, the less probable the transition to the other subsystem is.

f m s f f f m m s f f f f f f f f                        (A4.

  ) respectively, or undergo a transition to the other domain in the same zone x , thus transforming into the other particle with corresponding rate

  The objective guiding the search is to minimize the mean squared error between the breakthrough curve obtained by the DPMCS model for the small sub-domain and the mean curve obtained by the MODFLOW and MT3DMS solutions of the fr M fracture field realizations on the same small subdomain. The optimization procedure entails running the DPMCS model several times in correspondence of properly chosen sets of values for the inputs mean squared error with respect to the mean output of the MODFLOW and MT3DMS.

1 ; 5 , 5 ,

 155 Numerical experiments, consisting in the generation of many different realizations of a 2-D fracture field, are carried out to collect reference data to be used for determining the equivalent conducting properties for larger scales application. The size of the elementary domain, adopted to compute the upscaling parameters, is the right and left boundaries are given constant hydraulic head, whereas the upper and lower boundaries are set to noflow condition. To avoid boundary effects, the source is placed 5 cells, i.e. m from the upstream boundary and the breakthrough curve is computed 5 cells, i.e. m before the downstream boundary. Therefore, the actual distance along the horizontal direction is m 40 . A planar source is employed, whose length is set equal to m 40



  Algorithm search is performed to determine the optimal values of the parameters of the one-dimensional, DPMCS model. A constraint is set on Peclet's number to preserve the physics of the problem ( D uL Pe  where u is the velocity, L is a characteristic length and D is the dispersion coefficient): a substantial change in this parameter would, in fact, result in the passage from an essentially advective problem to a predominantly diffusive problem. Therefore, instead of optimizing both eq K and eq D , the equivalent conductivity and the reciprocal of the Peclet's number are encoded in the first two genes of the GA. The value of eq D can, then, be indirectly computed from these two parameters: is the discretization length in the horizontal direction (set equal to m x 1  

  of bits for their binary representation is 10. A population of 100 chromosomes evolving for 200 generations is employed.The GA search recursively solves the Monte Carlo scheme of Section 2in the first upstream compartment of a dual-permeability, one-

Figure 1 :

 1 Figure 1: Comparison of MODFLOW and MT3DS mean breakthrough curve (solid line) and Monte Carlo curve (dashed line) on the 50m  50m domain

  simplicity's sake, the daughter of the decay is not further considered in the transport analysis. Both upstream and downstream boundaries are given constant hydraulic head with total difference on the domain

1 by m 1 ,.

 11 and the top and bottom boundaries are set as no-flow. For the Monte Carlo simulation, a total of 6 10  M fracturon particles are injected in the first upstream compartment of a dual permeability, one-dimensional array of 19  The coarser spatial discretization employed for the Monte Carlo simulation is allowed by the equivalent parameters, which are valid at the scale of tens of meters. In 2, the breakthrough curve computed with the DPMCS model (dashed line) is compared with the reference one (solid line) obtained by averaging 10 simulations performed with MODFLOW and MT3DMS in correspondence of different fracture network realizations; the errorbars indicate the 1 ( = standard deviation) variability of the results obtained with different networks. Overall satisfactory agreement is achieved: the peak height and time are well captured, as

  1  R(no adsorption) is employed, whereas in the full-scale simulation the more realistic Pu 239 source taken from[START_REF] Cadini | Monte Carlo estimation of radionuclide release at a repository scale[END_REF] spreads approximately over a 1500 y time horizon; moreover, the slowing down the transport process. Simulations carried out on the same full-scale domain with instantaneous source and 1  R

Figure 2 :

 2 Figure 2: breakthrough curve at 200 m from the source: comparison between MODFLOW and MT3DMS (solid line) and DPMCS model (dashed line)

D,

  is the effective diffusion coefficient at the interface,  is a semi-empirical geometric factor and a is the characteristic half-width of the porous matrix. The terms in equations (4.a) and (4.b) are all referred to the relative volume of the corresponding subsystem, with the exception of s  , which is defined as the mass of solute per unit volume of bulk soil per unit time [Gerke and van Genuchten 1993a]. For the general form of equations (A1.a), (A1.b) and (A2) the interested reader should refer to [30][31].Then under the hypotheses of: i) one-dimensional domain, ii) no retardation ( of the matrix and the fractures system, respectively), constant dispersion coefficients it is possible to write equations (4.a) and (4.b) in terms of the concentrations referred to the total bulk volume (

  formally identical to (3.a) and (3.b); thus, by comparing the two systems of equations it is possible to write the following relationships for the forward and backward rates of the solutons and the fracturons: negative, equations (A5.a) and (A5.b) require an upper limit for the spatial discretization cell: forward and backward rates, the exchange rates can be expressed as:

  , resides, during its random walk, in compartment n at time k . At the end of the M simulated random walks of the contaminant particles, the values accumulated in the counters allow estimating the time-dependent probabilities of cell occupation by a soluton or a

	each compartment	n		, 1		2	,	....,	N	x		1	and each discrete time	k		, 1	2	,	....,	N	t	. During the
	simulation, a one is accumulated in the counters	Count		s	  k n ,	and	Count	f	  k n ,	if a soluton or
	fracturon, respectivelyfracturon, ) ( , t P s n and	n P	,	f	(	t	)		respectively:
																	Count	s	  k n ,	and	Count	f	  k n , , are associated to

is discretized in t N equally spaced time instants, with time step t  ; two different counters,

Table 1 : Stochastic parameters for the generation of the fracture fields

 1


G Ripamonti

G Lonati

P Baraldi

F Cadini

E Zio

Uncertainty propagation in a model for the estimation of the ground level concentration of dioxin/furans emitted from a waste gasification plant

Keywords: Uncertainty propagation, Hybrid probabilistic-possibilistic Monte Carlo, waste gasification plant, Environmental Impact Assessment

In this paper we compare two approaches for uncertainty propagation in a model for Environmental Impact Assessment (EIA). A purely Probabilistic (PMC) and a Hybrid probabilisticpossibilistic Monte Carlo (HMC) method are considered in their application for the estimation of the ground levels concentration of dioxin/furans emitted from a waste gasification plant. Under the condition of insufficient information for calibrating the estimation model parameters, HMC is shown to be a valid way for properly propagating parameters uncertainty to the model output, without adopting arbitrary and subjective assumptions on the input probability distribution functions. In this sense, HMC could improve the transparency of the EIA procedures with positive effects on the communicability and credibility of its findings.

INTRODUCTION

According to the current EU regulations (Directive 85/337/EEC as amended by Directives 97/11/EC and 2003/35/EC), an Environmental Impact Assessment (EIA) is required for public and private projects likely to have significant impacts on the environment by virtue inter alia of their nature, size or location. The scope of an EIA is the identification, description and assessment of a project's direct and indirect effects on human beings, fauna and flora, soil, water, air, climate, landscape, material assets and cultural heritage. The EIA is a part of the decisional process for the approval of a proposed project. The outcomes of the EIA help orient decision makers to the priorities, thus ensuring that environmental considerations are taken into account in an effective way. Furthermore, the EIA serves an important procedural role in the overall decision-making process by promoting public information and facilitating participation. This entails that an EIA be carried out taking into account that its findings will be handled not only by technical experts, but also communicated to decision makers and stakeholders of various nature, including the local population. Since in most common settings the proposed project does not meet popular consensus a priori, communicability of EIA results to the involved stakeholders plays a major and challenging role. To this purpose the EIA procedure must be open and transparent, reviewable and leading to robust and reproducible outcomes for informed decision-making and confident public awareness.

With reference to the proposal of new facilities responsible of atmospheric emissions (in the present work we consider a waste thermal treatment plant as an example), the EIA procedure consists of three main steps (Figure 1): a) the characterization of the source: estimation of the atmospheric emissions resulting from the operation of the plant, in terms of both the nature and the quantities of pollutants; b) the estimation of the atmospheric pollutant concentration levels in the area close to the proposed new installation, by means of proper atmospheric dispersion models; c) the assessment of the pollutant concentration levels at the receptors, resulting from the superposition of the estimated concentration to the existing background concentration levels.

The first objective of the EIA is then the demonstration that the resulting concentration levels at the receptors comply with both long-term (i.e.: annual average) and short-term (i.e.: hourly and daily averages) air quality limits for regulated pollutants. Secondly, as those facilities can be also responsible for the emission of both carcinogenic and non-carcinogenic trace pollutants, as metals or dioxins and furans (PCDD/Fs) in the case of waste thermal treatment plants, additional evaluations may also include human health risk assessment. This latter assessment is characterized by two main features: i) a considerable complexity of the mechanisms involved and hence large related uncertainties, ii) a general attitude of "aversion to risk", with implications on the way uncertainties are treated [START_REF] Guyonnet | Soft methods for treating uncertainties: applications in the field of environmental risks[END_REF].

The uncertainties that affect all the three steps of the EIA procedure may be originated from randomness of the physical processes involved or from lack of precise knowledge of them. In the end, regardless of the approaches and models adopted to estimate the plant emissions and their dispersion in the atmosphere, the final ambient concentrations estimates can be seriously affected by the uncertainty of the values of the source term parameters and by the variability of meteorological parameters.

On the contrary, the above mentioned attitude of "aversion to risk" often leads to answering to the question whether a given "acceptable" threshold may be exceeded, based on limited scenario calculations. In common practice, this translates in ignoring input parameter uncertainties to adopt a direct, deterministic approach in which conservative assumptions on the scenarios and parameter values are made, e.g.: i) full-load plant operation and ii) concentrations of the emitted pollutants equal to the regulatory limits. This approach leads to often excessively conservative results, not providing a realistic and objective description of the plant impact on the environment [START_REF] Lonati | Probabilistic health risk assessment of carcinogenic emissions from a MSW gasification plant[END_REF] To overcome this problem, the EIA procedure can be framed within a probabilistic approach for i) describing the uncertainties of the input parameters and ii) propagating them to the models outputs, i.e.: the atmospheric pollutant concentrations. In the probabilistic approach, the uncertain variables are described by means of probability distribution functions (PDFs), which are then propagated through the impact assessment model by Monte-Carlo (MC) simulation according to a pure Probabilistic Monte Carlo (PMC) method [START_REF] Schuhmacher | The use of Monte-Carlo simulation techniques for risk assessment: study of a municipal waste incinerator[END_REF][START_REF] Sonnemann | Framework for the uncertainty assessment in The Impact Pathway Analysis with an application on a local scale in Spain[END_REF][START_REF] Lonati | Health risk analysis of PCDD/F emissions from MSW incineration: comparison of probabilistic and deterministic approaches[END_REF]). Nevertheless, even in the PMC approach uncertainty and variability in the atmospheric dispersion model are usually neglected.

In general, the construction by statistical analysis of the empirical PDFs representing the uncertainty in the input parameters may not be an easy task due to the scarce availability of data [START_REF] Ferson | Different methods are needed to propagate ignorance and variability[END_REF]Ginzburg 1996, Baudrit et al. 2006). Within a probabilistic framework, the analyst traditionally solves this task forcing a subjective PDF on the scarce available data, within a Bayesan (subjective) view of probabilities. This presents two major drawbacks:

-doubts may be raised on the consistency and even conservativeness of the forced PDFs, in light of the incomplete/imprecise information at hand. The risk is to introduce artificial information that is in fact not there; -the subjectivity in the analysis is often neither declared nor justified, causing a loss of transparency in the procedure.

These drawbacks are of particular relevance, considering the EIA role as decision aiding and informative tool. A lack of transparency in the procedure can lead to a loss of credibility by the decision makers and public.

To cope with this situation of limited information, while avoiding to introduce subjective assumptions, a number of alternative representations of uncertainty have been proposed (e.g.: fuzzy set theory, evidence theory, possibility theory and interval analysis) and applied in different fields of risk assessment [START_REF] Klir | Fuzzy Sets and Fuzzy Logic: Theory and Applications[END_REF][START_REF] Aven | Some considerations on the treatment of uncertainties in risk assessment for practical decision making[END_REF]. Among them, possibility theory has received growing attention because of its representation power and its relative mathematical simplicity. In recent past years some applications of these alternative approaches have been exploited also in the environmental risk assessment field (Guyonnet et al. 2003[START_REF] Kentel | 2D Monte Carlo versus 2D fuzzy Monte Carlo health risk assessment[END_REF][START_REF] Li | An integrated fuzzy-stochastic modelling approach for risk assessment of groundwater contamination[END_REF][START_REF] Kumar | Partitioning total variance in risk assessment: Application to a municipal solid waste incinerator[END_REF].

In practice, different input parameters of the same model may be characterized by different amounts of information availability, thus possibly giving rise to both probabilistic and possibilistic uncertainty representations. In this regard, Hybrid probabilistic-possibilistic Monte Carlo methods (HMC) have been recently proposed for combining both representation frameworks in the same uncertainty propagation analysis (Baudrit et al. 2006[START_REF] Baraldi | A Combined Monte Carlo and Possibilistic Approach to Uncertainty Propagation in Event Tree Analysis[END_REF][START_REF] Flage | Possibility-probability transformation in comparing different approaches to the treatment of epistemic uncertainties in a fault tree analysis[END_REF].

In this work, the uncertainty representation and propagation problems are considered with reference to the first two steps of the EIA procedure (i.e.: emission estimation and atmospheric dispersion).

For ease of illustration, we make reference to a case study concerning the project of a new waste gasification plant. The plant, designed for a daily waste throughput of 900 metric tons, has been recently proposed for realization in the South-East of the city of Milan, Italy and is subject to the requirement of EIA for its potential impact on air quality.

The focus of the analysis is on dioxins and furans (PCDD/Fs) since the thermal treatment of waste is a well-known source of PCDD/Fs to the atmosphere [START_REF] Domingo | Health risk assessment of PCDD/PCDF exposure for the population living in the vicinity of a municipal waste incinerator[END_REF][START_REF] Schuhmacher | Long-term study of environmental levels of dioxins and furans in the vicinity of a municipal solid waste incinerator[END_REF]) and these pollutants have a relevant role in health risk assessment due to their carcinogenic potential and to their persistency in the environment, once released into the atmosphere [START_REF] Travis | Human exposure to dioxin[END_REF]Hattemer-Frey 1991, Mukerjee 1998).

The final objective of the work is the assessment of the contribution to PCDD/Fs long-term concentration (i.e.: annual average) in the ambient air in the surroundings of the proposed plant due to the flue gas emissions released at the stack of the plant. This information, other than a direct information on the environmental impact of plant emissions, is also the basic input for the subsequent health risk assessment step.

The results obtained through a HMC method for the assessment of PCDD/Fs ground-level concentrations are compared to those obtained by a standard PMC approach, and both the transparency of the EIA procedure and the communicability of results are discussed.

The manuscript is organized as follows: in Section 2 the equations used to model plant emissions and ambient air concentrations are presented; the characterization of the input parameters uncertainty for the development of HMC and PMC methods is illustrated in Section 3; Section 4 shows the output of the two propagation methods and discusses their methodological implications; finally, the advantages of the HMC method are drawn in the Conclusions. For completeness, the Appendix Section A1 provides the basics of possibility theory and Section A2 presents the pure PMC method for uncertainty propagation and its evolution to the Hybrid probabilistic-possibilistic Monte Carlo method under conditions of scarce data availability.

POLLUTANT AIR CONCENTRATIONS MODEL

The contribution to the ambient air concentration C air (fg m -3 ) on an annual average basis can be computed by multiplying an emission term Q (ng s -1 ), which quantifies the mass flow rate of PCDD/Fs released into the atmosphere, and an atmospheric dispersion factor DF (fg m -3 /ng s -1 ), which quantifies the ambient concentration per unit of mass flow rate:

DF Q C air   (1)
Since the atmospheric dilution of the emission occurs at different extents depending on the distance from the plant stack, C air and DF are actually space-dependent and take different values, C air (x,y) DF(x,y), at the receptors points located in the study domain, usually described by a Cartesian grid centered on the plant stack.

Equation (1) follows a simplified approach for the assessment of C air since it assumes that the PCDD/Fs are emitted in the gas-phase only: this approach, although still providing realistic results, reduces the complexity of the calculation model thus enabling the application of the uncertainty propagation methods.

The emitted PCDD/Fs mass flow rate Q in equation ( 1) has been computed on the basis of the plant daily throughput of waste, P (Mg waste day -1 ), the PCDD/Fs concentration in the emitted flue gas, C D (ng m -3 ), and the specific gas production V F (m 3 Mg waste -1 ):

24 3600     F D V C P Q (2)
As usual for this kind of evaluations (Van den [START_REF] Van Den Berg | Toxic Equivalency Factors (TEFs) for PCBs, PCDDs, PCDFs for humans and wildlife[END_REF], US EPA 2005), the PCDD/Fs concentration is expressed in terms of equivalent toxicity mass per unit gas volume at normal conditions (0° C, 101.3 kPa); coherently, the specific gas production refers to the same temperature and pressure conditions.

The dispersion term DF in equation ( 1 The DF values to be used in equation ( 1) for the estimation of the PCDD/Fs annual average concentration have been obtained for the 1681 nodes of a Cartesian grid centered on the plant (250 m cell spacing) through one-year long model simulations based on the hourly time series of locally measured meteorological data.

UNCERTAINTY REPRESENTATION AND PROPAGATION

The input parameters of the emission model (equation ( 2)) have been characterized as follows:

 the plant daily throughput P is a constant parameter whose value is assigned during project design; for the plant of this case study, P has been set to the value of 900 Mg waste day -1 ;  the PCDD/Fs concentration C D is an uncertain parameter whose value varies during normal gasification operation due to the fluctuations of the process parameters and to the heterogeneous and variable composition of the fed waste;

 the specific gas production V F is an uncertain parameter. Although its value is usually set in the plant design phase, actually V F can present variations during plant operation caused by fluctuations of the energy content of the fed waste as a consequence of its heterogeneous and variable composition.

Due to the rather limited applications of the waste gasification process, extended time series datasets containing C D and V F values collected during operation of similar plants are not available.

However, different studies have recently investigated the pollutant emission from waste gasification plants [START_REF] Klein | Gasification: An Alternative Process for Energy Recovery and Disposal of Municipal Solid Wastes[END_REF][START_REF] Yamada | Thermoselect Waste Gasification and Reforming Process[END_REF][START_REF] Porteous | Why energy from waste incineration is an essential component of environmentally responsible waste management[END_REF], Arena 2012, University of California 2009).

From these studies, 35 C D and four V F values have been derived.

The availability of a statistically significant set of C D average values, has allowed to represent the uncertainty on this model parameter by means of a probability distribution. A Kolmogorov-Smirnov's test considering lognormal, Weibull, Beta, logistic PDFs has been performed in order to properly choose the PDF best representative of the literature data. Figure 2 shows the selected Beta-PDF of parameters α = 0.36 and β = 1.32 and maximum value 0.07 ng m -3 , that is 30% less than the regulatory limit of 0.1 ng m -3 (European Union, 2000).

Conversely, for V F the available information (four literature values and the case study design value) is very scarce. Thus, two different alternative representations of the uncertainty affecting V F have been considered:  a subjective triangular PDF based on analyst judgment. The PDF range has been chosen to be [3360,6670] (m 3 Mg waste -1 ), which corresponds to the minimum and maximum values of the four literature values, and its mode has been chosen to be 5420 m 3 Mg waste -1 , which corresponds to the current case study design value (Figure 3);  a triangular possibility distribution with the same range of the subjective PDF and most likely value set equal to the case study design value (Figure 3). An example of the different wind conditions observed in the study area is given by the annual wind roses reported in Figure 4 for four years out of the ten considered.

For illustration purposes, only one grid point has been considered for the uncertainty propagation, chosen as the receptor most impacted by the plant emission (x M ,y M ), i.e. the grid node characterized by the highest average of the 10 yearly DF values.

Due to the small dataset available, two alternative representations of the DF(x M ,y M ) uncertainty have been considered:

 a trapezoidal PDF based on analyst judgement with its minimum equal to 0, its core ranging between 5.69•10 -3 and 4.19•10 -2 fg m -3 /ng s -1 , which corresponds to the minimum and maximum values of the ten estimated DF values, and its maximum equal to 0.21 fg m -3 /ng s -1 , which is the DF obtained for the receptor point assuming the worst day conditions observed during the 10year period for atmospheric dispersion all year long (Figure 5);

 a trapezoidal possibility distribution with the same features of the subjective PDF (Figure 5).

RESULTS

With respect to the uncertainty propagation, the PMC method has been applied to the case in which As expected, the cumulative distribution of the model output obtained by the PMC is within the belief and plausibility functions obtained by the hybrid approach.

Furthermore, the representation of the uncertainty on C air provided by the PMC method appears to be more concise and easy to be interpreted than the one provided by the HMC method. This is a consequence of the different representation of the uncertainty on the input parameters VF and DF(x M ,y M ): while the HMC method considers a set of cumulative distributions, the PMC forces all the uncertainty on V F and DF(x M ,y M ) to be represented by single distributions. However, given the scarcity of information available for the two parameters in this case study, the use of a specific PDF seems unjustified. Within the EIA procedure, for communication purposes, it is more effective to lump all the information contained in the obtained cumulative distribution of C air into a single value, such as for example the β percentile. To this aim, setting a degree of confidence β = 0.95, the PMC method provides a single value, C air (x M ,y M ) 95th = 2.80•10 -1 fg•m -3 , representative of the 95 th percentile of the distribution of the PCDD/Fs annual average at the most impacted receptor point in the area. On the contrary, with the HMC the 95 th percentile is an uncertain quantity itself, whose true value lies in the interval [1.15•10 -2 , 5.65•10 -1 ] fg•m -3 , i.e. the 95 th percentile of the Belief and Plausibility functions. Table 2 reports the values of different percentiles for both PMC and HMC methods: PMC values are always in the range of HMC values and the distance between the plausibility and belief functions increases from lower to higher percentiles. Notice that for any percentile considered the gap between the plausibility and the belief functions is in the range of one/two orders of magnitude:

the uncertainty on the VF and DF parameters appears to significantly affect the estimated ground level concentrations. In this respect, the capability of the HMC method to separately process the contribution of the different uncertainties gives rise to results that clearly show the effects of the lack of knowledge on the input parameters. This is a desirable quality leading to more informative and transparent outputs, also in the light of the subsequent calculations for health risk assessment.

It is worth noticing that the concentration estimates at the selected receptor for the most and the least favourable year for atmospheric dispersion obtained by means of the traditional deterministic approach are 0.03 fg•m -3 and 0.24 fg•m -3 , respectively, thus varying by almost one order of magnitude. This means that the traditional deterministic approach based on just one single year simulation, as is done in the common practice, can lead to either too precautionary or non-conservative estimates based on the arbitrary choice of the analyst about the reference year to be used in the model simulation. In this respect, a preliminary assessment of the representativeness of the year chosen for the weather model simulation or multi-year modelling are recommended for deterministic approach calculations.

Finally, notice that both PMC and HMC methods have been applied under the simplifying assumption of independence of the parameters C D and V F , although they are expected to be somehow correlated. Future developments of this work may thus include investigating the effects of the dependence between these parameters on the uncertainty propagation results, for example resorting to the approach suggested in [START_REF] Pedroni | Empirical comparison of methods for the hierarchical propagation of hybrid uncertainty in risk assessment, in presence of dependences[END_REF].

and  

) ( sup A P A   (a.2)
where inf and sup are with respect to all probability measures in P. Hence the necessity measure is interpreted as a lower level for the probability and the possibility measure is interpreted as an upper limit. Referring to subjective probabilities, the bounds reflect that the analyst is not able or willing to precisely assign his/her probability, and the bounds are the best he/she can do given the information available; in other words, he or she can only describe a subset of P which contains his/her probability [START_REF] Dubois | 337/EEC on the assessment of the effects of certain public and private projects[END_REF].

A.2. PROBABILISTIC-POSSIBILISTIC UNCERTAINTY REPRESENTATION AND PROPAGATION THROUGH A MODEL

Let us consider a model whose output

  n Y Y Y f Z ,..., , 2 1  is a function of n uncertain variables n i Y i , , 1 ,  
. The uncertainties of the first k variables can be represented by probability distributions

  y p i Y
properly derived from data, whereas for the remaining k n  the scarcity of information provide a weaker statistical base for a specific probability assignment, making the construction of PDFs a critical operation. Two different approaches are here presented to represent and propagate such uncertainties a. Pure probabilistic Monte-Carlo method (PMC)

The uncertainties of all n input variables are treated in terms of probability distributions. Probability representations of the k n  variables can be assigned i.e.: based on expert judgment. For example, let us considered a typical case in which the available information on a model parameter x is only that its values are located somewhere between a value x min and a value x max. In this case, a uni-

form probability distribution     min max / 1 x x x p   ,   x [x min ,
x max ] is typically assumed to represent the uncertainty on x. This approach appeals to: i) Laplace principle of insufficient reason according to which all that is equally plausible is equally probable and to ii) the maximum entropy approach [START_REF] Baraldi | Methods of Uncertainty Analysis in prognostics[END_REF] . However, doubts on the consistency of this uncertainty representation have been raised (Baudrit et al. 2006) as it seems that the insufficient knowledge may justify choices of specific functional probability distributions, like the uniform, but it should somehow account for the full set of possible probability distributions on X = |x min , x max |, so that the probability of value   x is allowed to take any values in |0,1| [START_REF] Baraldi | Methods of Uncertainty Analysis in prognostics[END_REF].

The PMC method then applies a single loop Monte Carlo simulation to propagate the uncertainties of the n variables through the model. The operative steps of the propagation procedure are the following: 1 Sample the

n k Y Y   , , 1   .
Different methods have been developed to derive possibility distributions from the available information (Baudrit & Dubois 2006, Dubois 1993). Often the analyst only knows that an uncertain variable can take values in a given range [a, b] wherein the most likely value is c. Triangular possibility function with the range |a, b| taken as base and c taken as vertex can be typically used to describe this information. It has been shown that the family of probability distributions defined by such possibility distribution with contains all the probability distributions with support [a, b] and mode c (Baudrit & Dubois 2006). In the HMC method, the propagation of the uncertainties expressed in terms of both PDFs and possibility curves is performed by combining the Monte Carlo technique with fuzzy interval analysis (Baudrit et al. 2006). The operative steps of the propagation procedure are the following: 1 sample the (Baudrit et al. 2006). Thus, one way to estimate the total uncertainty on   Y f is to provide a confidence interval at a given level of confidence, taking the lower and upper bounds from

th i  realization   i k i y y , , 1  of the random variable vector   k Y Y , , 1  2 select
    z Y f Bel z F , ) ( ) (     ,     z Y f Pl z F , ) ( ) (    
    z Y f Pl , ) (    and   ) (Y f Bel   z ,  
, respectively (Baudrit et al. 2006). Notice, however, that it is not possible to directly obtain information on the mean and standard deviation of the output distribution from the plausibility and belief distributions.

  ) has been computed by simulating the atmospheric transport and dispersion of the emitted PCDD/Fs in the surroundings of the plant by means of the ISCST3 model (Industrial Source Complex -Short-term version 3). Formerly a US EPA regulatory model for air dispersion modeling, this model is still regarded as a screening model for the assessment of the impact on air quality of atmospheric emissions from industrial sources (US EPA, 1995). Local meteorological conditions (wind speed and direction, ambient air temperature, atmospheric stability), source features (stack height, flue gas speed and temperature, pollutant mass flow rate) and geographical features of the area (terrain elevation, land use) are the basic input data for the dispersion model. Output of the ISCST3 model is the 2D field of the contribution to the groundlevel concentration data in the study domain around the plant. Since in this model the output con-centrations are linearly proportional to the emitted flow rate, the dispersion term DF can be evaluated running the model for a unit mass flow rate (1 ng s -1 ).

  Notice that, although the possibility and probability distributions used to represent the uncertainty on V F have similar shapes, they convey very different information: the possibility distribution summarizes a set of cumulative distributions, bounded by the so-called Necessity and Possibility functions (see Appendix a.1), whereas the subjective PDF defines just one specific cumulative distribution of that set.On the other hand, model predictions of the DF values in equation (1) are also uncertain quantities affected by i) the natural variability in the input parameters (such as meteorological variables and source features), ii) measurement errors and iii) modeling errors due to the difficulties in capturing the atmospheric behaviors[START_REF] Sax | A case study for assessing uncertainty in local-scale regulatory air quality modeling applications[END_REF] Isakov 2003, Rao 2005). Neglecting for simplicity the model uncertainties, it turns out that the DF values are mainly affected by the variability in the meteorological input parameters, which largely exceeds the effect of the variability of flue gas temperature and speed. In fact, the temperature is typically kept at a rather constant value during plant operation, whereas the flue gas speed plays a negligible role on the atmospheric behavior of the emitted plume. Here, the DF uncertainty has been estimated indirectly, by taking into account 10 different years of local meteorological data and by separately running the dispersion model for each input data set: thus, 10 yearly DF values have been estimated for each grid node in the computation domain.

  the uncertainty on C D , V F and DF(x M ,y M ) are represented by the PDFs, whereas the HMC method considers the PDF for C D and the possibility distributions for V F and DF(x M ,y M ). The number m of MC realizations of the parameters C D , V F and DF(x M ,y M ) used in the PMC method has been set to 1000. The HMC procedure, instead, has been run considering m = 1000 realizations of the parameter C D and, for each of these realizations, 21 α-cut values (range 0-1, step 0.05) for the possibilistic variables V F and DF(x M ,y M ).The cumulative distributions of the model output variable, C air (x M ,y M ) obtained by applying the PMC and HMC uncertainty propagation methods are presented in Figure6. The PMC method provides a single cumulative distribution (continuous line), whereas the HMC model provides the belief (lower dotted curve) and the plausibility (upper dashed curve) of the set A=(-,C air ), which can be interpreted as the boundary cumulative distributions of C air at the receptor of coordinates (x M , y M ).

Figure 6

 6 Figure 6 also shows that the HMC method allows to process separately the uncertainty on C D , represented by a probability distribution, from the uncertainty on V F and DF, represented by possibility distributions. These two contributions are explicitly visible in the results: in fact, the uncertainty on C D affects the slope of the belief and plausibility functions, while the distance between the belief and plausibility functions reflects the lack of knowledge of the parameters V F and DF. On the other hand, in the PMC method the contributions of the input parameters uncertainties on the model output are merged, giving rise to the slope of the output cumulative distribution.

3,

  Repeat steps 1 and 2 for many times. Then, derive the cumulative distribution   on the output uncertainty. Setting a percentile  , that corresponds to the desired degree of confidence, a concise output value of Z can be provided. The output realizations i found in step 2 can be directly used to compute other lumped indicators of the output distribution, such as its mean and standard deviation values.b. Hybrid probabilistic-possibilistic Monte-Carlo method (HMC)In the Hybrid probabilistic-possibilistic method, the uncertainties of the k n  variables are represented in terms of possibility distributions  

  a possibility value α [0,1] and the corresponding α -cuts of the possibility distributions e: the intervals of values of the possibilistic variables   step 2. and repeat for another α -cut; the fuzzy random realization (fuzzy interval) each α -cut 5 return to step 1. to generate a new realization of the random variables. The procedure is repeated for m i , , 1   : at the end of the procedure m realizations of fuzzy intervals are obtained, i.e.: m different realizations of possibility and necessity measures can be combined to obtain the A , this technique computes the belief and plausibility as the average of the possibility measures associated with each output fuzzy interval. The likelihood of the value   Y f passing a given threshold z can then be computed by considering the belief and the plausibility of the set as bounding, average cumulative distributions

CONCLUSION

In this work the uncertainties in the first steps for the atmospheric emission and dispersion estimation of an Environmental Impact Assessment (EIA) procedure have been treated by both a Hybrid probabilistic and possibilistic Monte Carlo method (HMC) and a pure Probabilistic Monte-Carlo method (PMC). Calculations have been performed for the case study of the project of a new waste gasification plant, focusing on dioxin/furan emissions that are a primary concern for human health due to their carcinogenic potential. The expected impact of the plant emissions has been quantified in terms of the contribution to the annual average concentration of dioxins and furans in ambient air at the most impacted receptor in the neighbourhood of the plant. The comparison between the results obtained by the HMC and PMC methods points out that the former is more effective in propagating the input uncertainties through the model when little information is available for some input parameters. In particular, i) the information provided by the output of the HMC method is more consistent with that (little) available for the input parameters and ii) uncertainty is processed more "transparently" than with PMC methods, avoiding arbitrary and subjective assumptions by the analyst on the input probability distribution functions and separating the contributions to the output uncertainty due to the probabilistic and possibilistic input parameters.

This research has been conceived as a preliminary study to assess the potential for the implementation of HMC methods to the complete EIA procedure. In this regard, the satisfactory outcomes of this first analysis foster future works towards the extension of the HMC method to the remaining steps of an EIA procedure, i.e. the evaluation of the human exposure through the different impact pathways and the subsequent risk assessment.

Finally, it is worth noticing that the principal users of the EIA findings are decision-makers often responsible for communicating the results to the population living in the area of the planned installation. With respect to environmental and health-related issues, there is in general social aversion to accept information expressed in terms of probability. Therefore, future studies will also have to investigate the way of post-processing the results of a whole EIA procedure developed within an hybrid probabilistic-possibilistic framework to make them communicable to the stakeholders and to make them more easily understandable. To this purpose, the abilities of the HMC method to avoid arbitrary assumptions and to provide results that explicitly report both probabilistic and possibilistic uncertainties could be useful for a transparent and clear post-process. for some y , it means that the outcome y is possible, i.e.: is just unsurprising, normal, usual [START_REF] Dubois | 337/EEC on the assessment of the effects of certain public and private projects[END_REF]. This is a much weaker statement than when probability is 1. The possibility function gives rise to probability bounds, upper and lower probabilities, referred to as necessity and possibility measures ( N ,  ). The possibility of an event

and the necessity measure   A N is defined by:

Let P(  ) be a family of probability distributions such that for all events