A Closed-form Solution for Source-term Emission of Xenon Isotopes from Underground Nuclear Explosions

Abstract

Isotopic ratios of radioactive xenons sampled in the subsurface and atmosphere can be used to detect underground nuclear explosions (UNEs) and civilian nuclear reactors. Disparities in the half-lives of the radioactive decay chains are principally responsible for time-dependent concentrations of xenon isotopes. Contrasting timescales, combined with modern detection capabilities, make the xenon isotopic family a desirable surrogate for UNE detection. However, without including the physical details of post-detonation cavity changes that affect radioxenon evolution and subsurface transport, a UNE is treated as an idealized system that is both closed and well mixed for estimating xenon isotopic ratios and their correlations so that the spatially dependent behavior of xenon production, cavity leakage, and transport are overlooked. In this paper, we developed a multi-compartment model with radioactive decay and interactions between compartments. The model does not require the detailed domain geometry and parameterization that is normally needed by high-fidelity computer simulations, but can represent nuclide evolution within a compartment and migration among compartments under certain conditions. The closed-form solution to all nuclides in the series 131–136 is derived using analytical singular-value decomposition. The solution is further used to express xenon ratios as functions of time and compartment position.

1 Introduction

Isotopic ratios and ratio correlations of radioactive xenon isotopes $^{131\text {m}}$Xe, $^{133\text {m}}$Xe, $^{133}$Xe, and $^{135}$Xe were proposed as possible indicators for detecting and discriminating underground nuclear explosions (Bowyer et al. 1998; Saey and De Geer 2005; Kalinowski and Pistner 2006; Kalinowski et al. 2010; Saey et al. 2010; Kalinowski 2011; Galan et al. 2018). Because of their suitable half-lives, these xenon isotopes have reasonable signal strength with low background concentrations (Saey et al. 2010; Sloan et al. 2016; Bowyer 2020). In addition to radioactive xenons, stable xenons produced from the decay and ingrowth reactions have also been used for recovering underground-nuclear-explosion (UNE) source-term activities (Aregbe et al. 1997; Cassata et al. 2014). Due to the complexity of radioactive decay and ingrowth, the xenon source term from a UNE was approximated in early multiphase transport models using some fraction of the total cumulative yields without considering delayed and time-dependent behavior (Carrigan et al. 1996; Sun and Carrigan 2014). In contrast to using simplified source terms, full decay and ingrowth networks were modeled assuming radionuclide evolution in a closed and well-mixed (i.e., batch-mode) detonation cavity (Kalinowski et al. 2010, 2011; Yamba et al. 2016). However, more realistic and detailed models considering the post-detonation state of the cavity and its effect on xenon isotopic evolution do not retain the idealized batch mode (Kalinowski and Liao 2014; De Geer 2013; Carrigan et al. 2020b) and closed-system assumptions (Carrigan et al. 2020b) of these isotopic evolutionary models. Similarly, the analytical models developed here do not depend on the batch-mode or closed-system assumptions by following the model development of Carrigan et al. (2020b). In the cavity-melt partitioning model of Carrigan et al. (2020b), xenon production and transport from a UNE are described using partial differential equations (PDEs) governing multiphase flow and transport and ordinary differential equations (ODEs) governing radioactive decay and ingrowth (Carrigan and Sun 2014; Sun and Carrigan 2014, 2016; Lowrey et al. 2015; Sun et al. 2015; Jordan et al. 2015; Sloan et al. 2016; Bourret et al. 2019; Harp et al. 2019; Carrigan et al. 2020a). To mitigate the difficulty and high computational expense of solving the coupled PDEs and ODEs of transport and the complex decay/ingrowth networks necessary in evaluating details of the linkage between the physical and chemical evolution of the cavity, closed-form solutions may be developed. Such an approach may be critical to rapidly developing predictive source terms for interpreting xenon-ratio observations.

In the context of our cavity-melt partitioning model (Carrigan et al. 2020b), the detonation cavity following a UNE is composed of several interacting subsystems or compartments (e.g., cavity, rock-melt puddle, and host rock, as shown in Fig. 1). A compartment is not expressed by a physical volume or a specific geometry, but rather by a mass of a well-mixed radioactive material that behaves uniformly and homogeneously in the compartment. Each compartment may exchange materials (radioactive and stable isotopes) with other compartments. Time-dependent material exchanges simulated in the compartments are (1) thermally induced condensation and rainout of refractory xenon precursors from cavity to rock melt; (2) back flow of radioxenons from rock melt to cavity gas; and (3) leakage of cavity gases into a fractured containment zone.

The radioactive decay and ingrowth of series 131, 133, and 135 in a single-batch reactor have been modeled numerically (Kalinowski et al. 2010) and analytically (Sun et al. 2015; Sloan et al. 2016; Yamba et al. 2018) by solving the first-order ODEs. The resulting solution of the first-order reaction networks (Fig. 2a, c, and e) has been further used to calculate a multi-isotope ratio chart (MIRC, Kalinowski et al. 2010). In reality, this single-batch model is unable to address nuclide redistribution among subsystems or its effect on xenon isotopic evolution in the gas phase of the post-detonation cavity (De Geer 2013; Carrigan et al. 2020b). The initial fraction of radionuclide inventory is partitioned between the cavity and rock-melt phase constituting the puddle (Carrigan et al. 2020b). The rainout rate of the refractories from the condensing vapor significantly influences the radioxenon signature. To include the effects of rainout and back diffusion, Pannecoucke (2017) and Pili et al. (2017) developed a numerical model of two compartments (cavity and puddle). Carrigan et al. (2020b) further developed a semi-analytical model involving a numerical solution for the partitioning and an analytical solution to radioactive decay and ingrowth.

Although numerical methods are available for integrating coupled ODEs with a linear interaction between the cavity and puddle, the stiffness of ODEs owing to greatly differing decay rates may cause convergence problems or require extra computational effort to obtain the solution. For this reason, closed-form solutions are always useful. Since Bateman (1910) first derived an analytical solution of the first-order sequential reactions without considering branching (multiple products from a single reactant) and converging (multiple reactants for a single product) connections, many analytical solutions have been developed for sequential decay/ingrowth in single-batch systems (e.g., Cetnar 2006; Yuan and Kernan 2007; Yuan 2010; Zhou et al. 2015) and in studying radionuclide transport in the subsurface (e.g., Pigford et al. 1980; van Genuchten 1985; Clement et al., 1998; Sun et al. 1999; Slodic̆ka and Baláz̆ová 2008). As shown in Fig. 2, branching reactions are involved in series 131–135 and converging reactions are contained in all six series. The assumption of absolutely sequential first-order reactions has limited the application of the Bateman family of solutions to decay/ingrowth series in the England and Rider (1994) database. More generalized formulation of analytical singular-value decomposition has been derived for a wide range of decay and ingrowth networks (Sun et al. 2008, 2012; Sloan et al. 2016) and the exact solutions were further coupled with numerical solutions of multi-phase transport (Carrigan et al. 2020b).

In this paper, we develop a multi-compartment model of isotopic evolution and migration that includes cavity, puddle, and host-rock subsystems, in which subsystem volumes and geometry remain unknown or uncertain. The conceptual model is shown in Fig. 1. Some fraction of the total radionuclide inventory is partitioned initially between the cavity and melt puddle. Subsequently, rapid cooling and condensation of refractory precursors occur with migration of the condensing isotopes into the melt phase from the gas/vapor phase according to their condensation temperatures (Note: The example provided in Fig. 3b of Carrigan et al. (2020b) was primarily for illustrative purposes. More detailed analyses of the post-detonation thermochemical state are required for a more accurate assessment of the rainout phase of inventory partitioning). The volatile xenon isotopes produced in the puddle may possibly diffuse back to the cavity, although no efficient mechanism for back migration from the puddle has been identified (Carrigan et al. 2020b). Together with xenon isotopes produced in the cavity, xenon isotopes are further driven into host rock.

Fig. 1

Conceptual model of a multi-compartment system. Xenon isotopes and their precursors are initially deposited in the cavity and puddle. In, Sn, Sb, Te-m, and Te precipitate from the cavity to puddle according to their condensation temperatures and a prescribed rainout rate. The rainout of each isotope is assumed to be linear with its cavity concentration when cavity temperature drops to its condensation threshold (Carrigan et al. 2020b). Xenon isotopes generated by chain reactions in the puddle may potentially diffuse back to the cavity. The accumulated xenon isotopes from both back diffusion and local chain reactions in the cavity will further transport to the host rock and atmosphere. R, D, and S in the cavity represent rainout, back diffusion, and seepage, respectively. The interaction between the puddle and host rock is neglected due to low puddle permeability

We derived analytical solutions of all radionuclide concentrations (in Fig. 2) in three compartments (cavity, puddle, and host rock) with interactions (rainout, back diffusion, and seepage). The derived solutions were further used to (1) formulate xenon-ratio solutions for studying the MIRC under various conditions and (2) connect to numerical models of multiphase flow and transport.

Fig. 2

Radioactive decay/ingrowth networks of series a 131, b 132, c 133, d 134, e 135, and f 136. Numbers on arrows denote branching factor (%). Red and green icons represent radioactive and stable xenon isotopes. The decay networks are based on England and Rider (1994)

2 Analytical Solution Development

In this section, we present ODEs of the multi-compartment system coupled with radionuclide decay/ingrowth networks and derive closed-form solutions for given network structures. The subsurface domain is composed of three compartments (gas-filled cavity, melt puddle, and host rock) that are connected by inter-compartment transport (precursor rainout, xenon back diffusion, and seepage from gas-filled cavity to host rock).

2.1 Multi-compartment Model

The mass change of all radionuclides in a given series in the three-compartment system described in Fig. 1 can be expressed as

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}\left[ \begin{array}{l} {\mathbf{c}}_{\mathrm{c}} \\ {\mathbf{c}}_{\mathrm{p}} \\ {\mathbf{c}}_{\mathrm{h}} \\ \end{array}\right] =\left[ \begin{array}{ccc} {\mathbf{A}}-{\mathbf{R}}-{\mathbf{S}} &{} {\mathbf{D}} &{} {\mathbf{O}} \\ {\mathbf{R}} &{} {\mathbf{A}}-{\mathbf{D}} &{} {\mathbf{O}} \\ {\mathbf{S}} &{} {\mathbf{O}} &{} {\mathbf{A}} \\ \end{array}\right] \left[ \begin{array}{l} \mathbf{c}_{\mathrm{c}} \\ {\mathbf{c}}_{\mathrm{p}} \\ {\mathbf{c}}_{\mathrm{h}} \\ \end{array}\right] , \end{aligned}$$

(1)

where ${\mathbf{c}}_{\mathrm{c}}$, ${\mathbf{c}}_{\mathrm{p}}$, and ${\mathbf{c}}_{\mathrm{h}}$ are the mass vectors in cavity, puddle, and host rock, respectively. ${\mathbf{A}}$ is a matrix of the first-order decay and ingrowth rates, ${\mathbf{R}}$ is the rainout-rate matrix from cavity to puddle, ${\mathbf{D}}$ is the back-diffusion coefficient matrix, and ${\mathbf{S}}$ is the seepage-rate matrix. ${\mathbf{O}}$ denotes $n \times n$ zero matrices and n is the number of radionuclides in the series.

2.2 First-order Reactions

A typical system of coupled first-order reactions is described using ODEs. For example, the rate equations of the 131 decay chain (Fig. 2a) in the cavity can be expressed as

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}\left[ \begin{array}{l} c_{\mathrm{c,1}} \\ c_{\mathrm{c,2}} \\ c_{\mathrm{c,3}} \\ c_{\mathrm{c,4}} \\ c_{\mathrm{c,5}} \\ c_{\mathrm{c,6}} \\ c_{\mathrm{c,7}} \\ c_{\mathrm{c,8}} \\ \end{array}\right] =\left[ \begin{array}{cccccccc} -k_1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{}0 &{}0\\ \alpha _1 k_1 &{} -k_2 &{} 0 &{} 0 &{} 0 &{} 0 &{}0 &{}0\\ 0 &{}\alpha _2 k_2 &{} -k_3 &{} 0 &{} 0 &{} 0 &{}0 &{}0\\ 0 &{} 0 &{} \alpha _3 k_3 &{}-k_4 &{} 0 &{} 0 &{}0 &{}0\\ 0 &{} 0 &{} \alpha _4 k_3 &{}\alpha _5 k_4 &{}-k_5 &{} 0 &{}0 &{}0\\ 0 &{} 0 &{} 0 &{}\alpha _6 k_4 &{}\alpha _7 k_5 &{}-k_6 &{} 0 &{}0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{}\alpha _8 k_6 &{}-k_7 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{}\alpha _8 k_6 &{} \alpha _{10}k_7 &{} -k_8\\ \end{array}\right] \left[ \begin{array}{l} c_{\mathrm{c,1}} \\ c_{\mathrm{c,2}} \\ c_{\mathrm{c,3}} \\ c_{\mathrm{c,4}} \\ c_{\mathrm{c,5}} \\ c_{\mathrm{c,6}} \\ c_{\mathrm{c,7}} \\ c_{\mathrm{c,8}} \\ \end{array}\right] = A~{\mathbf{c}}_{\mathrm{c}}, \end{aligned}$$

(2)

where $c_{c,i},~i=1,\ldots ,8$, denote the mass of $^{131}$In, $^{131}$Sn, $^{131}$Sb, $^{131\text {m}}$Te, $^{131}$Te, $^{131}$I, $^{131\text {m}}$Xe, and $^{131}$Xe in cavity, $k_i,~i=1,\ldots ,8$, are their decay rates, and $\alpha _i,~i=1,\ldots ,10$, are branching factors (see Fig. 11). Similarly, the decay and ingrowth ODEs are formulated based on England and Rider’s (1994) database for all series (131–136) in all three compartments.

2.3 Source-term Activities

Immediately following a detonation, radionuclides are partitioned between a gas/vapor phase in the cavity and rock-melt phase in the puddle. During the subsequent rapid cooling of the vapor phase, the vaporized iodine precursors (In, Sn, Sb, Te-m, and Te) condense out, mixing with the melt phase in the puddle compartment, while some small fraction of the xenon isotopes produced in the puddle may potentially be transported back to the cavity (Carrigan et al. 2020b). However, no significant short-term transport mechanism beyond chemical diffusion across the surface of the puddle into the gas phase has been identified (Carrigan et al. 2020b). Reintroducing xenon from the puddle back into the gas phase by diffusion represents a negligible to small contribution to the total gas available for release into the containment zone (Carrigan et al. 2020b). The precursor rainout is conceptualized as a thermally dependent linear process assuming Newtonian cooling in the cavity (Olsen 1967; Carrigan et al. 2020b) or using user-provided temperature profiles. In general, the rainout and back diffusion are expressed as

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}\left[ \begin{array}{l} {\mathbf{c}}_{\mathrm{c}} \\ {\mathbf{c}}_{\mathrm{p}} \\ \end{array}\right] =\left[ \begin{array}{cc} -{\mathbf{R}} &{} {\mathbf{D}} \\ {\mathbf{R}} &{}-{\mathbf{D}} \\ \end{array}\right] \left[ \begin{array}{l} {\mathbf{c}}_{\mathrm{c}} \\ {\mathbf{c}}_{\mathrm{p}} \\ \end{array}\right] , \end{aligned}$$

(3)

where ${\mathbf{R}}$ and ${\mathbf{D}}$ are diagonal matrices. For example, the rainout and back-diffusion matrices of series 131 (Fig. 2a) can be specified as

$$\begin{aligned} \begin{array}{lllllllll} &{} ^{131}\text {In} &{} ^{131}\text {Sn} &{} ^{131}\text {Sb} &{} ^{131\text {m}}\text {Te} &{} ^{131}\text {Te} &{} ^{131}\text {I} &{} ^{131\text {m}}\text {Xe} &{} ^{131}\text {Xe} \\ {\mathbf{R}}={\mathbf{diag}} &{} ({{\mathcal {H}}}(T_1)r_1, &{} {{\mathcal {H}}}(T_2)r_2, &{} {{\mathcal {H}}}(T_3)r_3, &{} {{\mathcal {H}}}(T_4)r_4, &{} {{\mathcal {H}}}(T_4)r_5, &{} 0, &{} 0, &{} 0)\\ {\mathbf{D}}={\mathbf{diag}} &{} (0, &{} 0, &{} 0, &{} 0, &{} 0, &{} 0, &{} d_7, &{} d_8),\\ \end{array} ~~~~ \end{aligned}$$

(4)

where $T_i,~i=1,~2,~3,~4$ are condensation temperatures of In, Sn, Sb, and Te (Te-m), respectively, $r_i,~i=1,\ldots , 5$, are rainout rates, $d_i,~i=7,8$ are back diffusion rates, and ${{\mathcal {H}}}(T_i)$ is the Heaviside step function with ${{\mathcal {H}}}(T_i)=0,~ \forall ~T > T_i$ and ${{\mathcal {H}}}(T_i)=1,~\forall ~T \le T_i$.

It also assumed that xenon isotopes produced in the cavity and diffused from puddle to 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. The mass change is expressed as

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}\left[ \begin{array}{l} {\mathbf{c}}_{\mathrm{c}} \\ {\mathbf{c}}_{\mathrm{h}} \\ \end{array}\right] =\left[ \begin{array}{cc} -{\mathbf{S}} &{} {\mathbf{O}} \\ {\mathbf{O}} &{} {\mathbf{S}} \\ \end{array}\right] \left[ \begin{array}{l} {\mathbf{c}}_{\mathrm{c}} \\ {\mathbf{c}}_{\mathrm{c}} \\ \end{array}\right] , \end{aligned}$$

(5)

where ${\mathbf{S}}$ is the diagonal matrix representing the seepage rate of xenon isotopes from cavity to host rock. For series 131, ${\mathbf{S}}$ in Eq. (5) is specified as

$$\begin{aligned} \begin{array}{lllllllll} &{} ^{131}\text {In}&{} ^{131}\text {Sn} &{} ^{131}\text {Sb} &{} ^{131\text {m}}\text {Te} &{} ^{131}\text {Te} &{} ^{131}\text {I} &{} ^{131\text {m}}\text {Xe} &{} ^{131}\text {Xe} \\ {\mathbf{S}}={\mathbf{diag}} &{} (0, &{} 0, &{} 0, &{} 0, &{} 0, &{} 0, &{} s_7, &{} s_8),\\ \end{array} \end{aligned}$$

(6)

where $s_i,~i=7,8,$ are seepage rates of $^{131\text {m}}$Xe and $^{131}$Xe.

2.4 Solution Methods

Vectors ${\mathbf{c}}_{\mathrm{c}}$, ${\mathbf{c}}_{\mathrm{p}}$, and ${\mathbf{c}}_{\mathrm{h}}$ are arranged in ${\mathbf {c}}$ to make the system matrix a lower triangular matrix,

$$\begin{aligned} \frac{\mathrm {d}{\mathbf {c}}}{\mathrm {d}t} = \mathscr {A}~{\mathbf {c}}, \end{aligned}$$

(7)

where ${\mathbf {c}}$ is rearranged mass vector and $\mathscr {A}$ is lumped matrix. For example, ${\mathbf {c}}$ and $\mathscr {A}$ can be expressed specifically for series 131 as shown in Appendix 1 ($\mathscr {A}_{131}$).

The matrix $\mathscr {A}$ in Eq. (7) can be further diagonalized as

$$\begin{aligned} \mathscr {A}={{\mathcal {S}}}\Lambda {{\mathcal {S}}}^{-1} \end{aligned}$$

(8)

(Sun et al. 2012, 2015; Sloan et al. 2016), where $\Lambda$ is an eigenvalue matrix of $\mathscr {A}$, ${{\mathcal {S}}}$ is the matrix whose columns are the eigenvectors of $\mathscr {A}$, and ${{\mathcal {S}}}^{-1}$ is the inverse matrix of ${{\mathcal {S}}}$.

The matrix of eigenvectors and its inverse are given analytically as

$$\begin{aligned} S_{ij} = \sum _{\zeta = 1}^{b_{ij}}\left\{ \alpha _{\zeta } \left[ -{{\mathcal {F}}}(j,i)\right] \prod _{l=1}^{n_{ij}-1}{{\mathcal {F}}}\left[ m(l),j\right] \right\} , \end{aligned}$$

(9)

and

$$\begin{aligned} S_{ij}^{-1} = \sum _{\zeta = 1}^{b_{ij}}\left\{ \alpha _\zeta \prod _{l=1}^{n_{ij}}{{\mathcal {F}}}\left[ m(l),i\right] \right\} , \end{aligned}$$

(10)

where

$$\begin{aligned} {{\mathcal {F}}}(i,j)=\frac{\lambda _i}{\lambda _i-\lambda _j}, \end{aligned}$$

i is the current species index and j an ancestor of species i, m(l) is the species index of the l-th ancestor of i, $n_{ij}$ is the generation number from species j to i, $\zeta$ is the branch-number index, and $b_{ij}$ is the number of branches that connect species i and j (Sun et al. 2012). $\lambda _i,~\forall i=1,2,\ldots ,N$ are the eigenvalues of $\mathscr {A}$ and are exactly the same values of diagonal elements of matrix $\mathscr {A}$. N is the length of the vector ${\mathbf {c}}$.

Substituting Eq. (8) in Eq. (7) and defining ${\mathbf{a}}={{\mathcal {S}}}^{-1}{\mathbf {c}}$,

$$\begin{aligned} \frac{\mathrm {d}{\mathbf{a}}}{\mathrm {d}t} = \Lambda ~{\mathbf{a}}. \end{aligned}$$

(11)

Each ODE in Eq. (11) is independent of the others, but with the first-order decay rate $\lambda$. In other words, the coupled ODEs (Eq. 7) in terms of ${\mathbf {c}}$ by ingrowth terms are decomposed into N ODEs with the same formulation,

$$\begin{aligned} \frac{\mathrm {d}a_i}{\mathrm {d}t} = \lambda _i~a_i, \quad \forall i=1,2,\ldots , N, \end{aligned}$$

(12)

where $a_i$ is the transformed concentration of $c_i$ and the solution is

$$\begin{aligned} a_i = a^0_i \exp (\lambda _i t), \quad \forall i=1,2,\ldots , N. \end{aligned}$$

(13)

Using ${\mathbf {c}}= {{\mathcal {S}}}~{\mathbf{a}}$, the solution is

$$\begin{aligned} {\mathbf {c}}(t) = {{\mathcal {S}}} \left[ \begin{array}{cccc} \exp (\lambda _1t) &{} 0&{} \cdots &{} 0\\ 0 &{} \exp (\lambda _2t) &{}\cdots &{} 0 \\ \vdots &{} \vdots &{} \ddots &{} \vdots \\ 0 &{} \cdots &{}0 &{} \exp (\lambda _nt) \\ \end{array}\right] {{\mathcal {S}}}^{-1}~{\mathbf {c}}^0, \end{aligned}$$

(14)

where ${\mathbf {c}}^0$ is the vector of initial concentrations.

For a given temperature history in the cavity (e.g., Fig. 3a, see Carrigan et al. 2020b), the solution (14) can be applied in 5 periods ($\textcircled{~i},~i=1,\ldots ,5$) that are calculated inversely using the temperature history. For example, if a simple Newtonian cooling is assumed (Olsen 1967), the rainout time is expressed

$$\begin{aligned} t = \frac{1}{\beta }\ln \frac{T_0-T_\infty }{T_{\mathrm{c}}-T_\infty } \end{aligned}$$

(15)

where $T_0$ [$^{\circ }$C] and $T_\infty$ [$^{\circ }$C] are the initial cavity temperature and ambient temperature, $T_{\mathrm{c}}$ [$^{\circ }$C] is the condensation temperature, and $\beta$ is the temperature decay factor. In each period, matrices of ${{\mathcal {S}}}$, ${{\mathcal {S}}}^-$, and $\Lambda$ are formulated specifically due to the temperature-dependent ${{\mathcal {H}}}(T_i)$ in (4). Then, the solution in the entire simulation time can be concatenated as

$$\begin{aligned} {\mathbf{c}}\big |_0^t = \left[ {\mathbf{c}}\big |_0^{t(\text {Sn})}~~ {\mathbf{c}}\big |_{t(\text {Sn})}^{t(\text {In})}~~ {\mathbf{c}}\big |_{t(\text {In})}^{t(\text {Sb})}~~ {\mathbf{c}}\big |_{t(\text {Sb})}^{t(\text {Te})}~~ {\mathbf{c}}\big |_{t(\text {Te})}^t\right] . \end{aligned}$$

(16)

The rainout time of precipitating isotopes for the temperature history (Carrigan et al. 2020b) is compared with their half-lives in Fig. 3b (9.375 × 10$^{-7} < t_{1/2}<$ 3.26 days). If t(Te) $\ll$ min($t_{1/2}$), the single-period solution (14) is applicable without considering ${{\mathcal {H}}}$. If t(Sn) $\gg$ max($t_{1/2}$), the effect of rainout becomes negligible. Otherwise, the multi-period solution (16) is used.

Fig. 3

a Plot of condensation time versus cavity temperature. b Comparison between half-lives of all precipitating isotopes and their condensation time for the illustrative temperature history given by Carrigan et al. (2020b)

2.5 Implementation and Simulations

The above solution was implemented and applied to series 131–136 with the independent yields, half-lives, and branching factors listed in Table 1 (England and Rider 1994). The matrices of series 131–136 with the compartment interactions were then derived for the sake of completeness. The matrix $\mathscr {A}_{131}$ of series 131 is provided as an example in Appendix A. The network of series 131 represented by $\mathscr {A}_{131}$ is partitioned into independent subsystems for formulating ${{\mathcal {S}}}$ and ${{\mathcal {S}}}^{-1}$. The technical details are given in Sun et al. (2015) and Sloan et al. (2016).

Closed-form solutions (Eq. 14) with derived ${{\mathcal {S}}}$ and ${{\mathcal {S}}}^{-1}$ were implemented in MATLAB (MathWorks, 2000). Six sets of inputs, including decay rates, branching factors, rainout rates, back-diffusion coefficients, seepage rates, and temperature history in the cavity are needed to obtain the solutions. Decay rates and branching factors are known from England and Rider (1994), while the remaining four inputs are UNE-specific, with optimized values that can possibly be used to recover the UNE history. The solutions of radionuclide mass in term of time and compartments were saved as outputs. Among those outputs, the mass (or radioactivity) fluxes of xenon isotopes from the cavity to host rock were the source term of numerical models for simulating subsequent multiphase flow and transport.

3 Results and Analyses

In this section, we simulated the evolution of all 43 nuclides (listed in Table 1) in the cavity, melt puddle, and host-rock compartments (see Fig. 1), calculated xenon fluxes from the cavity to host rock, and calculated xenon isotopic ratios for discriminating nuclear explosions and civilian nuclear applications.

3.1 Evolution of Radionuclide Inventories

When $r_i = 1\times 10^{-3}$ [s$^{-1}$], $d_i = 1\times 10^{-5}$ [s$^{-1}$], $s_i=1\times 10^{-7}$ [s$^{-1}$], initial inventory fraction in the puddle is 0.75 for 1 kt U235 fission with the temperature history as shown in Fig. 3a, and the xenon inventories in the cavity, puddle, and host rock were simulated as shown in Fig. 4. As shown in Fig. 4c, due to the differences in decay networks and half-lives, xenon isotopes appear in host rock at different times. This explains why it is rare to get all four xenon signals ($^{131\text {m}}$Xe, $^{133\text {m}}$Xe, $^{133}$Xe, and $^{135}$Xe) from a single gas sampling to generate a MIRC (Kalinowski et al. 2010). In this synthetic case, the common time window to detect the four xenons was about 12 days. The schedule of onsite inspection (gas sampling) should be made considering the detection window of desired xenon isotopes.

Inventories of xenon isotopes and their precursors in series 131 and 135 in the cavity and puddle, for example, are shown in Fig. 5. Unlike their precursors, xenon isotopes appear late, as delayed signatures. Both radioactive and stable xenons rise to peak and decay to baseline in the cavity and puddle, owing to diffusion from puddle to cavity ($D > 0$) and seepage from cavity to host rock ($S >0$). When $D = 0$ and $S = 0$, the stable xenons reach their respective asymptotic levels in the cavity and puddle. Similar figures can also be plotted and analyzed for other series to show the evolution of radionuclide activities in the three compartments.

Fig. 4

Inventories of xenon isotopes in the a cavity, b melt puddle, c host rock. d Detection time windows of $^{\mathrm{131m}}$Xe, $^{\mathrm{133m}}$Xe, $^{\mathrm{133}}$Xe, and $^{\mathrm{135}}$Xe with a detection limit of 1 mBq m$^{-3}$. Note that the mass of $^{\mathrm{134m}}$Xe and $^{\mathrm{135m}}$Xe is too low to be shown in c. Subplots a, b, and c share the legend in c

Fig. 5

Inventory evolution of xenon isotopes and their precursors in series 131 and 135 in the cavity and puddle

If the post-shot temperature in the cavity drops extremely fast below t(Te) (i.e., 988 $^{\circ }$C), the single-step solution (14) with constant rainout rates provides the same answer as the solution (16) with temperature-dependent rainout rates. However, when the condensation time is comparable with the half-lives of precipitating isotopes (Fig. 3b), a discrepancy between these two models is observed in xenon inventories in host rock (Fig. 6a). The relative error can reach 3% for this synthetic case.

Fig. 6

a Comparison of xenon inventories in host rock simulated using temperature independent (solid lines) and dependent rainout rates (dashed lines). b Relative error versus time

3.2 Xenon Fluxes to Host Rock

To capture the delayed signature of xenon isotopes, it is necessary to consider decay and ingrowth networks in reactive transport models (Carrigan et al. 2016). When all precursors are included, the simulation becomes computationally expensive. For example, to simulate $^{\mathrm{131m}}$Xe, $^{\mathrm{133m}}$Xe, $^{\mathrm{133}}$Xe, and $^{\mathrm{135}}$Xe, at least 22 species have to be considered in multiphase flow and transport models (Sun and Carrigan 2016; Sloan et al. 2016). To couple all six series with transport, 43 species must be considered. In general, those xenon precursors that are mainly contained in the cavity and puddle do not participate in the gas transport in host rock and atmosphere. If xenon fluxes are well defined by source-term activities (rainout, back diffusion, and seepage), transport in the host rock can be simplified as xenon transport with a total of 10 species (bold-faced species in Table 1) for all 6 series.

Figure 7 shows the source-term fluxes of all xenon isotopes from the cavity to host rock as functions of time. In this case, the fluxes of $^{\mathrm{131m}}$Xe, $^{\mathrm{133m}}$Xe, $^{\mathrm{133}}$Xe, and $^{\mathrm{135}}$Xe, that are used in the MIRC calculation peak at 9.67, 1.98, 3.04, and 0.375 days, respectively. Because $^{\mathrm{134m}}$I, the single precursor of $^{\mathrm{134m}}$Xe gets no reaction yield (see Fig. 2d) and only 2% of $^{\mathrm{134m}}$I decays to produce $^{\mathrm{134m}}$Xe, the $^{\mathrm{134m}}$Xe flux to host rock becomes negligible. The flux of the stable xenons fall off due to the reduced productivity of precursors in both the cavity and puddle, while the fluxes of radioactive xenons decline even faster with additional decay.

Fig. 7

Source-term fluxes of a all xenon isotopes and of b four xenons for developing the MIRC. Note that $^{\mathrm{134m}}$Xe flux is negligible

3.3 Isotopic Ratios and Event Time

Isotopic ratios of xenons and daughter species have been used for dating nuclear events inversely (De Geer 2012, 2013; Yamba et al. 2016). It has been stated that source-term partitioning, such as rainout from the cavity to the melt puddle, affects isotopic ratios and zero-time (time of detonation) estimates. The closed-form solution derived in this manuscript can be used to address the dependency of isotopic ratios on source-term activities and time: (1) isotopic ratios can be simulated using the forward model for known (or assumed) source-term activities as a function of time; (2) the zero time of a UNE can be estimated inversely for known or assumed source-term activities and measured isotopic ratios; and (3) source-term activities can be recovered using measured isotopic ratios and known event times.

As an example, we quantified source-term activities using measurements of the $^{\mathrm{131m}}$Xe/$^{\mathrm{133}}$Xe ratio from atmospheric samples detected in April 2013, first in Japan and several days later in Russia in April 2013 associated with the February 2013 declared DPRK UNE at the Punggye-ri Nuclear Test Site, with a known detonation time (Ringbom et al. 2014). Assuming a fast-decaying temperature in the cavity, the single-period solution (14) is used for interpreting the DPRK 2013 data. The optimization began with an educated guess on the range for each parameter (including the initial inventory fraction in the melt puddle, R, D, and S, defined in Eqs. (4) and (6). Using Latin hypercube (McKay et al. 1979; Tong 2005), 2000 parameter sets were generated within the defined parameter space. The $^{\mathrm{131m}}$Xe/$^{\mathrm{133}}$Xe ratio was calculated and compared with observed data at the time of measurement for all 2000 parameter sets. System parameters were identified by minimizing the sum of the absolute difference between modeled and measured values,

$$\begin{aligned} \begin{aligned}&\underset{\mathbf{X}}{\text {Minimize}}&\mathscr {F}\left( {\mathbf{X}}\right)&= \sum _{i=1}^{n_{\mathrm{E}}} \bigg | {{\mathcal {R}}}\left( {\mathbf{X}}, t_i \right) -\hat{{\mathcal {R}}}(t_i)\bigg | \\&\text {Subject to}&{\mathbf{X}} \in \Omega , \\ \end{aligned} \end{aligned}$$

(17)

where ${\mathbf{X}}$ is the vector of system parameters (i.e., the initial inventory fraction in the puddle, rainout rates, back-diffusion coefficients, and seepage rates) to be optimized, $\Omega$ represents the parameter space, ${{\mathcal {R}}}\left( {\mathbf{X}}\right)$ and $\hat{{\mathcal {R}}}$ are modeled and measured isotopic ratios (in terms of ${\mathbf{X}}$), and $n_{\mathrm{E}}$ is the number of measurements. The best match between model result and data was sorted from all sampling-based simulations and used to verify if it is close to the global minimum optimized using different methods such as the shuffled complex evolution (SCE) program developed by Duan et al. (1992).

The optimized results indicate that (a) 12.18% of inventory was initially deposited in the melt puddle, (b) precursors In, Sn, Sb, Te-m, and Te precipitated rapidly from cavity to melt puddle with rates of 2.256 × 10$^{-5}$, 7.521 × 10$^{-4}$, 7.400 × 10$^{-4}$, 6.236 ×  10$^{-5}$, 7.847 × 10$^{-4}$ [s$^{-1}$], respectively, (c) back-diffusion rates of Xe and Xe-m from the melt puddle to cavity 4.927 × 10$^{-7}$ [s$^{-1}$] is much lower than the precipitation rates of their precursors, and (d) the seepage rate from the cavity to host rock, 6.839 × 10$^{-7}$ [s$^{-1}$], is slightly higher than the back-diffusion rate. Although the SCE method is an efficient metaheuristic approach, it does not necessarily guarantee the true history of source-term activities due to limited data sets and a high-dimensional parameter space. Additional measurements of other isotopic ratios (including stable xenons) may improve our understanding of source-term activities. When the initial inventory fraction was assumed to be 60% (Borg 1975) in the melt puddle, the rainout rates were optimized as 5.856 × 10$^{-4}$, 2.129 × 10$^{-3}$, 2.3286 × 10$^{-3}$, 5.033 × 10$^{-5}$, 5.6390 × 10$^{-5}$ [s$^{-1}$] for In, Sn, Sb, Te-m, and Te, respectively. The corresponding back-diffusion coefficient and seepage rate are 3.739 × 10$^{-7}$ and 1.590 × 10$^{-6}$ [s$^{-1}$].

The MIRC (Kalinowski et al. 2010) was simulated using the optimized source-term activity parameters. The dashed blue line in Fig. 8a shows the fully fractionated (instantaneously produced or independent yield) and unfractionated (solid blue curve, England and Rider 1994) evolutionary correlations of $^{133\text {m}}$Xe/$^{131\text {m}}$Xe and $^{135}$Xe/$^{133}$Xe. Similar to Carrigan et al. (2020b), the correlation curves in the cavity, puddle, and host rock deviate from the ideal England and Rider curve. The MIRC depends on sample locations and source-term activities. Figure 8b shows time versus the $^{\mathrm{131m}}$Xe/$^{\mathrm{133}}$Xe ratios in the cavity, melt puddle, and host rock. As described by Carrigan et al. (2016), Japanese (three data sets) and Russian (two data sets) observations correspond to signals from the host rock (tunnel) and cavity, respectively, where the bulk of gases released in each case were emitted over short periods, as occurs when venting the air in a tunnel. The objective function (17) was specified as

$$\begin{aligned} \mathscr {F}\left( {\mathbf{X}}\right) = \sum _{i=1}^3 \bigg | {{\mathcal {R}}}_{\mathrm{h}}\left( {\mathbf{X}}, t_i \right) -\hat{{\mathcal {R}}}_{\mathrm{h}}(t_i)\bigg | +\sum _{i=1}^2 \bigg | {{\mathcal {R}}}_{\mathrm{c}}\left( {\mathbf{X}}, t_i \right) -\hat{{\mathcal {R}}}_{\mathrm{c}}(t_i)\bigg |, \end{aligned}$$

where the subscripts c and h represent the cavity and host-rock compartments. The best fit ($\mathscr {F}$ = 0.2087) shows both ratios in the cavity (magenta) and host rock (cyan) agree with the data.

Fig. 8

Compartmentally well-mixed evolution paths of isotopic ratios in cavity, puddle, and host rock. a MIRCs in the cavity, puddle, and host rock. The solid lines and the dashed curves represent the fully fractionated (independent yield) and unfractionated ratio correlations, respectively. Black squares are the correlation data of civilian nuclear applications (Topin et al. 2020). The dashed green line is established for distinguishing civilian applications from UNEs (Kalinowski et al. 2010). The time dimension is represented by dotted/dashed blue lines. b Time versus $^{\mathrm{131m}}$Xe/$^{\mathrm{133}}$Xe. Black circles and green squares represent Japanese and Russian observations (Ringbom et al. 2014), respectively

When source-term activities remain uncertain, the zero-time of a UNE can be estimated using uncertainty quantification. As shown in Fig. 9, the range of uncertainty of zero-time is a function of input uncertainties from source-term activity parameters. If the Japanese observation (with an averaged ratio of 0.1467) is assumed to be venting a tunnel following gas buildup from host rock, the time range was between 44 and 63 days according to the uncertainty quantification, while the true measurement was between 55 and 56 days. If the Russian data (1.32) is from the cavity directly, the uncertainty of zero-time was between 17 and 78 days. On the other hand, because true measurement was at 60 days, the signal was likely from the cavity. Although Japanese and Russian observations that are located in cyan and magenta zones (Fig. 9) indicate the possible story of source-term activities, the uncertainty of event time in terms of $^{\mathrm{131m}}$Xe/$^{\mathrm{133}}$Xe ratio can be reduced when additional data become available.

Fig. 9

Estimation of UNE detonation time based on $^{\mathrm{131m}}$Xe/$^{\mathrm{133}}$Xe observation. The red rectangle (in the cyan area) and white box (in the magenta area) represent Japanese and Russian data, respectively. The slight color difference in the cyan area indicates overlapped uncertainties

3.4 Xenon Source Term in Subsurface Transport

With the estimated xenon activities in host rock, the transport from host rock to ground surface may be simplified using uncoupled PDEs. For example, one-dimensional gas transport with constant and homogeneous velocity and dispersion is expressed as

$$\begin{aligned} \frac{\partial c}{\partial t} = {{\mathcal {D}}}\frac{\partial ^2c}{\partial x^2}-v\frac{\partial c}{\partial x} -\lambda c \end{aligned}$$

(18)

with boundary conditions (in host rock)

$$\begin{aligned} v c_0 = v c +{{\mathcal {D}}}\frac{\partial c}{\partial x}, \quad x=0, \end{aligned}$$

(19)

and

$$\begin{aligned} c = \frac{\partial c}{\partial x} =0, \quad x=\infty , \end{aligned}$$

(20)

and initial condition

$$\begin{aligned} c = 0, \quad x\in [0~\infty ], \quad t=0. \end{aligned}$$

(21)

In Eq. (18), c [Bq m$^{-3}$] is the radioactivity concentration, t [d] is time, ${{\mathcal {D}}}$ [m$^2$ d$^{-1}$] is the dispersion coefficient, v [m d$^{-1}$] is the flow velocity, x is distance, and $\lambda$ [d$^{-1}$] is decay rate.

The analytical solution for Eq. (18) with constant boundary condition $c_0$ was derived by Cleary and Ungs (1978) as

$$\begin{aligned} c(x,t)= & {} \frac{c_0v^2}{4\lambda {{\mathcal {D}}}} \left\{ 2\exp \left[ \frac{xv}{{\mathcal {D}}}-\lambda t \right] \cdot \mathrm {erfc} \left( \frac{x+vt}{2\sqrt{{{\mathcal {D}}}t}} \right) \right. \nonumber \\&+\left( \frac{u}{v}-1\right) \exp \left[ \frac{x}{2{{\mathcal {D}}}}\left( v-u\right) \right] \cdot \mathrm {erfc} \left( \frac{x-ut}{2\sqrt{{{\mathcal {D}}}t}} \right) \nonumber \\&- \left. \left( \frac{u}{v}+1\right) \exp \left[ \frac{x}{2{{\mathcal {D}}}}\left( v+u\right) \right] \cdot \mathrm {erfc} \left( \frac{x+ut}{2\sqrt{{{\mathcal {D}}}t}} \right) \right\} , \end{aligned}$$

(22)

where

$$\begin{aligned}u=\sqrt{v^2+4\lambda {{\mathcal {D}}}}.\end{aligned}$$

Because Eq. (18) is a linear partial differential equation, the principle of superposition can be used to simulate concentrations in the time-dependent source-boundary conditions. Then, the solution is expressed as

$$\begin{aligned} c(x,t) = c_1 {{\mathcal {A}}}(x,t) + \sum _{i=1}^{n_{\mathrm{T}}}{{\mathcal {H}}}(t-t_i)(c_{i+1}-c_i){{\mathcal {A}}}(x,t-t_i), \end{aligned}$$

(23)

where $c_i,~i=1,2,\ldots , n_{\mathrm{T}}$ are step-wise boundary concentrations, ${{\mathcal {A}}}$ represents the analytical solution of Eq. (22), and ${{\mathcal {H}}}(t)$ is the Heaviside step function with ${{\mathcal {H}}}(t)=0,~ t\le 0$ and ${{\mathcal {H}}}(t)=1,~t>0$.

Given ${{\mathcal {D}}} = 3.82$ m$^2$ d$^{-1}$, $v = 7.64$ m d$^{-1}$, $x = 100$ m, and boundary conditions (Fig. 10a), the radioactivities of xenon isotopes are predicted and compared with numerical results as shown in Fig. 10b. Although the analytical and numerical solutions were almost identical, the CPU time used for calculating analytical solutions was about 0.15% the time used for the numerical (NUFT) model (Nitao 1998; Hao et al. 2012). Because the NUFT model included xenon isotopes only, the CPU time was about 27% of the NUFT model, considering full decay and ingrowth networks. Although the analytical solution (Eq. 23) is limited to transport with constant and uniform velocity and dispersivity, numerical models (e.g., NUFT models) of transport can be used in place of ${{\mathcal {A}}}$ for simulating transport in realistic systems.

Fig. 10

Prediction of xenon radioactivities at ground surface. a Radioactivities in host rock near the cavity. b Radioactivities at ground surface calculated using the closed-form solution (solid) and NUFT model (dashed)

4 Conclusions and Discussion

We derived a closed-form solution to a multi-compartment system coupled with complex radioactive decay/ingrowth networks. The system was composed of three subsystems: the cavity, melt puddle, and host rock. The mass transport between compartments (subsystems) was conceptualized using rainout, back diffusion, and seepage. We extended previously published solutions from three standard chains (131, 133, and 135) in a well-mixed cell to six chains (131–136) in multi-compartment and interactive systems (the cavity, puddle, and host rock). The solution serves as a potential tool to address (1) radioactive decay and growth, (2) mass transport between subsystems, and (3) discrimination between UNEs and civilian nuclear applications.

It was assumed that precursors of iodine isotopes were well contained in the cavity and puddle and only xenon isotopes transport to the host rock and atmosphere. Time-dependent xenon concentrations in the cavity and host rock that are calculated using the derived solution can be considered as boundary conditions for non-isothermal and multiphase transport in the subsurface and atmosphere. In other words, iodine precursors may be eliminated in PDEs of subsurface transport, and the corresponding simulations become more efficient. As an example, the analytical solution of single-species transport in one dimension with a constant boundary condition (Eqs. 19 and 20, also see Cleary and Ungs 1978) was extended to predict xenon transport with time-dependent boundary conditions generated using the closed-form solution of multi-compartment source-term activities (Eq. 23). The same coupling scheme is applicable to high-fidelity transport models using numerical codes such as NUFT (Nitao 1998; Hao et al. 2012).

Uncertainties are always involved in monitored data from UNEs. It is difficult, if even possible, to quantify all physical processes with inadequate data. The lumped compartment model, with a reduced number of parameters, is a preferred approach to addressing the mass migration and full-scale radioactive decay and ingrowth. Although this paper does not focus on interpreting a particular UNE, the solution presented can also be used to study xenon isotopic data for understanding source-term activities and to examine hypotheses about event timing. Although we demonstrated the application of our solution for plotting a 4-isotopic MIRC, the solution can be used to calculate other ratios of xenon isotopes including stable xenons. The zero-time of UNEs has been mainly estimated using $^{140}$La/$^{140}$Ba produced from a sequential decay and ingrowth chain (De Geer 2012, 2013; Yamba et al. 2016). To incorporate other available isotopes produced from complex decay/ingrowth networks (with branching and converging reactions) for the zero-time estimate, our method can be applied to obtain closed-form solutions to any decay network described by England and Rider (1994). The derived solution is limited to the assumption that rainout rates are first order, but temperature dependent as demonstrated by Carrigan et al. (2020b). When the post-shot temperature in cavity drops extremely fast compared to half-lives of precipitating isotopes, the rainout can be approximated using a constant (temperature independent) rate. When the temperature drops much slower than radioactive decay, the impact of rainout becomes negligible. For all scenarios between those extreme conditions, temperature-dependent rainout model is preferred. We anticipate that the first-order and temperature-dependent rainout model can be an adequate solution for approximating and verifying high-fidelity models.

Appendix 1: State Variables and System Matrix of Series 131

The multi-compartment structure of series 131 was specified as shown in Fig. 11 with the corresponding vector of state variables ${\mathbf {c}}$ and system matrix $\mathscr {A}$ expressed (Table 1).

Fig. 11

Radioactive-decay chain of series 131 (Fig. 2a) with mass exchange among cavity, puddle, and host rock. $\alpha$ denotes branching factors. $r_i,~i=1,\ldots , 5$, $d_i,~i=7,~8$, and $s_i,~i=7,~8$ are, respectively, rainout, diffusion, and seepage rates. Only precursors of iodine isotopes are involved in the rainout from the cavity to melt puddle and only xenon isotopes can diffuse from puddle to cavity and seep from cavity to host rock. Letters ‘c’, ‘h’, and ‘p’ in circles represent cavity, puddle, and host rock, respectively. The number represents species as defined in Table 1

Table 1 Independent yield for U235f 100 fissions, half-life time ($t_{1/2}$), and branching factor ($\alpha$) (England and Rider 1994)

(24)

The decomposed reaction matrix $\Lambda$ in (11) is expressed,

(25)