The Study of Radioactive Fallout Source of Low-Equivalent Nuclear Bursts Based on Nuclear Cloud Simulation Using the CFD-DPM

Abstract

The activity-height distribution of radioactive particles in the stabilization cloud of a nuclear burst plays a crucial role in the radioactive fallout prediction model, serving as the source for transport, diffusion, and dose rate calculation modules. A gas-particle multiphase flow solver was developed using the OpenFOAM Computational Fluid Dynamics (CFD) library and discrete phase method (DPM) library under a two-way coupling regime to simulate the U.S. standard atmosphere of 1976 with good stability. The accuracy of the numerical model was verified through low-equivalent nuclear weapons tests, including RANGER-Able and BUSTER-JANGLE-Sugar, depicting reasonable spatio-temporal changes in cloud profiles. The initialization module of the Defense Land Fallout Interpretative Code (DELFIC) and activity-size distribution, which considered fractionation, were employed for nuclear fireball and radioactive particle initialization. Simulations indicated that the activity-height distribution of the stabilization cloud mainly concentrated on the lower third of air burst cloud caps, while settling near the burst center for surface or near-surface bursts. This study has confirmed the effectiveness of the gas-particle flow solver based on the CFD-DPM method in simulating low-equivalent nuclear clouds and enriching research on radioactive fallout prediction models.

1. Introduction

The residual radiation from nuclear bursts is a long-term effect that can cause lasting damage to the environment and human health over several years [1,2]. Predicting radioactive contamination has been a continuous concern since the birth of nuclear weapons. In the design of the radioactive fallout prediction model, the nuclear cloud and radioactive particles are critical source items that significantly impact the accuracy of the prediction model [3]. However, there is limited research available on the nuclear radioactive particle cloud in terms of a source item for fallout prediction models. Therefore, studying the nuclear radioactive particle cloud holds increasing research significance in civil nuclear defense and military operations.

Some studies have indicated that the accuracy of dose rate prediction is highly dependent on the distribution of radioactive particles in the stabilization cloud [4]. The historical nuclear weapons test image recognition provides sequential height data for the nuclear cloud rise, but there may be discrepancies between the visual range and actual distribution of radioactive particles, leading to deviations in observed fallout data compared to predictive results. The DELFIC (Defense Land Fallout Interpretative Code) model treats horizontal and vertical activity distribution as Gaussian [5], while the KDFOC3 (K-Division Defense Nuclear Agency Fallout Code) model describes the relationship between activity, particle size, and height as a “double-triangle” distribution [6]. Despite approximations suggesting that radioactive activity is concentrated in the lower third of the cloud cap [7], these models do not accurately reflect the spatial distribution of radioactive particles and activity.

After the comprehensive ban on nuclear tests, the numerical simulation using Computational Fluid Dynamics (CFD) has become a crucial method for studying nuclear clouds. When the nuclear fireball reaches pressure equilibrium with the atmosphere, this is referred to as the initial time [5]. Prior to the initial time, the reflected shock wave plays a significant role in the formation of the nuclear cloud. Some studies have utilized high Mach number models to simulate the process of the nuclear cloud rise in near-surface bursts [2,8,9]. After the initial time, only the rise of the fireball under buoyancy drive by low Mach number models should be considered [4]. Regardless of whether shock waves are involved, both models have yielded satisfactory predictions. Therefore, the simulation accuracy of spatio-temporal height of nuclear clouds is primarily attributed to modules for fireball formation and atmospheric environments. Some scholars utilize a single-phase flow model to conduct numerical simulations of nuclear cloud rise under the Eulerian coordinate system, without taking into account the cloud-borne radioactive particles. In this approach, the nuclear cloud is represented by high-temperature gas masses through a temperature or density field [10]. To obtain the spatio-temporal distribution of radioactive particles in the nuclear cloud, it is necessary to consider the multiphase flow method. Based on this, most researchers employ the discrete phase method (DPM) or discrete element method (DEM) to simulate weapon debris and soil engulfed by the fireball, coupled through one or two-way regimes with the gas phase in order to realize gas-particle cloud rise simulations [8,11]. Despite the widespread use of gas-particle coupling methods for describing explosive dust clouds, researchers have rarely studied in detail the spatial distribution of activity within stabilization clouds. Currently, most scholars primarily rely on semi-empirical theoretical models [7,12], or simple settlement models that ignore vertical turbulent mixing in order to approximately depict the activity distribution within stabilization clouds [5]. The increase in equivalent of burst generally leads to a higher stabilization height of nuclear clouds, disregarding atmospheric influences and focusing solely on the equivalent and height of the burst [13]. When nuclear clouds pass over the tropopause, they encounter more complex meteorological conditions with variable wind shear. To simplify these conditions for preliminary research purposes, nuclear bursts with an equivalent of less than 1 kiloton are categorized as low-equivalent and have a relatively lower stabilization height near the tropopause. At present, low-equivalent charges are commonly found in tactical nuclear weapons, such as the B61 family [14].

This paper utilized the Eulerian–Lagrangian multiphase flow solver based on the OpenFOAM framework to model the formation and evolution of nuclear clouds resulting from low-equivalent bursts. The atmospheric static equilibrium equations and “equilibrium fireball” model of DELFIC were integrated into the numerical framework to establish the initial conditions of stable atmospheric stratification and an initial nuclear fireball. The accuracy of the solver was validated using data from nuclear weapons tests, considering air bursts and surface or near-surface bursts, while also taking into account the distribution of radioactive activity among particles. OpenFOAM is an open-source CFD library widely employed in multiphase flow and heat transfer applications [15].

2. Materials and Methods

2.1. Datasets and Processing

2.1.1. Datasets of Historical U.S. Nuclear Weapons Tests

In this study, the validation of nuclear cloud simulations for low-equivalent nuclear bursts is conducted using datasets obtained from U.S. historical nuclear weapons tests, as documented in published test reports [16]. The time dependence of cloud upper (bottom)-edge height datasets from the RANGER-Able and BUSTER-JANGLE-Sugar nuclear weapons tests is utilized to validate air bursts and surface or near-surface bursts, respectively. The parameters of these tests are presented in

Table 1. The parameters of the nuclear weapons tests used in this paper.

Name	Burst Center Altitude (ASL *, m)	Height of Burst (AGL *, m)	Equivalent (kt)	Stabilization Cloud Upper (Bottom)-Edge Height (m)
RANGER-Able	957	323	1	4224 (2841)
BUSTER-JANGLE-Sugar	1285	1	1.2	4572 (3353)

* ASL, Above sea level; AGL, Above ground level.

.

2.1.2. Data Processing of U.S. Standard Atmosphere

The rise of a nuclear cloud is closely linked to factors, including burst equivalent, height of burst, burst center altitude, and meteorological parameters [17]. This paper utilizes the U.S. standard atmosphere parameters published by the U.S. National Oceanic and Atmospheric Administration (NOAA) in 1976 [18] to validate the accuracy of the model of atmospheric static equilibrium by considering changes in temperature, pressure, and density with altitude, as shown in Table S1 from the Supplementary Materials.

2.2. Researching Methodology

The flowchart in Figure 1 illustrates the methodology of this paper. Initially, the numerical model of the fluid phase of the nuclear cloud is established within the Eulerian framework. Subsequently, the numerical equation for the stable atmosphere is proposed based on a static equilibrium model, and the initial parameters for the nuclear fireball are obtained from the DELFIC-Initialization module. Following this, a trajectory tracking model for radioactive particles is developed within a Lagrangian framework, taking into account two-way coupling to simulate the behavior of radioactive particle-laden clouds. The rise time sequence of the nuclear cloud is then validated against historical test data to ensure the model’s reliability. Finally, activity-height distributions of the stabilization cloud are derived from number-size and activity-size distributions used in our numerical modeling framework.

2.3. Numerical Methods and Solver Development

The rise of a nuclear cloud is primarily driven by the difference in buoyancy between the high-temperature (low-density) nuclear fireball and the surrounding air. Once the temperature and density of the nuclear cloud reach equilibrium with the atmosphere, stabilization is considered to have been achieved. The buoyant rise process of a nuclear cloud is influenced by unsteady, compressible fluid mechanisms and turbulent mixing [17], with radioactive particles being carried by gas flow representing the burst product of noble gas and air [19]. Therefore, numerical modeling of nuclear cloud rise resembles low-Mach number heat transfer problems within the multiphase flow regime.

2.3.1. Governing Equations of Gas Flow (Eulerian Method)

In this study, the gas flow numerical method relies on the Finite Volume Method (FVM), which utilizes conservative integral forms for the governing equations of Eulerian phase [20]. The FVM is flexible due to direct discretization in physical space without transformation between physical and computational coordinate systems [21]. The governing equations are integrated over finite volumes within subdivided domains, with application of Gauss’ theorem to transform volume integrals of convection and diffusion terms into surface integrals, while the time derivative term is semi-discretized [22]. This work is currently based on solving the Reynolds Averaged Navier–Stokes (RANS) equations [22,23], where statistical averaging relies on time-averaged variables [5], as shown in Equations (1) and (2):

$${\displaystyle \frac{\partial \rho }{\partial t}}+\nabla \cdot \rho \mathbf{U}=0$$

(1)

$${\displaystyle \frac{\partial \rho \mathbf{U}}{\partial t}}+\nabla \cdot \left(\rho \mathbf{U}\mathbf{U}\right)=-\nabla p+\nabla \cdot \tau +\rho \mathbf{g}-S_p$$

(2)

where $\rho$, $\mathbf{U}$, $p$, and $\tau$ are density, velocity vector, pressure, and viscous stress tensor of the fluid phase, respectively. $t$ and $\mathbf{g}$ are the discretization of time and acceleration of gravity, respectively. $S_p$ is the source term. The pressure-based solver requires the state equation of the ideal gas to construct the pressure Poisson equation that is expressed by Equation (3):

$$p=\rho RT$$

(3)

where $R=291.0186\mathrm{J}/(\mathrm{k}\mathrm{g}\cdot \mathrm{K})$ and $T$ are the gas constant of air and temperature, respectively. The coupling of pressure–velocity–density of the continuous fluid is not necessary in solving the energy equation, which is given as

$${\displaystyle \frac{\partial \left(\rho h\right)}{\partial t}}+\nabla \cdot (\rho \mathbf{U}h)=\nabla \cdot (k\nabla T)+{\displaystyle \frac{\partial p}{\partial t}}$$

(4)

where $h$ and $k$ are specific enthalpy and thermal conductivity, respectively. The closure of the RANS equations requires additional equations for the unknown Reynolds stress components, as the set of RANS equations is not closed [24]. In this study, the standard k−ε turbulence model [22,23] is employed to close the RANS equations, which are presented in Equations (5) and (6):

$${\displaystyle \frac{\partial }{\partial t}}(\rho k)+\nabla \cdot (\rho \mathbf{U}k)-{\nabla }^2(\rho D_{k}k)=\rho G-\left({\displaystyle \frac{2}{3}}\rho \nabla \cdot \mathbf{U}k\right)-\left(\rho {\displaystyle \frac{\epsilon }{k}}k\right)$$

(5)

$${\displaystyle \frac{\partial }{\partial t}}(\rho \epsilon )+\nabla \cdot (\rho \mathbf{U}\epsilon )-{\nabla }^2(\rho D_{\epsilon }\epsilon )=C_1\rho G{\displaystyle \frac{\epsilon }{k}}-\left(\left({\displaystyle \frac{2}{3}}{C}_1-C_{3,RDT}\right)\rho \nabla \cdot \mathbf{U}\epsilon \right)-\left(C_2\rho {\displaystyle \frac{{\epsilon }^2}{k}}\right)$$

(6)

The standard k−ε is a typical two-equation turbulence model; Equations (5) and (6) are the governing equations for turbulent kinetic energy k and turbulent energy dissipation rate ε, respectively, where $G$, $D_k$, and $D_{\epsilon }$ are turbulent kinetic energy production rate, effective diffusivity for k, and effective diffusivity for ε, respectively. $C_1$ = 1.44, $C_2$ = 1.9, and $C_{3,RDT}$ = 0 are model coefficients [25,26].

2.3.2. Governing Equations of Particles Flow (Lagrangian Method)

According to some scholars [27,28], the primary forces acting on particles in a gas-solid multiphase system are the pressure gradient force, drag force, and gravity. Furthermore, the motion of individual particles is described by applying Newton’s second law in Lagrangian modeling, as shown in Equation (7).

$$m_i{\displaystyle \frac{d{\mathbf{U}}_i^p}{dt}}={\mathbf{F}}_{\nabla p}+{\mathbf{F}}_D+{\mathbf{F}}_g$$

(7)

where $m_i$ and ${\mathbf{U}}_i^p$ are the mass and velocities of particle i, respectively; ${\mathbf{F}}_{\nabla p}$, ${\mathbf{F}}_D$ and ${\mathbf{F}}_g$ are the pressure gradient force, drag force, and gravity acting on particle i, respectively. The pressure gradient force of the gas phase acting on particle i is given by

$${\mathbf{F}}_{\nabla p}=-V_i\nabla p$$

(8)

where $V_i$ denotes the volume of the spherical particle i, and $\nabla p$ is the pressure gradient. The drag force, as prescribed by Putnam [29] and Amsden et al. [30], is given as

$${\mathbf{F}}_D={\displaystyle \frac{3}{4}}{\displaystyle \frac{{\mu }_f{\mathrm{C}}_D{Re}_p}{{\rho }_p{d}_p^2}}$$

(9)

where ${\mu }_f$, $C_D$, and ${Re}_p$ are the dynamic viscosity of fluid at the cell occupying particle i, drag coefficient, and particle Reynolds number, respectively; ${\rho }_p$ and $d_p$ are the density and diameter of particle i, respectively. Based on Schiller and Naumann [31], the drag ${\mathrm{C}}_D$ is defined in Equations (10) and (11). ${Re}_p$ is expressed by Equation (12), which works for spherical particles and is suitable for high Reynolds numbers [27,32], as follows:

$$C_D={\displaystyle \frac{24}{{Re}_p}}\left(1+{\displaystyle \frac{1}{6}}{Re}_p^{{\displaystyle \frac{2}{3}}}\right)\mathrm{i}\mathrm{f}{Re}_p\le 1000$$

(10)

$$C_D=0.424\mathrm{if}{Re}_p>1000$$

(11)

$${Re}_p={\displaystyle \frac{\rho \left|{\mathbf{U}}_{\mathrm{r}\mathrm{e}\mathrm{l}}\right|d_p}{{\mu }_c}}$$

(12)

where $\rho$ and $\left|{\mathbf{U}}_{\mathrm{r}\mathrm{e}\mathrm{l}}\right|$ are the density of fluid at the cell occupying particle i and relative velocity between particle i and the fluid, respectively.

2.3.3. Numerical Solution Method of CFD-DPM

Because of the difference between Eulerian and Lagrangian methods, distinct numerical solutions for physical modeling need corresponding mathematical schemes [33].

In the aspect of the Eulerian system, the buoyantPimpleFoam solver in the heat-transfer library of OpenFOAM resources is utilized to solve the RANS equations for nuclear cloud rise problems [34,35,36]. This solver effectively handles the buoyancy term and low-Mach number flows with variable density, employing either the pressure-based PIMPLE or PISO algorithm [15,37,38]. A specific transformation, Equation (13), is used to manage the pressure gradient and buoyancy term, leading to a unified treatment in Equation (14) for improved robustness of the solver [39].

$$\left\{\begin{array}{l}{p}_{rgh}=p-\rho \mathbf{g}\cdot \mathbf{h}\\ \nabla p_{rgh}=\nabla p-\mathbf{g}\cdot \mathbf{h}\nabla \rho -\rho \mathbf{g}\\ -\nabla p+\rho \mathbf{g}=-\mathbf{g}\cdot \mathbf{h}\nabla \rho -\nabla p_{rgh}\end{array}\right.$$

(13)

$${\displaystyle \frac{\partial \rho \mathbf{U}}{\partial t}}+\nabla \cdot (\rho \mathbf{U}\mathbf{U})=-\mathbf{g}\cdot \mathbf{h}\nabla \rho -\nabla p_{rgh}+\nabla \cdot \tau -S_p$$

(14)

where the hydrodynamic pressure $p_{rgh}$ is used in the gradient term, $\mathbf{h}$ is the vector of the grid cell body center, and a new pressure gradient-buoyancy term appears on the RHS.

In terms of the Lagrangian system, the Discrete Parcel Method (DPM) is commonly utilized in gas-solid multiphase flow to alleviate computational burden. This method employs the concept of a particle group (parcel), where particles are organized into parcels based on their physical properties such as density and diameter. Each parcel is designed to ensure that particles within it share the same properties and experience common forces and motion history [27,40,41], resulting in low memory occupation [42]. Therefore, this work utilizes the discrete phase solver DPMFoam from the Lagrangian library of OpenFOAM resources.

Based on the classification of interaction types and volume fractions of particles [43], the multiphase flow solver is implemented with a two-way coupling scheme to account for momentum transfer from the particles to the fluid phase in the particle-laden cloud, based on Newton’s third law, which cannot be overlooked [44,45] in this study. Hence, the source term “$S_p$” is utilized to facilitate the two-way coupling, which is computed as a volumetric fluid-particle interaction force, expressed as follows:

$$S_p={\displaystyle \frac{{\sum }_{i=1}^{N_p} \left({\mathbf{F}}_{d,i}+{\mathbf{F}}_{\nabla p,i}\right)}{\rho V_{\mathrm{cell}}}}$$

(15)

where $V_{\mathrm{cell}}$ is the volume of a computational cell, and $N_p$ is the number of particles located in the cell in which the particles reside [33].

Furthermore, the impact of particle–particle interaction has been disregarded due to the primary focus on the transport and diffusion of radioactive particles, which is the influence of fluid on particles [17]. The flowchart illustrating the Eulerian–Lagrangian algorithm is presented in Figure 2, with the basic process outlined as follows:

• First, the Eulerian scheme reads the initial and boundary conditions as well as the initial information of particles for calculating the fluid phase interaction with the initial particles, and then transmits the momentum information to the Lagrangian system;
• Second, the gas forces (drag, buoyancy, etc.) acting on particles are obtained by the Lagrangian scheme, and then calculates the discrete phase (particle) momentum and position driven by the carrier gas;
• Third, the interaction force information for particles advances to the next time step and is adopted to renew the particle momentum and position of the Lagrangian scheme;
• Next, the latest particle information is communicated to the Eulerian scheme and is used to calculate the particle–fluid momentum interaction for the next time step of the Eulerian calculation;
• Finally, the loop between the Eulerian–Lagrangian framework continues until the termination time set by the customers, the particle positions are stored in the Eulerian grid cells of each time step and comprise the trajectories of particles.

2.4. Initialization of the Numerical Model

2.4.1. Construction of Stable Atmosphere Stratification

Meteorology significantly impacts the development of nuclear cloud rise [5], especially during the late stage of cloud rise. The vertical acceleration of cloud rise is influenced by atmospheric temperature layers, which results in a decelerating process. Additionally, the relatively low lapse rate of atmospheric temperature near the tropopause [46] leads to a large negative acceleration that strongly inhibits the rise of nuclear clouds.

The stabilization cloud of the low-equivalent burst generally remains within the troposphere, but when the equivalent passes over 1 kiloton, there is a significant chance that the cloud cap will enter the stratosphere. As the cloud rises, it is crucial to consider the impact of changes in ambient temperature, pressure, and density on its force [9]. Therefore, numerical discretization must be strictly carried out in accordance with the atmospheric thermal structure. Some scholars have utilized the U.S. standard atmosphere stratification published by NOAA in 1976 [31] as simulation input for initial ambient atmosphere [9,47], due to the universality of atmospheric temperature lapse rate and the exponential decay law of pressure.

The assumption for atmospheric static equilibrium can be summarized in the following three points [48]:

• First, the atmosphere is adiabatic;
• Second, the acceleration of the vertical motion of the atmosphere is far less than the acceleration of gravity, leading to the neglect of the acceleration term in the momentum equation [49];
• Third, the viscous force is neglected.

Therefore, the nuclear cloud rise can be considered to occur in a state of atmospheric static equilibrium, where the gravity of the atmosphere is balanced with the vertical pressure gradient force. Supposing that a block of air mass with a specific cross-sectional area and thickness, as shown in Figure S1 from the Supplementary Materials. The force equilibrium is given in Equation (16), its deformation is expressed in gradient form of the pressure term in Equation (17):

$$-dpdA-\rho \mathbf{g}dAdz=0$$

(16)

$$\nabla p=-\rho \mathbf{g}$$

(17)

where z is the height from the surface of the Earth, $dA$ and $dz$ are cross-sectional area and thickness of the air block, and other variables have the same meanings as in Section 2.3.1.

The assumption of static equilibrium in the troposphere implies that there is no change in enthalpy, which is equivalent to an adiabatic atmosphere. Therefore, the LHS of Equation (4) equals zero based on Assumption 1, and the momentum Equation (2) can be approximated as steady, without convection and Reynolds stress issues corresponding to Assumptions 2 and 3. Only the buoyancy and heat flux divergence terms are considered in the RANS equation, as shown in Equation (18). At this point, the RANS equation is equivalent to atmospheric static equilibrium.

$$\left\{\begin{array}{c}\nabla p=-\rho \mathbf{g},\\ \nabla \cdot [k\nabla T]=0,\end{array}\right.$$

(18)

Based on this, the U.S. standard atmosphere defines temperature as a linear function of height. Traditionally, up to 86 km, the function is expressed as a series of seven successive linear equations corresponding to seven layers denoted by the subscript b (0 to 6), but only up to the tropopause will be calculated in this work, as shown in Equation (19) [18].

$$T=T_b+L_b\cdot \left(Z-Z_b\right),$$

(19)

$T_b$ = 288.15 K is temperature of the first layer (b = 0), $L_b$ = $dT/dZ$ is the temperature gradient, identical to $-6.5$ $\mathrm{K}/\mathrm{k}\mathrm{m}$ (b = 0) or 0 $\mathrm{K}/\mathrm{k}\mathrm{m}$ (b = 1), corresponding to the troposphere and isothermal layers. $Z_b$ is the initial altitude of corresponding layer.

Consequently, the two equations of pressure versus height can be obtained from the numerical integration of Equation (18) [18]. The equation for the troposphere layer is given by Equation (20):

$$p=P_b\cdot {\left[{\displaystyle \frac{T_b}{T_b+L_b\cdot (Z-Z_b)}}\right]}^{{\displaystyle \frac{g_0}{R\cdot L_b}}}$$

(20)

and for the isothermal layer, it is given by Equation (21):

$$p=P_b\cdot \mathrm{e}\mathrm{x}\mathrm{p}\left[{\displaystyle \frac{-g_0\cdot \left(Z-Z_b\right)}{R\cdot T_b}}\right]$$

(21)

where $P_b$ is the reference-value, $P_b=\mathrm{101,325}\mathrm{p}\mathrm{a}$ (b = 0) or $P_b=1.278\times {10}^5$ $\mathrm{p}\mathrm{a}$ (b = 1). $g_0$ = 9.81 m/s² is the acceleration of gravity at sea level.

2.4.2. Construction of the Equilibrium Fireball of Nuclear Explosions

In previous studies, Lundquist et al. utilized the WRF-based atmospheric modeling framework to simulate the behavior of radioactive particles in the nuclear cloud from the “Grable” test [4]. The energy proportion of the “ideal hot bubble” was set at 1/6 by Lundquist et al., while Arthur et al. chose 1/3 equivalent [50]. Taking into account the fireball contacts the ground, using a proportion of 1/6 leads to a more realistic simulation. Hao et al. employed an empirical method to set the proportion of nuclear cloud energy at 35% equivalent for non-contact situations [11]. In the DELFIC initialization module, 45% of the nuclear explosion energy is allocated to heating the air, soil, and water within the fireball [5]. Glasstone and Dolan noted that approximately 35% to 45% of the energy from nuclear explosions in the atmosphere is received as thermal energy in both the visible and infrared portions of the spectrum [7]. Therefore, it is considered reasonable to assume that the fireball produced by nuclear weapons explosions contains an energy proportion of 30% to 45%. Additionally, the DELFIC provides a comprehensive calculation program for determining initial fireball parameters such as size and temperature increment. Hence, the DELFIC code has been selected to describe the initial fireball in this study.

The initial fireball is regarded as an oblate sphere with an eccentricity of 0.75 in the two-dimensional computational domain [5]. In order to define the initial fireball in the OpenFOAM framework, the oblate sphere fireball is transformed into a circular shape while keeping equal area to guarantee the consistent proportion of energy conversion. The equivalent circle radius is given by

$$R_0=\sqrt{H_{c,i}\cdot R_{c,i}}$$

(22)

where $R_0$, $H_{c,i}$, and $R_{c,i}$ are equivalent radius of the circle, vertical radius and horizontal radius of the oblate sphere, respectively.

In order to confirm the type of phase interaction, it is necessary to consider the volume fraction that reaches the maximum value at the initial fireball moment during the simulation. Therefore, the mass of fallout and the volume of the initial fireball for air bursts and surface or near-surface bursts are needed for calculation, as shown in Equations (23)–(25):

$$m_s=k_{\Lambda }{W}^{3/3.4}{(180-\Lambda )}^2(360+\Lambda ),0\le \Lambda \le 180$$

(23)

$$m_s=90.7,\Lambda >180$$

(24)

where $m_s$, $W$, and $\Lambda$ are the mass of fallout, equivalent, and scaled height of the burst, respectively. $k_{\Lambda }$ is a constant equal to 0.07741. Equation (23) is for the surface or near-surface bursts, and Equation (24) is for air bursts.

$$V_i=\left(m_{a,i}+m_{w,i}\right)R_a{T}_i^{*}/P$$

(25)

where V_i, R_a, $T_i^{*}$, and P are the initial fireball volume, ideal gas law constant of air, virtual temperature in the cloud at initial time, and ambient pressure. m_a,i and m_w,i are the mass of gas and water of the initial cloud, respectively.

Hence, the volume fraction of the radioactive particles can be obtained via Equation (26):

$$F=m_s/({\rho }_p\cdot V_i)$$

(26)

where F and ${\rho }_p$ are the volume fraction and density of particles, respectively. According to Equation (26), the volume fractions of radioactive particles range from 0.12 to 0.19% (4.7 × 10⁻⁶ to 8.8 × 10⁻⁵%) for the surface or near-surface bursts (air bursts) of equivalents ranging from 0.05 to 1 kiloton, which indicates that the interaction regime between particles and turbulence has satisfied the judgement for two-way coupling [43].

2.4.3. Construction of Radioactive Particles

The radioactive particles in the nuclear cloud are formed through the condensation of vaporized debris from the nuclear weapon and surrounding substances, generated by the tremendous temperature and entrainment during the initial stage of a nuclear burst. This process also involves the incorporation of fission products, leading to the radioactive property of these particles [51]. The size of these radioactive particles within the initial fireball can vary from sub-micrometers to millimeters, resulting in varying levels of radioactivity based on observed results. To achieve a reasonable distribution of activity, it is essential to determine the number-size distribution of these radioactive particles. Generally, the lognormal distribution is widely used to represent this number-size distribution [1,3]. Therefore, we have chosen to use a lognormal distribution for depicting the number-size distribution of radioactive particles. The number of particles $N\left({\delta }_i\right)$ within a diameter range (${\delta }_i{,\delta }_{i+1}$) coincides with a lognormal distribution, as described by Equation (27) [5]:

$$N({\delta }_i)={\displaystyle \frac{1}{\sqrt{2\pi }\ln \left(S\right)}}{\int }_{{\delta }_i}^{{\delta }_{i+1}}{\displaystyle \frac{1}{\delta }}\exp [-{\displaystyle \frac{1}{2}}({\displaystyle \frac{\mathrm{l}\mathrm{n}\delta -\mathrm{l}\mathrm{n}{\delta }_{50}}{\ln S}}{)}^2]d\delta$$

(27)

where ${\delta }_{50}$ and $S$ represent the geometric median diameter and geometric standard deviation, respectively. In the DELFIC model, the parameters ${\delta }_{50}$ = 0.15 μm and $S$ = 2.0 are adopted for air bursts, while ${\delta }_{50}$ = 0.407 μm and $S$ = 4.0 are employed for surface or near-surface bursts [5].

The radioactive particles are assumed to have a spherical shape and a constant density distribution, with a density of 2650 kg/m³ in this study according to the range of fallout particle densities (2600–2700 kg/m³) [3]. The weapon debris is primarily condensed in the equilibrium fireball at the initial time, similar to the observed radioactive particles from pure air bursts with diameters ranging from 0.1 to 4 μm by Nathans et al. [52]. Soil particles, observed by Spriggs et al. at the NTS site [53], are initially located on the surface with diameters ranging from 4 to 2800 μm. In this initialization, we randomly and evenly place particles representing weapon debris with four median diameters in the range of 0.1 to 4 μm within the fireball described in Section 2.4.2, while soil particles represented by six median diameters in the range of 4 to 2800 μm are placed on the surface.

In consideration of the high concentration of the radioactive particles in the nuclear cloud, the maximum number of parcels in the manual-injection type of OpenFOAM Lagrangian library is limited to 1 × 10⁶ for each particle group with uniform-size distribution. Therefore, the initial particle number concentration is about 1.64 × 10⁹ (7.19 × 10⁴) cm⁻³ based on Li et al. [3] for surface or near-surface bursts (air bursts), which is adopted to establish number-size distributions in conformity with the lognormal regarding the median particle diameters listed in

Table 2. The initialization of the number-size distributions of particles, indicated using the number concentration of initial particles in the surface or near-surface bursts (air bursts), that corresponding to the two values shown outside and inside the brackets, respectively.

Median Particle Diameter/μm	Number Concentration of Particles/cm−3	Median Particle Diameter/μm	Number Concentration of Particles/cm−3
0.1	2.56 × 108 (2.00 × 104)	250	4.23 × 105 (0)
0.5	6.62 × 108 (4.88 × 104)	500	2.74 × 103 (0)
1	3.00 × 108 (2.74 × 103)	1000	221.56 (0)
4	3.43 × 108 (222.79)	1500	12.10 (0)
50	8.11 × 107 (0.078)	2800	2.43 (0)

for ten size intervals.

Due to differences in solidification temperatures between carrier materials (i.e., weapon debris and soil) and fission products, fractionation occurs during particle formation. When the solidification temperature of the fission products is higher than that of the carrier materials, the fission products will condense with the carrier while it is still molten or condensing, resulting in a uniform distribution of radioactivity over the entire condensed particle, known as volume distribution. Conversely, when the fission products condense after the carrier particles have solidified, radioactivity is uniformly distributed over the particle’s surface, known as surface distribution. In fact, most radioactive particles exhibit characteristics of radioactivity distribution that fall between volume and surface distributions [54]. One advantage of the lognormal distribution is its ability to easily present the surface and volume distributions of particle radioactivity in a bimodal form, calculated as the weight average with coefficients $C_1$ and $C_2$ [55]. The surface and volume distributions conform to the lognormal as the number-size distribution. However, the logmean diameter $\mathrm{l}\mathrm{n}{\delta }_{50}$ of the surface and volume distributions can be obtained by substituting the $\mathrm{l}\mathrm{n}{\delta }_{50}$ in Equation (28) with $\mathrm{l}\mathrm{n}{\delta s}_{50}=\mathrm{l}\mathrm{n}{\delta }_{50}+2(\mathrm{l}\mathrm{n}S{)}^2$ and $\mathrm{l}\mathrm{n}{\delta v}_{50}=\mathrm{l}\mathrm{n}{\delta }_{50}+3(\mathrm{l}\mathrm{n}S{)}^2$, respectively. For air bursts, the geometric median diameter of particle activity is calculated to be 0.39 μm for surface distribution and 0.63 μm for volume distribution; for surface or near-surface bursts, it is calculated to be 19 μm for surface distribution and 130 μm for volume distribution. The geometric standard deviation is also S, resulting in a bimodal distribution form of activity-size given as

$$A\left({\delta }_i\right)=C_{1}As\left({\delta }_i\right)+C_{2}Av\left({\delta }_i\right)$$

(28)

where $A\left({\delta }_i\right)$ represents the proportion of particle activity in the diameter range (${\delta }_i{,\delta }_{i+1}$). $As\left({\delta }_i\right)$ and $Av\left({\delta }_i\right)$ are the surface and volume distributions of particle activity in the diameter range (${\delta }_i{,\delta }_{i+1}$), respectively. In this work, Bridgman’s conclusions are adopted, in which $C_1=0.32$ and $C_2=0.68$ [56]. Therefore, the radioactive particles’ activity-size distributions conform to the bimodal lognormal regarding the median particle diameters listed in

Table 3. The proportion of total activity with respect to the median particle diameter (the two values outside and inside the brackets represent the results of the surface or near-surface bursts and air bursts, respectively).

Median Particle Diameter/μm	Activity Percentage $\mathit{A}\left({\mathit{\delta }}_{\mathit{i}}\right)$/%	Median Particle Diameter/μm	Activity Percentage $\mathit{A}\left({\mathit{\delta }}_{\mathit{i}}\right)$/%
0.1	0 (1.04)	250	36.41 (0)
0.5	0.14 (44.23)	500	11.11 (0)
1	0.41 (34.53)	1000	6.68 (0)
4	4.03 (19.92)	1500	2.20 (0)
50	36.35 (0.28)	2800	1.75 (0)

The sum of the activity of each particle size interval is 100% in theory, but in the surface or near-surface bursts, 0.92% of the activity remains in particles sizes outside the diameter range of 0.1 to 2800 μm.

for ten size intervals.

The analysis from Table 3 indicates that the predominant activity in air bursts is attributed to smaller weapon fragments within the initial fireball. Conversely, for surface or near-surface bursts, the majority of activity comes from larger particles generated by contaminated soil, with a concentration on particle diameters ranging from 4 to 2800 μm occupying 98.53%. To conserve computing resources, particles smaller than 1 μm will be disregarded in simulations of surface or near-surface bursts, while particles larger than 4 μm will be ignored for air burst simulations. The geometry diagrams of the initial fireball with weapon debris and soil particles are provided in Figure S2a,b in the Supplementary Materials for the RANGER-Able and BUSTER-JANGLE-Sugar tests, respectively.

3. Test Cases of the Gas-Particle Nuclear Cloud

3.1. Numerical Discretization of Computational Cases

3.1.1. Numerical Schemes

The temporal terms in the Eulerian governing equations are discretized using the Eulerian scheme, while the pressure gradient terms are discretized using the Gauss linear scheme. The Laplacian terms are discretized using the Gauss linear corrected scheme. The convective terms of the RANS and standard k−ε turbulence equations are discretized by applying the Gauss upwind scheme and Gauss linear scheme for the divergence terms. The density of the fluid equation is solved using a Preconditioned Conjugate Gradient method (PCG) solver with a preconditioner of Simplified Diagonal-based Incomplete Cholesky (DIC), and the hydrodynamic pressure $p_{rgh}$ is solved using the Geometric Algebraic Multi-Grid (GAMG) with a smooth-solver of the Gauss–Seidel method. When utilizing the PIMPLE algorithm with an out-loop number identical to 1, it is transformed into the PISO algorithm, and an inner-loop number of 10 for velocity–pressure coupling is sufficient to achieve convergence by setting residuals to 10⁻⁸ and 10⁻⁶ for hydrodynamic pressure and momentum-turbulence equations, respectively.

The Lagrangian governing equations are solved using a two-way coupling approach with Eulerian integration to ensure the statistical significance of the discrete phase movement [57]. To maintain stability in the coupling method between Eulerian and Lagrangian phases, the maximum Courant number is set to 0.3 at each time step, preventing discrete parcels from traveling across more than one cell per time-step.

3.1.2. Mesh and Boundary Conditions

In this case, the computational resolution of the fluid phase and discrete particles is conducted in a pseudo two-dimensional geometry measuring 16 × 16 × 0.05 km (X × Z × Y), with a mesh consisting of 320 × 320 × 1 cells. The boundary conditions for the atmospheric boundaries are set as type patch for the side faces (Y-Z) and top face (X-Y), while the front and back faces (X-Z) are designated as type axisymmetric to simulate the three-dimensional cloud within a pseudo two-dimensional framework. The initialization of the atmosphere in Section 2.4.1 (Dirichlet boundary condition) dictates that the fluid field be patched accordingly, with the bottom face (X-Y) representing ground surface being set as type wall with zero-gradient pressure and no-slip shear boundary conditions (Robin boundary condition). The geometry and boundary conditions of this case are shown in Figure S3 from the Supplementary Materials.

The OpenFOAM “block Mesh” tool is used to generate hexahedral structured grids with uniform cell size, and the “check Mesh” command is employed to check the quality of the mesh. The mesh quality report is shown in

Table 4. The results of mesh quality check.

Parameters	Value
Max Aspect Ratio	1
Avg. Non-orthogonality	0
Max Skewness	7.07493 × 10−14
Boundary Patches	4
Cells	102,400
Overall Domain Bounding Box (m)	(−8000 −25 0), (8000 25 16,000)

.

3.2. Validation of Atmospheric Stratification

The pressure-based solver is utilized with the Eulerian model to simulate the fluid phase of atmospheric flow field, which is initialized with a stable atmosphere stratification under the assumption of atmospheric static equilibrium. The time-averaged numerical atmospheric parameters (the first 500 s) are compared with the U.S. standard atmosphere to validate the ability of this numerical model in handling atmospheric stratification. The numerical simulation error statistics are presented in Figure 3, and it is observed that the numerical solutions of the stable atmosphere demonstrate good performance of atmospheric parameters consistent with the U.S. standard atmosphere of 1976.

3.3. Validation and Analysis of Nuclear Cloud Rising

The upper (bottom)-edge height of the computational fluid phase of the stabilization cloud for the RANGER-Able test is 4257 (2907) m, compared to the observed height of 4224 (2841) m, as shown in Figure 4a. This indicates that the spatio-temporal shift of the computational cloud using the scalar field of temperature is highly consistent with the visible nuclear cloud observed from aircraft and theodolite measurements. For the BUSTER-JANGLE-Sugar test, the upper (bottom)-edge height of the computational fluid phase of stabilization cloud is 4485 (3035) m versus the observed height of 4580 (3416) m, and the diameter of 1900 m, versus the observed diameter of 1532 m, as shown in Figure 4b. The time-varying height trend for DELFIC’s upper (bottom)-edge shows slower rise development and lower stabilization heights for both edges, compared to the observed. Therefore, the computational method of this work demonstrates more accurate performance with respect to the nuclear cloud rise for both the RANGER-Able and BUSTER-JANGLE-Sugar tests.

In the fluid phase of stabilization, nuclear cloud profiles of the RANGER-Able and BUSTER-JANGLE-Sugar tests are depicted by scalar and vector fields, clearly visible in the temperature and density contour maps in Figure 5a,b,d,e. This is further demonstrated by the velocity magnitude of vector fields presented in Figure 5c,f, respectively. The pressure of the nuclear cloud has reached equilibrium with the ambient atmosphere, resulting in a balance between the pressure gradient force and drag force, as shown in Figure S4b,e from the Supplementary Materials. The turbulent kinetic energy fields at the stabilization time point are given in Figure S4c,d, indicating more intense turbulent movement at the vortex-ring of the cloud cap due to occurring entrainment. From Figure 5 and Figure S4, it can be observed that scalar and vector fields for air bursts and surface or near-surface bursts exhibit almost identical physical characteristics for stabilization clouds. However, differences in activity-height distribution observed from data cannot be reflected in fluid fields. Therefore, further study on the spatio-temporal distribution of radioactive particles within the discrete phase is necessary.

For air bursts, the computational particle distributions in the stabilization cloud of four particle size groups were initialized in the initial fireball, as shown in Figure S5 from the Supplementary Materials. This indicates that the RANGER-Able test represents a non-stem nuclear cloud without soil entrainment. It is evident that soil particles do not enter the fireball and become contaminated with radioactive materials. Therefore, the nuclear cloud loaded with activity from the RANGER-Able test has been carried by condensation particles entirely from weapon debris.

For surface or near-surface bursts, the computational particle distributions in the stabilization cloud of weapon debris and soil were initialized in the initial fireball and on the surface with eight particle size groups, as shown in Figure S6 from the Supplementary Materials. This indicates that the BUSTER-JANGLE-Sugar test represents a cap-connects-stem nuclear cloud while soil particles have been rolled up into the fireball. The activity of the surface or near-surface burst cloud is carried by soil particles pulled up from the surface, containing induced activity and forming relatively large radioactive particles with volatile radionuclides condensed onto their surfaces. Contrary to common belief, particles of 50 μm in size remain higher in the stabilization cloud than smaller particles in surface or near-surface bursts, as depicted in Figure S6. It can be inferred that differences in temperature and density of the initial fireball generate a pressure gradient force near the surface, significant enough to draw soil particles into the cloud. Consequently, particles with a middle diameter range possess greater upward kinetic energy than lighter weapon debris particles. However, when particle size exceeds a certain range, settling occurs during the upward process before reaching the stabilization cloud. Therefore, it can be concluded that particles of middle diameter range around 50 μm determine the maximum distribution height of activity, while dust stem and fallout of the stabilization cloud are formed by larger-sized particles exceeding 50 μm.

In general, when the height of burst (HOB) of a nuclear burst exceeds 54.86 W^0.4 (where W is the equivalent of the nuclear weapons in kiloton), local fallout may not pose a significant threat. However, global fallout should still be considered. The rise of the nuclear cloud from the RANGER-Able test was analyzed using a time-varying Eulerian–Lagrangian multiphase flow field, as depicted in Figure S7 from the Supplementary Materials. It is evident that no radioactive particles settled during the rise of the RANGER-Able test cloud, consistent with expectations for free-fallout conditions. Similarly, analysis of the nuclear cloud rise from the BUSTER-JANGLE-Sugar test revealed settlement of radioactive particles during its rise, resulting in serious local fallout near the burst center, as shown in Figure 6.

Figure 6 and Figure S7 indicate that radioactive particles are primarily concentrated near the vortex ring due to the anisotropic pressure gradient force of the fireball and turbulent flow field. Therefore, it can be inferred that the activity mainly gathers in the vortex-ring, similar to Jodin’s ring distribution around the symmetry of the cloud in the cloud cap [58] and DELFIC’s Gaussian distribution in the vertical and horizontal directions of the cloud [56].

3.4. Analysis of Stabilization Nuclear Cloud Activity Distribution

3.4.1. Activity Calculation of Stabilization Cloud

One minute after the bursts of nuclear weapons with the fission equivalent $Q$ kt, the total activity of $\beta$-rays is calculated as [54]:

$$A_{\beta }=3.7\times {10}^{21}Q\left(Bq\right),$$

(29)

Similarly, one minute after the bursts of nuclear weapons with the fission equivalent $Q$ kt, the total activity of $\gamma$-rays is [54]:

$$A_{\gamma }=1.48\times {10}^{21}Q\left(Bq\right),$$

(30)

Way and Wigner [59] put forward a semi-empirical treatment of the decay rate of fission products post the bursts of nuclear weapons, which can be expressed as follows [54]:

$$A=A_0({\displaystyle \frac{t}{t_0}}{)}^n.$$

(31)

In Equation (31), A₀ and A are the total activity sums of the $\beta$-rays and $\gamma$-rays at time t₀ (min) and t (min), respectively. The exponent term n is a constant value depending on the span of time; usually, the mean value of the term n is taken as −1.2, so-called the “t^−1.2 law” [59].

Therefore, when t₀ and A₀ were taken as the values one minute after the bursts of Equation (31), the total activity of the BUSTER-JANGLE-Sugar test (RANGER-Able) with a fission equivalent of 1.2 (1) kiloton at the time point of post-burst 4.32 (3.95) min of the observed stabilization cloud can be calculated based on Equations (29)–(31). The theoretical activity values of the fission products in the stabilization clouds are 8.9482 × 10²⁰ (9.9635 × 10²⁰) Bq, respectively.

3.4.2. Validation of the Activity Distribution

The statistical results of the vertical spatial distribution of activity in the stabilization cloud are crucial for the subsequent calculation of transport and diffusion modules in fallout prediction. Initially, the vertical spatial distribution of radioactive particles at a stabilization time is calculated through post-processing numerical computations. Additionally, the total activity of fission products can be determined using the method outlined in Section 3.4.1. Therefore, by combining the activity-size distribution of radioactive particles from Section 2.4.3, it is possible to obtain the activity-height distribution of the stabilization cloud in the vertical direction. Equation (32) can be utilized to acquire the activity-height distribution based on stratification using a “disk toss” method [5,12,19], with six layers being set for air bursts and surface or near-surface bursts according to AFPT (ARL Fallout Prediction Technique) model considerations [12].

$$A_j={\sum }_{i=1}^{i=n}{F}_i{\displaystyle \frac{X_{ij}}{X_i}}$$

(32)

where $n$ is the number of radioactive particle groups divided by size, $F_i$ represents the proportion of total activity in a particle group $i$ confirmed from Table 3, $X_i$ represents the total number of particles in a particle group $i$, $X_{ij}$ represents the number of particles in a particle group $i$ of the layer $j$, and $A_j$ represents the proportion of activity in the layer $j$.

Due to the calculation of the activity-height distribution being based on the proportion of particles rather than the actual number, it is advisable to minimize the number of particles in order to reduce calculation costs. Therefore, it is sufficient for the number distributions of particles to ensure that the interval with the least proportion (i.e., median particle diameter of 2800 μm) has an adequate representation of the spatial distribution statistical law within the stabilization cloud. As a result,

Table 5. Out of (In) the brackets: the computational number of particles (percentage of activity, %) at each layer of the stabilization cloud of the RANGER-Able test with regard to the median particle diameter (μm). The total particle number sums to 4 × 10⁴. The total percentage of activity sums to 99.72% based on Table 4, and the total activity equals 8.5207 × 10²⁰ according to Equations (24)–(26) as a computational activity value.

	Median Particle Diameter (μm)	P1	P2	P3	P4	Total
Layers	Height (AGL, m)	0.1 μm	0.5 μm	1 μm	4 μm
1	3919~4508	1891 (0.176)	4468 (7.2705)	248 (5.6191)	19 (3.0523)	6626 (16.1179)
2	3329~3919	2239 (0.2084)	5535 (9.0068)	304 (6.8879)	25 (4.0161)	8103 (20.1192)
3	2739~3329	7041 (0.6555)	17,178 (27.9527)	972 (22.0231)	80 (12.8516)	25,271 (63.4829)
4	0~2739	0 (0)	0 (0)	0 (0)	0 (0)	0 (0)
	Total	11,171 (1.04)	27,181 (44.23)	1524 (34.53)	124 (19.92)	4 × 104 (99.72)

presents the number (activity) spatial distribution of the stabilization cloud in the vertical direction for the RANGER-Able test, and

Table 6. Out of (In) the brackets: the computational number of particles (percentage of activity, %) at each layer of the stabilization cloud of the BUSTER-JANGLE-Sugar test with regard to the median particle diameter (μm). The total particle number sums to 7 × 10⁴. The total percentage of activity sums to 98.94% based on Table 4, and the total activity equals 7.693 × 10²⁰ according to Equations (24)–(26) as a computational activity value.

	Median Particle Diameter (μm)	P3	P4	P5	P6	P7	P8	P9	P10	Total
Layers	Height (AGL, m)	1 μm	4 μm	50 μm	250 μm	500 μm	1000 μm	1500 μm	2800 μm
1	3401~4081	0 (0)	0 (0)	1209 (4.39)	0 (0)	0 (0)	0 (0)	0 (0)	0 (0)	1209 (4.39)
2	2721~3401	682 (0.03)	690 (0.28)	919 (3.34)	694 (2.53)	314 (0.35)	194 (0.13)	21 (0.01)	1 (0)	3515 (6.66)
3	2041~2721	3081 (0.13)	3163 (1.27)	2435 (8.85)	672 (2.45)	713 (0.79)	398 (0.27)	46 (0.02)	0 (0)	10,508 (13.78)
4	1361~2041	6237 (0.26)	6147 (2.48)	1859 (6.76)	358 (1.3)	143 (0.16)	140 (0.09)	24 (0.01)	1 (0)	14,909 (11.06)
5	681~1361	0 (0)	0 (0)	169 (0.61)	0 (0)	0 (0)	46 (0.03)	13 (0.01)	0 (0)	228 (0.65)
6	10~681	0 (0)	0 (0)	132 (0.48)	0 (0)	0 (0)	24 (0.02)	38 (0.02)	0 (0)	194 (0.51)
7	0~10	0 (0)	0 (0)	3277 (11.91)	8276 (30.13)	8830 (9.81)	9198 (6.14)	4858 (2.14)	4998 (1.75)	39,437 (61.89)
	Total	10,000 (0.41)	10,000 (4.03)	10,000 (36.35)	10,000 (36.41)	10,000 (11.11)	10,000 (6.68)	5000 (2.2)	5000 (1.75)	7 × 104 (98.94)

does so for the BUSTER-JANGLE-Sugar test.

In the simulation of air bursts, the majority of activity is concentrated in the lower third of the cloud cap of the RANGER-Able test, as shown in Figure 7a. However, the total activity distribution in the AFPT model [12] is evenly divided into three layers with a percentage of 33.1% each, as presented in Figure 7b. This discrepancy is probably due to the fact that the activity-size distribution in AFPT is based on an inversion of sedimentation measurement results from two air burst tests, which do not accord with a lognormal distribution. In contrast, DELFIC [5] and Glasstone and Dolan [7] have taken into account fractionation caused by condensation at varying rates of refractory and volatile radionuclides [19], resulting in an activity distribution between surface and volume distributions for radioactive particles in practice.

On the contrary, DELFIC’s and Glasstone-Dolan’s activity distributions primarily focus on the local fallout resulting from surface or near-surface bursts, with limited consideration given to the global fallout yield from air bursts occurring above the free-fallout height. For example, in the model proposed by Glasstone and Dolan, radioactive particles with a log-median particle radius ranging from 20 to 200 μm account for 87% of the total activity, while particles with a log-median particle radius less than (over) 20 (200) μm contribute only 12 (1)% of the activity. Therefore, these distribution models are not applicable to air bursts.

In the simulation of surface or near-surface bursts, the majority of activity of 62.4% is concentrated on the surface, with approximately 25% and 12% remaining in the cloud cap and dust stem respectively for the BUSTER-JANGLE-Sugar test, as depicted in Figure 7c. However, the activity-height distribution of AFPT indicates that most of the activity remains in the cloud cap, with a significant portion staying in the dust stem. While the percentage of activity caused by settlement near the burst center is minimal, as shown in Figure 7d. This distribution differs significantly from the simulation result of the BUSTER-JANGLE-Sugar test. The activity distribution model proposed by Glasstone and Dolan shows relatively higher agreement with the simulation result, with discrepancies mainly occurring in the lower dust stem layer and percentage of settlement around 10%, as illustrated in Figure 7e. Similar to the conclusions drawn from the RANGER-Able test, the discrepancies between AFPT model and simulation results are attributed to errors resulting from the dose rate measurement inversion of AFPT, which is based on the SUNBEAM-Small Boy test does not adhere to a lognormal distribution pattern.

4. Conclusions and Discussion

In this study, we have developed a numerical simulation method for nuclear cloud rise based on the low Mach number model of gas-particle multiphase flow within the framework of CFD-DPM. The fluid phase is described by the RANS equations of the Eulerian system, while the discrete phase is represented by the DPM of the Lagrangian system in accordance with Newton’s second law. The coupled model utilizes an Eulerian–Lagrangian framework with a two-way coupling regime. The buoyantPimpleFoam solver and DPMFoam solver were utilized to construct the coupled multiphase solver in OpenFOAM codes, and sufficient numerical stability was achieved when using the PISO algorithm (the out-loop number of PIMPLE equals 1). To initialize the stratification of the stable atmosphere, the atmospheric static equilibrium equations were implemented. In comparison with the U.S. standard atmosphere of 1976, the computational atmospheric parameters (temperature, density, pressure, and vertical velocity component) within an altitude range of 16 km exhibit high consistency, with standard deviations of 0.0435, 632.943, 0.0028, and 0.079, respectively. However, this model has weaknesses in simulating wind direction at different altitudes due to high calculation costs associated with the 3D model and sensitivity of numerical algorithms to boundary conditions which can lead to divergence when initial velocity values are applied at lateral boundaries based on meteorological wind data from observations.

The DELFIC-Initialization module was utilized for the initialization process to create a thermal source with specific temperature, initial height and velocity, taking into consideration the applicability of low Mach number in this numerical algorithm. Additionally, lognormal distributions of particle number-size and activity-size were employed to simulate radioactive particles formed in the initial fireball and lofted soil, with fractionation described using surface and volume distribution of activity. Subsequent nuclear cloud simulations of two representative types of low-equivalent nuclear weapons tests demonstrated reasonable similarity to observations in terms of cloud rise rate, stabilization height and diameter. The cloud profile of the RANGER-Able test (BUSTER-JANGLE-Sugar) represents a no-stem (cap connects stem) nuclear cloud from air bursts (surface or near-surface bursts), consistent with the criterion of free-fallout height. The turbulent flow field in the cloud, as indicated by the turbulent kinetic energy and velocity magnitude of the vector field, suggests the presence of violent vortex movement during cloud rise. It can be confirmed that the primary cause of vortex formation is attributed to turbulent mixing resulting from the pressure gradient force originating from a thermal source, leading to the entrainment of cold and high-density air into the hot and low-density fireball until pressure equilibrium is reached.

In the simulation of the RANGER-Able test, the activity-height distribution of the stabilization cloud was found to be discrepant with the trisection of the AFPT, whereas activity mainly concentrated in the lower third of the cloud cap. In the simulation of the BUSTER-JANGLE-Sugar test, the activity-height distribution of stabilization cloud was found to be significantly discrepant with the AFPT, but it showed relatively higher agreement with Glasstone and Dolan’s model. Therefore, this numerical model demonstrates reasonable accuracy in predicting the activity-height distribution of stabilization clouds as a source item for fallout prediction models of low-equivalent range (less than 1 kiloton) in types of air bursts and surface or near-surface bursts. Therefore, the existing distributions of activity-height and activity-size have limitations when derived from dose rate measurements without considering the fractionation of radioactive particles.

DELFIC is widely used and highly regarded in predicting radioactive fallout from nuclear bursts, but as a one-dimensional model, it lacks a detailed description of the actual distribution of radioactivity within the cloud. This research has incorporated a semi-empirical equilibrium fireball model based on observations of nuclear clouds and activity-size distribution from DELFIC, achieving good verification results compared to other source models. However, there are still limitations in simulating large-equivalent nuclear bursts and complex meteorological conditions. Future work will consider larger equivalent nuclear bursts while entering the stratosphere, as well as shallow-buried underground nuclear bursts, which have received relatively little public research attention but will be increasingly important in the future.