- Email: [email protected]
Contaminant transport in bidimensional porous media via biased Monte Carlo simulation
Pergamon
NucL Energy, Vol. 25, No. 16, pp. 1301-1316, 1998 © 1998 Elsevier Science Ltd. All fights reserved PII: S0306-4549(98)00015-2 Printed in Great Britain 03064549/98 $19.00 + 0.00 Ann.
C O N T A M I N A N T T R A N S P O R T IN B I D I M E N S I O N A L P O R O U S M E D I A VIA B I A S E D M O N T E C A R L O SIMULATION M. MARSEGUERRA* and E. ZIO Department of Nuclear Engineering, Polytechnic of Milan, Via Ponzio 34/3, 20133 Milano, Italy (Received 23 December 1997)
A b s t r a c t ~ n e of the most vulnerable spots of nuclear power systems is the hazard posed by the radioactive wastes at the end of the fuel cycle. Currently, the most widely pursued end solution consists of disposing the radioactive wastes in repositories buried in deep geologic formations of presumed high stability. Licensing of such repositories is conditional on the assessment of the reliability of such confinement for quite extensive periods of time ranging from 10 000 to 1 million years. Over such long time scales the repositories may leak, thus releasing contaminants to the surrounding medium. The released contaminant may reach the groundwater and be transported for substantial distances, eventually returning to the biosphere, where it poses a hazard for human health and the environment. In the present paper, we propose a phenomenological model based on the Kolmogorov-Dmitriev theory of branching stochastic processes in which a distributed source releases contaminant which is transported over a bidimensional medium whose physical properties are stochastic in space. The proposed model is implemented within a Monte Carlo scheme which exploits variance reduction techniques. © 1998 Elsevier Science Ltd. All rights reserved
1. INTRODUCTION
In the field of radioactive waste disposal, the technical evaluation of the environmental protection activities and the performance assessment of the devised solutions have to rely, to a considerable extent, on results obtained by predictive computer codes which in turn are based on mathematical models of the system. In general, the demonstration of the safety of these systems is not a straightforward task because of the long-term potential hazard posed by the radioactive waste. Due to this long time frame (thousands or even millions of years) over which safety compliance is sought, it is evident that there can be no direct experimental proof of disposal performance in terms of safety and reliability. Inevitably, then, the demonstration of compliance with the specific regulations is based on a theoretical evaluation of the disposal system. *Author for correspondence. Fax: 0039-2-2399-6309; e-mail: [email protected] 1301
1302
M. Marseguerra and E. Zio
In this context, one of the predominant phenomena which needs to be properly simulated is the transport of contaminants in groundwater. Indeed, the intention of placing these disposal systems in deep geologic formations is such as to make groundwater a fundamental vector of transport for the return of the contaminant to the biosphere. Computational modelling of the relevant processes influencing groundwater transport involves the consideration of a variety of physical and chemical phenomena. As an open system, groundwater may exchange mass and energy with the neighbouring soil, atmosphere and surface waters. The exchange occurs through physical-chemical processes of adsorption, ion-exchange, infiltration, evaporation, inflow, outflow, etc. Human activities can, also, strongly affect the characteristics of groundwater, and its chemical and biological constituents, by introducing substances which are often harmful. If the level of such constituents is beyond a given standard, then the groundwater is said to be contaminated or polluted. Mathematical models have been developed which can deal with the complex physical structures and irregular geometric shapes of aquifer systems and can describe the complex physical-chemical-biological phenomena that typically occur in the porous media hosting the groundwater flow. The main approaches used to model the contaminant transport in groundwater can be categorized as: 1. classical advection-dispersion theory; 2. Monte Carlo simulation; 3. transport theory. The classical advection-dispersion theory constitutes the most common approach to modelling contaminant transport and is based on a mass balance over a control volume and the assumption of validity of Fick's law to describe the mechanisms of hydrodynamic dispersion. The resulting partial differential equation needs to be coupled to the equation describing the flow of groundwater in the medium to determine the groundwater velocity field, i.e. the advective velocity. A common way for taking into account phenomena of adsorption-desorption in the host rock consists of reducing the dispersion coefficient and the advective velocity by means of a retardation factor on the basis of a linearity assumption between the amount of contaminant dissolved in the fluid and that trapped in the matrix of the host rock. A major drawback of this model is that it introduces a local parameter, the dispersivity, which does not satisfactorily account for the experimental data since it actually turns out to be an integral parameter dependent on the whole system under investigation. Indeed, many field tracer tests, covering various spatial scales, performed to determine the values of dispersivity for different media, give indications that the dispersivity values depend on the scale of the test. In general, the difference among values obtained in field and laboratory tests is attributed to the effects of heterogeneities on the macroscopic flow field. Since most heterogeneities in geological materials occur at a scale larger than that which can be included in borehole samples, dispersivity values from tests on small samples give a property of the medium at a scale of insufficient size for general use in prediction of dispersion in the field. In other words, the estimated values of dispersivity depend on the distance or time of tracer propagation (the so-called scale effect). For example, Gelhar et al. (1992) reviewed the data for 59 different field sites and found values of longitudinal dispersivity ranging over several orders of magnitude. This fact reveals the inconsistency of considering this parameter as a physical property of the medium only.
Contaminant transport in bidimensional porous media
1303
Monte Carlo methods have been applied to problems of contaminant transport in groundwater (Prickett et al., 1981). A large number of histories of particle motion are generated by random sampling from appropriate probability distribution functions. Each particle is engaged in two kinds of movements: one is advection, represented by the movement of the tracer particle jointly with the flow field; the other is dispersion, which may be seen as a random fluctuation superimposed on the average movement towards the zones with smaller contaminant concentrations. The particles are then followed in their movements according to equations describing their motion and other physical-chemical-biological processes such as adsorption and desorption. More recently, a new model of contaminant transport in fractured media has been developed by Williams (1992), based on an analogy with neutron transport. This model interprets the points of branching of the fractures as collision points in correspondence of which a pseudo-scattering reaction occurs, and the average distance between successive branching points is treated as an effective mean free path. The balance equation for the mass concentration results in an integro~lifferential transport equation for its expected value, the randomness being implictly contained in the kernel as the distribution of the lengths and orientations of the fractures. This approach has the important advantage that the many similarities with the time-dependent neutron transport equation render several of the existing, validated numerical codes almost directly applicable. However, this theoretically sound approach might bear some limitations in that it becomes burdensome when many transformations (physical and chemical) of the contaminant substances need to be considered. In a previous paper (Marseguerra and Zio, 1997), the present authors have proposed a new model for the description of contaminant transport in groundwater systems, based on the Kolmogorov-Dmitriev theory of branching stochastic processes (Kolmogorov and Dmitriev, 1947). The approach is similar to that of Antonopoulos-Domis et al. (1995) and Lee and Lee (1995) in that it is based on a discretization of the space variable and on consideration of the individual transitions that the contaminant particles may undergo according to specified probabilities. According to this model, several kinds of particles corresponding to the different contaminant locations and states are introduced and, then, suitable systems of coupled ordinary differential equations for the number of particles of each kind are derived. The problem is, thus, tackled in a very direct and straightforward manner and the analysis is significantly simplified by working in terms of probability generating functions (PGFS). A first advantage of the approach is that it has a flexible structure which is well suited for dealing with multidimensional geometries and also for modelling various phenomena such as chemical reactions and radioactive decay. Besides the expected value of the particles' populations, it can also yield equations for the variances-covariances and higher moments. Moreover, the adsorption~lesorption processes may be easily taken into account, so that the level of radioactive contaminant that remains trapped in the host medium and, therefore, the total activation at a site may be estimated. An obvious limitation of the analytical approach to the solution lies in the fact that the larger the variety of contaminant particles considered, the higher the order of the system of coupled differential equations to be solved: this might lead to numerical problems in the solution. For these reasons, in this paper we attempt to inherit the Kolmogorov-Dmitriev theory of branching stochastic processes within a Monte Carlo scheme of simulation. The great
1304
M. Marseguerra and E. Zio
advantage is that with the current power of modern computer systems, there are basically very few limitations on the number of particle kinds and interactions that one can simulate. In addition, the proper use of variance reduction techniques can further increase the efficiency of the calculation in terms of accuracy. Moreover, the additional flexibility introduced by the Monte Carlo scheme is such as to allow the introduction of elements of stochasticity in the properties of the media. This approach falls very well along the lines of stochastic hydrogeology, where the properties of the media are considered as stochastic processes in order to account for all the occurring heterogeneities (de Marsily, 1986). The paper begins, in Section 2, with a short outline of the concepts underlying the application of the Kolmogorov-Dmitriev theory of branching processes to the contaminant transport in groundwater (Marseguerra and Zio, 1997): In Section 3 we present the bidimensional stochastic model for contaminant transport and all the possible events and processes that are considered. Section 4 is devoted to the actual Monte Carlo simulation whose results are reported in Section 5. A biasing technique tailored for our model is also described in Section 4. Some remarks on the approach and some indications on future research on the subject conclude our work.
2. THE KOLMOGOROV-DMITRIEV MODEL In this section we shall summarize the basic features of the Kolmogorov-Dmitriev model as applied to the transport of contaminants in groundwater. The interested reader should refer to Marseguerra and Zio (1997) for further details. Let us consider the problem of contaminant transport within a medium which, in principle, may be mono-, bi- or tridimensional. The spatial domain along which the transport occurs is subdivided in Nz discrete zones, z = 1, 2 . . . . . Nz. In each zone we wish to determine the amount of contaminant dissolved in the water and that which is adsorbed in the host medium. To this aim we introduce two categories of particles, namely the solutons, which are the particles of contaminant dissolved in groundwater, and the trappons, which are the particles of contaminant trapped in the host medium. We assume that each soluton can, in principle, give rise to several types of trappons. This allows us to account for different types of trapping rock materials characterized by different retardation mechanisms such as physical adsorption, chemical adsorption, direct incorporation into the rock structure (mineralization) or precipitation caused by enhanced concentrations at the rock surface. If we denote by Sz, z = 1, 2 . . . . . Nz, the number of solutons in zone z, called z-solutons, and by t/, z = 1,2 . . . . . Nz, i = 1, 2 . . . . . Art, the number of trappons of the ith kind in zone z, called iz-trappons, we are dealing with m = N z ( N t + 1) particle kinds. Obviously, the generalization to several kinds of solutons is straightforward. Within a time interval dr, one soluton can undergo an adsorption transition in zone z, thus transforming into a trappon of one of the N t different kinds, in zone z. We denote by b~,z the rate at which this transition occurs. Alternatively, the soluton can travel to a neighbouring zone z' during dr, with transition rate bs,z~z,, z' = 1, 2 . . . . . Nz, z' # z. Concerning one i-trappon in zone z, it can only 'desorb' to a soluton with reaction rate d. The trappons, then, show many similarities, in their adsorbing-desorbing behaviour, with the neutron precursors of reactor physics.
Contaminant transport in bidimensional porous media
1305
Finally, the contaminant might undergo a radioactive decay with rate 2. For simplicity's sake we assume that the daughter of the decay is neither radioactive nor contaminant, so that it is of no further interest in the analysis. Obviously, this assumption may be easily removed. By working in the P G F domain and exploiting the properties of the PGFs, we arrive at the following system of Nz(Nt + 1) coupled differential equations for the expected values S(z*, 31lsz, 0) and T"(z*, rl 1,z, 0) of the numbers of z*-solutons and t~z*-trappons at time 3, given that we started from one soluton in z at time r = 0:
O--S(z*,rll,2,0)=03
/
b~z.+ Ebs,z--,~,+A S(z*,rllsz, O) Zt~Z*
Nt
+ Z bs,e~z.S(z', 311,z, O) + E d ' T i ( z * , vii,z, O) z'#z*
z* =
1,2 . . . . . Nz
i=l
(1) 0
z" • , 311,z, O) r~ (z*, 311,z, 0) -- b,z.S(z
(d z" + A)~(z*, 3[lsz, O)
z * = 1,2 . . . . . Nz
~ = 1,2 ..... Nt (2)
These equations appeal to our physical intuitions as they represent the balance between production and destruction processes. In equation (1), the first addendum at RHS accounts for all processes that make the solutions disappear from z*, by transformation in a trappon, by travel to a different zone or by radioactive decay; the second term is a production term of solutons that, starting from another zone z' ~ z*, reach z* at time 3; the last term is also a production term representing desorption, at time 3, of the contaminant trapped in zone z*. A similar interpretation may be given to the terms in equation (2). The initial conditions for the above equations are
S(z*,Oll,z,O)=Szz. ~i(z*, 011,~, 0 ) = 0 Let us now consider a case in which a Poissonian source of variable strength s(r) emits solutons in the location z during the time interval (0, 3s). In this case, the following differential equation for the evolution in time of the expected values of the soluton and trappon populations is obtained:
i + z,,#z. Ss(z*, 31z, O, 3s)= s(r),~z.z - t FNt b'z* ~ bs,z.--,z,, +2 I Ss(z*, 31z, O, 3s) Nt + Zdir~(z i=l
*, r[2, O, 3s) + Z bs,z,-->z*Ss(zt, l"12, O, 3s) z'# z*
z* = 1,2 . . . . . N~
(3)
1306 0 __
0r
M. Marseguerra and E. Zio
T s1"*(z* , rlz, 0, rs)
=
b ,~.(l"* _~z._o,~
_* , rlz, 0, r~) - (d r + 2)T](z*, rlz, 0, rs)
z * = 1,2 . . . . . /*=1,2 .....
Nz Nt (4)
The physical meaning of the various terms of destruction and production follows immediately. The main advantage of the Kolmogorov-Dmitriev approach is that it allows the computing of the higher moments of the relevant distributions. In particular, it turns out that in this case of a Poissonian source the process remains rigorously Poissonian so that the variances of the populations equal the respective mean values.
3. A BIDIMENSIONAL STOCHASTIC M O D E L FOR CONTAMINANT TRANSPORT Let us consider an area in which flow of groundwater, and the associated contaminant transport, are occurring. In our bidimensional model, this region of space is discretized in (NZrNZj) rectangular pixels, where NZi and NZj are the number of discretization steps in the horizontal and vertical directions, respectively. The location of the source of contamination in this region can, then, be identified with corresponding elements of coordinates (Is, J~). The migration of the contaminant from the source to the surrounding zones is modelled by means of eight transition rates (Fig. 1): f and b are the probability densities of the contaminant making a forward and backward jump to the adjacent neighbouring pixel in the horizontal direction, respectively; u and d are the probability densities of the contaminant making an upward and downward jump to the adjacent neighbouring pixel in the vertical direction, respectively; fu = x/f 2 + u2 (forward-upward),
bu
\
U ,L
\/ b
f
".
/ bd
/
fa
,p
d
\ fd
Fig. 1. Possible moves of the contaminant.
Contaminant transport in bidimensional porous media
1307
f d = v/f 2 + d 2 (forward-downward), bu = ~ u2 (backward-upward), and bd --- V ~ + d 2 (backward-downward) are probability densities for diagonal moves. By a proper assignment of the values of these parameters one may also consider any kind of subdomain D within the above defined rectangular domain, simply by setting to zero the transition probabilities outside D. During its travel within the region, the contaminant may undergo physical-chemical processes of adsorption and desorption which typically slow down its transport. These phenomena can be modelled quantitatively by means of two additional transition rates, ads, for the adsorption process, and des, for the desorption process. The rates of movement f, b, u, d are characteristic of the flow regime which is ultimately determined by the physical properties of the host medium; similarly, the processes of adsorption and desorption are intrinsic characteristics of the solid matrix. Possible heterogeneities of the medium may be accounted for by considering these parameters as random variables in space, with assigned distributions. Finally, the contaminant may also decay, if radioactive. The evolution of the decay products, if still hazardous, can very well be followed simply by introducing new kinds of particles with their characteristic rates of movement within, and interaction with, the host matrix. A mission time, Tmiss, is defined as the period of interest for the groundwater transport of contaminants, and a given number NF of final locations of interest are selected within the region as observation points for the final amount of contaminant concentration at Tmiss. These locations could represent drinking water wells, public living areas, discharge points and the like.
4. MONTE CARLO SIMULATION: A NUMERICAL EXAMPLE We consider a rectangular region of space subdivided into NI.Nj = 40.20 = 800 identical, rectangular pixels. The bounds of this region are considered absorbing. A source of contaminant is located in the southwest corner (1,1), whereas ArE = 5 observation wells are located in the upper region of the east bound at locations (40,20), (40,18), (40,16), (40,14), (40,12), respectively (Fig. 2). To account for spatial heterogeneity, the transition rates f, b, u, d, ads, des are assumed to be all log-normally distributed with percentile parameters given in Table 1. The Monte Carlo trials are divided in bins and at the beginning of each bin a different set of values of these parameters is sampled from their respective distributions. For what concerns the nature of the contaminant, we have considered a radioactive element which decays with a time constant of 20 000 years in a non-hazardous daughter whose transport poses no additional risk and, thus, needs not be modelled further. The transport of the contaminant is followed over a period of 105 years, which is a typical regulatory time for high level radioactive waste repositories. Both analogue and biased simulations have been considered. For the biased Monte Carlo simulation a tailored technique has been proposed which favours contaminant moves directed towards one of the assigned observation locations. At the beginning of each trial, one of the NF observation locations is selected at random as the target pixel (IF, Jr) of contaminant accumulation for that trial.When the contaminant particle enters a defined area of interest around the target pixel we favour those transitions which bring the contaminant closer to the target pixel. The area of interest around the pixel (IF, Jr) is
1308
M. Marseguerra and E. Zio
J~
...............................
4"~0 4~ 4~,i6 4~,4
h 1
40
I
Fig. 2. The region in which contaminant transport occurs. Table 1. 5th and 95th percentiles of the log-normal distributions of the parameters f, b, u, d, ads, des Parameter f b u d ads des
P5
P95
5x10 -3 5x10 -4 5x10 -5 5.5x 10-4 5x 10-3 2x 10-4
8x10-3 8×10 -4
8x10 -5 8.5x 10 -4
8x 10-3 3x 10-4
a circle of radius r centred at the pixel: it is only within this region that biasing takes place (Fig. 3). By so doing, we bias only histories which already point in the right direction, towards (IF, JF), thus avoiding forcing indiscriminately all the histories, even those directed in wrong directions for which the biasing would require great amounts of time and significant reductions of the weights eventually collected. Furthermore, in order to have a certain control on the weights, we gradually increase the biasing within the region of influence by employing a smooth function 0(/) which starts from 00 on the circumference and increases towards the centre (IF, JF), up to a final value OF. More specifically, the current position (Ic, Jc) of the contaminant is compared to (IF, JF) and the transition rates of all those moves that bring the contaminant closer to the target are suitably increased. The amplification factor here adopted is: 0(0= l+a(0e-~
l
where I is the distance from the current to the final location and
(5)
Contaminant transport in bidimensional porous media
1309
J
/
Fig. 3. Circle of influence for biasing. ot = (OF -- 1).e~ _
-_____._0__0 (r
-,o,(V)
(6)
l = ~//(Ic -- IF) 2 Jr- ( J c -- J F ) 2
Fig. 4 shows a typical behaviour of 0(/). It is important to stress that the values of the functions ot and/~ in equation (6) should be chosen coherently with the choice of r, so as to guarantee a 'balanced' biasing. In other words, if r contains too many pixels, the successive multiplications of the weight by 0 can result in a vanishing final weight, a situation which should be avoided. In favouring certain moves, we obviously disfavour others which, when sampled, largely increase the weight. The presence of a small number of histories carrying large weights can significantly deteriorate the statistics of the results at later times. To circumvent this problem, and obtain statistically significant results throughout the mission time without increasing the number of trials beyond reasonable limits, we coupled our biasing technique with the splitting technique (Hughes et al., 1958): when the weights of the unfavoured histories exceed pre-cstablished thresholds, the history is split into two (or more) subhistories from that time point on, each one having half (or less) of the initial weight. The splitting thus generates a number of subhistories which distribute their contributions to reinforce the statistics.
5. RESULTS First an analogue simulation with l06 histories was performed, requiring a computer time of 1.72 min on an Alpha Station 250 4/266. Figure 5 reports the time evolution of the contaminant concentration in correspondence of the five target locations. As can be seen,
1310
M. Marseguerra and E. Zio 10
i
Amplification factor for biasing i
i
i
,
10
15
20
25
30
35
40
I
Fig. 4. Behaviour of the amplification factor 0 vs distance from the target pixel. the event of the contaminant travelling from the source in the south-west corner to the target point (40,20) at the opposite north-east corner is rather rare, whereas more probable events are those concerning the contaminant going to more central targets on the east bound, such as (40,14) and (40,12), for example. Contaminant concentration at the target locations: analo case 1.2 i 10.4 i(40,12)
0.8
o
0.6
0.4
0.2
1
2
3
4
5 time(y)
6
7
8
9
10 x 104
Fig. 5. Time evolution of the contaminant concentration at the target locations: analogue simulation.
Contaminant transport in bidimensional porous media
1311
This behaviour is confirmed, as expected, by the biased simulation. Again, 106 histories were performed, requiring 6.33 min of computer time on the same machine. The radius of the circle of influence of each target pixel was chosen to be equal to 40, with a starting bias amplification factor 00 --- 2 on the circumference rising up to a value OF = 10 at the pixel. Splitting is set up in such a way as to occur when the weight carried by the contaminant particle in the current history goes beyond a threshold value of 1.5. This has generated 2 026246 additional subhistories which have allowed us to keep the magnitude of the collected weights to an average of 4.75x 10 -4 with a standard deviation of 5.30x 10-6. This ensures good statistical results, which are shown in Fig. 6: the greater smoothness of the curves is an indication of better statistical accuracy as compared to the analogue case. The efficiency of the biasing obviously depends on how rare the events of interest are: if they are reasonably probable, then biasing is likely to be unnecessary; on the other hand, if they are quite rare, then, biasing can lead to a significant improvement in the accuracy of the estimates. In our case, the move probabilities are such that arrival at the target pixels is less and less probable as we go up from the target in (40,12) to that in the northeastern corner (40,20). In Table 2 we compare the analogue and biased cases on the basis of the usual figure of merit given by the inverse of the product of the variance times the computer time required. In the present case, each Monte Carlo game has five different target pixels, so that the figure 1 The of merit of the computation is actually a vector M with five components Mr,-= (~)itc," computer time tci appearing in the ith component of M results from the assumption that the total computer time is equally shared among all target pixels. For each time instant, the standard deviation (trc)i appearing in the ith component of M gives a measure of the accuracy of the estimate of the contaminant concentration, C, in the ith target location. For t * = 2 • l04 years, the vector of figures of merit M is reported in the table, each component being referred to a different pixel. It can be seen that the biasing is not efficient Contaminant concentration at the target locations: biased case
x 10 4 1.2 ~
~
T
T
r
r
r
T
6
7
8
9
(40,12) 1
0.8
0.6
0.4
0.2
1
2
4
5
time(y)
10 x 10~
Fig. 6. Time evolution of the contaminant concentration at the target locations: biased simulation.
1312
M. Marseguerra and E. Zio
Table 2. Figures ofmerit for the estimates ofconcentration in the various pixels at t" = 2 x 104 years; tc, = computer time (min) for the ith target pixel (total computer time divided by the number of target pixels); (trc)i = standard deviation of the concentration in the ith target pixel at t*; M i = I Target (IF, JF) (40,12) (40,14) (40,16) (40,18) (40,20)
(trc)i analogue (to, = 0.34 min)
(ac)i biased (tc~ = 1.27 min)
Mi analogue (t~, = 0.34 rain)
Mi biased (t~, = 1.27 rain)
3.03× 10 - 6
1.81 × 10-6
2.13x 10-6 1.40× 10-6
3.09x 10 - 6
2.17× 1011 6.40× 10I1 1.48× 1012 3.52× 1012 0
2.41x 10I1 8.25x 1011 1.25x 1013 5.31X 1013
2.52x 10-6 1.22x 10 -7 4.95x 10-s
8.56x 10 - 7
cx~
3.23x 1014
for the relatively probable event of arrival in (40,12) but it becomes more and more efficient, with respect to the analogue simulation, as we go towards the more distant targets. In particular, no contribution is carried by the analogue simulation to the farthest target (40,20), thus giving rise to a null figure of merit. As expected, then, the efficiency of the biasing depends on the probability of arriving at the target pixels and, therefore, on the values of the transition rates of movement between adjacent pixels. The less probable are the movements leading in the direction of the targets, the more efficient is the biasing. A global comparison between the analogue and the biased computations can be done in terms of the average of the components of the vector M. In the present case we have Manalog =
1.17 × 10 j2
Mbiased = 7.78 X
1013
which correspond to a mean gain of 66 in favour of the biased simulation. The higher accuracy achieved by the biasing is also shown in the comparison of Figs 7 and 8, which report the distribution of arrival times at the target location (40,12) in the analogue and biased case, respectively, and associated error bands. For this quantity, Table 3 reports the comparison in terms of the standard figure of merit. At each time, the accuracy of the results is quantified in terms of the standard deviation % of the distribution of arrival times, p, in (40, 12). For t = 105 years, the figure of merit M = ~ shows I P~ that, in this case, the biased simulat'on pe r fo r n l s almost six times better than the analogue one. In the above example, in correspondence of t = Tmiss the analogue histories carry contributions to the distribution only in the target pixel (40,12), whereas no contribution Table 3. Comparison between the analogue and biased simulations: to, = computer time (min) for the target pixel (40,12) (total computer time divided by the number of target pixels); trq = standard deviation of the probability of arrival in (40,12) at Tmiss= 105 years; M =
Analogue Biased
tc, (min)
ap
M
0.34 1.27
9.95 × 10 - 7 2.17× 10 -7
2.94× 1012 1.67× 1013
Contaminant transport in bidimensional porous media
13 l 3
Distribution of arrival times at the location (40,12): analog case i i + + i i
10 -s
1.8.
/
1.6 ................................
t
...................................
......
1.4
1.2
1 0" 0.8
0.6
0.4
0.2
0
1
2
3
4
7
8
9
10
time(y)
X 10 4
Fig. 7. Distribution of arrival times at location (40, 12) analogue simulation. arrives in the other targets. The latter situation corresponds to a value of zero in the associated components of the vector of figure of merit. On the contrary, the biased histories are such as to provide contributions to all target pixels, i.e. to all components of the vector of figure of merit. To further investigate this aspect, we have considered a second case study concerning the transport of Tc-99 (2 = 3.24 x 10 -6 years -1) in the same spatial region. The properties
Distribution of arrival times at the location (40,12): biased case ! ! ! ~ ! ! +
x 10 -s
1.8
i
:
i
:
i
i
1.4 . . . . . . . . . . . . .
i ................
:. . . . . . . . . . .
1.2 ..............
: : .................
:
i :.......
:
: .........
:. . . . . .
i : .......
:. . . . . .
.......
: 1 ............................
!
....
: . . . .
.........
.........
!
:..............................
: ..........
i ..........
3
6
7
0"0. 8
0.6
0.4
0.2
6'
1
2
4
5 time(y)
8
9
10 x 10 4
Fig. 8. Distribution of arrival times at location (40,12): biased simulation.
M . Marseguerra and E. Zio
1314 x 10-r
Contaminant concentration at (40,12): analog case
6
5
4
o3
2
1
I!111II II Illilllllllll
I:111IIItIIl:lllllll[ll: 0
li
~
f
2
3
4
. ~ 5 time(y)
9
10 x 104
Fig. 9. Time evolution of the contaminant concentration in (40,12): analogue simulation.
of the host medium are not varied except for the percentiles of the distribution of the downward rate, d, which are an order of magnitude larger than those of Table 1. This renders the event of contaminant accumulation in the north-eastern target points much less probable. No modification was made to the parameters for the biased simulation. Figures 9 and 10 report the time evolution of the contaminant concentration in the target location (40,12) as resulting from the analogue and biased simulations, respectively. x 10-r
Contaminant concentration at (40,12): biased case
v
i i
3 ........ !....... i .........
0
0
2
:~ .........
i..........
3
4
i........
......
6
time(y)
i ..........
7
i .........
-
! .......
9
10 x 104
Fig. 10. Time evolution of the contaminant concentration in (40,12): biased simulation.
Contaminant transport in bidimensional porous media
1315
Table 4. Comparison between the analogue and biased simulations: to, = computer time (min) for the target pixel (40,12) (total computer time divided by the number of target pixeh); a e = standard deviation of the probability of arrival in (40,12) at Tmiss = 105 years; M = to,
Analogue Biased
(min)
1.20× 102 7.7
trp
M
8.29x 10-8 4.62× 10-8
1.21 x 1012 6.08× 1013
It can be noted how the arrival of contaminant at this location is almost three orders of magnitude less probable than in the first case examined. The analogue case required 108 histories and a computer time of 974min to provide results of any statistical significance, whereas 106 histories were performed for the biased case which, however, benefited from the contributions of 7 896 086 split subhistories. Overall, the computer time required in the biased case was 38.46 min and 191 248 weights arrived at a destination in one of the target points, each weight bringing, on average, a contribution of 1.25x 10-6 with a standard deviation of 6.81 x 10-8. Table 4 reports the details of the comparison between the analogue and biased cases: the values of the figure of merit M show that in this case of very rare events the biased simulation provides a performance which is 50 times better than that of the analogue simulation.
6. CONCLUSIONS Many environmental management and technical evaluation activities rely on the extensive use of predictive computer models to describe the phenomena occurring over the time frames required by the regulations. Natural groundwater transport is one such phenomenon which is of particular interest because groundwater is a principal vector through which hazardous contaminant can circulate to the biosphere and, thus, pose a significant risk to the health of the human population and to the safety of the environment. Several approaches to modelling groundwater contaminant transport have been proposed in the past. In this paper we have focused on a relatively new model, previously proposed by the authors, which is based on the theory of branching stochastic processes of Kolmogorov and Dmitriev. The flexible structure offered by the model allows one to describe explicitly many of the non-linear phenomena which actually occur in practice. The degree of detail which can be reached in the analytical model is, however, limited by the high order of the system of coupled differential equations to be solved. To circumvent this problem, in this paper we have proposed to structure the Kolmogorov and Dmitriev approach within a pure Monte Carlo framework. The great advantage is that with the current power of modern computer systems, there are basically very few limitations on the number of particle kinds and interactions that one can simulate. In addition, the proper use of variance reduction techniques can further increase the efficiency of the calculation in terms of accuracy. In this context, a biasing technique has been introduced, which forces the contaminant to travel towards the accumulation points of interest within the region under study. This is coupled with a splitting technique to ensure a proper distribution of the contributions to the statistics of the results throughout the
1316
M. Marseguerra and E. Zio
mission time. The efficiency of this technique has been proven through a sample application on a bidimensional region. Finally,the additional flexibility introduced by the Monte Carlo scheme has allowed us to accommodate elements of stochasticity in the physical properties of the media to account for possibly existing heterogeneities.
REFERENCES
Antonopoulos-Domis, M., Clouvas, A. and Marseguerra, M. (1995) On the compartmental modeling of cesium migration in soils. Nucl. Sci. & Eng. 121,461. Gelhar, L. W., Welty, C. and Rehfeldt, K. R. (1992) A critical review of data on fieldscale dispersion in aquifers. Water Resour. Res. 28(7), 1955. Hughes, D. J., Sanders, J. E. and Horowitz, J. (1958) Physics and mathematics. Progress in Nuclear Energy 2, 315-368. Kolmogorov, A. N. and Dmitriev, N. A. (1947) Dokl Akad. Nank 5553 56, 7. Lee, Y. and Lee, K. J. (1995) Nuclide transport of decay chain in the fractured rock medium: a model using continuous time Markov process. Ann. Nucl. Energy 22(2), 71. Marseguerra, M. and Zio, E. (1997) Modelling the transport of contaminants in groundwater as a branching stochastic process. Ann. Nucl. Energy 24(8), 625-644. Marsily, G. de (1986) Quantitative Hydrogeology. Academic Press, San Diego, California. Prickett, T. A., Naymik, T. G. and Lonnguist, C. G. (1981) A 'random-walk' solute transport model for selected groundwater quality evaluations. ISWS/BUL-65/81, Bulletin 65, State of Illinois, Illinois Department of Energy and Natural Resources, Champaign. Williams, M. M. R. (1992) A new model for describing the transport of radionuclides through fractured rock. Ann. Nucl. Energy 19(10), 791.
NucL Energy, Vol. 25, No. 16, pp. 1301-1316, 1998 © 1998 Elsevier Science Ltd. All fights reserved PII: S0306-4549(98)00015-2 Printed in Great Britain 03064549/98 $19.00 + 0.00 Ann.
C O N T A M I N A N T T R A N S P O R T IN B I D I M E N S I O N A L P O R O U S M E D I A VIA B I A S E D M O N T E C A R L O SIMULATION M. MARSEGUERRA* and E. ZIO Department of Nuclear Engineering, Polytechnic of Milan, Via Ponzio 34/3, 20133 Milano, Italy (Received 23 December 1997)
A b s t r a c t ~ n e of the most vulnerable spots of nuclear power systems is the hazard posed by the radioactive wastes at the end of the fuel cycle. Currently, the most widely pursued end solution consists of disposing the radioactive wastes in repositories buried in deep geologic formations of presumed high stability. Licensing of such repositories is conditional on the assessment of the reliability of such confinement for quite extensive periods of time ranging from 10 000 to 1 million years. Over such long time scales the repositories may leak, thus releasing contaminants to the surrounding medium. The released contaminant may reach the groundwater and be transported for substantial distances, eventually returning to the biosphere, where it poses a hazard for human health and the environment. In the present paper, we propose a phenomenological model based on the Kolmogorov-Dmitriev theory of branching stochastic processes in which a distributed source releases contaminant which is transported over a bidimensional medium whose physical properties are stochastic in space. The proposed model is implemented within a Monte Carlo scheme which exploits variance reduction techniques. © 1998 Elsevier Science Ltd. All rights reserved
1. INTRODUCTION
In the field of radioactive waste disposal, the technical evaluation of the environmental protection activities and the performance assessment of the devised solutions have to rely, to a considerable extent, on results obtained by predictive computer codes which in turn are based on mathematical models of the system. In general, the demonstration of the safety of these systems is not a straightforward task because of the long-term potential hazard posed by the radioactive waste. Due to this long time frame (thousands or even millions of years) over which safety compliance is sought, it is evident that there can be no direct experimental proof of disposal performance in terms of safety and reliability. Inevitably, then, the demonstration of compliance with the specific regulations is based on a theoretical evaluation of the disposal system. *Author for correspondence. Fax: 0039-2-2399-6309; e-mail: [email protected] 1301
1302
M. Marseguerra and E. Zio
In this context, one of the predominant phenomena which needs to be properly simulated is the transport of contaminants in groundwater. Indeed, the intention of placing these disposal systems in deep geologic formations is such as to make groundwater a fundamental vector of transport for the return of the contaminant to the biosphere. Computational modelling of the relevant processes influencing groundwater transport involves the consideration of a variety of physical and chemical phenomena. As an open system, groundwater may exchange mass and energy with the neighbouring soil, atmosphere and surface waters. The exchange occurs through physical-chemical processes of adsorption, ion-exchange, infiltration, evaporation, inflow, outflow, etc. Human activities can, also, strongly affect the characteristics of groundwater, and its chemical and biological constituents, by introducing substances which are often harmful. If the level of such constituents is beyond a given standard, then the groundwater is said to be contaminated or polluted. Mathematical models have been developed which can deal with the complex physical structures and irregular geometric shapes of aquifer systems and can describe the complex physical-chemical-biological phenomena that typically occur in the porous media hosting the groundwater flow. The main approaches used to model the contaminant transport in groundwater can be categorized as: 1. classical advection-dispersion theory; 2. Monte Carlo simulation; 3. transport theory. The classical advection-dispersion theory constitutes the most common approach to modelling contaminant transport and is based on a mass balance over a control volume and the assumption of validity of Fick's law to describe the mechanisms of hydrodynamic dispersion. The resulting partial differential equation needs to be coupled to the equation describing the flow of groundwater in the medium to determine the groundwater velocity field, i.e. the advective velocity. A common way for taking into account phenomena of adsorption-desorption in the host rock consists of reducing the dispersion coefficient and the advective velocity by means of a retardation factor on the basis of a linearity assumption between the amount of contaminant dissolved in the fluid and that trapped in the matrix of the host rock. A major drawback of this model is that it introduces a local parameter, the dispersivity, which does not satisfactorily account for the experimental data since it actually turns out to be an integral parameter dependent on the whole system under investigation. Indeed, many field tracer tests, covering various spatial scales, performed to determine the values of dispersivity for different media, give indications that the dispersivity values depend on the scale of the test. In general, the difference among values obtained in field and laboratory tests is attributed to the effects of heterogeneities on the macroscopic flow field. Since most heterogeneities in geological materials occur at a scale larger than that which can be included in borehole samples, dispersivity values from tests on small samples give a property of the medium at a scale of insufficient size for general use in prediction of dispersion in the field. In other words, the estimated values of dispersivity depend on the distance or time of tracer propagation (the so-called scale effect). For example, Gelhar et al. (1992) reviewed the data for 59 different field sites and found values of longitudinal dispersivity ranging over several orders of magnitude. This fact reveals the inconsistency of considering this parameter as a physical property of the medium only.
Contaminant transport in bidimensional porous media
1303
Monte Carlo methods have been applied to problems of contaminant transport in groundwater (Prickett et al., 1981). A large number of histories of particle motion are generated by random sampling from appropriate probability distribution functions. Each particle is engaged in two kinds of movements: one is advection, represented by the movement of the tracer particle jointly with the flow field; the other is dispersion, which may be seen as a random fluctuation superimposed on the average movement towards the zones with smaller contaminant concentrations. The particles are then followed in their movements according to equations describing their motion and other physical-chemical-biological processes such as adsorption and desorption. More recently, a new model of contaminant transport in fractured media has been developed by Williams (1992), based on an analogy with neutron transport. This model interprets the points of branching of the fractures as collision points in correspondence of which a pseudo-scattering reaction occurs, and the average distance between successive branching points is treated as an effective mean free path. The balance equation for the mass concentration results in an integro~lifferential transport equation for its expected value, the randomness being implictly contained in the kernel as the distribution of the lengths and orientations of the fractures. This approach has the important advantage that the many similarities with the time-dependent neutron transport equation render several of the existing, validated numerical codes almost directly applicable. However, this theoretically sound approach might bear some limitations in that it becomes burdensome when many transformations (physical and chemical) of the contaminant substances need to be considered. In a previous paper (Marseguerra and Zio, 1997), the present authors have proposed a new model for the description of contaminant transport in groundwater systems, based on the Kolmogorov-Dmitriev theory of branching stochastic processes (Kolmogorov and Dmitriev, 1947). The approach is similar to that of Antonopoulos-Domis et al. (1995) and Lee and Lee (1995) in that it is based on a discretization of the space variable and on consideration of the individual transitions that the contaminant particles may undergo according to specified probabilities. According to this model, several kinds of particles corresponding to the different contaminant locations and states are introduced and, then, suitable systems of coupled ordinary differential equations for the number of particles of each kind are derived. The problem is, thus, tackled in a very direct and straightforward manner and the analysis is significantly simplified by working in terms of probability generating functions (PGFS). A first advantage of the approach is that it has a flexible structure which is well suited for dealing with multidimensional geometries and also for modelling various phenomena such as chemical reactions and radioactive decay. Besides the expected value of the particles' populations, it can also yield equations for the variances-covariances and higher moments. Moreover, the adsorption~lesorption processes may be easily taken into account, so that the level of radioactive contaminant that remains trapped in the host medium and, therefore, the total activation at a site may be estimated. An obvious limitation of the analytical approach to the solution lies in the fact that the larger the variety of contaminant particles considered, the higher the order of the system of coupled differential equations to be solved: this might lead to numerical problems in the solution. For these reasons, in this paper we attempt to inherit the Kolmogorov-Dmitriev theory of branching stochastic processes within a Monte Carlo scheme of simulation. The great
1304
M. Marseguerra and E. Zio
advantage is that with the current power of modern computer systems, there are basically very few limitations on the number of particle kinds and interactions that one can simulate. In addition, the proper use of variance reduction techniques can further increase the efficiency of the calculation in terms of accuracy. Moreover, the additional flexibility introduced by the Monte Carlo scheme is such as to allow the introduction of elements of stochasticity in the properties of the media. This approach falls very well along the lines of stochastic hydrogeology, where the properties of the media are considered as stochastic processes in order to account for all the occurring heterogeneities (de Marsily, 1986). The paper begins, in Section 2, with a short outline of the concepts underlying the application of the Kolmogorov-Dmitriev theory of branching processes to the contaminant transport in groundwater (Marseguerra and Zio, 1997): In Section 3 we present the bidimensional stochastic model for contaminant transport and all the possible events and processes that are considered. Section 4 is devoted to the actual Monte Carlo simulation whose results are reported in Section 5. A biasing technique tailored for our model is also described in Section 4. Some remarks on the approach and some indications on future research on the subject conclude our work.
2. THE KOLMOGOROV-DMITRIEV MODEL In this section we shall summarize the basic features of the Kolmogorov-Dmitriev model as applied to the transport of contaminants in groundwater. The interested reader should refer to Marseguerra and Zio (1997) for further details. Let us consider the problem of contaminant transport within a medium which, in principle, may be mono-, bi- or tridimensional. The spatial domain along which the transport occurs is subdivided in Nz discrete zones, z = 1, 2 . . . . . Nz. In each zone we wish to determine the amount of contaminant dissolved in the water and that which is adsorbed in the host medium. To this aim we introduce two categories of particles, namely the solutons, which are the particles of contaminant dissolved in groundwater, and the trappons, which are the particles of contaminant trapped in the host medium. We assume that each soluton can, in principle, give rise to several types of trappons. This allows us to account for different types of trapping rock materials characterized by different retardation mechanisms such as physical adsorption, chemical adsorption, direct incorporation into the rock structure (mineralization) or precipitation caused by enhanced concentrations at the rock surface. If we denote by Sz, z = 1, 2 . . . . . Nz, the number of solutons in zone z, called z-solutons, and by t/, z = 1,2 . . . . . Nz, i = 1, 2 . . . . . Art, the number of trappons of the ith kind in zone z, called iz-trappons, we are dealing with m = N z ( N t + 1) particle kinds. Obviously, the generalization to several kinds of solutons is straightforward. Within a time interval dr, one soluton can undergo an adsorption transition in zone z, thus transforming into a trappon of one of the N t different kinds, in zone z. We denote by b~,z the rate at which this transition occurs. Alternatively, the soluton can travel to a neighbouring zone z' during dr, with transition rate bs,z~z,, z' = 1, 2 . . . . . Nz, z' # z. Concerning one i-trappon in zone z, it can only 'desorb' to a soluton with reaction rate d. The trappons, then, show many similarities, in their adsorbing-desorbing behaviour, with the neutron precursors of reactor physics.
Contaminant transport in bidimensional porous media
1305
Finally, the contaminant might undergo a radioactive decay with rate 2. For simplicity's sake we assume that the daughter of the decay is neither radioactive nor contaminant, so that it is of no further interest in the analysis. Obviously, this assumption may be easily removed. By working in the P G F domain and exploiting the properties of the PGFs, we arrive at the following system of Nz(Nt + 1) coupled differential equations for the expected values S(z*, 31lsz, 0) and T"(z*, rl 1,z, 0) of the numbers of z*-solutons and t~z*-trappons at time 3, given that we started from one soluton in z at time r = 0:
O--S(z*,rll,2,0)=03
/
b~z.+ Ebs,z--,~,+A S(z*,rllsz, O) Zt~Z*
Nt
+ Z bs,e~z.S(z', 311,z, O) + E d ' T i ( z * , vii,z, O) z'#z*
z* =
1,2 . . . . . Nz
i=l
(1) 0
z" • , 311,z, O) r~ (z*, 311,z, 0) -- b,z.S(z
(d z" + A)~(z*, 3[lsz, O)
z * = 1,2 . . . . . Nz
~ = 1,2 ..... Nt (2)
These equations appeal to our physical intuitions as they represent the balance between production and destruction processes. In equation (1), the first addendum at RHS accounts for all processes that make the solutions disappear from z*, by transformation in a trappon, by travel to a different zone or by radioactive decay; the second term is a production term of solutons that, starting from another zone z' ~ z*, reach z* at time 3; the last term is also a production term representing desorption, at time 3, of the contaminant trapped in zone z*. A similar interpretation may be given to the terms in equation (2). The initial conditions for the above equations are
S(z*,Oll,z,O)=Szz. ~i(z*, 011,~, 0 ) = 0 Let us now consider a case in which a Poissonian source of variable strength s(r) emits solutons in the location z during the time interval (0, 3s). In this case, the following differential equation for the evolution in time of the expected values of the soluton and trappon populations is obtained:
i + z,,#z. Ss(z*, 31z, O, 3s)= s(r),~z.z - t FNt b'z* ~ bs,z.--,z,, +2 I Ss(z*, 31z, O, 3s) Nt + Zdir~(z i=l
*, r[2, O, 3s) + Z bs,z,-->z*Ss(zt, l"12, O, 3s) z'# z*
z* = 1,2 . . . . . N~
(3)
1306 0 __
0r
M. Marseguerra and E. Zio
T s1"*(z* , rlz, 0, rs)
=
b ,~.(l"* _~z._o,~
_* , rlz, 0, r~) - (d r + 2)T](z*, rlz, 0, rs)
z * = 1,2 . . . . . /*=1,2 .....
Nz Nt (4)
The physical meaning of the various terms of destruction and production follows immediately. The main advantage of the Kolmogorov-Dmitriev approach is that it allows the computing of the higher moments of the relevant distributions. In particular, it turns out that in this case of a Poissonian source the process remains rigorously Poissonian so that the variances of the populations equal the respective mean values.
3. A BIDIMENSIONAL STOCHASTIC M O D E L FOR CONTAMINANT TRANSPORT Let us consider an area in which flow of groundwater, and the associated contaminant transport, are occurring. In our bidimensional model, this region of space is discretized in (NZrNZj) rectangular pixels, where NZi and NZj are the number of discretization steps in the horizontal and vertical directions, respectively. The location of the source of contamination in this region can, then, be identified with corresponding elements of coordinates (Is, J~). The migration of the contaminant from the source to the surrounding zones is modelled by means of eight transition rates (Fig. 1): f and b are the probability densities of the contaminant making a forward and backward jump to the adjacent neighbouring pixel in the horizontal direction, respectively; u and d are the probability densities of the contaminant making an upward and downward jump to the adjacent neighbouring pixel in the vertical direction, respectively; fu = x/f 2 + u2 (forward-upward),
bu
\
U ,L
\/ b
f
".
/ bd
/
fa
,p
d
\ fd
Fig. 1. Possible moves of the contaminant.
Contaminant transport in bidimensional porous media
1307
f d = v/f 2 + d 2 (forward-downward), bu = ~ u2 (backward-upward), and bd --- V ~ + d 2 (backward-downward) are probability densities for diagonal moves. By a proper assignment of the values of these parameters one may also consider any kind of subdomain D within the above defined rectangular domain, simply by setting to zero the transition probabilities outside D. During its travel within the region, the contaminant may undergo physical-chemical processes of adsorption and desorption which typically slow down its transport. These phenomena can be modelled quantitatively by means of two additional transition rates, ads, for the adsorption process, and des, for the desorption process. The rates of movement f, b, u, d are characteristic of the flow regime which is ultimately determined by the physical properties of the host medium; similarly, the processes of adsorption and desorption are intrinsic characteristics of the solid matrix. Possible heterogeneities of the medium may be accounted for by considering these parameters as random variables in space, with assigned distributions. Finally, the contaminant may also decay, if radioactive. The evolution of the decay products, if still hazardous, can very well be followed simply by introducing new kinds of particles with their characteristic rates of movement within, and interaction with, the host matrix. A mission time, Tmiss, is defined as the period of interest for the groundwater transport of contaminants, and a given number NF of final locations of interest are selected within the region as observation points for the final amount of contaminant concentration at Tmiss. These locations could represent drinking water wells, public living areas, discharge points and the like.
4. MONTE CARLO SIMULATION: A NUMERICAL EXAMPLE We consider a rectangular region of space subdivided into NI.Nj = 40.20 = 800 identical, rectangular pixels. The bounds of this region are considered absorbing. A source of contaminant is located in the southwest corner (1,1), whereas ArE = 5 observation wells are located in the upper region of the east bound at locations (40,20), (40,18), (40,16), (40,14), (40,12), respectively (Fig. 2). To account for spatial heterogeneity, the transition rates f, b, u, d, ads, des are assumed to be all log-normally distributed with percentile parameters given in Table 1. The Monte Carlo trials are divided in bins and at the beginning of each bin a different set of values of these parameters is sampled from their respective distributions. For what concerns the nature of the contaminant, we have considered a radioactive element which decays with a time constant of 20 000 years in a non-hazardous daughter whose transport poses no additional risk and, thus, needs not be modelled further. The transport of the contaminant is followed over a period of 105 years, which is a typical regulatory time for high level radioactive waste repositories. Both analogue and biased simulations have been considered. For the biased Monte Carlo simulation a tailored technique has been proposed which favours contaminant moves directed towards one of the assigned observation locations. At the beginning of each trial, one of the NF observation locations is selected at random as the target pixel (IF, Jr) of contaminant accumulation for that trial.When the contaminant particle enters a defined area of interest around the target pixel we favour those transitions which bring the contaminant closer to the target pixel. The area of interest around the pixel (IF, Jr) is
1308
M. Marseguerra and E. Zio
J~
...............................
4"~0 4~ 4~,i6 4~,4
h 1
40
I
Fig. 2. The region in which contaminant transport occurs. Table 1. 5th and 95th percentiles of the log-normal distributions of the parameters f, b, u, d, ads, des Parameter f b u d ads des
P5
P95
5x10 -3 5x10 -4 5x10 -5 5.5x 10-4 5x 10-3 2x 10-4
8x10-3 8×10 -4
8x10 -5 8.5x 10 -4
8x 10-3 3x 10-4
a circle of radius r centred at the pixel: it is only within this region that biasing takes place (Fig. 3). By so doing, we bias only histories which already point in the right direction, towards (IF, JF), thus avoiding forcing indiscriminately all the histories, even those directed in wrong directions for which the biasing would require great amounts of time and significant reductions of the weights eventually collected. Furthermore, in order to have a certain control on the weights, we gradually increase the biasing within the region of influence by employing a smooth function 0(/) which starts from 00 on the circumference and increases towards the centre (IF, JF), up to a final value OF. More specifically, the current position (Ic, Jc) of the contaminant is compared to (IF, JF) and the transition rates of all those moves that bring the contaminant closer to the target are suitably increased. The amplification factor here adopted is: 0(0= l+a(0e-~
l
where I is the distance from the current to the final location and
(5)
Contaminant transport in bidimensional porous media
1309
J
/
Fig. 3. Circle of influence for biasing. ot = (OF -- 1).e~ _
-_____._0__0 (r
-,o,(V)
(6)
l = ~//(Ic -- IF) 2 Jr- ( J c -- J F ) 2
Fig. 4 shows a typical behaviour of 0(/). It is important to stress that the values of the functions ot and/~ in equation (6) should be chosen coherently with the choice of r, so as to guarantee a 'balanced' biasing. In other words, if r contains too many pixels, the successive multiplications of the weight by 0 can result in a vanishing final weight, a situation which should be avoided. In favouring certain moves, we obviously disfavour others which, when sampled, largely increase the weight. The presence of a small number of histories carrying large weights can significantly deteriorate the statistics of the results at later times. To circumvent this problem, and obtain statistically significant results throughout the mission time without increasing the number of trials beyond reasonable limits, we coupled our biasing technique with the splitting technique (Hughes et al., 1958): when the weights of the unfavoured histories exceed pre-cstablished thresholds, the history is split into two (or more) subhistories from that time point on, each one having half (or less) of the initial weight. The splitting thus generates a number of subhistories which distribute their contributions to reinforce the statistics.
5. RESULTS First an analogue simulation with l06 histories was performed, requiring a computer time of 1.72 min on an Alpha Station 250 4/266. Figure 5 reports the time evolution of the contaminant concentration in correspondence of the five target locations. As can be seen,
1310
M. Marseguerra and E. Zio 10
i
Amplification factor for biasing i
i
i
,
10
15
20
25
30
35
40
I
Fig. 4. Behaviour of the amplification factor 0 vs distance from the target pixel. the event of the contaminant travelling from the source in the south-west corner to the target point (40,20) at the opposite north-east corner is rather rare, whereas more probable events are those concerning the contaminant going to more central targets on the east bound, such as (40,14) and (40,12), for example. Contaminant concentration at the target locations: analo case 1.2 i 10.4 i(40,12)
0.8
o
0.6
0.4
0.2
1
2
3
4
5 time(y)
6
7
8
9
10 x 104
Fig. 5. Time evolution of the contaminant concentration at the target locations: analogue simulation.
Contaminant transport in bidimensional porous media
1311
This behaviour is confirmed, as expected, by the biased simulation. Again, 106 histories were performed, requiring 6.33 min of computer time on the same machine. The radius of the circle of influence of each target pixel was chosen to be equal to 40, with a starting bias amplification factor 00 --- 2 on the circumference rising up to a value OF = 10 at the pixel. Splitting is set up in such a way as to occur when the weight carried by the contaminant particle in the current history goes beyond a threshold value of 1.5. This has generated 2 026246 additional subhistories which have allowed us to keep the magnitude of the collected weights to an average of 4.75x 10 -4 with a standard deviation of 5.30x 10-6. This ensures good statistical results, which are shown in Fig. 6: the greater smoothness of the curves is an indication of better statistical accuracy as compared to the analogue case. The efficiency of the biasing obviously depends on how rare the events of interest are: if they are reasonably probable, then biasing is likely to be unnecessary; on the other hand, if they are quite rare, then, biasing can lead to a significant improvement in the accuracy of the estimates. In our case, the move probabilities are such that arrival at the target pixels is less and less probable as we go up from the target in (40,12) to that in the northeastern corner (40,20). In Table 2 we compare the analogue and biased cases on the basis of the usual figure of merit given by the inverse of the product of the variance times the computer time required. In the present case, each Monte Carlo game has five different target pixels, so that the figure 1 The of merit of the computation is actually a vector M with five components Mr,-= (~)itc," computer time tci appearing in the ith component of M results from the assumption that the total computer time is equally shared among all target pixels. For each time instant, the standard deviation (trc)i appearing in the ith component of M gives a measure of the accuracy of the estimate of the contaminant concentration, C, in the ith target location. For t * = 2 • l04 years, the vector of figures of merit M is reported in the table, each component being referred to a different pixel. It can be seen that the biasing is not efficient Contaminant concentration at the target locations: biased case
x 10 4 1.2 ~
~
T
T
r
r
r
T
6
7
8
9
(40,12) 1
0.8
0.6
0.4
0.2
1
2
4
5
time(y)
10 x 10~
Fig. 6. Time evolution of the contaminant concentration at the target locations: biased simulation.
1312
M. Marseguerra and E. Zio
Table 2. Figures ofmerit for the estimates ofconcentration in the various pixels at t" = 2 x 104 years; tc, = computer time (min) for the ith target pixel (total computer time divided by the number of target pixels); (trc)i = standard deviation of the concentration in the ith target pixel at t*; M i = I Target (IF, JF) (40,12) (40,14) (40,16) (40,18) (40,20)
(trc)i analogue (to, = 0.34 min)
(ac)i biased (tc~ = 1.27 min)
Mi analogue (t~, = 0.34 rain)
Mi biased (t~, = 1.27 rain)
3.03× 10 - 6
1.81 × 10-6
2.13x 10-6 1.40× 10-6
3.09x 10 - 6
2.17× 1011 6.40× 10I1 1.48× 1012 3.52× 1012 0
2.41x 10I1 8.25x 1011 1.25x 1013 5.31X 1013
2.52x 10-6 1.22x 10 -7 4.95x 10-s
8.56x 10 - 7
cx~
3.23x 1014
for the relatively probable event of arrival in (40,12) but it becomes more and more efficient, with respect to the analogue simulation, as we go towards the more distant targets. In particular, no contribution is carried by the analogue simulation to the farthest target (40,20), thus giving rise to a null figure of merit. As expected, then, the efficiency of the biasing depends on the probability of arriving at the target pixels and, therefore, on the values of the transition rates of movement between adjacent pixels. The less probable are the movements leading in the direction of the targets, the more efficient is the biasing. A global comparison between the analogue and the biased computations can be done in terms of the average of the components of the vector M. In the present case we have Manalog =
1.17 × 10 j2
Mbiased = 7.78 X
1013
which correspond to a mean gain of 66 in favour of the biased simulation. The higher accuracy achieved by the biasing is also shown in the comparison of Figs 7 and 8, which report the distribution of arrival times at the target location (40,12) in the analogue and biased case, respectively, and associated error bands. For this quantity, Table 3 reports the comparison in terms of the standard figure of merit. At each time, the accuracy of the results is quantified in terms of the standard deviation % of the distribution of arrival times, p, in (40, 12). For t = 105 years, the figure of merit M = ~ shows I P~ that, in this case, the biased simulat'on pe r fo r n l s almost six times better than the analogue one. In the above example, in correspondence of t = Tmiss the analogue histories carry contributions to the distribution only in the target pixel (40,12), whereas no contribution Table 3. Comparison between the analogue and biased simulations: to, = computer time (min) for the target pixel (40,12) (total computer time divided by the number of target pixels); trq = standard deviation of the probability of arrival in (40,12) at Tmiss= 105 years; M =
Analogue Biased
tc, (min)
ap
M
0.34 1.27
9.95 × 10 - 7 2.17× 10 -7
2.94× 1012 1.67× 1013
Contaminant transport in bidimensional porous media
13 l 3
Distribution of arrival times at the location (40,12): analog case i i + + i i
10 -s
1.8.
/
1.6 ................................
t
...................................
......
1.4
1.2
1 0" 0.8
0.6
0.4
0.2
0
1
2
3
4
7
8
9
10
time(y)
X 10 4
Fig. 7. Distribution of arrival times at location (40, 12) analogue simulation. arrives in the other targets. The latter situation corresponds to a value of zero in the associated components of the vector of figure of merit. On the contrary, the biased histories are such as to provide contributions to all target pixels, i.e. to all components of the vector of figure of merit. To further investigate this aspect, we have considered a second case study concerning the transport of Tc-99 (2 = 3.24 x 10 -6 years -1) in the same spatial region. The properties
Distribution of arrival times at the location (40,12): biased case ! ! ! ~ ! ! +
x 10 -s
1.8
i
:
i
:
i
i
1.4 . . . . . . . . . . . . .
i ................
:. . . . . . . . . . .
1.2 ..............
: : .................
:
i :.......
:
: .........
:. . . . . .
i : .......
:. . . . . .
.......
: 1 ............................
!
....
: . . . .
.........
.........
!
:..............................
: ..........
i ..........
3
6
7
0"0. 8
0.6
0.4
0.2
6'
1
2
4
5 time(y)
8
9
10 x 10 4
Fig. 8. Distribution of arrival times at location (40,12): biased simulation.
M . Marseguerra and E. Zio
1314 x 10-r
Contaminant concentration at (40,12): analog case
6
5
4
o3
2
1
I!111II II Illilllllllll
I:111IIItIIl:lllllll[ll: 0
li
~
f
2
3
4
. ~ 5 time(y)
9
10 x 104
Fig. 9. Time evolution of the contaminant concentration in (40,12): analogue simulation.
of the host medium are not varied except for the percentiles of the distribution of the downward rate, d, which are an order of magnitude larger than those of Table 1. This renders the event of contaminant accumulation in the north-eastern target points much less probable. No modification was made to the parameters for the biased simulation. Figures 9 and 10 report the time evolution of the contaminant concentration in the target location (40,12) as resulting from the analogue and biased simulations, respectively. x 10-r
Contaminant concentration at (40,12): biased case
v
i i
3 ........ !....... i .........
0
0
2
:~ .........
i..........
3
4
i........
......
6
time(y)
i ..........
7
i .........
-
! .......
9
10 x 104
Fig. 10. Time evolution of the contaminant concentration in (40,12): biased simulation.
Contaminant transport in bidimensional porous media
1315
Table 4. Comparison between the analogue and biased simulations: to, = computer time (min) for the target pixel (40,12) (total computer time divided by the number of target pixeh); a e = standard deviation of the probability of arrival in (40,12) at Tmiss = 105 years; M = to,
Analogue Biased
(min)
1.20× 102 7.7
trp
M
8.29x 10-8 4.62× 10-8
1.21 x 1012 6.08× 1013
It can be noted how the arrival of contaminant at this location is almost three orders of magnitude less probable than in the first case examined. The analogue case required 108 histories and a computer time of 974min to provide results of any statistical significance, whereas 106 histories were performed for the biased case which, however, benefited from the contributions of 7 896 086 split subhistories. Overall, the computer time required in the biased case was 38.46 min and 191 248 weights arrived at a destination in one of the target points, each weight bringing, on average, a contribution of 1.25x 10-6 with a standard deviation of 6.81 x 10-8. Table 4 reports the details of the comparison between the analogue and biased cases: the values of the figure of merit M show that in this case of very rare events the biased simulation provides a performance which is 50 times better than that of the analogue simulation.
6. CONCLUSIONS Many environmental management and technical evaluation activities rely on the extensive use of predictive computer models to describe the phenomena occurring over the time frames required by the regulations. Natural groundwater transport is one such phenomenon which is of particular interest because groundwater is a principal vector through which hazardous contaminant can circulate to the biosphere and, thus, pose a significant risk to the health of the human population and to the safety of the environment. Several approaches to modelling groundwater contaminant transport have been proposed in the past. In this paper we have focused on a relatively new model, previously proposed by the authors, which is based on the theory of branching stochastic processes of Kolmogorov and Dmitriev. The flexible structure offered by the model allows one to describe explicitly many of the non-linear phenomena which actually occur in practice. The degree of detail which can be reached in the analytical model is, however, limited by the high order of the system of coupled differential equations to be solved. To circumvent this problem, in this paper we have proposed to structure the Kolmogorov and Dmitriev approach within a pure Monte Carlo framework. The great advantage is that with the current power of modern computer systems, there are basically very few limitations on the number of particle kinds and interactions that one can simulate. In addition, the proper use of variance reduction techniques can further increase the efficiency of the calculation in terms of accuracy. In this context, a biasing technique has been introduced, which forces the contaminant to travel towards the accumulation points of interest within the region under study. This is coupled with a splitting technique to ensure a proper distribution of the contributions to the statistics of the results throughout the
1316
M. Marseguerra and E. Zio
mission time. The efficiency of this technique has been proven through a sample application on a bidimensional region. Finally,the additional flexibility introduced by the Monte Carlo scheme has allowed us to accommodate elements of stochasticity in the physical properties of the media to account for possibly existing heterogeneities.
REFERENCES
Antonopoulos-Domis, M., Clouvas, A. and Marseguerra, M. (1995) On the compartmental modeling of cesium migration in soils. Nucl. Sci. & Eng. 121,461. Gelhar, L. W., Welty, C. and Rehfeldt, K. R. (1992) A critical review of data on fieldscale dispersion in aquifers. Water Resour. Res. 28(7), 1955. Hughes, D. J., Sanders, J. E. and Horowitz, J. (1958) Physics and mathematics. Progress in Nuclear Energy 2, 315-368. Kolmogorov, A. N. and Dmitriev, N. A. (1947) Dokl Akad. Nank 5553 56, 7. Lee, Y. and Lee, K. J. (1995) Nuclide transport of decay chain in the fractured rock medium: a model using continuous time Markov process. Ann. Nucl. Energy 22(2), 71. Marseguerra, M. and Zio, E. (1997) Modelling the transport of contaminants in groundwater as a branching stochastic process. Ann. Nucl. Energy 24(8), 625-644. Marsily, G. de (1986) Quantitative Hydrogeology. Academic Press, San Diego, California. Prickett, T. A., Naymik, T. G. and Lonnguist, C. G. (1981) A 'random-walk' solute transport model for selected groundwater quality evaluations. ISWS/BUL-65/81, Bulletin 65, State of Illinois, Illinois Department of Energy and Natural Resources, Champaign. Williams, M. M. R. (1992) A new model for describing the transport of radionuclides through fractured rock. Ann. Nucl. Energy 19(10), 791.