On the sensitivity of a marine dispersion model to parameters describing the transfers of radionuclides between the liquid and solid phases

On the sensitivity of a marine dispersion model to parameters describing the transfers of radionuclides between the liquid and solid phases R. Peria ´n˜ez  Departamento de Fı ´sica Aplicada I, E.U. Ingenierı ´aT´ecnica Agrı ´cola, Universidad de Sevilla. Ctra. Utrera km 1, 41013 Sevilla, Spain Abstract The sensitivity of a marine dispersion model for non-conservative radionuclides, previously developed and validated for the English Channel, to parameters describing the exchanges between the liquid and solid phases (suspended matter and bottom sediments) has been studied using a Monte Carlo method. A probability distribution is assigned to each parameter. They are sampled to obtain a set of model parameters and a model run is carried out. This process is repeated to obtain a distribution of model outputs. Partial correlation coefficients are calculated to assess the relative influence of each parameter on model output. Errors are also assigned to model results. Three situations are studied: an instantaneous release of radionuclides, a continuous release and the case of a contaminated sediment behaving as a long-term source of radionuclides. Calculations have also been carried out for two radionuclides with different geochemical behaviour: 137 Cs and 239,240 Pu. The results indicate that all parameters are relevant, depending on the phase we are interested in obtaining the result and on the source term (instantaneous, continuous or due to sediments). However, parameters that are, in general, more influential are kinetic rates, mixing depth in the sediment and mean radius of suspended and sediment particles. This suggests that including several particle sizes in future radionuclide dispersion models could lead to an improvement in model results. Differences have also been found with respect to the relevance of some parameters depending on the geochemical behaviour of the radionuclide.  Tel.: +34-954-486474; fax: +34-954-486436. E-mail address: [email protected] (R. Peria ´n˜ez). Keywords: Monte Carlo; Modelling; Sensitivity; Radionuclides; Partial correlations; Model uncertainty 1. Introduction Numerical models have been widely used to assess the effect of radioactivity releases to the environment. Model predictions may be used as the basis for important decisions on issues such as emergency response and countermeasures, waste management, and environmental remediation and restoration. The consequences of such decisions may be significant in terms of human, ecological or economic costs (Thiessen et al., 1999). Thus, it is essential to evaluate the reliability of the model predictions. Model uncertainties have been typically investigated through simple sensitivity studies in which each of the parameters in the model was varied in turn, leaving the others fixed at their nominal values (Peria ´n˜ez et al., 1996; Tappin et al., 1997; Margvelashvily et al., 1999). Indeed, sensitivity of a model to a parameter that has a high natural variability or about which little is known can increase the uncertainty associated with model predictions. Now, the Monte Carlo method has proven to be generally a more suitable approach in carrying out a quantitative analysis of sensitivity and propagation of uncertainties in radionuclide transport codes (Whicker et al., 1999). It is also possible to calculate partial correlation coefficients (r part ) for relationships between model outputs and values of the parameters selected for analysis. A partial correlation estimates the linear correlation between an output variable and a parameter after removing the effects of the other parameters. A rank of model parameters can be obtained according to their partial correlation coefficients, which gives an indication of their relative influence on model output (Whicker et al., 1999). The effort for obtaining site-specific values for model parameters, for a given model application, should be focused on the most influential ones. A parameter with a high r part may deserve more study to improve confidence in it and perhaps to reduce its uncertainty. In contrast, research on a parameter with a small r part may not be useful because of the low relationship between its value and the output. The objective of this work consists of studying the sensitivity of a marine dispersion model for non-conservative radionuclides previously developed and validated for the English Channel (where radionuclides are released from Cap de La Hague nuclear fuel reprocessing plant). In particular, only the model sensitivity to parameters governing the exchanges of radionuclides between the liquid and solid phases (suspended matter and bottom sediments) will be studied. Model sensitivity to parameters related to tidal propagation and sediment dynamics is not included in this work. Model sensitivity will be studied in three typical situations, which correspond to the different source terms that can generally be found: (1) the case of a hypothetical instantaneous release of radionuclides from the reprocessing plant, (2) the case of a continuous release and (3) the case in which there is no external input of radionuclides but the bottom sediments are initially contaminated, so that they behave as a long-term source of previously released waste radionuclides. Errors will be assigned to model predictions and r part s will be calculated for the three situations studied. This will allow establishing a ranking of partial correlation coefficients to analyze the relative influence of each parameter. Calculations are made for two radionuclides with different geochemical behaviour: 137 Cs and 239,240 Pu. This type of sensitivity analysis has been carried out in box models. For instance, Breshears et al. (1992) present such a study for an agroecosystem model that represents processes affecting milk contamination after a fallout deposition of radionuclides. Carroll and Harms (1999) studied uncertainty in a model that simulates the dispersion of radionuclides in a shallow Arctic Bay. The model uses partition coefficients to describe the distribution of radionuclides between water and suspended matter particles. Probability distributions were assigned to suspended matter concentration and distribution coefficients k d s (parameters needed to calculate the partition coefficient). These distributions were sampled and probability distributions were then obtained for the partition coefficient. The dispersion of radionuclides was then computed for the partition coefficient probability distribution mode and for the lowest and the highest values. Thus, in this work, probability distributions are used for the partition coefficients, although the sensitivity analysis is carried out in the simple way mentioned above. This is the first time, to the author’s knowledge, that a detailed Monte Carlo-based sensitivity study is carried out on a marine dispersion model for non-conservative radionuclides. The model and method used for the study are described in the next section. Following this, results are presented and discussed. 2. Materials and methods 2.1. The model The model used for the study has been described and tested before (Peria ´n˜ez, 2003), so few details are given here. It is a long-term dispersion model for non-con- servative radionuclides. The model is based upon residual (averaged) circulation obtained from a previously validated tidal model of the Channel (Peria ´n˜ez and Reguera, 1999). The averaged suspended matter concentrations over the Channel have also been obtained from a sediment model which includes advection/diffusion of particles, settling, deposition and erosion of the sediment. Details can be seen in Peria ´n˜ez (2000). The model domain, showing the residual circulation for average wind speed and direction (southwest 6 m/s), is presented in Fig. 1. The position of Cap de La Hague nuclear fuel reprocessing plant is also indicated in the map. The measured suspended matter concentrations in the Channel (Eisma and Kalf, 1987) range from some 3 ppm along the French shore to 5 ppm along the British shore, although concentrations up to 10 ppm can be found in some locations. The model results are in agreement with these values (Peria ´n˜ez, 2000). Since the Channel is a highly dynamic environment in terms of tides, the sediments are essentially composed of coarse material (Boust, 1999). Mud deposits can only be found along both shores in the areas of weaker tidal currents (amplitudes of the order of 0.5 m/s or less). The model considers that radionuclides can be present in three phases: solution, suspended matter and active bottom sediments (particles with a diameter <62.5 lm). It is considered that the exchange of radionuclides between the liquid and solid phases is governed by a single reversible reaction. Thus, transfer of radionuclides from water to the solid phase is governed by a coefficient k 1 and the inverse process by a coefficient k 2 , which are denoted kinetic transfer coefficients (dimensions [T] 1 ). The adsorption of radionuclides depends on the surface of particles per water volume unit, denoted as the exchange surface (Peria ´n˜ez, 1999, 2000, 2002a): k1¼v1ðSmþSsÞð1Þ where S m and S s are the exchange surfaces, per unit water volume (dimensions [L] 1 ), for suspended matter and bottom sediments, respectively, and v 1 is a parameter with the dimensions of a velocity denoted as the exchange velocity. Assuming spherical particles and a step function for the size distribution of suspended particles, it can be easily obtained (Peria ´n˜ez, 1999, 2000, 2002a) that: Sm¼3m qRð2Þ Ss¼3Lf / RH ð3Þ where mis the suspended matter concentration, qis the suspended particle density, Ris the mean radius of suspended matter and active sediment particles, Lis the average mixing depth (the distance to which the dissolved phase penetrates the sediment), fgives the fraction of active sediment, His the water depth and /is a 0 5 10 15 20 25 30 35 40 45 50 010 20 30 40 50 60 70 10 cm/s ENGLAND FRANCE Cotenin Dover Strait La Hague Fig. 1. Residual currents in the Channel for average wind conditions. Each unit in the xand yaxes is 5000 m. correction factor that takes into account that not all the surface of the sediment particles is in contact with water since it will be partially hidden by other particles. This formulation has been successfully applied in the simulation of the dispersion of non-conservative radionuclides in estuarine (Peria ´n˜ez et al., 1996, 2002a) and marine environments (Peria ´n˜ez, 1999, 2000). Three partial differential equations are then obtained (Peria ´n˜ez, 2003), whose solutions give the time evolution of specific activities in water, suspended matter and active sediment. In the case of water, C d (Bq/m 3 ), is: @Cd @t¼ðadv þdifÞk1Cdþk2Csmþk2 AsLqsf/ Hð4Þ where C s and A s are specific activity (Bq/kg) in suspended matter and the active bottom sediment, respectively. Since the model is two-dimensional, it is considered that C d and C s have a uniform distribution over the water column. The sediment bulk density, q s , is expressed in kg/m 3 and (adv þdif) means advective plus diffusive transport of radionuclides. Kinetic coefficient k 1 is obtained from Eqs. (1)–(3); thus it can be seen that k 1 depends on the suspended matter concentration and the distribution of active sediments over the model domain. On the other hand, k 2 is considered constant since the desorption process does not depend upon the available particle surface, as is the case with adsorption. The equation for the temporal evolution of specific activity in suspended matter particles is: @ðmCsÞ @t¼ðadv þdifÞþk1Cdk2CsmþEDmð5Þ where k 1 is now given by the first term of Eq. (1). Symbol ED m is the erosion– deposition term. Finally, the equation for the temporal evolution of specific activity in the active bottom sediment is: @As @t¼k1 CdH Lqsfk2As/þEDsð6Þ where k 1 is given by the second term of Eq. (1). The form of the erosion–deposition terms for suspended matter and bottom sediments, ED m and ED s , can be seen in Peria ´n˜ez (2003). For the case study (3), redissolution from a contaminated sediment, the kinetic model is replaced by a two-step model consisting of two consecutive reversible reactions. The first, governed by the same coefficients k 1 and k 2 , is a reversible adsorption or isotopic exchange with non-specific sites on particle surfaces. The second and slower reaction, governed by rates k 3 and k 4 (forward and backward, respectively), is a reversible transfer of radionuclides to more specific sites in the particles. The second process may be a slow diffusion of ions into pore and interlattice spacings, inner complex formation or a transformation, such an oxidation. Thus, in the two-step model, the solid phase is divided into two phases, which are denoted as non-specific and specific. This has been done since very recent experiments (see for instance Ciffroy et al., 2001) have shown that a two-step model is more adequate than a one-step model to simulate both the uptake and release kinetics. On the other hand, such a two-step model has been implemented in a plutonium dispersion model of the Irish Sea (Peria ´n˜ez, 2002b). Results indicate that there are no appreciable differences between both models when simulating the dispersion of radionuclides introduced from an external source of radionuclides over relatively short time scales (several months), but a two-step model should be used if the source is due to redissolution from contaminated sediments since in this case a one-step model produces a too fast decrease in sediment activity if compared with experimental evidences (see Peria ´n˜ez, 2002b). Thus, for cases (1) and (2), a one-step model is used since it is computationally cheaper, and the two-step model is used for case (3). To solve the equations, a spatial and temporal discretization is carried out: the Channel was divided into 3750 grid cells (forming a matrix 75 50). The grid extends from 4.0vW to 1.5vE and from 48.3vto 51.0vN (see Fig. 1). The grid cell size is Dx¼Dy¼5000 m (xand ymeasured eastward and northward, respectively) and time step is Dt¼3600 s. The equations are solved explicitly using second order accuracy finite difference schemes. 2.2. Method for the sensitivity analysis The sensitivity analysis has been limited to parameters governing the exchanges of radionuclides between the liquid and solid phases: v 1 ,k 2 ,q,q s ,R,L,/and f,as well as k 3 and k 4 in the case of redissolution. Other parameters concerning tidal propagation or sediment dynamics have not been included in this study since we are mainly interested in the description of the transfers between phases as commented above. A Gauss probability distribution is assigned to each parameter, which is defined by its average value and standard deviation. A random sampling of each distribution is carried out using a Monte Carlo method, so that a set of model parameters is obtained, and the model is run. This process is repeated 200 times, so that a distribution of model results is obtained. The set of model results allows assigning uncertainties to model output and calculating partial correlation coefficients. Partial correlation coefficients are calculated in the following way: if a set of variables ðx1;x2;...;xpÞis given, the r part between variables 1 and 2, for instance, is calculated obtaining the correlation coefficient in the regression: e1;3;4;...;p¼Ae2;3;4;...;pþBð7Þ where e1;3;4;...;pand e2;3;4;...;pare the residues of the multiple regression of variables 1 and 2, respectively, with respect to the rest of variables 3;...;p. The correlation coefficient for regression Eq. (7) is the partial correlation coefficient between variables 1 and 2. In the three cases analyzed, the model can produce an extensive amount of information: distribution maps of radionuclide concentrations for each phase and at any desired time, and/or temporal evolution of concentration in each phase and for each desired point of the model domain. Considering that 200 simulations are carried out for each case, the volume of data obtained would be rather difficult to handle. Thus, to simplify the model output, we have fixed our attention in a given target grid box and at a defined instant of time after the beginning of the simulation (t¼0). This target means a sensitive point in the domain, where a local population (human or not) could receive an external dose due to contamination of water and/or sediments or an internal dose due to ingestion of contaminated material, for instance. The instant of time at which results are obtained is the time required for the peak activity to travel from La Hague to the target grid box. Thus, the model would provide information on maximum activity levels reached in that point (in water, suspended matter and bottom sediments) so as to determine, for instance, if any kind of remedial action should be taken. For simulations carried out with 137 Cs, we have defined grid box (50,20) as the target and the results are obtained 40 days after the beginning of the release from La Hague or the redissolution from the sediment. In the case of 239,240 Pu, target grid box is (30,26). This grid box is closer to La Hague since due to the high reactivity of plutonium, this radionuclide moves slower than Cs. 3. Results and discussion The values given to the parameters are presented in Table 1 for the 137 Cs application. Each value, which is the nominal value previously used in this model (Peri- a ´n˜ez, 2003), is taken as the average value for each Gauss distribution. The errors given in Table 1 (arbitrarily taken as 20% of the nominal values) correspond to the standard deviation of each probability distribution. The probability distribution of each parameter is sampled using a Monte Carlo method over a 3rdistance around the average value. Details on the selection of the nominal values for the parameters can be seen in Peria ´n˜ez (2003) and will not be repeated here. Initial conditions for the simulations are the following. In cases 1 and 2, the model is started from zero concentrations in water, suspended matter and bottom sediments over all the domain. In the case of an instantaneous release, an input of Table 1 Nominal value and error assigned to each parameter in the model for the application to 137 Cs Symbol Description Value qParticle density 2600 520 kg=m3 q s Sediment bulk density 900 180 kg=m3 v 1 Exchange velocity ð2:10:4Þ108m=s k 2 Kinetic coefficient ð1:16 0:23Þ105s1 LSediment mixing depth 0:10 0:02 m /Correction factor 0:10 0:02 RParticle mean radius 15 3lm fActive sediment fraction 0:10 0:02 k 3 Kinetic coefficient ð1:16 0:23Þ107s1 k 4 Kinetic coefficient ð1:16 0:23Þ108s1 5:01012 Bq of dissolved 137 Cs is assumed at t¼0 and at grid box (26,24), where wastes from La Hague are released. In the case of a continuous release, 5:0 1010 Bq are introduced each time step. In study 3, (redissolution) an initial specific activity of 1:0105Bq=kg was considered in the active fraction of the bottom sediment in grid box (26,24) (in the specific-site sediment phase), which means an inventory of 2:25 1013 Bq in the entire compartment. Zero concentrations are assumed for water, suspended matter and sediment (unless in the mentioned point) over all the domain. As an example, the distribution of dissolved 137 Cs over the model domain 40 days after the instantaneous release from La Hague is presented in Fig. 2. The position of the target grid box is also shown in the map. It can be seen that the patch of radionuclides essentially moves along the French shore and that the target grid box is effectively in the area of the activity peak at this instant of time. In the case of a continuous release, a banded structure showing higher concentrations along the French shore is now extending from La Hague to Dover (see Peria ´n˜ez, 2003). In the case of redissolution, concentrations at La Hague in water and suspended matter initially increase due to the input from the sediment and after decrease due to advective/diffusive processes. A patch extending along the French shore is then observed. The sediments are contaminated as the patch travels above them. The frequency histograms of results have been obtained for C d ,C s and A s in the 137 Cs instantaneous release experiment as an example. Such histograms have been fitted to Gauss distributions. The results are presented in Table 2, where the average value and width of the distributions are given together with the correlation coefficient of the numerical fitting. It can be seen that, in general, the fitting of the histograms to the distributions is good: Gauss distributions are obtained for the model results from the probability distributions of parameters. Uncertainty can be assigned to model results now. Thus, it can be concluded that for the target grid 10 20 30 40 50 60 70 10 20 30 40 50 France England La Hague Fig. 2. Distribution of dissolved 137 Cs (Bq/m 3 ) 40 days after an instantaneous release at La Hague. The position of the target grid box is indicated by the triangle. box 40 days after the discharge, 137 Cs concentrations are 35:41:5Bq=m3,5:3 2:3Bq=kg and 11 4Bq=kg for water, suspended matter and active bottom sediment, respectively (errors are 1r). The model results are more precise in the case of the dissolved phase (relative error 4%), while uncertainties are larger for the suspended matter and bottom sediments (relative errors 43% and 36%, respectively). This method for estimating model uncertainties could be applied to any other modelling application. Absolute values of the partial correlation coefficients for the three cases studied and for each parameter are presented in Fig. 3 in the case of 137 Cs. The most influential parameters are arbitrarily defined as those having r part above 0.5 (Breshears et al., 1992). Moderately influential parameters are defined as those having r part between 0.3 and 0.5. It can be seen that all parameters are relevant in the description of the transfers between the liquid and solid phases, depending on the phase in which we are interested and on the source term. For instance, if we are interested in obtaining specific activity in suspended matter particles (due to any reason), it is clear that particle density, q, must be precisely specified for the three situations investigated (rpart >0:77). This parameter, however, is not relevant if our interest is focused in water or bottom sediments (rpart <0:16 in all cases). Sediment bulk density, q s , is relevant for A s in the three situations (rpart >0:79). It is also relevant for the dissolved phase in the case of redissolution (rpart ¼0:84), but is only moderately relevant (rpart ¼0:47) for suspended matter in this case. It may be due to the fact that transfer to suspended matter from the sediment is indirect, the dissolved phase being an intermediate step. In cases of instantaneous or continuous releases, q s is not relevant either for the suspended matter or for the dissolved phase (rpart <0:08). The exchange velocity v 1 and the kinetic coefficient k 2 are always relevant parameters, for all situations and phases (rpart >0:59). The only exception is v 1 in the dissolved phase for the redissolution experiment, where the influence is moderate (rpart ¼0:42). It seems that the activity content in water is essentially more dominated by the transfer from the sediment to water (desorption is dominating) at La Hague than by the subsequent transfer to suspended matter and bottom sediments that is produced as the patch of redissolved radionuclides travels along the Channel. The mixing depth in the sediment, L, is relevant for the dissolved phase in the three cases (rpart >0:75) since the kinetic coefficient k 1 depends on L(see Eq. (3)) Table 2 Average value ( xx), width (w) and correlation coefficient for the numerical fitting of model results to Gauss distributions for each phase in the 137 Cs instantaneous release experiment. The width of the distribution is w¼2r Phase  xxwr 2 C d (Bq/m 3 ) 35.39 3.04 0.949 C s (Bq/kg) 5.28 4.60 0.829 A s (Bq/kg) 11.12 7.97 0.680