Modelling of U, Th, Ra and 137Cs radionuclides behaviour in rivers. Comparison with field observations

Applied Mathematical Modelling 25 (2000) 57±77

www.elsevier.nl/locate/apm

Modelling of U, Th, Ra and 137Cs radionuclides behaviour in rivers. Comparison with ®eld observations M.J. Rodriguez-Alvarez a, F. Sanchez

b,*

a b

Dpto. Matem atica Aplicada, Univ. Polit ecnica de Valencia, Camino de Vera s/n, E-46020 Valencia, Spain Instituto de Fisica Corpuscular (I.F.I.C.), C.S.I.C.-Universidad de Valencia, Edi®cio Institutos de Paterna, Aptdo. 2085, E-46071 Valencia, Spain Received 30 November 1999; received in revised form 10 April 2000; accepted 30 May 2000

Abstract The aim of this work is the description and modelling of the major processes which govern the dynamic behaviour of radionuclides in rivers. A multi-compartmental box model of the river has been developed in order to consider not only hydrodynamic equations but also other complex phenomena as adsorption/desorption of radionuclides to sediments. The model has been applied to the Jucar river (Spain) for U, Th, Ra and 137 Cs radionuclides. The results of the model are compared with ®eld measurements for the case of naturally occurring radionuclides. Ó 2000 Elsevier Science Inc. All rights reserved. Keywords: River modelling; Environmental radionuclides; Simulation nuclear plant accident

1. Introduction Naturally occurring radionuclides in surface waters are not transported by the water in a simple way. Other complex phenomena account for the global behaviour of such radionuclides, mainly due to processes which intercommunicate sediment, suspended matter and water phase. The aim of this work is to develop a model that describes the dynamic behaviour of these radionuclides in rivers. Previously [1±3], uranium, thorium and radium activity measurements were done by the authors along the Jucar river (one of the main rivers in the eastern Spain, 500 km long). This will allow us to verify the model results by comparing with these available data. Moreover, distribution factors calculated with these measurements [1,2] are needed to properly develop our model. Beside the Jucar river is located the Cofrentes Nuclear Power Plant, the Power Plant outlet being placed just in the Jucar river. However, no uranium, thorium or radium release have been detected during these radionuclides measurements on water, sediment and suspended matter samples collected very near the Cofrentes Nuclear Plant outlet [1±3]. Thus, it can be assumed that the radioactivity levels of these isotopes will depend mainly on the geological substratum and

*

Corresponding author. Tel.: +34-96-398-3501; fax: +34-96-398-3488. E-mail address: ®lomeno@i®c.uv.es (F. Sanchez).

0307-904X/00/$ - see front matter Ó 2000 Elsevier Science Inc. All rights reserved. PII: S 0 3 0 7 - 9 0 4 X ( 0 0 ) 0 0 0 3 5 - 4

58

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

dynamics of the river. Finally we will apply our model to a nearly space-time point source of 137 Cs (1:3  1010 Bq) in water due to a hypothetical Cofrentes Nuclear Plant accident. The parameters which will determine the water depth all along the river are the climatology together with the water outlets, drains, irrigation canals and dams located along the river. We have as input data the river ¯ow and the water depth for several gauging stations placed at the river. In our case when trying to obtain the equations that describe our system we do not take into account the periodic tides. This phenomenon mainly determines velocities and depth of river water [4] in the mouth proximity. However for our system (Jucar river ¯ows into Mediterranean Sea), these tides are negligible in the river mouth vicinity. 2. The model The main di€erence between our model and those developed by other authors [5±8], is that we do not expect to have man-made radionuclides discharge e€ects. Thus, we will only consider naturally occurring radionuclides (uranium, thorium and radium) originating in river sediments itself. In this way, terms which are normally neglected in the other models (which could be called in these models as ``background'' terms) become dominant terms in our model, while the negligible terms for our model will be those named as external contributions to radionuclides activity. We will apply our model to the Jucar river (from Alarc on dam up to the river mouth on the Mediterranean Sea). This implies (Fig. 1) about 70% of the whole river length. Several simpli®cations are done in our model. As is usual in this kind of modelling [8±10], we adopt a onedimensional description of the river, by integration of the equations over the river cross-section averaged values. In this case the equations are simpler, the number of parameters needed by the model being smaller. The river was divided into nine sections in correspondence with the nine gauging stations located along the river. These sections were sub-divided into cells each 5.5 km long, i.e., the spatial resolution of our model will be 5.5 km. A mean water depth is assumed for the dams located on the river. Mean water depth for dams is obtained by dividing total water dam volume by its total surface. As water velocity in dams will be nearly zero, we will assume that sedimentation processes of suspended matter will be predominant inside such dams. Consequently in our model, suspended matter concentration will be nearly zero inside dams. Starting from the conservation ¯uid equations [8], we will try to describe water, sediments and radionuclides motion by considering that there is matter exchange between adjoining cells. Water motion is described by means of hydraulic and hydrodynamic variables of the river, in the socalled hydrodynamic submodel. Starting from data obtained by the hydrodynamic submodel, we can calculate suspended matter and sediments concentrations for each individual cell and their variation in time (sediment submodel). Finally, radionuclides activities in water, sediments and suspended matter and their time variation along the river can be modelled by using data coming from the hydrodynamic and sediment submodels. 2.1. The hydrodynamic submodel For this submodel, we need as input the data that describe the river geometry: river width and altitude of the river bed over the sea level for each one of the cells we have divided the river. The equation of continuity we have to solve is oq ~ ˆ Q; ‡ div…qU† ot

…1†

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

59

Fig. 1. Jucar river. Schematic representation.

together with the equation of motion ~ ~ o…qU† ~
grad†…qU† ~ ˆ ÿg grad‰q…Z ‡ D†Š ÿ gqU ~ jUj ; ‡ …U …2† ot C2Z ~ is the water where div and grad stand for divergence and gradient, respectively, and (Fig. 2) U ~ velocity, jUj the water velocity modulus (m/s), Z the water depth averaged over cross- and lengthsection (m), q the cell water volume (m3 ), Q the net tributary discharge (m3 /s), g the acceleration of gravity (9.8 m/s2 ), C the Chezy coecient ([C 2 ] ˆ m/s2 ), t the time (s) and D is the altitude of the river bed over sea level (m). We have veri®ed with the data of the gauging stations placed along the Jucar river that river ¯ow and water depth are nearly constant in time (on a six-month basis), if we except very short and sporadic torrential rain periods. In this sense, as we have not the e€ect of periodic tides, we have solved our equations without considering the time dependence of these equations. That is, if we impose to our model stationary solutions from the hydrodynamic point of view, Eq. (1) becomes ~ ˆ Q: div…qU†

…3†

60

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

Fig. 2. Schematic representation of river geometry.

~ has only As we have adopted a one-dimensional description of the river we assume that U component in the x-direction (i.e., we adopt x-direction as the longitudinal coordinate axis di~ ~ ˆ Ux~i ˆ U~i†. For simplicity, in the following we will write U instead of jUj. rected downstream, U With this assumption, U represents water velocity in the motion direction …x† averaged over river cross-section (m/s) and Eq. (3) becomes o … AZDxU † ˆ Q; …4† ox where A is the river width averaged over length-section and Dx is the cell length (Fig. 2). Eq. (4) can be solved if we consider two adjacent cells (i and i ‡ 1) by introducing ®nite differences Ai‡1 Zi‡1 DxUi‡1 ÿ Ai Zi DxUi ˆQ Dx

…5†

and we ®nally obtain Ai‡1 Zi‡1 Ui‡1 ÿ Ai Zi Ui ˆ Q:

…6†

If we simplify Eq. (2), we obtain (assuming again one-dimensional case in the x-direction) U

o…qU † oq…Z ‡ D† U2 ‡g ‡ gq 2 ˆ 0: ox ox C Z

…7†

If we consider two adjacent cells (i and i ‡ 1) by introducing ®nite di€erences we can rewrite Eq. (7)       qi‡1 Ui‡1 ÿ qi Ui qi‡1 Zi‡1 ÿ qi Zi qi‡1 Di‡1 ÿ qi Di U2 U ‡g ‡g ‡ gq 2 ˆ 0; …8† Dx Dx Dx C Z where qi ˆ Ai Zi Dx:

…9†

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

61

U , Z and q have been calculated by averaging over two adjacent cells: Uˆ

Ui‡1 ‡ Ui ; 2

…10†

qˆ

qi‡1 ‡ qi ; 2

…11†

Zˆ

Zi‡1 ‡ Zi : 2

…12†

Eqs. (6) and (8) conform a two-equation system with two unknowns Ui‡1 and Zi‡1 (cell i ‡ 1 velocity and water depth, respectively). We use as initial values of velocity and water depth those obtained from the nine gauging stations placed along the river. That is, the gauging stations data will be the boundary conditions for our model. Concerning Chezy ``friction'' coecient (C), there are described in the literature several methods to obtain it. Some authors calculate Chezy coecient by making use of expressions [11±14] that can become very complicated or by using some parameters which have to be evaluated directly from river observations (soil type, river slope, kind of river bed vegetation, etc.) [15±17]. However, there is another way [5,18,19] in which Chezy coecient is obtained from the model calibration. This last approximation is the one we have followed in this work. 2.2. Application of the hydrodynamic submodel to Jucar river We have applied the hydrodynamic submodel to the Jucar river. Our hydrodynamic submodel predicts velocities and water depth along the river for each one of the 63 cells (5.5 km long) we have subdivided the river. The validity of the model is determined by its ability or not to predict velocities and water depth on gauging stations without dams, as for the gauging stations located on dam exit the water ¯ow is regulated arti®cially. We have used for Chezy coecient a value of 0.049 m1=2 /s for the main part of modelled river length (56 cells, i.e., about 310 km), while for the last 40 km a value of 0.029 m1=2 /s has been assumed. These values are compatible with those used by other authors [5,11±19]. It is interesting to point out that the last stretch of the Jucar river is canalized arti®cially, which implies that the hydrodynamic friction coecient (Chezy coecient) should be lower in this part of the river, as is the case. In Fig. 3, we show velocities and water depth calculated with the hydrodynamic submodel together with the values recorded at the nine gauging stations. It can be seen that the model computed and recorded values agree quite well. For gauging stations numbers 2, 5, 6 and 7, a bigger di€erence is observed between predicted and recorded values. These stations are located on dams exit and, as it was explained before, the water ¯ow is regulated arti®cially and it is not possible to predict these values by means of river modelling. 2.3. The sediment submodel After the capability of our hydrodynamic submodel to predict the observed values of water depth and velocity of the river are veri®ed, we will try to describe sediment behaviour by studying the processes in which they are involved. These processes that describe water/sediments interaction are very important as ®nal radionuclide behaviour is mainly controlled by adsorption/desorption and lixiviation between water and sediment particles. Such processes depend not only of physical properties (sediment size, water velocity, water depth, etc.) but also on chemical characteristics (pH, conductivity, temperature, etc.) of water and sediment. Our sediment submodel will try to predict the amount of suspended matter and sediment and their

62

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

Fig. 3. (a) Water velocities calculated with hydrodynamic submodel along the Jucar river (), together with recorded

velocities (H); (b) water depth calculated with hydrodynamic submodel along the Jucar river (), together with recorded water depth (I).

time variation for each one of the cells in to which we have divided the river. The equations are based on the conservation of suspended matter and sediment. For suspended matter we can write o…qS† ~ ˆ QS ‡ … ÿ SED ‡ RES† q ; ‡ div…qUS† ot Z and for sediments

…13†

oB ˆ SED ÿ RES; …14† ot where Q, q, t and U have the same meaning as in Eq. (1), and S is the suspended matter (g/m3 ), B the bottom sediments (g/m2 ), SED the sedimentation ¯ux (g/m2 s) and RES is the resuspension ¯ux (g/m2 s). To evaluate the sedimentation and resuspension ¯uxes we have assumed that these ¯uxes are proportional to the suspended matter and bottom sediments amount, respectively [20±22]: SED ˆ K1 S;

…15†

RES ˆ K2 B;

…16†

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

63

where K1 and K2 are the sedimentation and resuspension constants, respectively. There are different ways to evaluate these parameters [23±25]. We have assumed that K1 is inversely proportional to water velocity while K2 is assumed to be directly proportional to water velocity. Moreover, we have imposed some restrictions to our sediment submodel in order to be as close as possible to the real situation one: 1. Suspended matter and sediments concentration should never reach negative values. If we obtain a negative concentration value the algorithm forces the sedimentation (resuspension) process to stop. 2. K2 has to take into account the natural ``resistance'' of suspended matter to increase inde®nitely. Some authors [22,25] explain this resistance in terms of the increase of the collisions among dissolved particles which tend to diminish their concentration in water. Thus, we can write K1 and K2 as: a …17† K1 ˆ ; U   S K2 ˆ bU 1 ÿ ; …18† SS where a and b are adjustable parameters and SS is the saturation concentration value of suspended matter. Moreover, better results are obtained if we assume a linear dependence on velocity (U) of K1 and K2 instead of high order dependence. Initial values of suspended matter …S† and bottom sediment …B† for each one of the cells that the river is subdivided are unknown a priori. Nevertheless, these values can be taken from averaged measurements of the speci®c river to be modelled. In the next section, we will give these initial values when applying the sediment submodel to Jucar river. To solve Eqs. (13) and (14), we follow the Lax±Wendrof two-step method [26]: 1. To solve the di€erential equations starting from the initial values. 2. We do a Dt time step and we solve the di€erential equations in time. The initial conditions at this stage are those calculated in point 1. 3. We obtain for each individual cell the new variable values, by doing an increase of Dx in space coordinates. 4. Steps 2 and 3 are repeated to cover the whole time interval to be studied. Thus Eq. (13) can be rewritten as t t
ÿ q Ui‡1 Si‡1 ÿ qiÿ1 Uiÿ1 Siÿ1 qi t‡1 …Si ÿ Sitÿ1 † ‡ i‡1 ÿ QSit ‡ ÿ K1 Sit ‡ K2 Bti Ai Dx ˆ 0; 2Dt 2Dx

…19†

where i, i ‡ 1 and i ÿ 1 refer to adjacent cells (and spatially separated by Dx), while the labels t, t ‡ 1 …ˆ t ‡ Dt† and t ÿ 1 …ˆ t ÿ Dt† refer to the time considered. Eq. (14) that accounts for sediments behaviour can be written as Bt‡1 ÿ Bitÿ1 i ˆ ‡K1 Sit ÿ K2 Bti ; 2Dt where

…20†

Sit ˆ

Sitÿ1 ‡ Sit‡1 ; 2

…21†

Bti ˆ

Bitÿ1 ‡ Bt‡1 i ; 2

…22†

64

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

K1 ˆ

a ; Ui

  Sitÿ1 ; K2 ˆ bUi 1 ÿ SS

…23†

…24†

substituting Eqs. (21)±(24) in Eqs. (19) and (20), we have a two-equation system with two unand Sit‡1 that can be solved for each one of the cells by using as boundary conditions knowns: Bt‡1 i the recorded values at the gauging stations. We have to point out that in our model we consider Sitÿ1 to obtain the resuspension constant K2 , Eq. (24), instead of Sit . This approximation simulates mathematically the slowing down of the suspended matter ``production'' process in order to avoid it to increase inde®nitely. Finally, our model has been done in such a way that time steps Dt accomplish the Courant±Friedercih±Lewy rule [27] jU jDt 6 Dx:

…25†

2.4. Application of the sediment submodel to Jucar river We have applied the sediment submodel we have developed to Jucar river. As we did previously for the hydrodynamic submodel, our model predicts now the sediment and suspended matter behaviour for each one of the 63 cells into which we have divided the Jucar river. We have solved the equations assuming Dt ˆ 600 s. By using lower Dt values, we have not observed signi®cative di€erences on the model results. However, Dt cannot be arbitrarily big, as we must be to ful®l the Courant±Friedercih±Lewy rule (Eq. (25)). We have obtained K1 and K2 by ®tting averaged suspended matter and sediments concentrations predicted by our model with those measured previously [1,2,28]. We will use in our model a ˆ 5  10ÿ4 m2 /s2 and b ˆ 1  10ÿ5 mÿ1 (Eqs. (17) and (18)). Concerning the suspended matter saturation value (SS, Eq. (18)), we have adopted for SS the higher value recorded at the nine gauging stations located in the Jucar river during two years …1996; 1997†. This implies SS ˆ 100 g/m3 . As it was pointed out, initial values of suspended matter …S† and bottom sediment …B† for each one of the cells that the river is subdivided are unknown a priori. Nevertheless, we have assumed as initial suspended matter values, the mean value we have obtained from the Jucar river samples (S initial  30 g/m3 , for all cells). However, we have not any measurement for sediment bottom concentration …B†. We have found that if we assume di€erent values for initial bottom sediment concentration only a change on the ®tted a and b values is observed. For this reason, we have assumed a value of Binitial ˆ 100 g/m2 , for all cells. This value is in agreement with those normally used in the literature [7]. Despite using conservation matter equations (Eqs. (13) and (14)) we obtain small ¯uctuations, nearly periodic in time, for suspended matter and sediments concentration averaged over the 63 cells (Fig. 4). These small ¯uctuations are mainly due to: 1. Resuspension/sedimentation processes which are considered by our equations. 2. Tributary discharges of water imply little contribution/leakage of suspended matter to our system. We have veri®ed that the results of our model depend slightly on the initial values assigned to suspended matter and bottom sediment concentrations. This could be explained if we assume that the sedimentation/resuspension processes are mainly governed by the hydrodynamic conditions of the system. However, as it can be expected, the solutions of our equations are

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

65

Fig. 4. Model prediction for the temporary evolution of suspended matter and sediment concentration spatially

averaged over the 63 cells in to which the Jucar river is subdivided.

strongly a€ected if we change the resuspension/sedimentation parameters (K1 and K2 , Eqs. (15) and (16)). In this sense, when changes on K1 and K2 are done, bottom sediment concentrations predicted by our model are more a€ected than calculated suspended matter concentrations. In fact, in the limit K1 ˆ K2 ˆ 0, bottom sediment concentrations predicted by the model do not depend on time (Eq. (20)). However, this is not the case for suspended matter concentrations (Eq. (19)). The predicted spatial distribution of suspended matter and bottom sediment for t ˆ 123 h (5 days) and t ˆ 7:5 days can be seen in Fig. 5. These times will have some interest when studying activities behaviour along the Jucar river (see the next section). It is interesting to notice that suspended matter concentration can vary, for an individual cell, up to 100% during this interval time, while bottom sediment variation within a speci®c cell is relatively much smaller. As it was already pointed out, the predicted suspended matter concentration are zero on dams while bottom sediment concentrations are constant on these sites. Along the Jucar river there are four dams: cells 7 (``Picazo'' dam), 31 and 32 (``Embarcaderos'' dam), 38 (``Millares'' dam) and 46, 47, 48, 49 (``Tous'' dam).

66

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

Fig. 5. (a) Model prediction of suspended matter concentration for the 63 cells in which the Jucar river has been divided;

(b) same as in (a) for bottom sediment concentration. Two di€erent t values have been considered.

2.5. The radionuclides submodel The last step on our modelling of radionuclides behaviour on rivers is the radionuclides submodel. Radionuclides submodel will describe the behaviour (in space and time) of radionuclides concentration in water, suspended matter and bottom sediment. Our model should include transport, adsorption and leaching of radionuclides. As the mean life of radionuclides studied is much bigger than the interval time considered in our model, no decay of these radionuclides will be considered in our equations. The main assumptions in our model are the following: 1. Water/sediment interaction is the sum of sediments leaching by water and sediments adsorption from water, which implies radionuclides transfer among water () suspended matter () bottom sediment. 2. Adsorption process is considered proportional to water activity of the radionuclide considered while leaching process is assumed to be proportional to suspended matter activity. This implies that leaching process of bottom sediments is negligible. 3. The proportionality constant for the above processes (we will assume the same value for both processes as they are inverse processes) will depend on the radionuclide considered, although very small variations are expected for the same chemical element. We call this coecient Kd in m3 /g units.

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

67

If equilibrium is reached, Kd will be equal to activities ratio between suspended matter and water. That is, if we have a situation of equilibrium for ionic exchange, according to Kd de®nition we can write AS ˆ Kd AW ;

…26†

where AW is the water activity (mBq/m3 ) for the radionuclide considered and AS is the suspended matter activity (mBq/g) for the radionuclide considered. We are going to calculate the number of atoms, and consequently the activity, that ``passes'' from one phase to another. This calculation is very similar to that done in [29]. If the number of the radionuclides in water is bigger than those expected in an ionic equilibrium situation, then N atoms of the radionuclide considered with associated activity DAW (mBq/m3 ) will be ``transferred'' to suspended matter, in such a way that for the new situation we can write AS ‡

DAW ˆ Kd …AW ÿ DAW †; S

DAW ˆ

S …Kd AW ÿ AS †; 1 ‡ Kd S

…27† …28†

where we have followed the same notation as previously. Eq. (28) is also valid when the transfer takes place from suspended matter to water (in this case DAW will be negative). Then, we can write the equations that govern radionuclides in water as follows:     o o DAW q S qAS S q ˆ Kd AW ; …qAW † ‡ …qUAW † ˆ ÿ …29† ÿ ot ox s 1 ‡ Kd S 1 ‡ Kd S s s where s (in s) is the characteristic time for the ionic exchange between dissolved phase and suspended matter phase. Similarly, for suspended matter we can write     o o S qAS S qAW …qSAS † ‡ …SqUAS † ˆ ÿ ‡ Kd ot ox 1 ‡ Kd S 1 ‡ Kd S s s q …30† ‡ …RES AB ÿ SED AS † ; Z where Z is the river water depth, Eq. (2), and AB is the bottom sediment activity (mBq/g) due to the radionuclide of interest. For bottom sediments, the equation that governs radionuclide activity behaviour is o …BAB † ˆ SED AS ÿ RES AB : ot

…31†

To solve this system of three di€erential equations (29)±(31), we will follow the same procedure described before for the sediment submodel (Lax±Wendrof two-step method [26]). In this way, Eq. (29) transforms into
qi‡1 Ui‡1 AtWi‡1 ÿ qiÿ1 Uiÿ1 AtWiÿ1 qi ÿ t‡1 tÿ1 ‡ AWi ÿ AW i 2Dt 2Dx     t qAtSi qAWi Sit Sit ˆ ÿ Kd : t t 1 ‡ Kd Si s 1 ‡ Kd Si s

…32†

68

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

For suspended matter we solve Eq. (30) as follows: t t qi‡1 Ui‡1 AtSi‡1 Si‡1 ÿ qiÿ1 Uiÿ1 AtSiÿ1 Siÿ1 qi t‡1 tÿ1 tÿ1 S ÿ A S † ‡ …At‡1 Si i 2Dt Si i 2Dx   t   t t
ÿ qASi qAWi Si Sit ˆÿ ‡ Kd ‡ ÿ K1 Sit AtSi ‡ K2 Bti AtBi Ai Dx t t 1 ‡ Kd Si s 1 ‡ Kd Si s

…33†

and for bottom sediments we solve Eq. (31) t‡1 tÿ1 tÿ1 Bt‡1 i ABi ÿ Bi ABi ˆ ‡K1 Sit AtSi ÿ K2 Bti AtBi ; 2Dt

…34†

where AtSi ˆ

AStÿ1 ‡ At‡1 Si i ; 2

…35†

AtBi ˆ

t‡1 Atÿ1 Bi ‡ ABi ; 2

…36†

AtWi ˆ

AWtÿ1 ‡ At‡1 Wi i : 2

…37†

t‡1 t‡1 We have a system of three equations with three unknowns (At‡1 Wi , ASi and ABi ) that can be solved by assuming as initial boundary conditions experimental values of activities in water, suspended matter and bottom sediments obtained along the system of interest.

3. Application of the model to Jucar river Finally, once we have modelled the major processes which govern the dynamic behaviour of radionuclides in rivers we have applied our model to the Jucar river. Uranium, thorium and radium activity measurements were done previously by the authors in 10 points along the Jucar river [1,2,28]. This will allow us to verify the model results by comparing with these experimental data. Moreover, we will use in our equations distribution factors calculated with these measurements [1,2]. As it was done previously for the sediment submodel, we have solved the equations assuming Dt ˆ 600 s and Dx ˆ 5:5 km. From the di€erent uranium, thorium and radium radionuclides we have measured activity data in Jucar river, we have applied our model to 238 U, 230 Th and 226 Ra isotopes. We have chosen these radionuclides as they are the most abundant in nature for the chemical element considered (in the case of uranium and thorium) while 226 Ra is the only radium isotope we had experimental values for the three phases involved in the model (water, suspended matter and bottom sediments samples) [1,2,28]. By assuming as initial boundary conditions experimental values of activities in water, suspended matter and bottom sediments obtained by the authors along Jucar river, we have allowed the system to evolve by means of the equations that de®ne the model. We have found, similarly as it was for the sediment submodel, that the results of the model vary slightly if the initial values vary within certain ``limits''. To study this e€ect we have used two di€erent initial conditions:

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

Table 1

Distribution factors (Kd ) for ID. 1 2 3 4 5 6 7 8 9 10

238

U,

232

Th and

Kd (238 U) (1/kg) 3

…0:9  0:3†  10 …1:1  0:3†  103 …0:52  0:04†  103 …0:37  0:07†  103 …1:45  0:09†  103 …1:56  0:09†  103 …1:6  0:3†  103 …1:04  0:12†  103 …0:77  0:04†  103 …0:19  0:02†  103

226

69

Ra along the Jucar river (obtained from [1,2]) Kd (232 Th) (1/kg) 5

…1:4  0:2†  10 …3:1  0:7†  105 …1:2  0:2†  105 …1:1  0:3†  105 …2:3  0:5†  105 …1:8  0:3†  105 …4:2  0:8†  105 …2:0  0:4†  105 …1:4  0:2†  105 …1:5  0:3†  105

Kd (226 Ra) (1/kg) …18:4  1:7†  103 …8:4  0:9†  103 …24  2†  103 …9:3  1:2†  103 …6:3  0:7†  103 …9:9  1:8†  103 …24:2  1:9†  103 …21:4  1:7†  103 …11:1  1:4†  103 …3:3  0:5†  103

1. For each radionuclide considered, its initial activity is assumed to be the mean value of the experimental activity (10 sampling points) we have measured along the Jucar river (i.e., the same value for the 63 cells into which we have subdivided the river). 2. By taking the measured activity as initial activity values for the 10 cells in which we have experimental values. For the rest (53 cells) the initial activity value is assigned by interpolation from the experimental one values. After about ®ve iterations in time (50 min) the results of the model do not depend on which one, from the above initial conditions, has been used. We have used for Kd (Eq. (26)) the values of Table 1, which were obtained previously [1,2]. In Fig. 6, we show the model prediction for the temporary evolution of 238 U, 232 Th and 226 Ra activities in water, suspended matter and bottom sediment, respectively. The activity values are space averaged over the 63 cells into which the Jucar river is subdivided. The averaged activity values oscillate slightly around its initial mean value. This behaviour is quite similar to what we have found before for the averaged suspended matter and bottom sediments concentrations (Fig. 4). Although not shown in Fig. 6, the same behaviour is predicted for the rest of the phases. Nevertheless, this similarity on predicted temporary evolution for sediments and radionuclides should be expected as the equations that govern the radionuclides evolution are quite similar to those used for sediments. We show in Figs. 7±9 the predicted spatial distribution for 238 U, 232 Th and 226 Ra activities for t ˆ 123 h (5 days). For this t value, mean di€erences between predicted and experimental values are lower than 20% for all the radionuclides considered. It is very important to point out that these status of minimum di€erence are repeated periodically over a period of about 7 days, while each minimum di€erence condition lasts for about 6 h. Moreover, these situations of minimum di€erence happen at the same time for the three radionuclides considered. These results induce us to investigate about the possibility that the status of minimum di€erence be associated to speci®c suspended matter (bottom sediments) situations. We have veri®ed that di€erences between suspended matter concentrations for status of minimum di€erence which are also lower than 20%, as it was the case for radionuclides activities. The same situation occurs for concentrations of bottom sediments, whose di€erences are lower than 20% during the situations of minimum di€erence. It should be remarked that the agreement between the activities predicted by our model and the experimental values (H in Figs. 7±9), which are assumed as initial values …t ˆ 0† for our model, is quite good. We show also in Figs. 7±9 the model predictions for a random t value (about 2 days after the minimum di€erence situation discussed above). We have found that for a speci®c cell, predicted variations in time can be much higher than those we found for the mean value activities

70

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

Fig. 6. (a) Model prediction for the temporary evolution of 238 U activity in water; (b) same as (a) for 232 Th in suspended

matter; (c) same as (a) for 226 Ra in bottom sediments. The activity values are spatially averaged over the 63 cells in which the Jucar river is subdivided.

averaged over the 63 cells (Fig. 6). These variations in the activity for a speci®c cell at two di€erent t values can reach up to 100% (Figs. 7±9). Finally, the agreement between predicted and experimental activities will allow us to apply our model to a nearly space-time point source of 137 Cs in water, due to a hypothetical Cofrentes Nuclear Plant accident. 4. Modelling of

137

Cs release from Cofrentes Nuclear Power Plant

As a direct application of our model we have modelled a hypothetic accidental release of 137 Cs in Jucar river water by Cofrentes Nuclear Power Plant. Cofrentes Nuclear Plant is located at about 190 km from the Jucar river mouth (Fig. 1), the Power Plant outlet being placed just in the river. The 137 Cs activity amount (space-time point source) we have simulated as accidentally ejected by the Power Plant is 1:3  1010 Bq in aqueous phase. This quantity seems to be quite

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

71

238 U activity in suspended matter; (b) same as (a) for bottom sediments; (c) same as (a) for water. The (H) represent the experimental values, which are assumed as initial values for the model. Two di€erent t values have been considered.

Fig. 7. (a) Model prediction for the spatial distribution of

realistic [9]. We have assumed Kd ˆ 14  10ÿ3 m3 gÿ1 for 137 Cs [9], as no experimental values for such radionuclide are available for Jucar river. Our model does not take into account meteorological e€ects, i.e., strong rains during a 3±4 day period. We show in Fig. 10(a) the predicted spatial distribution for 137 Cs activity for t ˆ 6 h, t ˆ 8 h (for which 137 Cs reaches Embarcaderos dam) and t ˆ 59 days (for which 137 Cs starts to exit from Embarcaderos dam). Notice that we assume the whole 137 Cs release from Nuclear Power Plant occurs at t ˆ 0. It is interesting to point out that two months after the 137 Cs release, the radionuclide contamination has moved forward only 20 km from its original ejection point, due to the Embarcaderos dam e€ect on retarding water velocity. Another important e€ect of Embarcaderos

72

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

Fig. 8. (a) Model prediction for the spatial distribution of 232 Th activity in suspended matter; (b) same as (a) for bottom

sediments; (c) same as (a) for water. The (H) represent the experimental values, which are assumed as initial values for the model. Two di€erent t values have been considered.

dam is the 137 Cs activity dilution (103 mBq/l, which implies a factor of 102 below the initial values) due to mixing processes inside the dam. We see, Fig. 10(a), that our model predicts a behaviour quite similar for water, bottom sediment and suspended matter. In Fig. 10(b), we show the predicted spatial distribution for 137 Cs activity for t ˆ 60 days (137 Cs reaches Millares dam), t ˆ 84 days (137 Cs starts to exit from Millares dam) and t ˆ 85 days, just when the 137 Cs contamination arrives at Tous dam. Again it is interesting to point out that due to Millares dam the ``heading face'' of the polluted wave is delayed for a period of about 24 days (see Fig. 10(b)).

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

73

Fig. 9. (a) Model prediction for the spatial distribution of 226 Ra activity in suspended matter; (b) same as (a) for bottom

sediments; (c) same as (a) for water. The (H) represent the experimental values, which are assumed as initial values for the model. Two di€erent t values have been considered.

In Fig. 11, we show the model prediction for the spatial distribution of 137 Cs activity at four di€erent times. It is interesting to notice that the polluted heading face needs about 83 days to pass through Tous dam (from day 85±86 up to day 168). Nevertheless, it only takes one day to go from Tous dam to the Jucar river mouth (cell number 63). Finally, we can conclude that after 170 days the 137 Cs contamination arrives at the mouth of the river. Two di€erent zones can be observed, for times lower than that needed by the contaminated heading face to arrive at river mouth (Fig. 11): 1. Before Tous dam (cell 46). In this zone 137 Cs radioactivity levels are of the order of 6  102 mBq/l in water (10 Bq/g for bottom sediments and suspended matter).

74

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

137 Cs activity in suspended matter, bottom sediment and water for three di€erent evolution times related with Embarcaderos dam; (b) same as (a) but for three di€erent times related with Millares and Tous dams (see text).

Fig. 10. (a) Model prediction for the spatial distribution of

2. After Tous dam (cell 50). In this zone 137 Cs radioactivity levels are much lower (1 mBq/l in water while for bottom sediments and suspended matter we have about 7  10ÿ3 Bq/g). Consequently, the main part of 137 Cs radioactivity dilution takes place at Tous dam. We show also in Fig. 11 the expected spatial distribution of 137 Cs activity about 7 months after the simulated accident. It can be seen that at this time (t ˆ 208 days), 137 Cs contamination for cells close to the Cofrentes Nuclear Plant has become nearly zero. This is because we have not considered any migration term in the equations of our model. Nevertheless, 137 Cs radioactivity still remains at high levels after the ®rst dam, located at cells 31±33 (Embarcaderos dam). This fact shows us that well after the release of 137 Cs, dams could act not only as the main responsible for

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

Fig. 11. Model prediction for the spatial distribution of

75

137

Cs for four di€erent times after a hypothetic accident from Cofrentes Nuclear Power Plant. The polluted heading face arrives at the Jucar river mouth after 169 days.

the dilution of 137 Cs but also as continual re-emitters of such element, becoming for a very long time secondary pollution sources.

5. Conclusions The main purpose of this work was the description and modelling of the major processes which govern the dynamic behaviour of radionuclides in rivers. To do this, we have developed a multicompartmental box model of the river in order to consider not only hydrodynamic equations but also other complex phenomena as adsorption/desorption of radionuclides to sediments. We have developed the model by simulation of three complementary aspects of river behaviour: Hydrodynamic, sediment and radionuclide submodels.

76

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

The model has been applied to the Jucar river (Spain) for 238 U, 230 Th, 226 Ra and 137 Cs radionuclides. We have considered three di€erent phases: water, suspended matter and bottom sediment the river being divided into 63 cells each 5.5 km long. We have used as input data for our model together with dynamic variables (river ¯ow and the water depth for several gauging stations) the experimental values for the distribution factors (Kd ) of 238 U, 230 Th and 226 Ra measured previously [1,2]. We have found that our model reproduces quite well the experimental values for the activities of 238 U, 230 Th and 226 Ra radionuclides, being the mean di€erence between predicted and experimental values lower than 20% for all the radionuclides considered. We have found that activities model predictions are nearly periodic in time, being associated to nearly periodic concentration levels of suspended matter and bottom sediments. Finally as a direct application of our model we have modelled a nearly space±time point source of 137 Cs (1:3  1010 Bq) in water due to an hypothetical Cofrentes Nuclear Power Plant accident. For this study no cesium migration e€ects have been considered and the Jucar river has been assumed as the only way for the transmission of 137 Cs. We can conclude that the heading face of released 137 Cs needs about 170 days to arrive at the Mediterranean sea, mainly due to the delaying e€ect of the dams located along Jucar river. Finally, we have observed that our model predicts that well after the release of 137 Cs, t ˆ 7 months, dams could act not only as the main responsible for the dilution of 137 Cs but also as continual re-emitters of such element, becoming for a very long time, secondary pollution sources. References [1] M.J. Rodrõguez-Alvarez, F. S anchez, The transfer of uranium from sediment to water along Jucar river (Spain), J. Radioanal. Nucl. Chem. 242 (2) (1999) 297±307. [2] F. S anchez, M.J. Rodrõguez-Alvarez, E€ect of pH, temperature, conductivity and sediment size on thorium and radium activities along Jucar river (Spain), J. Radioanal. Nucl. Chem. 242 (3) (1999) 671±681. [3] M.J. Rodrõguez-Alvarez, F. S anchez, E. Navarro, in: M. Garcõa-Le on, G. Madurga (Eds.), Proceedings of the Third International Summer School, Huelva, Spain, World Scienti®c, Singapore, 1994, p. 549. [4] D.T. Pugh, Tides, Surges and Mean Sea-Level. A Handbook for Engineers and Scientist, Wiley, Chichester, 1987. [5] D. Prandle, A numerical model of the southern North Sea and river Thames, Institute of Oceanographic Sciences, Report 4, IOS/R14, 1974. ~ez, J.M. Abril, M. Garcõa-Le [6] R. Peri an on, Modelling the dispersion of non-conservative radionuclides in tidal waters - part 1: conceptual and mathematical model, J. Environ. Radioactivity 31 (2) (1996) 127±141. ~ez, J.M. Abril, M. Garcõa-Le [7] R. Peri an on, A modelling study of 226 Ra dispersion in an estuarine system in southwest Spain, J. Environ. Radioactivity 24 (1994) 159±174. [8] J.S. Smitz, E. Everbecq, E. Frere, Modelling of radionuclides transport by sediments in waterways, Third Research Coordination Meeting Casaccia, Roma, 2±7 December 1995. [9] P. Benes, M. Cernik, O. Sl avic, Modelling of migration of 137 Cs accidentally released into a small river, J. Environ. Radioactivity 22 (1994) 279±293. [10] W. Li, Numerical simulation of cooling water Shantou thermal power plant, in: Lee, Cheung, Balkema, Rotterdam (Eds.), Environmental Hydraulics, 1991, pp. 1217±1223. [11] J.J. Dronkers, Tidal computations for rivers, coastal area and seas, in: Proceedings of the American Society of Civil Engineers, J. Hydraul. Div. 95 (HY1) (1969) 29±77. [12] D. Knight, Hydraulics of ¯ood channels, in: K. Beven, P. Carling (Eds.), Floods: Hydrological, Sedimentological and Geomorphological Implications, Wiley, Chichester, 1989, pp. 83±105. [13] Committee on Hydromechanics of the Hydraulics Division, Friction Factors in open Channels, Progress Report of the Task Force on Friction Factors in Open Channels, Proceedings of the American Society of Civil Engineers, J. Hydraul. Div. 89 (HY2) (1963) 97±143. [14] J.F. Lyness, W.R.C. Myers, Comparisons between measured and numerically modelled, in: W.R. White, J. Waltts (Eds.), Unsteady Flows in a Compound Channel Using Di€erent Representations of Friction Slope, Second International Conference on River Flood Hydraulics, Wiley, Chichester, 1994, pp. 383±991.

M.J. Rodriguez-Alvarez, F. Sanchez / Appl. Math. Modelling 25 (2000) 57±77

77

[15] H.B. Herbich, S. Shulits, Large-scale roughness in open-channel ¯ow, in: Proceedings of the American Society of Civil Engineers, J. Hydraul. Div. 90 (HY6) (1964) pp. 203±230. [16] V.T. Chow, Open Channel Hydraulics, McGraw-Hill, New York, 1959. [17] V.T. Chow, L.W. Mays, D.R. Maidment, Applied hydrology McGraw-Hill series in water resources and environmental engineering, 1998. [18] L.K. Ghosh, S.C. Patel, M. Arora, Numerical modelling of ¯ow and wind induced circulation of hot water discharged into a reservoir, in: Lee, Cheung, Balkema, Rotterdam (Eds.), Environmental Hydraulics, 1991, pp. 883±888. [19] F. Karim, Bed con®guration and hydraulic resistance in alluvial-channel ¯ows, J. Hydraul. Eng 121 (1) (1995) 15±25. [20] O. Starosolszky (Eds.), Applied Surface Hydrology,Water Resources Publications, Budapest, 1987. [21] J.C.J. Nihoul, B.M. Jamart (Eds.), Three-Dimensional Models of Marine and Estuarine Dynamics, Elsevier, Amsterdam, 1987. [22] R.J. Garde, K.G.R. Raju, Mechanics of Sediment Transportation and Alluvial, Stream Problems, edited by Wiley Eastern Limited, 1995. [23] G. Govers, Empirical relationships for the transport capacity of overland ¯ow, in: D.E. Walling, A. Yair, S. Berkowicz, (Eds.), Erosion, Transport and Deposition Processes. Proceedings of the Jerusalem Workshop, March±April 1987, IAHS 189 (1990) 45±63. [24] F. Engelung, J. Fredsoe, in: V. Chow (Ed.), Advances in Hydroscience, Academic Press, Vol. 13, 1982, 187±214. [25] P.D. Komar, Flow-competence evaluations of the hydraulic parameters of ¯oods: an assessment of the technique, in: K. Beven, P. Carling (Eds.), Floods: Hydrological, Sedimentological and Geomorphological Implications, Wiley, Chichester, 1989, pp. 107±134. [26] W.H. Press et al., Numerical Recipes the Art of Scienti®c Computing, Cambridge University Press, Cambridge, 1986. [27] D. Prandle, A modelling study of the mixing of 137 Cs in the seas of the European continental shelf, Philos. Trans. R. Soc. London, Ser. A 310 (1984) 407±436. [28] M.J. Rodrõguez-Alvarez, Medidas de Uranio,Torio y Radio en Muestras Ambientales por Espectrometrõa alfa. Aplicaci on a un Modelo Unidimensional de Transporte de Radion uclidos en el Rõo J ucar, Ph.D. Thesis, Universidad de Valencia, Spain (in Spanish), 1998. [29] J.M. Abril, M. Garcõa-Le on, A 2D 4-phases marine dispersion model for non-conservative radionuclides part 1: conceptual and computational model, J. Environ. Radioactivity 20 (1993) 71±88.