The Mathematical Modelling Of Radionuclide Transport In The Subsurface Environment Around The Chornobyl Nuclear Power Plant

The Shelter constructed above the destroyed Unit-4 of the Chor-nobyl Nuclear Power Plant contains 20 MCi of nuclear fuel. More than 1000 m 3 of intermediate-level radioactive water is disposed in its basement. The purpose of this work is to simulate migration of the radionuclides 90 Sr and 137 Cs from the Shelter into the subsurface environment to evaluate their migration rate and migration paths. A mathematical model accounting for the coupled transport of water and radionuclides in variably saturated media is used. The lack of suitable experimental and eld studies excluded the possibility of complete validating the model. Results of the simulations will be useful for future studies of the environmental impact of the Shelter.


Introduction
A Shelter containment structure constructed upon the ruins of the destroyed reactor building of the Chornobyl Unit-4 was completed in November 1986.This containment was built in very high radioactive elds in the short time span of six months by using unique and remote construction methods.However, these remote methods for pouring fresh concrete did not always permit good observation and view of work, and in a number of cases high radioactive elds made welding key containment components impossible.As a result of this, openings and breaks in the walls and the roof of the Shelter are estimated to total about 1200 m 2 .
By current estimates, the Shelter contains from 135 to 180 tons of nuclear fuel with an activity of up to 20 000 000 Ci 1].Within the Shelter the largest part of the nuclear fuel from the completely destroyed reactor core is contained in fuel-containing materials (FCM), and the rest { in fuel assemblies and their fragments.The FCM can be represented as lava-like materials with a glass silicate matrix containing many inclusions of various uranium compounds, nely dispersed high-level active fuel dust, and aerosols.Experimental observations show a gradual increase in the migration rate of radionuclides from the FCM caused by their permanent interaction with the water and atmosphere inside the Shelter.Contamination of the environment by radioactive substances from the Shelter may occur in two principal ways: (1) release of fuel dust into the atmosphere through the openings in the Shelter if there is a collapse of internal structures within the Shelter, and (2) migration of radionuclides through the subsurface together with the water which is inside the Shelter.
According to previous research 2,3] the Shelter contains more than 3000 m 3 of intermediate-level and low-level radioactive water.More than 1000 m 3 of intermediate-level radioactive water is disposed in the basement of Unit-4.Most of the contaminated water is located in the turbo-generator hall and rooms on the rst level in the northern part of the Shelter 4]  Shelter, water pathways and locations within the Shelter.North-south cross-section along Axis 47. the Shelter as: precipitation, water used for dust suppression, moisture condensation onto building structures, and water leakage from the service lines of the 3rd unit of the Chornobyl NPP.Most precipitation onto the roof penetrates through cracks to the interior of the Shelter.The volume of this precipitation is approximately 1000 m 3 per year.For dust suppression, 15 m 3 of water are sprayed twice a month within the Shelter.Water enters the Shelter through the above mentioned openings and ows downwards into the lower regions of the reactor building where most of the fuel is located.Upon movement through the FCM and interaction with them, the water is contaminated by radionuclides.The contaminated water penetrates into the subsurface environment.The total activity of the intermediate-level radioactive water lies within the range of 2 10 12 to 9 10 13 Bq/l.Cesium isotopes constitute the main portion of the water activity; the cesium activity varies from 1.6 10 5 to 5.5 10 7 Bq/l.The total 90 Sr activity in the water was determined to be 3.6 10 3 to 1.1 10 6 Bq/l.
The Chornobyl NPP is located on the right bank of the Prypiat river.The NPP site is a at river bench.The aquifer system is composed of ne, medium and sometimes coarse alluvial sands with interlayers of sand loam and loam.Medium sands prevail.According to eld studies the hydraulic conductivity of the alluvial sands ranges from 6.0 to 15.0 m day ?1 .These alluvial sands the thickness of which varies from 25 to 40 m are underlain by clay marls of the Kyiv Eocene suite.The marls have very low permeability and act as an aquitard.The in ltration is estimated between 80 mm and 200 mm per year 5].
The purpose of this work is to simulate migration of 90 Sr and 137 Cs through geologic media due to di erent scenarios of water leakage from the Shelter.

Model description 2.1. Water ow equations
Mathematical description of water transport in a variably saturated subsurface implies assumptions which are commonly done for this type of problems.It is assumed that uid transport through the soil occurs in response to pressure gradients and gravitational body forces following Darcy's ow equations.The porous medium is considered to be rigid.Pressure gradients in the gas phase are so negligible that gas pressure remains e ectively constant at the atmospheric pressure.
Neglecting hysteresis and temperature gradients, the movement of moisture through a saturated-unsaturated porous medium is described by @ @t ( S ) = @ @x i k r k ij @p @x j + g @z @x j where t is the time (T ); x i are the cartesian spatial coordinates (L); is the porosity of the medium; is the uid density (ML ?3 ); k ij are components of the intrinsic permeability tensor (L 2 ); k r is the relative permeability (0 6 k r 6 1); is the uid dynamic viscosity (ML ?1 T ?1 ); p is the uid pressure (ML ?1 T ?2 ); g is the gravitational constant (LT ?2 ); z is the elevation measured from a reference datum (L); f represents sources or sinks of uid in the system (ML ?3 T ?1 ); and S is the volumetric saturation.Initial conditions for equation (2.1) consist of initial values of pressure, speci ed by p(x; 0) = '(x); where x is a spatial coordinate vector and ' is a prescribed function.
Boundary conditions may be Dirichlet (type 1): p(x; t) = g 1 (x; t) on ? 1 ; or Neumann (type 2): k r k ij @p @x j + g @z @x j n i = g 2 (x; t) on ? 2 ; where g 1 is the prescribed pressure on the boundary ? 1 ; n i is the ith component of the unit vector normal to the boundary ? 2 , directed outward; and g 2 is the prescribed outward uid ux normal to ? 2 .? 1 and ? 2 comprise the entire boundary of the domain.Usually ? 1 is the part of boundary at which surface-water bodies are intercepted by aquifers; ? 2 are the boundaries across which uid passes at a speci c rate (impervious boundaries, areas of in ltration or evapotranspiration, points of injection and withdrawal, etc.).

Transport equations
Let c l represent the volumetric radionuclide activity in the aqueous phase (ML ?3 ); c s { the volumetric exchangeable sorbet activity in the solid phase (ML ?3 ); and c f { the volumetric radionuclide activity which is xed in the mineral lattice (ML ?3 ).Assuming that the ion-exchange reaction transferring activity between aqueous and solid exchangeable phases has a suciently short timescale, we will consider this reaction as instantaneous which is described by a linear equilibrium isotherm.Transfers of activity to and from the xed phase are modelled by rst-order rate constants sf and fs , respectively, as shown in gure 2 6].
. Schematic representation of the kinetic sorption model.The leaching of radionuclides from the fuel particles is modelled by a rst-order equation: @c p @t = ?(p + )c p ; c p (x; 0) = c 0 p (x); where c p is the radionuclide activity in the fuel particles per soil solid volume (ML ?3 ); p is the rst-order constant of radionuclide leaching from fuel particles (T ?1 ); is the radionuclide decay constant (T ?1 ); c 0 p is the initial radionuclide activity in the fuel particles (ML ?3 ).
Ignoring chemical di usion in the solid phase, an equation for adjectivedispersive transport of radionuclide activity in the exchangeable phase in variably saturated media may be written in the form: @c @t = @ @x i SD ij @c l @x j ?@ @x i (v i c l ) ? c + (1 ?T )( p c p + fs c f ?sf c s ); (2.2) where c is the radionuclide activity per total volume (ML ?3 ).
D ij is the dispersion tensor de ned by where L and T are the longitudinal and transverse dispersiveness (L), respectively; ij is the Kronecker delta; j v j is the magnitude of the Darcy velocity; D 0 is the molecular di usion coe cient (L 2 T ?1 ); and is the tortuosity.The total radionuclide activity is de ned by: where T is the total porosity of the medium; and volumetric phase activities are interrelated through the solid-aqueous partition coe cient k d , according to the relationship where s is the soil solid density (ML ?3 ).Initial conditions for equation (2.2) consist of the speci cation of initial radionuclide activity in the aqueous phase: c l (x; 0) = c 0 (x): Boundary conditions for in ow boundaries are Dirichlet (type 1): c l = c 1 on ? 3 ; or Cauchy (type 3): where c 1 is the prescribed radionuclide activity on the boundaries ? 3 ; n i is the unit vector normal to boundaries, directed inward; v b is the prescribed uid ux normal to the boundaries ? 4 ; and c b is the prescribed radionuclide activity associated with the prescribed uid ux.
At out ow boundaries the assumption is often made that transport across the boundaries, according to 7], occurs by advection only and the Neumann (type 2) boundary condition is speci ed by @c l @x j = 0 on ? 5 ; where ? 3 + ? 4 + ? 5 = ?is the boundary of the domain.
The transport of activity in the xed phase is described by: @c f @t = sf c s ?fs c f ?c f ; where c 0 f is the initial activity in the xed phase.

Simulation results
Numerical simulations of the 90 Sr and 137 Cs migration from the Shelter into and through the subsurface environment were conducted for a twodimensional vertical cross-section selected along the streamlines of the groundwater ow.These streamlines were obtained from eld observations of groundwater level distribution in the vicinity of the Chornobyl NPP 5].The geology of this cross-section is shown in gure 3.
The MSTS code 8,9] adapted to implement the model described above was used for the numerical simulations.The ow domain was discretized by a non-uniform rectangle mesh, the nodes of which were spaced at 1 m { 100 m intervals in the horizontal direction, and from 0.2 to 3 m in the vertical direction.Time steps of 30 days were used.Values for the hydraulic and species properties used in the simulations are given in table  The prescribed values of the hydraulic head were speci ed along the left and right boundaries, the top of the Shelter's foundation, and the boundaries on which surface water bodies are intercepted by the aquifer.At the ground surface the liquid ow rate was equal to the in ltration rate or evapotranspiration rate.The annual amount of in ltration was 80 mm.The marls were considered as impervious.The water depth was assumed to be 2 m in the turbo-generator hall, and 1.6 m in rooms on the rst level in the northern part of the Shelter.Seasonal uctuations of the water level in the canal varied between 109.5 m and 110.2 m; in Azbuchin Lake { between 106.0 m and 106.8 m; and for the hydraulic head at the left boundary { between 109.5 m and 110.0 m 5].These seasonal uctuations of the water level in surface water bodies and the variability over time of precipitation and evapotranspiration were taken into account in the simulations.
Before the Chornobyl accident a special drainage system operated to prevent under ooding of the Chornobyl NPP building's foundation.After the Chornobyl accident this drainage system was damaged and shut down.As a result of that a groundwater level increase was observed until 1994.The nowadays groundwater level below the Shelter varies from 109.5 m to 110.5 m depending on seasonal groundwater uctuations and precipitation.The groundwater increase after the Chornobyl accident was not taken into account in the simulations.The initial condition for the liquid phase is illustrated in gure 4.
A zero radionuclide ux condition was employed at all the boundaries, # # # excluding the Shelter basement and the boundaries that contact surface water bodies.The 90 Sr activity of water from the basement was 2:3 10 6 Bq/l, and 2800 Bq/l in Azbuchin Lake 5].Similarly, 137 Cs activity in the basement water and Azbuchin Lake was 7.1 10 7 and 84 Bq/l, respectively.From 1986 to 1995, the radionuclide activity of canal water was determined from the annual average activity of water samples from the Cooling pond.
The radionuclide activity of canal water after 1995 was speci ed as the radionuclide activity of canal water in 1995.For the area surrounding the Chornobyl NPP, land surface contamination by 137 Cs in 1995 was from 3.8 10 ?3 to 1.9 10 ?4 Ci per kg of soil 2]. 90Sr land surface contamination was equal to one half of 137 Cs and surface contamination.Just after the accident, exchangeable 90 Sr was 5% of the total 90 Sr activity and exchangeable 137 Cs was 10% of the total 137 Cs activity.The subsurface was assumed to be uncontaminated before the accident.The date of the Chornobyl accident was taken as the starting date for the simulations.
There are no reliable data on water leakage beyond the bounds of the Shelter.However, analysis of water pathways and water locations within the Shelter shows that there is a possibility of water leakage from the turbogenerator hall, northern rooms of the Shelter, and through the cascade wall.The cascade wall was lled with concrete, boxed radioactive materials collected from the roof and surroundings, building debris, and structural metals.So, it has a lot of openings and breaks through which part of precipitation onto the cascade wall penetrates into the subsurface.Therefore, in order to estimate migration rates of radionuclides from the Shelter through the subsurface, numerical simulations were conducted for the following cases: Case 1. Water leakage occurs from rooms on the rst level in the northern part of the Shelter.Case 2. In addition to case 1, one half of the precipitation onto the cascade wall penetrates into the outside of the Shelter.Case 3. Water from the turbo-generator hall, rooms in the northern part of the Shelter, and a half of the precipitation onto the cascade wall penetrates into the subsurface environment.The e ects of various water leakage scenarios on radionuclide migration from the Shelter through the subsurface are compared with case 2. Figures 5  and 6 show the predicted activity distribution of 90 Sr and 137 Cs radionuclides in the aqueous phase in 1995 and 2045 for case 2.
The in uence of di erent water leakage into the outside of the Shelter on migration of 90 Sr through the subsurface is shown in gures 7-9.Analysis of Comparisons of predicted activity-depth pro les of 90 Sr and 137 Cs radionuclides in the aqueous phase with the measurements in 1995 are given in gures 10 and 11.These predicted activity-depth pro les were selected at point x = 1075 m which approximately corresponds to the location of wells 1-g, 2-g and 4-g in the northern part of the Chornobyl NPP building area.The annual mean values of radionuclide activity obtained on measurements in wells 1-g, 2-g, 4-g and 6-g 2,3] are denoted by circles.A triangle denotes the radionuclide activity detected in well C10 in 1995   and denote mean values of radionuclide activity detected in well C10 and in wells 1-g, 2-g, 4-g and 6-g, respectively.
is located approximately 150 m east of wells 1-g and 4-g.

Conclusions
This work has concentrated on modelling 90 Sr and 137 Cs migration from the Shelter into and through the subsurface environment.A two-dimensional mathematical model accounting for the coupled transport of water and radionuclides in variably saturated porous media was used.The hydraulic conductivity, the longitudinal and transverse dispersiveness are sensitive parameters for radionuclide transport modelling.The key parameters for radionuclide exchange between soil water and soil matrix are the distribution coe cient or k d , and rst-order rate constants sf and fs .Due to the limited availability of eld data, these parameters were assigned by experts.As a result of that only a qualitative radionuclide transport modelling was achieved.However, modelling results presented in this work illustrate the basic e ects of di erent possible water leakage into the outside of the Shelter on the migration rate of radionuclides through the subsurface.The simulation results showed that, so far, the groundwater contamination detected in the observation wells is only de ned by the surface land contamination.Much more quantitative testing must be conducted to evaluate the environmental impact of the Shelter on the subsurface contamination.

Figure 10 .Figure 11 .
Figure 10.Comparison of predicted and measured activity-depth pro les of 90 Sr in the aqueous phase in 1995.Modelgenerated pro les were obtained for (|) case 2 and(---) case 3. 4 and denote mean values of radionuclide activity detected in well C10 and in wells 1-g, 2-g, 4-g and 6-g, respectively.
(see gure 1).It is also reasonable to consider such major water sources within Water Figure 1.Schematic representation of the Chornobyl Unit-4

Table 1 .
1. Values of the rst-rate constants sf , fs were taken similar to the experimentally obtained values of these parameters for the upper soil layer of the left bank oodplain of the Prypiat river.Parameters sf and fs vanished for the Values of parameters used in numerical simulations.