Convection-Diffusion Model for Radon Migration in a Three-Dimensional
Confined Space in Turbulent Conditions

Convection and diffusion are the main factors affecting radon migration. In this paper, a coupled diffusion-convection radon migration model is presented taking into account turbulence effects. In particular, the migration of radon is simulated in the framework of the k-ε turbulence model. The model equations are solved in a complex 3D domain by the finite element method (FEM). Special attention is paid to the case study about radon migration in an abandoned air defense shelter (AADS). The results show that air convection in a confined space has a great influence on the radon migration and the radon concentration is inversely proportional to the wind speed.

Abstract: Convection and diffusion are the main factors affecting radon migration. In this paper, a coupled diffusion-convection radon migration model is presented taking into account turbulence effects. In particular, the migration of radon is simulated in the framework of the k-ε turbulence model. The model equations are solved in a complex 3D domain by the finite element method (FEM). Special attention is paid to the case study about radon migration in an abandoned air defense shelter (AADS). The results show that air convection in a confined space has a great influence on the radon migration and the radon concentration is inversely proportional to the wind speed. Time difference between the two-point measurement time (h) k Turbulent kinetic ε Turbulence dissipation rate P k Turbulen ε t kinetic (k) due to the average velocity gradient E Radon exhalation rate (Bq/m 3 /s) C 1 ,C 2 Concentration of radon at two adjacent moments (Bq/m 3

) A
Opening area of radon-collecting hood (m 2 )

Introduction
Pollutants in the confined space are more likely to accumulate, causing harm to people health [1,2]. Radon and its daughters are special among these pollutants because of their radioactivity causing lung cancer [3][4][5]. Exposure to radon gas and its daughters has long been recognized as a potential health hazard, and long-term inhalation of high radon concentrations may cause lung cancer [6]. Therefore, knowledge of the radon migration in the confined space is of great significance for predicting radon concentration, which provides guides for reducing the harm of radon. In previous studies, common methods are one-dimensional or two-dimensional analytical solution methods [7,8]. However, these methods are not applicable to the three-dimensional case with complex boundary conditions. At present, there are a few researches on this issue. Rabi et al. [9] simulated the radon concentration distribution in an experimental model room by FORTRAN software. However, the model is relatively simple, only considering the distribution of radon in laminar flow. Collignan et al. [10] and Akbari et al. [11] simulated the indoor radon migration based on convective-diffusion without considering the effect of radon decay term. When the radon migrates in the confined space with high air velocity or complex geometric structure, the turbulence flow must be considered for the numerical model [12,13]. In this paper, we established a three-dimensional model of radon migration in confined space with consideration of turbulence flow and radon diffusion-convection process. A case study of radon migration in an AADS was numerically analyzed by the model.

Radon Migration Model
The main theory of radon migration is mainly a combination of the diffusion and convection flow of radon gas [14][15][16]. The radon migrating in the confined space can be simulated using the following equation [10,15]: where C is the radon concentration distribution (Bq/m 3 ), D is the diffusion coefficient of radon in the air (1.05 Â 10 -5 m 2 /s), U is the velocity (m/s), λ is the decay constant of radon (2.1 Â 10 -6 /s), F is radon source term (Bq/m 3 /s).
Assuming that radon source term in the air is equal to 0, the steady-state radon migration model in confined space may be expressed as: The key to numerically solving Eq. (2) is to determine the second term, ∇(UC) which represents the radon convective velocity is closely related to the distribution of the radon concentration in space, and they are coupled with each other. Because air flow in the confined space is highly susceptible to the turbulent flow, it is necessary to use the turbulence model to obtain the convective velocity in the confined space, so as to establish a coupling relationship between the convective velocity and radon diffusion.

Air Turbulence Model in Confined Space
In this paper, the air turbulent flow is simulated using the k-ε turbulence model. The k-ε turbulence model is a two-equation model to simulate the air turbulent flow. The model includes two parameters: k (turbulent kinetic) and ε (turbulence dissipation rate) [17,18]. In general, the air turbulent flow in the confined space is assumed to be small and incompressible so that it can be simulated by solving two variables (velocity and length scale).

Numerical Approach
The FEM is used to numerically calculate the steady-state radon migration model and air turbulence model in the confined space. The FEM is an efficient numerical calculation method, which includes three processes: mesh, element analysis and solving approximate variational equations [21]. In the FEM, the weak form of the three-dimensional steady-state radon migration model is expressed as: where Ω is the integration domain, in this paper, the whole confined space.
If the radon exhalation rate (the Neumann boundary) is used as a boundary condition, according to the Green formula, Eq. (7) is rewritten as [22]: where g is the radon exhalation rate, and Γ is the boundary of the radon exhalation. Eq. (8) is the weak form of the radon migration model in three-dimensional space.
It is generally believed that free tetrahedral mesh element (FTME) is a good choice for threedimensional FEM, and its interpolation function can be expressed as: where a, b, c and d are constants. Assuming that the values of the four vertices of the FTME are C 1 , C 2 , C 3 and C 4 , respectively, so: From Eqs. (9) and (10), it can be obtained as: where N 1 (x, y, z), N 2 (x, y, z), N 3 (x, y, z), N 2 (x, y, z) are vertex coordinate functions. Eq. (14) can be expressed as a matrix form, . Such matrix forms facilitate programming and computeraided computing. Firstly, the element stiffness matrix is used to integrate the Eq. (8) on the grid element.
Secondly, the whole model is calculated by using the balance principle. Thirdly, use constraints to solve linear equations.
The conventional CFD method was applied to numerical calculation of gas turbulence model, which is not discussed in this paper. Schematic diagram of the model proposed in this paper can be seen in Fig. 1. First, the CFD method is used to calculate the convective velocity of the confined space, which is substituted into the steady-state radon migration model. Then, the partial differential equation of the steady-state radon migration model is numerically calculated to obtain the distribution of radon concentrations in the confined space by finite element method.

Model Validation
In this research, a room of dimensions 3.01 m Â 3.01 m Â 3 m (Fig. 2) reported by Chauhan et al. [23] Fig. 3, and the related results proposed by the literature are shown in Fig. 4. Noted that the model for numerical calculation in the literature is time-dependent (transient model), while our model is a steady-state model for predicting radon concentration for an infinitely long time mathematically. This difference determines that calculation results in this paper are larger than those in the literature, and it is not difficult to find from Fig. 4 that the calculated radon concentration in the literature gradually approaches the results in this paper with the increase of time. The verification shows that the method proposed in this paper is effective and reliable.

Case Study
AADS is a typical confined space with high radon concentration, which is potentially harmful to people health. In this section, we use the aforementioned model to analyze radon migration in an AADS. The AADS  The initial radon concentration is determined by the radon exhalation rate of the surfaces of the 4 lanes and storage room in the AADS. The radon exhalation rate is measured by the local static method, as shown in Fig. 6. The local static method measures the radon exhalation rate by the increase of radon concentration in unit time by a sealed container [24]. The radon exhalation rate can be obtained by where E is the radon exhalation rate (Bq/m 3 ·s), C 1 and C 2 are radon concentrations at two times (Bq/m 3 ), A is the area of the radon exhalation (m 2 ), Δt is the time difference between two times (h), V is the volume of radon-collecting hood (m 3 ). The radon exhalation rates of the AADS are listed in Tab. 1.
The inlet velocity can be calculated by the ventilation rate. The ventilation rate of the shelter is 10 h -1 under the condition of natural ventilation. Suppose the AADS has one inlet surface, the amount of inlet air is equal to that of outlet air and there is no air flow from the lining, ground and walls in the AADS, the velocity of the inlet boundary is expressed as [9]: where λ V is the ventilation rate (h -1 ), A vent is the ventilation area (3.57 m 2 ), V cave is the volume of the AADS (338.50 m 3 ). The pressure of the outlet boundary is 0 Pa.

Results and Discussion
Figs. 7 and 8 present numerical calculation results of wind speed and radon concentration of the AADS at z = 1 m. The entrance wind speed of the No. 1 lane is small, and the distribution is more uniform, showing the state of the laminar flow. The radon concentration of the No. 1 lane is also uniformly distributed. In the storage room, due to the complex geometry, the air flow in the storage room is in a turbulent state. The higher concentration of radon from the No. 1 lane follows the air flow into the No. 2 lane, forming a bend-shaped migration path in the storage room, while the radon concentration in the other areas of the storage room is small. Due to the rapid change of the geometry of the air flow channel, the air flow of the No. 2 lane is changed from a relatively smooth layer to a turbulence distribution. The wind speed at the inside of the corner is small, while the outside is large. The distribution of radon concentration is basically consistent with that of the wind speed in the AADS, which means that radon migration in confined space is greatly affected by air convection.      [25]. The radon concentration at the entrance of the confined space is the highest and the minimum at about 1/5 of the longitudinal length, that is close to 2.18 m. The reason is that the diffusion length of radon, defined as the distance it diffuses during its mean lifetime, is 2.18 m in air. Beyond this distance, the radon diffusing from the entrance is neglected, therefore, the radon concentration is the lowest at 2.18 m away from the entrance.   Fig. 12, it can be seen that if the diffusion term is removed from the radon migration model, the numerical results of the radon concentration will be inconclusive, that is because the velocity of the wall in the k-ε turbulence model is 0, when the radon exhalation rate is taken as the boundary condition, exhalation radon cannot participate in convective action. From Fig. 13, if the convection term is removed   from the radon migration model, the radon concentration in the confined space shows a significant gradient change from the precipitated center because of the radon migration only depending on the diffusion of the concentration gradient.

Conclusions
In this paper, a three-dimensional model is proposed to study radon migration in confined space and numerically calculated by the FEM. In the model, radon convection migration is simulated by the k-ε turbulence model, while the Fick law is used for radon diffusion migration, and the two are coupled to reflect the interaction between air convection and diffusion during radon migration. The validation study shows that the proposed model is effective and reliable. The model was used to analyze radon migration in an AADS, which is a typical confined space. Results revealed that radon migration in confined space is greatly affected by air convection, and the radon concentration is affected by the turbulence flow and inversely proportional to the wind speed. The radon concentration at the entrance of the confined space is the highest and the minimum at about 2.18 m away from the entrance. By comparing with the simulation results of removing diffusion and convection term from the model, when the radon exhalation rate is taken as the boundary condition, it will cause the calculation results to deviate greatly, affecting the calculation accuracy.