- Email: [email protected]
Monte Carlo estimation of radionuclide release at a repository scale
Annals of Nuclear Energy 37 (2010) 861–866
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Monte Carlo estimation of radionuclide release at a repository scale F. Cadini a,*, J. De Sanctis a, T. Girotti a, E. Zio a, A. Luce b, A. Taglioni b a b
Dipartimento di Energia, Politecnico di Milano, Via Ponzio 34/3, I-20133 Milan, Italy ENEA CR Saluggia, Via Crescentino 41, 13040 Saluggia (VC), Italy
a r t i c l e
i n f o
Article history: Received 9 September 2009 Received in revised form 18 February 2010 Accepted 19 February 2010 Available online 3 April 2010 Keywords: Radioactive waste repository Radionuclide migration Compartment model Monte Carlo simulation
a b s t r a c t Prediction of radionuclide release is a central issue in the performance assessment of nuclear waste repositories. This requires modeling of the radionuclide migration processes through the repository barriers, accounting for the related uncertainties. The present paper illustrates a Monte Carlo simulationbased compartment model in which detailed, local-scale modeling feeds a global-scale analysis of the repository, at reasonable computational expenses. An application to a realistic case study is presented to verify the feasibility of the approach. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction The objective of the engineered and natural barriers of a radioactive waste repository is to prevent the release of radionuclides and retard their migration to the groundwater, and eventually to the biosphere (IAEA, 1999a). The assessment of the performance of these barriers, both in the repository design and operative phases, is a complex task which entails (Helton and Sallaberry, 2009; Yim and Simonson, 2000): (i) the identification of the scenarios that influence the repository behavior throughout its lifetime; (ii) the estimation of their probabilities of occurrence; (iii) the estimation of their consequences associated to the release of the radionuclides, typically in terms of the expected dose received by the defined critical group of people; and (iv) the evaluation of the uncertainties associated to the aforementioned estimates. In this framework, the quantitative analysis of the processes of radionuclide migration across the barriers of the repository to the main intake paths, plays a fundamental role. Eventually, the analysis is to be made at the repository scale, accounting also for macroscopic effects such as those due to the spatial arrangements of the waste containers in the repository structures (Kawasaki et al., 2005; Kawasaki and Ahn, 2006) and the heterogeneity of the media (Tsujimoto and Ahn, 2005; Williams, 1992; Williams, 1993). For example, the assumption of radionuclides point release, often made for ease of computation, may lead to inaccuracies in the estimation of the repository safety performance, since there might be * Corresponding author. E-mail address: [email protected] (F. Cadini). 0306-4549/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2010.02.011
non-linear interactions among the radionuclides released by different canisters in the same groundwater stream and the radionuclide travel time across the barriers of the repository may vary by orders of magnitude. In practice, this can lead to significant computational expenditures if the analysis is performed by extension to the repository scale of the detailed flow and mass transport codes used at the single waste package scale. Also, various national Nuclear Regulatory Agencies are issuing new licensing standards for radioactive waste repositories, which require the estimation of the expected dose to some defined critical groups and the associated uncertainty of both aleatory (stochasticity of the future system behavior) and epistemic (lack of knowledge of the model parameter values) types (Helton and Sallaberry, 2009). This also can lead to quite substantial computational efforts with the flow and transport codes classically used. An approach which seems suitable for effectively handling the above issues is that of proceeding to a compartmentalized modeling of the radionuclide migration at the different scales of the media of the repository. Detailed modeling would be applied to the smaller domain to characterize the migration process in a single compartment; the results of this modeling would serve as input of the model of the repository area compartments, based on a coarser discretization and, thus, computationally less intensive. Simplifying Markovian hypotheses are often introduced to model the stochasticity of the transfer process across the compartments (Kawasaki et al., 2005; Kawasaki and Ahn, 2006; Lee and Lee, 1995; Xu et al., 2007). In this paper, instead, we resort to Monte Carlo simulation for modeling particle migration in a compartmentalized repository; the simulation scheme allows modeling realistic features of system behavior without encountering
862
F. Cadini et al. / Annals of Nuclear Energy 37 (2010) 861–866
the difficulties (or simplifications thereof) brought by the need of finding analytical or numerical solutions, and naturally accounting for the aleatory uncertainty related to the stochasticity of the migration process. The proposed approach is illustrated through its application to a near surface repository design concept studied by ENEA, the Italian New technologies, Energy and Environment Agency (Marseguerra et al., 2001). The results are compared to those obtained by application of the detailed model to the full repository scale. The paper is organized as follows. In Section 2, the formulation of the process of radionuclide migration within a compartmentalized medium and the basic principles of Monte Carlo simulation are recalled. In Section 3, the Monte Carlo simulation model of radionuclide migration is applied to the ENEA case study (Kawasaki et al., 2005; Marseguerra et al., 2001). Finally, some conclusions on the advantages and limitations of the proposed approach are provided in Section 4.
2. Compartment model of radionuclide migration For simplicity, but with no loss of generality, let us assume that the radionuclide migration occurs within a two-dimensional medium discretized into N = Nx Ny compartments (Fig. 1), where Nx and Ny are the number of compartments along the x and y axes, respectively. An additional compartment N + 1 is added to represent the environment in which the domain is embedded. For simplicity of illustration, and with no loss in the characterization of the migration process, radioactive decay is here neglected; its inclusion within the Monte Carlo simulation modeling schemes poses no particular additional challenge. Let us also suppose that the radionuclide migration across the compartments can be modeled by a continuous-time Markov process. The state variable X(t) represents the position state of a single radionuclide migrating within the compartment matrix, i.e., X(t) = n implies that the radionuclide is in compartment n at time t, n = 1, 2, . . . , N + 1 (IAEA, 1999b). The stochastic process of radionuclide migration described by the state variable X(t) is completely described by the row vector
PðtÞ ¼ ½P1 ðtÞP2 ðtÞ PNþ1 ðtÞ
ð1Þ
where Pn(t) is the probability that the radionuclide is in state n at time t, i.e., X(t) = n. Accounting for the Markov property, if the radionuclide is in state (compartment) i at time t (i.e., X(t) = i), the probability of reaching state j at time t + v does not depend on the states (compartments) X(u) visited by the radionuclide at times u prior to t (i.e., 0 6 u < t); in other words, given the present state (compart-
n=1
2
...
Ny
2
Ny +1
Ny + 2
...
2N y
3
...
...
...
...
...
...
...
...
...
...
...
...
nx = 1
... ...
N = Nx × N y
Nx
ny = 1
2
...
Ny
Fig. 1. An Nx Ny matrix of N compartments.
ment) X(t) occupied by the radionuclide, its future behavior is independent of the past (Zio, 2009):
P½Xðt þ v Þ ¼ jjXðtÞ ¼ i; XðuÞ ¼ xðuÞ; 0 6 u < t ¼ P½Xðt þ v Þ ¼ jjXðtÞ ¼ i
ð2Þ
The conditional probabilities
P½Xðt þ v Þ ¼ jjXðtÞ ¼ i i; j ¼ 1; 2; . . . ; N þ 1
ð3Þ
are called the transition probabilities of the Markov process. If the transition probabilities do not depend on the time instants t and t + v but only on the width of the separating time interval v, then the Markov process is said to be homogeneous:
P½Xðt þ v Þ ¼ jjXðtÞ ¼ i ¼ pij ðv Þ t; v > 0 and i; j ¼ 1; 2; . . . ; N þ 1
ð4Þ
Considering a time interval v = Dt sufficiently small that only one transition can occur and applying the Taylor expansion of Eq. (4), the one-step transition probability from compartment i to compartment j can be written as:
pij ðDtÞ ¼ P½Xðt þ DtÞ ¼ jjXðtÞ ¼ i ¼ kij Dt þ hðDtÞ
ð5Þ
where kij is the rate of the transition from state (compartment) i to state (compartment) j, and lim
Dt!0
hðDtÞ Dt
¼ 0.
In the continuous time domain (i.e., for an infinitesimal time interval Dt ? 0), the stochastic time Tij that the radionuclide resides in state (compartment) i before making a transition to state (compartment) j is, then, exponentially distributed with parameter kij and the migration process is probabilistically described by the following system of ordinary differential equations:
dP ¼ PðtÞ K dt
ð6Þ
where
2
Nþ1 P k1j 6 6 j¼2 6 K¼6 6 k21 4
k12
Nþ1 P j¼1; j–2
k2j
3 k1Nþ1 7 7 7 7 k2Nþ1 7 5
ð7Þ
is the transition rate matrix. In principle, Eq. (6) allows finding analytical solutions for the dynamics of the state probabilities P(t), provided the values of the compartment transition rates kij (7) are available. However, in many realistic cases, some of the above assumptions must be relaxed, e.g. to account for non-homogeneities in time and space. In such cases, the transition times between the compartments can no longer be described by exponential distributions and analytic solutions are often not available, thus rendering numerical approximation schemes mandatory. In this regard, Monte Carlo simulation can offer an effective procedure for estimating the probabilities Pn(t) that a radionuclide is in compartment n at time t, n = 1, 2, . . . , N + 1 and, consequently, the probability density function of the release in the environment pdfenv(t). The stochastic migration process of a large number M of radionuclides in the Nx Ny domain is simulated by repeatedly sampling the transitions of each individual radionuclide particle across all compartments, from the proper transition probability density functions. The random walk of the individual radionuclide is simulated either until it exits the domain to the N + 1 compartment ‘‘environment”, or until its lifetime crosses the time horizon T of the analysis. The time horizon T is subdivided in Nt time instants at a distance Dt to each other; a counter C(n,k) is associated to each compartment n = 1, 2, . . . , N + 1 and each discrete time instant k = 1, 2, . . . , Nt. During the simulation, a one is accumulated in
F. Cadini et al. / Annals of Nuclear Energy 37 (2010) 861–866
the counter C(n,k) every time a radionuclide particle resides, during its random walk, in compartment n at a time k. At the end of the M simulated random walks of the radionuclide particles, the values accumulated in the counters allow estimating the time-dependent probabilities of compartment occupation Pn(t) as:
Pn ðkÞ ffi
Cðn; kÞ ; M
ð8Þ
Analogously, the release probability density function can be estimated as:
pdfenv ðtÞ ffi
CðN þ 1; kÞ ; M Dt
kDt < t < ðk þ 1ÞDt
ð9Þ
In this framework, the transitions between the compartments are described by stochastic laws which should be derived by proper physical considerations: in the simple case of a Markovian migration process, the constant transition rates can be identified
863
by comparison to the deterministic parameters of classical advection dispersion models (Marseguerra and Zio, 2002; Xu et al., 2007); when such description does not yield sufficiently accurate results, then the compartment model described above must be fed by more realistic distributions estimated by means of detailed models applied at local scale (Kawasaki and Ahn, 2006), as it will be shown in the next Section.
3. Estimation of the release from a near surface repository In what follows, the Monte Carlo simulation-based compartment model illustrated in the previous Section is applied to estimate the release of radionuclides from a near surface repository of a design concept studied by ENEA (Marseguerra et al., 2001). The main containment structures of the disposal facility are the waste packages, the modules or containers, the cells and the
Fig. 2. Conceptual design of the waste package (left) and the module (right) (Marseguerra et al., 2001).
Fig. 3. The 5 6 8 array of modules in a cell (Marseguerra et al., 2001).
864
F. Cadini et al. / Annals of Nuclear Energy 37 (2010) 861–866
disposal units; these constitute a set of multiple barriers against water flow and radionuclide migration. The waste packages (Fig. 2, left) consist of steel drums containing the radioactive waste immobilized in a concrete solidified matrix. The diameter of the waste package is 0.791 m and its height is 1.1 m, for a total volumetric capacity of around 400 l. The module is a concrete boxshaped structure, covered and sealed with a concrete top cover, which contains six waste packages (Fig. 2, right); the empty spaces between the packages are filled by a special concrete-based filling, i.e. bentonite. The module external dimensions are: 3.05 m of
Fig. 4. Conceptual scheme of the proposed compartment model, S = source; CT = crossing time.
Fig. 5.
239
Pu release probability density function from a module (left) and
length, 2.09 m of width and 1.7 m of height, whereas the internal dimensions are: 2.75 m of length, 1.79 m of width and 1.37 m of height. The modules are arranged in 5 6 8 arrays within cells of concrete structure built below the natural ground level, but above the normal water table level (Fig. 3). Finally, the disposal unit is a concrete structure embedding a row of 6–10 cells, whereas the whole disposal facility is in general made up of several units, typically in parallel rows, each one representing an independent system which can be built and operated individually without interfering with the other units. Under reasonable simplifying hypotheses of (i) identical modules, (ii) a prevalent vertical mass transport, and (iii) a symmetry in the lateral diffusion spread, the estimation of the probability of release into the unsaturated zone below the repository can be reduced to a one-dimensional problem of estimating the release from a column of five identical modules (Fig. 4). Then, the proposed compartment representation of the repository column consists in a one-dimensional array in which a compartment corresponds to a module. Each compartment is provided with a radionuclide source term, whose time distribution is equal to that of the upstream module release. Physically in fact, the radionuclide released by the upstream module can be assumed as a line source for the downstream module. For those source particles, it seems appropriate to assume a distribution of the transition times to the next compartment (crossing times, CT) equal to the distribution of the times spent by the radionuclides in crossing the module, i.e., in reaching the bottom of the module. The estimation of the release out of the repository column is performed by the Monte Carlo simulation of M = 5500 5 = 27,500 particles of the radionuclide 239Pu, equally distributed among the compartments of the column and released according to the source time distribution. This particular radionuclide has been chosen since it represents a long term threat to the environment for its radioactivity (T1/2 = 2.411 104 year) and toxicity. Obviously, it is expected that its activity concentration in a near surface repository is kept below according to the IAEA safety requirements (IAEA, 1994, 1999b). The transition times of the particles between the compartments are sampled from the crossing time distribution. The mission time T = 1000 year is divided into Nt = 1000 time channels of width Dt = 1 year. For simplicity, and with no loss of generality, it is conservatively assumed that (i) the protective action of the concrete cell roof and ceiling, and the backfill layers can be neglected; (ii) the whole column is saturated and a constant water head of
239
Pu crossing time probability density function of a module (right).
865
F. Cadini et al. / Annals of Nuclear Energy 37 (2010) 861–866 -3
Pu239 release probability density [1/years]
2.5
x 10
Compartment model Detailed model
2
1.5
1
0.5
0
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Time [years] Fig. 6. Comparison between the 239Pu release probability density function from a five-modules column obtained by the compartment model and that obtained by the detailed MASCOT model of (Marseguerra et al., 2003).
0.15 m is applied at the top of the higher module, so that the water head at the top of the column h(z = 5 1.7 m) = 5 1.7 + 0.15 = 8.65 m; (iii) the water head at the bottom of the column is h(z = 0) = 0; (iv) each module bears a constant water head gradient DDhz ¼ 1:0176, where Dh = 8.65 m and Dz = 5 1.7 m = 8.5 m is the column height; (v) the (239Pu radioactive decay and its generation from the decay chains of other radionuclides are neglected, and (vi) the migration of 239Pu occurs under the hypothesis of linear isothermal (Marseguerra et al., 2001). Both the source term and the transition times are estimated by means of the detailed model Monte Carlo Analysis of Subsurface Contaminant Transport (MASCOT), developed by the LASAR group of the Department of Energy of the Politecnico di Milano (http:// lasar.cesnef.polimi.it) (Marseguerra et al., 2003). The approach is based on the Kolmogorov–Dimitriev theory of branching stochastic processes (Kolmogorov and Dimitriev, 1947). The corresponding model is evaluated by Monte Carlo simulation, where a large number of particles of 239Pu are followed in their migration through the module media according to the appropriate probability distribution functions describing their transport, based on the physical– chemical properties of the radionuclide 239Pu and the hosting media (Marseguerra et al., 2001). The flexibility of the proposed model allows taking into account complex geometries and heterogeneities in the physical–chemical properties of the media crossed. The resulting source term and transition times distribution estimates are reported in Fig. 5, left and right respectively. The total release probability from the column of five modules is shown in Fig. 6 (crosses). The five peaks represent the releases of the five modules delayed by the downstream compartment crossing. The first peak is the unaltered contribution of the last module of the column, i.e., its release. The macroscopic effects of the array of modules are (i) progressive delays of the releases in the unsaturated zone below the repository and (ii) progressive spread of the releases due to the diffusion process. The comparison in Fig. 6 between the releases estimated by the compartment model (crosses) and those estimated by the application of the detailed model to the whole column (solid line) shows a satisfactory agreement, with a total CPU time of 514 s for the compartment model (including the release and crossing time distribution estimation, 186 s and 325 s, respectively), which is about four times shorter than that offered by the MASCOT detailed model, i.e. 2143 s.
4. Conclusions The objective of the performance assessment of a radioactive waste repository is that of verifying the compliance of the expected doses to the critical groups to the limits imposed by law by the various National Regulatory Agencies, taking into account the related uncertainties. To this aim, a fundamental task is the modeling of the radionuclide migration through the repository barriers, typically due to water infiltration and consequent transport of the radionuclides out of the repository to the main intake paths. In this paper, a Monte Carlo simulation scheme has been employed to model the processes of radionuclide migration to the environment within a compartmentalized representation of the physical domain in which the transport occurs. The application of the method to a realistic repository case study has shown the power offered by the flexibility intrinsic in the stochastic simulation, which allows to handle, for example, the heterogeneity of the migration domain and the complex time-dependent distributions of the radionuclide sources and transition times. Furthermore, the compartment representation of the domain in which the radionuclide migration occurs allows the extension of the results of the detailed modeling of the physical transport processes occurring at the compartment scale to the global scale of the whole repository, at lower computational times. In the application presented, the individual radionuclide sources have been considered independent, i.e. there is no interference of the sources with the release of radionuclides. This is intuitively suitable for radionuclides that are released congruently or instantaneously from the waste matrix, or in general, for those cases where adjacent sources do not affect releases of radionuclides in the compartment of interest. If there is interference of sources/plumes with the release of radionuclides or the release is dependent on the concentration field in the medium, the base Monte Carlo approach is not suitable. On the other hand, the flexibility of the simulation approach allows further developments of the method towards the modelization of non-linear phenomena in the release processes (e.g. concerning leaching, diffusion, dissolution and solubility limits) and transport phenomena (e.g. concerning fractures, colloids and bio-transport), and the inclusion of uncertainty propagation analyses. These developments will be the subject of future research investigations.
866
F. Cadini et al. / Annals of Nuclear Energy 37 (2010) 861–866
Acknowledgements This work has been funded by the Ente per le Nuove Tecnologie, L’Energia e l’Ambiente (ENEA) within the Framework Program Agreement with the Italian Ministry of Economic Development for the research area N. 5.2.5.8 ‘‘NUOVO NUCLEARE DA FISSIONE”. References Helton, J.C., Sallaberry, J.C., 2009. Conceptual basis for the definition and calculation of expected dose in performance assessments for the proposed high-level radioactive waste repository at Yucca Mountain, Nevada. Reliability Engineering and System Safety 94, 677–698. IAEA, 1994. Classification of Radioactive Waste. IAEA Safety Series, 111-G-1.1, Vienna. IAEA, 1999a. Safety Assessment for Near Surface Disposal of Radioactive Waste. Safety Standards Series No. WS-G-1.1, IAEA eds., Austria. IAEA, 1999b. Near Surface Disposal of Radioactive Waste. IAEA Safety Standards Series – Requirements, WS-R-1, Vienna. Kawasaki, D., Ahn, J., 2006. Compartment Model for Particle Migration with Residence Time Distribution. IHLRWM 2006, April 30–May 4, Las Vegas, NV. Kawasaki, D., et al., 2005. Markov Chain Model Particle Migration at the Repository Scale. Global 2005, October 9–13, Tsukuba, Japan. Kolmogorov, A.N., Dmitriev, N.A., 1947. Stochastic branching processes. Doklady Akademi Nauk SSSR 56 (1), 7–10 (in Russian).
Lee, Y.M., Lee, K.J., 1995. Nuclide transport of decay chain in the fractured rock medium: a model using continuous time markov process. Annals of Nuclear Energy 22 (2), 71–84. Marseguerra, M., Zio, E., 2002. Basics of the Monte Carlo Method with Application to System Reliability. LiLoLe-Verlag GmbH, Hagen, Germany. Marseguerra, M. et al., 2003. Monte Carlo simulation of contaminant release from a radioactive waste deposit. Mathematics and Computers in Simulation 62, 421– 430. Marseguerra, M., et al., 2001. Sviluppo di un Modello Stocastico e sua Implementazione in un Codice Monte Carlo. Internal Report, Politecnico di Milano (in Italian). Tsujimoto, K., Ahn, J., 2005. Large-scale Simulation of Rock Fracture Network for Mass Transport in Geologic Repository by Earth Simulator. Global 2005, October 9–13, Tsukuba, Japan. Williams, M.M.R., 1992. A new model for describing transport of radionuclides through fractured rock. Annals of Nuclear Energy 19 (10–12), 791–824. Williams, M.M.R., 1993. Radionuclide transport in fractured rock a new model: application and discussion. Annals of Nuclear Energy 20 (4), 279–297. Xu, S. et al., 2007. Criteria for resolution-scales and parameterization of compartmental models of hydrological and ecological mass flows. Journal of Hydrology 335, 364–373. Yim, Man-Sung, Simonson, S.A., 2000. Performance assessment models for low level radioactive waste disposal facilities: a review. Progress in Nuclear Energy 36 (1), 1–38. Zio, E. (Ed.), 2009. Computational Methods for Reliability and Risk Analysis, Series on Quality, Reliability and Engineering Statistics, vol. 14. World Scientific Publishing Co., Pte. Ltd., Singapore.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Monte Carlo estimation of radionuclide release at a repository scale F. Cadini a,*, J. De Sanctis a, T. Girotti a, E. Zio a, A. Luce b, A. Taglioni b a b
Dipartimento di Energia, Politecnico di Milano, Via Ponzio 34/3, I-20133 Milan, Italy ENEA CR Saluggia, Via Crescentino 41, 13040 Saluggia (VC), Italy
a r t i c l e
i n f o
Article history: Received 9 September 2009 Received in revised form 18 February 2010 Accepted 19 February 2010 Available online 3 April 2010 Keywords: Radioactive waste repository Radionuclide migration Compartment model Monte Carlo simulation
a b s t r a c t Prediction of radionuclide release is a central issue in the performance assessment of nuclear waste repositories. This requires modeling of the radionuclide migration processes through the repository barriers, accounting for the related uncertainties. The present paper illustrates a Monte Carlo simulationbased compartment model in which detailed, local-scale modeling feeds a global-scale analysis of the repository, at reasonable computational expenses. An application to a realistic case study is presented to verify the feasibility of the approach. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction The objective of the engineered and natural barriers of a radioactive waste repository is to prevent the release of radionuclides and retard their migration to the groundwater, and eventually to the biosphere (IAEA, 1999a). The assessment of the performance of these barriers, both in the repository design and operative phases, is a complex task which entails (Helton and Sallaberry, 2009; Yim and Simonson, 2000): (i) the identification of the scenarios that influence the repository behavior throughout its lifetime; (ii) the estimation of their probabilities of occurrence; (iii) the estimation of their consequences associated to the release of the radionuclides, typically in terms of the expected dose received by the defined critical group of people; and (iv) the evaluation of the uncertainties associated to the aforementioned estimates. In this framework, the quantitative analysis of the processes of radionuclide migration across the barriers of the repository to the main intake paths, plays a fundamental role. Eventually, the analysis is to be made at the repository scale, accounting also for macroscopic effects such as those due to the spatial arrangements of the waste containers in the repository structures (Kawasaki et al., 2005; Kawasaki and Ahn, 2006) and the heterogeneity of the media (Tsujimoto and Ahn, 2005; Williams, 1992; Williams, 1993). For example, the assumption of radionuclides point release, often made for ease of computation, may lead to inaccuracies in the estimation of the repository safety performance, since there might be * Corresponding author. E-mail address: [email protected] (F. Cadini). 0306-4549/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2010.02.011
non-linear interactions among the radionuclides released by different canisters in the same groundwater stream and the radionuclide travel time across the barriers of the repository may vary by orders of magnitude. In practice, this can lead to significant computational expenditures if the analysis is performed by extension to the repository scale of the detailed flow and mass transport codes used at the single waste package scale. Also, various national Nuclear Regulatory Agencies are issuing new licensing standards for radioactive waste repositories, which require the estimation of the expected dose to some defined critical groups and the associated uncertainty of both aleatory (stochasticity of the future system behavior) and epistemic (lack of knowledge of the model parameter values) types (Helton and Sallaberry, 2009). This also can lead to quite substantial computational efforts with the flow and transport codes classically used. An approach which seems suitable for effectively handling the above issues is that of proceeding to a compartmentalized modeling of the radionuclide migration at the different scales of the media of the repository. Detailed modeling would be applied to the smaller domain to characterize the migration process in a single compartment; the results of this modeling would serve as input of the model of the repository area compartments, based on a coarser discretization and, thus, computationally less intensive. Simplifying Markovian hypotheses are often introduced to model the stochasticity of the transfer process across the compartments (Kawasaki et al., 2005; Kawasaki and Ahn, 2006; Lee and Lee, 1995; Xu et al., 2007). In this paper, instead, we resort to Monte Carlo simulation for modeling particle migration in a compartmentalized repository; the simulation scheme allows modeling realistic features of system behavior without encountering
862
F. Cadini et al. / Annals of Nuclear Energy 37 (2010) 861–866
the difficulties (or simplifications thereof) brought by the need of finding analytical or numerical solutions, and naturally accounting for the aleatory uncertainty related to the stochasticity of the migration process. The proposed approach is illustrated through its application to a near surface repository design concept studied by ENEA, the Italian New technologies, Energy and Environment Agency (Marseguerra et al., 2001). The results are compared to those obtained by application of the detailed model to the full repository scale. The paper is organized as follows. In Section 2, the formulation of the process of radionuclide migration within a compartmentalized medium and the basic principles of Monte Carlo simulation are recalled. In Section 3, the Monte Carlo simulation model of radionuclide migration is applied to the ENEA case study (Kawasaki et al., 2005; Marseguerra et al., 2001). Finally, some conclusions on the advantages and limitations of the proposed approach are provided in Section 4.
2. Compartment model of radionuclide migration For simplicity, but with no loss of generality, let us assume that the radionuclide migration occurs within a two-dimensional medium discretized into N = Nx Ny compartments (Fig. 1), where Nx and Ny are the number of compartments along the x and y axes, respectively. An additional compartment N + 1 is added to represent the environment in which the domain is embedded. For simplicity of illustration, and with no loss in the characterization of the migration process, radioactive decay is here neglected; its inclusion within the Monte Carlo simulation modeling schemes poses no particular additional challenge. Let us also suppose that the radionuclide migration across the compartments can be modeled by a continuous-time Markov process. The state variable X(t) represents the position state of a single radionuclide migrating within the compartment matrix, i.e., X(t) = n implies that the radionuclide is in compartment n at time t, n = 1, 2, . . . , N + 1 (IAEA, 1999b). The stochastic process of radionuclide migration described by the state variable X(t) is completely described by the row vector
PðtÞ ¼ ½P1 ðtÞP2 ðtÞ PNþ1 ðtÞ
ð1Þ
where Pn(t) is the probability that the radionuclide is in state n at time t, i.e., X(t) = n. Accounting for the Markov property, if the radionuclide is in state (compartment) i at time t (i.e., X(t) = i), the probability of reaching state j at time t + v does not depend on the states (compartments) X(u) visited by the radionuclide at times u prior to t (i.e., 0 6 u < t); in other words, given the present state (compart-
n=1
2
...
Ny
2
Ny +1
Ny + 2
...
2N y
3
...
...
...
...
...
...
...
...
...
...
...
...
nx = 1
... ...
N = Nx × N y
Nx
ny = 1
2
...
Ny
Fig. 1. An Nx Ny matrix of N compartments.
ment) X(t) occupied by the radionuclide, its future behavior is independent of the past (Zio, 2009):
P½Xðt þ v Þ ¼ jjXðtÞ ¼ i; XðuÞ ¼ xðuÞ; 0 6 u < t ¼ P½Xðt þ v Þ ¼ jjXðtÞ ¼ i
ð2Þ
The conditional probabilities
P½Xðt þ v Þ ¼ jjXðtÞ ¼ i i; j ¼ 1; 2; . . . ; N þ 1
ð3Þ
are called the transition probabilities of the Markov process. If the transition probabilities do not depend on the time instants t and t + v but only on the width of the separating time interval v, then the Markov process is said to be homogeneous:
P½Xðt þ v Þ ¼ jjXðtÞ ¼ i ¼ pij ðv Þ t; v > 0 and i; j ¼ 1; 2; . . . ; N þ 1
ð4Þ
Considering a time interval v = Dt sufficiently small that only one transition can occur and applying the Taylor expansion of Eq. (4), the one-step transition probability from compartment i to compartment j can be written as:
pij ðDtÞ ¼ P½Xðt þ DtÞ ¼ jjXðtÞ ¼ i ¼ kij Dt þ hðDtÞ
ð5Þ
where kij is the rate of the transition from state (compartment) i to state (compartment) j, and lim
Dt!0
hðDtÞ Dt
¼ 0.
In the continuous time domain (i.e., for an infinitesimal time interval Dt ? 0), the stochastic time Tij that the radionuclide resides in state (compartment) i before making a transition to state (compartment) j is, then, exponentially distributed with parameter kij and the migration process is probabilistically described by the following system of ordinary differential equations:
dP ¼ PðtÞ K dt
ð6Þ
where
2
Nþ1 P k1j 6 6 j¼2 6 K¼6 6 k21 4
k12
Nþ1 P j¼1; j–2
k2j
3 k1Nþ1 7 7 7 7 k2Nþ1 7 5
ð7Þ
is the transition rate matrix. In principle, Eq. (6) allows finding analytical solutions for the dynamics of the state probabilities P(t), provided the values of the compartment transition rates kij (7) are available. However, in many realistic cases, some of the above assumptions must be relaxed, e.g. to account for non-homogeneities in time and space. In such cases, the transition times between the compartments can no longer be described by exponential distributions and analytic solutions are often not available, thus rendering numerical approximation schemes mandatory. In this regard, Monte Carlo simulation can offer an effective procedure for estimating the probabilities Pn(t) that a radionuclide is in compartment n at time t, n = 1, 2, . . . , N + 1 and, consequently, the probability density function of the release in the environment pdfenv(t). The stochastic migration process of a large number M of radionuclides in the Nx Ny domain is simulated by repeatedly sampling the transitions of each individual radionuclide particle across all compartments, from the proper transition probability density functions. The random walk of the individual radionuclide is simulated either until it exits the domain to the N + 1 compartment ‘‘environment”, or until its lifetime crosses the time horizon T of the analysis. The time horizon T is subdivided in Nt time instants at a distance Dt to each other; a counter C(n,k) is associated to each compartment n = 1, 2, . . . , N + 1 and each discrete time instant k = 1, 2, . . . , Nt. During the simulation, a one is accumulated in
F. Cadini et al. / Annals of Nuclear Energy 37 (2010) 861–866
the counter C(n,k) every time a radionuclide particle resides, during its random walk, in compartment n at a time k. At the end of the M simulated random walks of the radionuclide particles, the values accumulated in the counters allow estimating the time-dependent probabilities of compartment occupation Pn(t) as:
Pn ðkÞ ffi
Cðn; kÞ ; M
ð8Þ
Analogously, the release probability density function can be estimated as:
pdfenv ðtÞ ffi
CðN þ 1; kÞ ; M Dt
kDt < t < ðk þ 1ÞDt
ð9Þ
In this framework, the transitions between the compartments are described by stochastic laws which should be derived by proper physical considerations: in the simple case of a Markovian migration process, the constant transition rates can be identified
863
by comparison to the deterministic parameters of classical advection dispersion models (Marseguerra and Zio, 2002; Xu et al., 2007); when such description does not yield sufficiently accurate results, then the compartment model described above must be fed by more realistic distributions estimated by means of detailed models applied at local scale (Kawasaki and Ahn, 2006), as it will be shown in the next Section.
3. Estimation of the release from a near surface repository In what follows, the Monte Carlo simulation-based compartment model illustrated in the previous Section is applied to estimate the release of radionuclides from a near surface repository of a design concept studied by ENEA (Marseguerra et al., 2001). The main containment structures of the disposal facility are the waste packages, the modules or containers, the cells and the
Fig. 2. Conceptual design of the waste package (left) and the module (right) (Marseguerra et al., 2001).
Fig. 3. The 5 6 8 array of modules in a cell (Marseguerra et al., 2001).
864
F. Cadini et al. / Annals of Nuclear Energy 37 (2010) 861–866
disposal units; these constitute a set of multiple barriers against water flow and radionuclide migration. The waste packages (Fig. 2, left) consist of steel drums containing the radioactive waste immobilized in a concrete solidified matrix. The diameter of the waste package is 0.791 m and its height is 1.1 m, for a total volumetric capacity of around 400 l. The module is a concrete boxshaped structure, covered and sealed with a concrete top cover, which contains six waste packages (Fig. 2, right); the empty spaces between the packages are filled by a special concrete-based filling, i.e. bentonite. The module external dimensions are: 3.05 m of
Fig. 4. Conceptual scheme of the proposed compartment model, S = source; CT = crossing time.
Fig. 5.
239
Pu release probability density function from a module (left) and
length, 2.09 m of width and 1.7 m of height, whereas the internal dimensions are: 2.75 m of length, 1.79 m of width and 1.37 m of height. The modules are arranged in 5 6 8 arrays within cells of concrete structure built below the natural ground level, but above the normal water table level (Fig. 3). Finally, the disposal unit is a concrete structure embedding a row of 6–10 cells, whereas the whole disposal facility is in general made up of several units, typically in parallel rows, each one representing an independent system which can be built and operated individually without interfering with the other units. Under reasonable simplifying hypotheses of (i) identical modules, (ii) a prevalent vertical mass transport, and (iii) a symmetry in the lateral diffusion spread, the estimation of the probability of release into the unsaturated zone below the repository can be reduced to a one-dimensional problem of estimating the release from a column of five identical modules (Fig. 4). Then, the proposed compartment representation of the repository column consists in a one-dimensional array in which a compartment corresponds to a module. Each compartment is provided with a radionuclide source term, whose time distribution is equal to that of the upstream module release. Physically in fact, the radionuclide released by the upstream module can be assumed as a line source for the downstream module. For those source particles, it seems appropriate to assume a distribution of the transition times to the next compartment (crossing times, CT) equal to the distribution of the times spent by the radionuclides in crossing the module, i.e., in reaching the bottom of the module. The estimation of the release out of the repository column is performed by the Monte Carlo simulation of M = 5500 5 = 27,500 particles of the radionuclide 239Pu, equally distributed among the compartments of the column and released according to the source time distribution. This particular radionuclide has been chosen since it represents a long term threat to the environment for its radioactivity (T1/2 = 2.411 104 year) and toxicity. Obviously, it is expected that its activity concentration in a near surface repository is kept below according to the IAEA safety requirements (IAEA, 1994, 1999b). The transition times of the particles between the compartments are sampled from the crossing time distribution. The mission time T = 1000 year is divided into Nt = 1000 time channels of width Dt = 1 year. For simplicity, and with no loss of generality, it is conservatively assumed that (i) the protective action of the concrete cell roof and ceiling, and the backfill layers can be neglected; (ii) the whole column is saturated and a constant water head of
239
Pu crossing time probability density function of a module (right).
865
F. Cadini et al. / Annals of Nuclear Energy 37 (2010) 861–866 -3
Pu239 release probability density [1/years]
2.5
x 10
Compartment model Detailed model
2
1.5
1
0.5
0
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Time [years] Fig. 6. Comparison between the 239Pu release probability density function from a five-modules column obtained by the compartment model and that obtained by the detailed MASCOT model of (Marseguerra et al., 2003).
0.15 m is applied at the top of the higher module, so that the water head at the top of the column h(z = 5 1.7 m) = 5 1.7 + 0.15 = 8.65 m; (iii) the water head at the bottom of the column is h(z = 0) = 0; (iv) each module bears a constant water head gradient DDhz ¼ 1:0176, where Dh = 8.65 m and Dz = 5 1.7 m = 8.5 m is the column height; (v) the (239Pu radioactive decay and its generation from the decay chains of other radionuclides are neglected, and (vi) the migration of 239Pu occurs under the hypothesis of linear isothermal (Marseguerra et al., 2001). Both the source term and the transition times are estimated by means of the detailed model Monte Carlo Analysis of Subsurface Contaminant Transport (MASCOT), developed by the LASAR group of the Department of Energy of the Politecnico di Milano (http:// lasar.cesnef.polimi.it) (Marseguerra et al., 2003). The approach is based on the Kolmogorov–Dimitriev theory of branching stochastic processes (Kolmogorov and Dimitriev, 1947). The corresponding model is evaluated by Monte Carlo simulation, where a large number of particles of 239Pu are followed in their migration through the module media according to the appropriate probability distribution functions describing their transport, based on the physical– chemical properties of the radionuclide 239Pu and the hosting media (Marseguerra et al., 2001). The flexibility of the proposed model allows taking into account complex geometries and heterogeneities in the physical–chemical properties of the media crossed. The resulting source term and transition times distribution estimates are reported in Fig. 5, left and right respectively. The total release probability from the column of five modules is shown in Fig. 6 (crosses). The five peaks represent the releases of the five modules delayed by the downstream compartment crossing. The first peak is the unaltered contribution of the last module of the column, i.e., its release. The macroscopic effects of the array of modules are (i) progressive delays of the releases in the unsaturated zone below the repository and (ii) progressive spread of the releases due to the diffusion process. The comparison in Fig. 6 between the releases estimated by the compartment model (crosses) and those estimated by the application of the detailed model to the whole column (solid line) shows a satisfactory agreement, with a total CPU time of 514 s for the compartment model (including the release and crossing time distribution estimation, 186 s and 325 s, respectively), which is about four times shorter than that offered by the MASCOT detailed model, i.e. 2143 s.
4. Conclusions The objective of the performance assessment of a radioactive waste repository is that of verifying the compliance of the expected doses to the critical groups to the limits imposed by law by the various National Regulatory Agencies, taking into account the related uncertainties. To this aim, a fundamental task is the modeling of the radionuclide migration through the repository barriers, typically due to water infiltration and consequent transport of the radionuclides out of the repository to the main intake paths. In this paper, a Monte Carlo simulation scheme has been employed to model the processes of radionuclide migration to the environment within a compartmentalized representation of the physical domain in which the transport occurs. The application of the method to a realistic repository case study has shown the power offered by the flexibility intrinsic in the stochastic simulation, which allows to handle, for example, the heterogeneity of the migration domain and the complex time-dependent distributions of the radionuclide sources and transition times. Furthermore, the compartment representation of the domain in which the radionuclide migration occurs allows the extension of the results of the detailed modeling of the physical transport processes occurring at the compartment scale to the global scale of the whole repository, at lower computational times. In the application presented, the individual radionuclide sources have been considered independent, i.e. there is no interference of the sources with the release of radionuclides. This is intuitively suitable for radionuclides that are released congruently or instantaneously from the waste matrix, or in general, for those cases where adjacent sources do not affect releases of radionuclides in the compartment of interest. If there is interference of sources/plumes with the release of radionuclides or the release is dependent on the concentration field in the medium, the base Monte Carlo approach is not suitable. On the other hand, the flexibility of the simulation approach allows further developments of the method towards the modelization of non-linear phenomena in the release processes (e.g. concerning leaching, diffusion, dissolution and solubility limits) and transport phenomena (e.g. concerning fractures, colloids and bio-transport), and the inclusion of uncertainty propagation analyses. These developments will be the subject of future research investigations.
866
F. Cadini et al. / Annals of Nuclear Energy 37 (2010) 861–866
Acknowledgements This work has been funded by the Ente per le Nuove Tecnologie, L’Energia e l’Ambiente (ENEA) within the Framework Program Agreement with the Italian Ministry of Economic Development for the research area N. 5.2.5.8 ‘‘NUOVO NUCLEARE DA FISSIONE”. References Helton, J.C., Sallaberry, J.C., 2009. Conceptual basis for the definition and calculation of expected dose in performance assessments for the proposed high-level radioactive waste repository at Yucca Mountain, Nevada. Reliability Engineering and System Safety 94, 677–698. IAEA, 1994. Classification of Radioactive Waste. IAEA Safety Series, 111-G-1.1, Vienna. IAEA, 1999a. Safety Assessment for Near Surface Disposal of Radioactive Waste. Safety Standards Series No. WS-G-1.1, IAEA eds., Austria. IAEA, 1999b. Near Surface Disposal of Radioactive Waste. IAEA Safety Standards Series – Requirements, WS-R-1, Vienna. Kawasaki, D., Ahn, J., 2006. Compartment Model for Particle Migration with Residence Time Distribution. IHLRWM 2006, April 30–May 4, Las Vegas, NV. Kawasaki, D., et al., 2005. Markov Chain Model Particle Migration at the Repository Scale. Global 2005, October 9–13, Tsukuba, Japan. Kolmogorov, A.N., Dmitriev, N.A., 1947. Stochastic branching processes. Doklady Akademi Nauk SSSR 56 (1), 7–10 (in Russian).
Lee, Y.M., Lee, K.J., 1995. Nuclide transport of decay chain in the fractured rock medium: a model using continuous time markov process. Annals of Nuclear Energy 22 (2), 71–84. Marseguerra, M., Zio, E., 2002. Basics of the Monte Carlo Method with Application to System Reliability. LiLoLe-Verlag GmbH, Hagen, Germany. Marseguerra, M. et al., 2003. Monte Carlo simulation of contaminant release from a radioactive waste deposit. Mathematics and Computers in Simulation 62, 421– 430. Marseguerra, M., et al., 2001. Sviluppo di un Modello Stocastico e sua Implementazione in un Codice Monte Carlo. Internal Report, Politecnico di Milano (in Italian). Tsujimoto, K., Ahn, J., 2005. Large-scale Simulation of Rock Fracture Network for Mass Transport in Geologic Repository by Earth Simulator. Global 2005, October 9–13, Tsukuba, Japan. Williams, M.M.R., 1992. A new model for describing transport of radionuclides through fractured rock. Annals of Nuclear Energy 19 (10–12), 791–824. Williams, M.M.R., 1993. Radionuclide transport in fractured rock a new model: application and discussion. Annals of Nuclear Energy 20 (4), 279–297. Xu, S. et al., 2007. Criteria for resolution-scales and parameterization of compartmental models of hydrological and ecological mass flows. Journal of Hydrology 335, 364–373. Yim, Man-Sung, Simonson, S.A., 2000. Performance assessment models for low level radioactive waste disposal facilities: a review. Progress in Nuclear Energy 36 (1), 1–38. Zio, E. (Ed.), 2009. Computational Methods for Reliability and Risk Analysis, Series on Quality, Reliability and Engineering Statistics, vol. 14. World Scientific Publishing Co., Pte. Ltd., Singapore.