Implications of Underground Nuclear Explosion Cavity Evolution for Radioxenon Isotopic Composition

Isotopic ratios of radioxenons sampled in the atmosphere or subsurface can be used to verify the occurrence of an underground nuclear explosion (UNE). Differences in the half-lives of radioactive xenon precursors and their decay-chain networks produce different time-dependent concentration profiles of xenon isotopes allowing isotopic ratios to be used for tracking UNE histories including estimating the time of detonation. In this study, we explore the potential effects of post-detonation cavity processes: precipitation of iodine precursors, gas seepage, and prompt venting on radioxenon isotopic evolution which influences UNE histories. Simplified analytical models and closed-form solutions yielding a potentially idealized radioactive decay/ingrowth chain in a closed and well-mixed system typically have limited application by not including the partitioning of the radionuclide inventory between a gas phase and rock melt created by the detonation and by ignoring gas transport from the cavity to host rock or ground surface. In reality, either subsurface transport or prompt release that is principally responsible for gas signatures violates the closed-system (or batch-mode) assumption. A closed-form solution representing time-dependent source-term activities is extended by considering the cavity partitioning process, slow seepage, and/or prompt release of gases from the cavity and applied to realistic systems.


Introduction
Monitoring for underground nuclear explosions (UNEs) uses atmospheric detection of radionuclide gases released by either seepage or prompt venting at the ground surface of the test site following a UNE (Saey & De Geer, 2005;Carrigan et al., 2016). Advances in the technology for measuring radioxenon concentrations have enhanced the detection of both UNE signatures and the global background of radioxenons produced by civilian facilities including nuclear reactors and medical isotope production plants (Le Petit et al., 2015;Achim et al., 2016;Gueibe et al., 2017;Haas et al., 2017;Hoffman & Berg, 2018). Isotopic ratios of radioxenons (e.g., ½131m 133m 133 135 Xe) have been assessed as indicators of UNEs (Bowyer et al., 1998;Saey & De Geer, 2005;Kalinowski & Pistner, 2006;Kalinowski et al., 2010;Saey et al., 2010;Kalinowski, 2011;Galan et al., 2018). Different half-lives of these xenon isotopes originating from three different decay networks can produce detectable signals (Saey et al., 2010;Sloan et al., 2016;Bowyer, 2020). The radioactivity of each xenon isotope is characterized by a unique dependence on time as a result of its decay chain and half-lives of precursors reaching peak activity at a distinctive time. In the absence of physical partitioning of the radioactive inventory in the postdetonation cavity, seepage, and prompt venting, radioxenon profiles will only reflect the mechanism of radioactive decay and ingrowth in a closed and well-mixed (i.e., batch-mode) system (De Geer, 2013;Kalinowski & Liao, 2014;Yamba et al., 2016). A realistic UNE system, in general, may be composed of a cavity and melt puddle ( Fig. 1) with possible post-detonation partitioning of radionuclides from the gas-filled cavity into the underlying puddle of molten rock (R), gas seepage into the fractured zone of containment (S), and/or prompt venting of larger volumes of cavity gas (V) either of which is potentially capable of producing detectable xenon signatures at ground surface or downwind. In the context of our cavity-melt partitioning model (Carrigan et al., 2020;Sun et al., 2021), the detonation system following a UNE is described by several interacting subsystems or compartments (e.g., cavity, rock-melt puddle, and host rock). Each compartment may exchange radionuclides with other compartments. Although the slow seepage from the cavity to host rock was considered by Sun et al. (2021) in the analytical solution of generalized systems, prompt venting during a specified period of time has not been included. In the fully coupled numerical model of multiphase transport and radionuclide decay and ingrowth Bourret et al., 2021), the transport of xenon isotopes produced from chain reactions in the cavity was accurately described. However, the signature of promptly vented xenons was approximated using the xenon evolution in the cavity directly . Since the prompt venting alters mass balance in the cavity, a full-scale mass exchange needs to be considered among all compartments. In this paper, we derived a closed-form solution to a full-scale model considering time-dependent mass exchanges among interactive compartments: (1) thermally induced condensation (rainout) of refractory iodine precursors from cavity to rock melt, (2) seepage of cavity gases into a fractured containment zone, and (3) prompt venting of cavity gases to ground surface. The back diffusion from melt puddle to cavity is not considered because of its low gas flux compared to other exchange processes.

Model and Solution Development
In this section, we propose ordinary differential equations (ODEs) of the multi-compartment system coupled with radionuclide decay/ingrowth networks and source-term activities, and derive closed-form solutions for given network structures (Fig. 2, see also Bourret et al., 2021). The study domain is composed of multiple compartments (gas-filled cavity, melt puddle, and host rock) connected by intercompartment mass exchange (precursor rainout, seepage from gas-filled cavity to host rock, and/or prompt venting from the cavity to ground surface).
The mass change of all radionuclides in each decay chain in the multi-compartment system described in Fig. 1 can be expressed as d dt where c c , c p , c h , and c v are the mass vectors in cavity, puddle, host rock, and prompt release, respectively. A is a matrix of the first-order decay and ingrowth rates, R is the rainout-rate matrix from cavity to puddle, S is the seepage-rate matrix, and Q is the prompt venting matrix. O denotes n Â n zero matrices and n is the number of radionuclides in the decay chain. The residual heat of a UNE is translated into a temperature distribution in the cavity and surrounding area. The post-shot temperature distribution can be estimated with a Hugoniot phenomenological equation of state model, as described by Butkovich (1974): where T [ C] is the temperature, T 0 [ C] is the initial temperature prior to detonation, R [m] is the radial distance from the working point (in a spherical coordinate system), Y [kt] is the energy-equivalent yield, and q [g cm À3 ] is the bulk density of the rock. A simple Newtonian cooling is assumed (Olsen, 1967) for calculating the rainout time as where T m [ C] is the post-shot temperature in the cavity, T c [ C] is the condensation temperature, b is the temperature decay factor, and t T=2 is the half-time of cavity temperature. Using (2) and assuming a cavity radius of 20 m, the cavity temperature can be calculated as shown in Fig. 3a for an uncompacted soil (1.2 g cm À3 ) and two denser rock types with three different rock densities. The corresponding rainout times of In, Sn, Sb, and Te-m/Te are shown in Fig. 3b.

Precipitation of Precursors
The post-detonation precipitation of the refractory radioactive inventory involves a time-dependent mass transfer from the vapor-filled cavity to the rock melt or puddle at the bottom of the cavity. This is a complex problem that is dependent on several parameters (e.g., cavity pressure, temperature, yield, cavity-gas composition, etc.) and for the purposes of this paper only a somewhat general and simplified model is considered here. It is assumed that during the rapid cooling of the post-detonation vapor phase as described in Sect. 2.1, the vaporized iodine precursors (In, Sn, Sb, Te-m, and Te) condense out, mixing with the melt phase in the puddle compartment (Pili et al., 2017;Carrigan et al., 2020). Since no significant short-term diffusion mechanism has been identified from the puddle to cavity (Carrigan et al., 2020), any back diffusion from the puddle to the cavity gas is neglected. The rainout is expressed as where R is a diagonal rainout-rate matrix. For where T i ; i ¼ 1; 2; 3; 4, are the assumed condensation temperatures of In, Sn, Sb, and Te-m/Te, respectively, r i ; i ¼ 1; . . .; 5, are rainout rates, and HðT i Þ is the Heaviside step function with

Seepage of Isotopic Xenon
It is also assumed that xenon isotopes produced in the cavity transport away to the host rock, and the mass loss from the cavity, due to excess pressure, is linearly proportional to the mass in the cavity, which is consistent with the ideal gas law. The mass change is expressed as where S is the diagonal matrix representing the seepage rates of xenon isotopes from cavity to host rock. For decay chain 131, S in Eq. (6) is specified as where s 1 and s 2 are seepage rates of 131m Xe and 131 Xe.

Prompt Venting of Isotopic Xenon
Prompt release (venting) driven by pressure gradient (between the cavity and ground surface) is described by one-dimensional advection and expressed in the form of first-order reaction where v [m s À1 ] is the velocity of venting flow, L [m] is the venting distance, q [s À1 ] is the mass loss rate due to prompt venting, and Pðt 0 ; t 1 Þ [-] is a boxcar function defined as The mass change due to the prompt venting is described in matrix format as d dt For decay chain 131, Q in Eq. (9) is specified as where t 0 and t 1 [s] are the starting and ending time of prompt venting.

Solution Development
Equation (1) is arranged as a lower triangular matrix where A is lumped matrix diagonalized as where K is an eigenvalue matrix of A, S is the matrix whose columns are the eigenvectors of A, and S À1 is the inverse matrix of S. The analytical formulation of S and S À1 is described in Sun et al. (2012Sun et al. ( , 2015Sun et al. ( , and 2021. Defining a ¼ S À1 c, Eq. (11) is expressed in terms of a as with a transformed initial condition a 0 ¼ S À1 c 0 Vol. 180, (2023) Implications of where c 0 is the vector of initial concentrations. Each ODE in Eq. (13) is independent of the others, but with the first-order decay rate k. In other words, the coupled ODEs (Eq. 11) in terms of c by ingrowth, rainout, seepage, and prompt venting, are decomposed into N ODEs with the same formulation, where a i is the transformed concentration of c i and the solution of a i is Using c ¼ S a, the solution of c is The closed-form solution can be verified for absolutely sequential reactions by comparing the Bateman equation (Bateman, 1910) and for an integrated system (with all source-term activities) by comparing a numerical solution using ode113s in MATLAB (MathWorks, 2000).

Results and Analyses
In this section, we simulate the evolution of xenon isotopes in the cavity, melt puddle, host-rock, and promptly vented gas, and xenon isotopic ratios for demonstrating the impact of source-term activities on xenon signatures. Assuming zero initial inventory in the puddle for 1 kt U235 fission, four models are conceptualized as listed in Table 1 for demonstrating the effect of source-term processes on xenon inventories and isotopic ratios of activities. A reference model (#0) representing the idealized case of isotopic evolution in a cavity determined entirely by the radioactive decay chain (Kalinowski et al., 2010;Kalinowski, 2011;Yamba et al., 2016) is used to compare models #1, 2, 3 for demonstrating the individual effect of rainout, seepage, and prompt venting, respectively. Model #4 includes all four processes (Table 1)

Evolution of Xenon Isotopes
Assuming the cavity temperature decaying from 3000 to 20 C with temperature half-life of 600 [s] and a prompt venting between 1 Â 10 5 [s] (1.16 [d]) and 1.864 Â 10 5 [s] (2.16 [d]), the xenon inventories of Model #4 in the cavity and promptly vented gas were simulated as shown in Fig. 4.
As shown in Fig. 4a-c, the xenon inventories in the cavity are clearly reduced due to the prompt venting between 10 5 and 1.864 Â 10 5 [s]. 135 Xe in vented gas peaks during the venting period due to its short half-life, while other three xenons peak after the one-day venting (Fig. 4d). A late-time venting may not contain detectable signature 135 Xe because of its short half-life. The inventory profiles shown in Fig. 4 can be used to calculate isotopic ratios of xenon radioactivities. Fig. 5a shows the Multi-Isotope Ratio Chart (MIRC, the correlation of 133m Xe/ 131m Xe and 135 Xe/ 133 Xe) in the cavity and host rock with individual and combined source activities. While Table 1 Deterministic models with multiple physical processes 1400 Y. Sun et al. Pure Appl. Geophys. rainout shifts the England and Rider curve (blue) to the right (see also Carrigan et al., 2020), seepage (cyan) and prompt venting (dashed magenta) move the curve back to the left on the plot. The timedependent ratio of 131m Xe to 133 Xe (in the cavity) shown in Fig. 5b indicates that rainout brings the England and Rider (ER) curve (Model #0, blue) down toward the x-axis while both seepage and prompt venting lift the ratio value above the blue curve under the ER condition. The comparison between observed data and Model #0 curve reveals if the system is closed with rainout or open with possible seepage or prompt venting. The observed data of 131m Xe/ 133 Xe ratio from atmospheric samples associated with the February 2013 declared DPRK UNE at the Punggye-ri Nuclear Test Site, with a known detonation time (Ringbom et al., 2014), may indicate the precursor precipitation occurred in a leaky system. Three Japanese data sets that are located below the ER curve (blue curve) and close to the 131m Xe/ 133 Xe ratio in host rock (green curve) reflect a possible seepage from the cavity to host rock since those samples were taken off site. Both 131m Xe/ 133 Xe ratios in cavity (red curve) and host rock (green curve) behave similarly after 0.2 d in this example. If those samples (circles) were taken from the cavity directly, the red curve might be used to interpret the early-time precipitation. Two Russian data points (green squares) located above the ER curve may indicate a possible venting from the cavity.

Effect of Rainout on Xenon Activity Ratios
In absence of seepage and prompt venting, the rainout effect on both MIRC and 131m Xe/ 133 Xe ratio in the cavity is simulated using Model #1 and shown in Fig. 6a, b for three different rainout rates [1 Â 10 À4 , 1 Â 10 À3 , 1 Â 10 À2 ] [s À1 ] under Newtonian cooling condition from 3000 to 20 C with temperature half-life of 600 [s]. As shown in Fig. 6a, a fast rainout moves the correlation to the right away from the ER curve (batch-mode closed without rainout). Similarly, the fast rainout results in a low ratio of 131m Xe/ 133 Xe (Fig. 6b). In addition to the rainout-rate effect, the temperature profile also affects MIRC and 131m Xe/ 133 Xe ratio. Fig. 6c, d show the MIRC and 131m Xe/ 133 Xe ratio for t T=2 ¼ [300 600 1200] [s] and r ¼1 Â 10 À3 [s À1 ]. All possible rainout scenarios (for given r and t T=2 ) move both MIRC and 131m Xe/ 133 Xe ratio to the right side of their ER curves.

Effect of Seepage and Prompt Venting on Xenon Activity Ratios
When s is assumed to be [1 Â 10 À6 , 1 Â 10 À5 , 1 Â 10 À4 ] [s À1 ] without considering rainout and prompt venting, the MIRC and 131m Xe/ 133 Xe ratio in the cavity are simulated using Model #2 as shown in Fig. 7a, b. A high seepage rate moves the MIRC correlation away from the ER curve towards the left and produces 131m Xe/ 133 Xe ratio higher than the curve subject to idealized batch-mode conditions. Similarly to seepage effect, the time of prompt venting (Model #3) deviates the MIRC (Fig. 7c) and 131m Xe/ 133 Xe ratio (Fig. 7d) curves away from the standard ER curves after the prompt venting starts. Although a different venting rate q alters xenon inventories in the cavity, it does not change the MIRC and 131m Xe/ 133 Xe ratio. The time scale of seepage and venting is referred to OTA (1989).

Discussion and Conclusion
For predicting time-dependent isotopic evolution, the idealized standard model (ER) assuming a batchmode closed cavity may not be appropriate when iodine precursors precipitate from the cavity to the melt puddle, xenon isotopes transport to host rock, or prompt venting occurs. We extended our closed-form solution to a multi-compartment system coupled with complex radioactive decay/ingrowth networks, temperature-dependent precursor precipitation, cavity gas seepage, and specified prompt venting. The standard model (#0 in Table 1) serves as a reference with deviations potentially indicating precipitation of refractories, xenon seepage, and/or prompt venting have occurred and if the cavity is closed. For example, data points of cavity samples located above the ER curve in 131m Xe/ 133 Xe(t) (Fig. 5b) indicate possible seepage or prompt venting while those data points of the cavity samples below the ER curve demonstrate an early-time precipitation of iodine precursors. Only early-time processes (i.e., precipitation, seepage, and prompt venting) under thermal conditions are considered in our model and analytical solution. Late-time driving force of gas transport (e.g., barometric pumping) in fractured rock is not considered. A prompt venting that can be physically defined by operational activities (e.g., starting and ending time), uncertain hydrodynamic and geological conditions (e.g., depth of burial, rock properties, venting pathways) is described by venting rate q with a boxcar function (Pðt 0 ; t 1 Þ in Eq. 8). Although the Newtonian cooling is used in this paper for describing the fast cooling processes (from 3000 to 988 C), the true temperature profile either from an onsite monitoring or hydrodynamic modeling can be inputted to calculate precursor precipitation schedule. Then, the difference between observed xenon activity ratios and those simulated using ER model may be possibly used to calculate yield number inversely and recover historic source-term activities. If the cavity temperature remains above 988 C for a long time (e.g., t T=2 [ 100 [d]), the precursors of iodine will have enough time to decay away for complete xenon Vol. 180, (2023) Implications of Underground Nuclear Explosion Cavity production in the cavity and the effect of rainout is eliminated. Then, the analytical solution with r ' 0 applies for seepage and prompt venting only. The generalized analytical solution can be used as a handy and robust tool for conducting sensitivity analyses and verifying high-fidelity computer codes with a low computational cost. In this paper, we have only considered precipitation of the refractory precursors of iodine. We have assumed most or all iodine remains in the cavity gas phase and does not precipitate into the underlying melt or onto hot cavity surfaces. Such an assumption is probably justifiable for gas seepage rather than for rapid venting. If rapid venting occurs, it is expected that isotopes of iodine may be released with radioxenon, which will have the effect of shifting the MIRC toward the ER ingrowth curve and away from the rainout curve. Future analytical models will attempt to quantify this effect of iodine loss during the period of rapid venting.