Numerical simulation of radioactive pollution diffusion in 3D stream

In an earlier study, we showed how Pulsed Neutron Activation (PNA) analysis could be implemented to estimate the flow speed in a one-dimensional pipeline [1]. The neutron capture by target stable nuclides, like Oxygen-16 in water, leads to radioactive Nitrogen-16 emitting high-energy gamma-rays, which travels with flowing water in a channel and can be detected and analyzed downstream in the pipe to get information about its average transport time. Then by having the traveled distance, i.e., the distance between the gamma source and the gamma detector, one can obtain the speed of the activity and, consequently, the flow speed and mass rate in the pipeline. This report shows how we can extend the technique to study pollution diffusion, including radioactive pollution, in a three-dimensional stream.

TRANSPORT EQUATION IN THREE DIMENSIONS

The transport of pollution in a three-dimensional stream with the diffusion coefficient D flowing with velocity v, when there is source or sink , can be described by the convection-diffusion equation

Ct=(DC)(Cv)+,(1)

where C is the pollution concentration. The transport equation (1) in a Cartesian coordinate system, where the fluid is moving with the constant velocity of v can be represented as

C(x,y,z,t)t=(Dx2C(x,y,z,t)x2+Dy2C(x,y,z,t)y2+Dz2C(x,y,z,t)z2)(2)
(vxC(x,y,z,t)x+vyC(x,y,z,t)y+vzC(x,y,z,t)z)(3)
+(x,y,z,t)(4)

where

  • C(x,y,z,t): the pollution concentration at location (x,y,z) and time t,

  • Dx,Dy,Dz: the horizontal and vertical components of diffusion coefficient. For an anisotropic medium, like rivers and lakes, Dx,Dy,Dz are unequal, while for an isotropic medium, the diffusion proceeds at the same rate in horizontal and vertical directions,

  • vx,vy,vz: the flow velocity components in the x, y, and z directions. They are constant in each spatial direction for the studied contaminated region,

  • (x,y,z,t): the pollution emission at location (x,y,z) and time t.

Considering the boundary conditions, the pollution concentration must converge to zero as x, y, or z tends to infinity

limx±C(x,y,z,t)0,(5)
limy±C(x,y,z,t)0,(6)
limz±C(x,y,z,t)0;(7)

and the initial condition,

C(x,y,z,t=0)=Mδ(xx0)δ(yy0)δ(zz0),(8)

where M is the released point-source pollution at location (x0,y0,z0) and time t=0. In a 1D model in xdirection, MMLyLz represents the surface mass flux or total mass per unit area, where LyLz is the area scale of neglected dimensions. Similarly, in a 2D model in the xy direction, MMLz represents linear mass flux or total mass per unit length, where Lz is the area scale of neglected dimension. The total mass of pollution at time t must be equal to the number of undecayed pollution at time t, i.e. M(t). So, the boundary and initial conditions are

V𝑑VC(x,y,z,t)=M(t).(9)

If the source of the contamination is radioactive material, the mass of undecayed material at time t is given by the radioactive decay law

M(t)=M0eλt,(10)

where λ=ln(2)T1/2 is the decay constant for the induced radioactive material. For an instantaneous point-source pollution

C(x,y,z,t)=M0(4πt)3/2(DxDyDz)1/2(11)
×exp((xx0vxt)24Dxt(yy0vyt)24Dyt(zz0vzt)24Dztln(2)T1/2t).(12)

NUMERICAL RESULTS

To study the pollution diffusion from a radioactive source in a river or lake, we develop a computational toolkit in Matlab to solve the transport equation and calculate pollution concentration as a function of time at different locations.

Instantaneous release of radioactive pollution in a 1D pedagogical stream:

For the first implementation, we study the diffusion of radioactive pollution from an instantaneous release in a 1D water canal flowing at a constant speed. In the image above, we show the pollution release in a one-dimensional water canal flowing in the x direction with a constant speed vx. The Matlab scripts listed in Appendix A, create 2D motion plots (in AVI and GIF formats) for the evolution of pollution concentration.

FIG. 2: Few screenshots of 2D motion plots obtained from the Matlab toolkit that calculates the pollution diffusion obtained from instantaneous release from radioactive neutron source 16N (with a half-life T1/2 = 7.13 s) as a function of time t and distance x, in a water canal flowing with a constant speed vx = 0.3 m/s and a diffusion coefficient

Instantaneous release of radioactive pollution in a 2D pedagogical stream:

In this section, we study the diffusion of radioactive pollution from an instantaneous release in a 2D water stream flowing at constant speeds vx and vy in the x and y directions. The Matlab scripts listed in Appendix B, create 2D motion plots (in AVI and GIF formats) for the evolution of pollution concentration from a release at point ( x0, y0).

FIG. 3: Few screenshots of 2D motion plots obtained from the Matlab toolkit that calculates the pollution diffusion obtained from instantaneous release from radioactive neutron source 16N (with a half-life T1/2 = 7.13 s) at (x0 = 0 m, y0 = 0.5 m) as a function of time t and distances x and y, in a water river flowing with a constant speeds vx = 0.3 m/s and vy = 0.0 m/s and diffusion coefficients Dx = Dy = 0.001.

Instantaneous release of radioactive pollution in a 3D pedagogical stream:

In this section, we study the diffusion of radioactive pollution from an instantaneous release in a 3D water stream flowing at constant speeds vx, vy, and vz in the x, y, and z directions. The Matlab scripts listed in Appendix B, create 3D motion plots (in AVI and GIF formats) for the evolution of pollution concentration from a release at point (x0, y0, z0).

FIG. 3: Few screenshots of 3D motion plots obtained from the Matlab toolkit that calculates the pollution diffusion obtained from instantaneous release from radioactive neutron source 16N (with a half-life T1/2 = 7.13 s) at (x0, y0, z0) as a function of time t and distances x and y, in a water river flowing with a constant speeds vx = 0.3 m/s, vy = 0.2 m/s, and vz = 0.1 m/s and diffusion coefficients Dx = 0.001 Dy = Dz = 0.01.