A Model Study of a Homogeneous Light-Water Thorium Reactor

This work presents a computational study of a 232Th-based homogeneous light-water reactor. Thorium reactors have been proposed as an alternative to the uranium fuel cycle since they exploit both the availability of thorium and its ability to afford fissile uranium isotopes by a sequence of neutron captures. Besides 233U, as a result of the neutron captures, a significant amount of 234U (36.3%) and 6.46% of 235U are formed in the reactor under study. More importantly, the proposed simulation points out the possibility of a continuous withdrawal of the uranium isotopes without compromising the criticality and the power output of the reactor. This withdrawal affords the fissile material for the startup of reactors other than the first one, which requires a one-time only limited amount of fissile material. The significant molar fraction of the 234U (0.17) in the extracted fuel does not pose a limitation on weapon proliferation, as a consequence of its high fission cross section for high-energy neutrons.


Introduction
Thorium reactors have been studied as an alternative to uranium-based reactors mainly because thorium is more widely available than uranium [1,2]. In a thorium reactor, the fissile material 233 U is obtained by the non-fissile 232 Th by the sequence of neutron (n) capture and β − decays Besides the availability of fuel, other advantages have been suggested, such as the reduced production of heavy actinides, in particular 239 Pu that carries the possibility of nuclear weapon proliferation. The fuel cycle of a thorium based reactor hinges on the fission of 233 U and does not contain any 238 U, the precursor of 239 Pu. Furthermore, the occurrence of 232 U by the reaction (n, 2n) on 233 U is thought to impede the usage of 233 U for nuclear weapons [3] because of the relatively short half-life time (t 1/2 ) of 232 U (68.9 yr) that entails the presence of its hard-gamma emitting daughter 208 Tl (t 1/2 = 3.053 min). It is not the purpose of this work to discuss the extensive literature on this topic; interested readers may refer to the excellent review by Ault [4]. Most existing suggestions for the configuration of thorium reactors use molten salts as medium and graphite as a moderator [5,6]. Furthermore, fast reactors simulations have been suggested, either based on molten salts [7,8], lead cooled [9], or water. Among the latter we have a boiling water reactor (BWR) study [10], and pressurized water reactor (PWR) studies based on mixed uranium-plutonium [11] or uranium-thorium [12] fuels. A reduction in power density leading to a reduced temperature of the fuel elements was predicted for a mixed oxide UO 2 -ThO 2 configuration compared to pure UO 2 [13]. Furthermore, a reduction of burnable neutron poisons was pointed out for Small Modular Reactors operating on UO 2 -ThO 2 [14]. However, to the author's knowledge, no study has been published concerning a light-water PWR based solely on thorium (except for the startup phase, requiring a fissile component). The fissile material is being introduced into the reactor under study only on the first two days, i.e., the startup phase. After that, the reactor produces enough fissile material to support the production of energy. More importantly, the proposed simulation points out the possibility of a continuous withdrawal of the uranium isotopes without compromising neither the criticality nor the power output of the reactor. Thus, the proposed reactor could be exploited as a uranium source for the existing uranium-based reactors in the case the uranium supply would become problematic.
After the few thorium based reactors that have been operating in the past were discontinued, a few programs exist, like the 300 MW Advanced Heavy Water Reactor planned in India, making use of an external feed of plutonium. The THOREX (Thorium-Uranium 233 Extraction) reprocessing cycle has been discussed in [15].
This work develops a model study for a homogeneous thorium reactor in light water with an external feed that needs to be operating only for the first two days after startup. The requirement for the startup of a thorium reactor implies the use of fissile material [16], since 232 Th is not itself fissile. To this purpose, both 235 U and 239 Pu have been used. Instead, this work explores the possibility of using the 233 U bred from a thorium reactor operating with continuous reprocessing that would also avoid the necessity of 233 Pa extraction [17].
Since the homogeneous reactor would breed a mixture of 233 U, 234 U, and 235 U, it would be possible in principle to extract the bred uranium isotopes to start up other reactors. Although the abundance and role of the bred uranium isotopes is prominent in thorium reactors, it has non been the subject of many investigations so far. This possibility is being investigated in this paper along with the feasibility of maintaining criticality in the presence of a significant sink for all uranium isotopes. The consequence of the composition of the withdrawn uranium on weapon proliferation is briefly examined. The collected fuel could thus be used to start a new thorium reactor without resorting to uranium enriched in the isotope 235 U. In particular, we also want to assess whether the amount of non-fissile 234 U in the bred fuel would impair criticality during the startup phase. The reactor would operate at a nearly stationary state with a steady input of 232 Th.
In the sequence of neutron captures, beta decays, and fissions taking place in the environment of a reactor core, there are many nuclides representing branching points, because alternative processes to neutron capture, such as the (n, 2n) reaction, are taking place. In this study, we only take into consideration the main sequence from 232 Th to 235 U. We include the production of 235 U because we deem it important to quantify its formation in view of the eventuality of nuclear proliferation applications.

Modeling Nuclear Processes in the Reactor
In order to determine the feasibility of a continuous withdrawal of uranium isotopes from a thorium-based reactor, it is necessary to have a detailed time evolution of the densities of all nuclides. The model has to include both sources and sinks of the uranium isotopes, along with the necessary control procedures for the neutron flux, aimed to keep the reactor critical and affording the desired power output. After establishing the form of all the rate coefficients, the equations for the time evolution can be formulated in terms of all the relevant densities. In particular, the neutron density (or equivalently the flux) needs to be evaluated by detailing all its source, sink, and transport terms. The neutron diffusion term D∆n (D is the diffusion coefficient and ∆ the Laplacian operator), describing the changes in neutron density as a function of time and space, can be simplified by approximating the neutron depletion in the reactor volume as homogeneous, and due to the neutron flux across the surface. In this way, the time derivative of the neutron density, determined by the neutron flux at the reactor boundary, can in turn be evaluated by solving the spatial eigenvalue equation for neutron density along with the Neumann boundary condition. This procedure allows the problem to be formulated as a system of ordinary differential equations.
Considering a system of nuclei and neutrons (with energy ε, velocity v n , in the center of mass reference frame, and volume V) with reduced mass m and giving a reaction with energy-dependent cross section σ, the corresponding rate coefficient η for the process involving thermal neutrons would have the form η = σv n , that is where the second equality is obtained by definingε = ε/kT (k being the Boltzmann constant and T the temperature). The quantities ρ(ε) and Λ represent the density of translational states and the DeBroglie wavelength, respectively, (h is the Planck constant). We take into consideration the following reactions and decays, 233 Th 233 Pa 234 Pa plus the elastic scattering of thermal neutrons by 232 Th and 1 H. All rate coefficients for β − decay are indicated with the letter λ with either an ordering subscript or the subscript "dn" for the emission of delayed neutrons. The rate coefficients for neutron-nuclei interactions are indicated by η, with the subscript f for "fission", c for "capture" by heavy nuclei, and x for the neutron capture by either the fission products or the control neutron absorber. We use the reference density of metallic thorium ρ 0 = 11.7 g cm −3 ≡ 3.04 × 10 28 m −3 to define the reduced densities as a ij = ρ ij /ρ 0 , ρ ij being the number density. We follow the convention used in the Manhattan project, indicating the reduced density with the first lowercase letter of the element symbol, its subscript being the last digit of the atomic number (i) and the last digit of the mass number (j) [18]. For example, the reduced density of 232 90 Th, 233 91 Pa, and 235 92 U are noted as t 02 , p 13 , and u 25 , respectively. The same indexes label the various sources and sinks s. We use for the reduced density of neutrons the same symbol used for neutrons, n, and the fission fragments in Equations (4)-(16) are indicated as FF. Since the rate coefficients in processes (4)-(16) can be calculated from known data, we may write a set of differential equations for the time evolution of the species, once the initial conditions for the number density for all species are specified. In the non-dimensional time unit λ 13 t = 2.9741 × 10 −7 s −1 t, the evolution of the homogeneous system for reactions (4)-(16) is given by the equationṡ with the dot indicating the derivative with respect to λ 13 t, and s ij = source ij /λ 13 ρ 0 , and the parameters given in Table 1. The variable h indicates the reduced density of 1 H and the variable f indicates the reduced number density of the fission fragments, while the same letter as a subscript indicates fission. The parameters χ 1 = 1.27 and χ 2 = 0.00624 represent the net number of the prompt and delayed neutrons emitted after fission, respectively [19]. In Equation (25) we summarize the effect of the neutron absorber 113 Cd (obeying the equationċ = −ω x cn + s c , the source term s c being subject to a feedback control to keep the power output at the desired level) with the neutron poison fission products y, exhibiting a high cross section for neutron capture (as 135 Xe and 149 Sm, having total a fission yield of 0.074 and obeying the equationẏ = −ω x yn + 0.074 Ω f n). Representing all neutron absorbers by the single variable x and the absorption cross section σ x = 3.12 × 10 4 b, and summing these two equations we obtain Equation (25).
In the conventional time and density units, the diffusion coefficient D (length 2 time −1 ) for neutrons may be expressed through the diffusion length L d as with The subscript el refers to the elastic scattering of neutrons. Defining L = v/λ 13 , and expressing L d in terms of the reduced parameters we have The reduced diffusion coefficient D length 2 thus is Neglecting the contribution to the neutron flux of both the source s n and the delayed neutrons in Equation (24), we have:ṅ = κn + D∆n.
Assuming constant coefficients, a solution of the form n = ψ t ψ r [20], and separating the variables, we obtain the time equationψ with solution the spatial eigenvalue equation with the Neumann boundary condition [21] dψ and solution The quantity L t in Equation (40) is the transport length given by The constant C was chosen to match the initial conditions for a homogeneous neutron number density n 0 in a spherical reactor of radius R n 0 4 3 that gives The complete solution is and the Neumann boundary condition (40) becomes If we approximate the time derivative of the neutron density in the whole volume in terms of the instantaneous flux at the boundary ϕ R , as 4π R R 2 ϕ R /V = 3ϕ R /R, and express ϕ R according to Equation (45) as where ∂ r ≡ ∂/∂r, the the time derivative of the neutron reduced density becomes D ξ 2 n, and, in non-dimensional units.
Solving Equation (46) for ξ, we may evaluate the expression (48). Within this approximation, Equation (24) can be writtenṅ and the system of Equations (17)-(26) becomes a system of ordinary differential equations.

Simulation Details
The actual value of Dξ 2 is obtained iteratively by solving Equation (46) for ξ at the end of a simulation, and use this value back in Equation (49) until self consistency is reached. To validate the approximation given by Equation (48), we report in Table 2 the fairly stable values of ξ, sin ξR/R, and Dξ 2 at different times during a simulation. For s 23 = s 24 = s 25 = 0, we can look for a stationary solution of the system of the Equations (17)-(26). The stationary reduced densities of the system (17)-(26) (denoted with a bar) are as followsn Defining the ratio z between the power density at time t and the required value for the power density ρ w , (the average recoverable energy per fission being E f = 170 MeV [19]), the condition z = 1 at ω c02 t 02 = const givess We thus have stationary values for all state variables (except x and s x ) in terms of the value of t 02 and ρ w . To obtainx, we equate the expression for κ in Equations (49) and (27) and solve for x. Since κ is also a function of x through the term ω l , that in turns depends on x because of L t , Ω c , Ω f , and Ω el , this operation must be done iteratively. Oncex is known, Equation (25) affords the value ofs x . We actually do not find positive solutions for bothx ands x unless R is very large (R > 15 m), but nevertheless the above expressions turn out to be useful approximations to a true stationary state.
The values of all state variables are given in Table 3 for a reactor with a stationary source of 232 Th s 02 = 2.03 × 10 −4 (≡ 22.13 kg m −3 yr −1 ); the corresponding burnup rates and fractional power are given in Table 4. The processed uranium with the composition given in Table 3 corresponds to the molar fractions 0.5724, 0.3629, and 0.06463 for 233 U, 234 U, and 235 U, respectively. In reality, a proper stationary-state does not exist for realistic values of R. However, a pseudo-steady state can be attained: i.e., a state where all the values of the state variables change on a time scale much longer than the typical reactor operation time. For example, 5 years after startup, we haveṫ 02 /t 02 = −1.44 × 10 −4 yr −1 ,u 23 /u 23 = −3.50 × 10 −4 yr −1 ,u 24 /u 24 = 2.39 × 10 −2 yr −1 , andu 25 /u 25 = 2.96 × 10 −2 yr −1 .
We now address the issue of providing the fissile material to start a reactor. We begin by observing that, by including a sink of uranium of the form s ij = −s fs02ûij /Ŝ, where 0 ≤ s f ≤ 1 andŜ =û 23 +û 24 +û 25 , the stationary point of the system (17)-(26) (denoted byâ) satisfies the nonlinear system of equationsû 23 with the corresponding expression for the stationary 232 Th source beinĝ Since we do not have the option to actually operate the sink at the true stationary point, the actual expression used for the sink is and the composition of the extracted uranium fuel is given by integration over time of the sink term The solution of the system of Equations (60)-(62) is obtained by iteration until the composition given by Equation (65) converges to a stable value. We thus have the option of running the reactor while withdrawing, for example by continuous processing, the fissile material necessary to start a new reactor, with fuel of composition given by Equation (65).

Results and Discussion
We begin detailing the results of this computational study [22] with a simulation of the reactor with a startup phase of 2.0 d, where each individual source of uranium is set to the value s ij = 38.92ū ijũij /S, matching the composition given by Equation (65) and totaling s 23 + s 24 + s 25 = 2.02 × 10 −2 ≡ 6.13 kg m −3 d −1 . At the end of the startup phase, the amount of uranium put into the reactor is 12.26 kg m −3 , and the uranium source is switched off. The ratio of the reduced density of the uranium species to the corresponding pseudo-stationary density is plotted in Figure 1 during the startup phase.
The corresponding reduced density of 232 Th equals its initial valuet 02 = 0.05, equivalent to 581.12 kg m −3 . This value corresponds to the concentration of 2.53 mol/L, within the limits of thorium nitrate solubility in water given in [23] as S Th = 2.615 + 0.010 (T − 25) mol/L for 5 ≤ T ≤ 60 • C. The reduced density of hydrogen, assumed to be constant, is equal to 1.07. All values of the state variables are given as ratios to the pseudo-stationary state densities with s ij = 0 at the power density ρ w = 50 MW m −3 . In the two-day startup run with ρ w = 50 MW m −3 , the maximum value of s x (2.72 × 10 −4 ) corresponds to 3.99 × 10 −2 kg m −3 d −1 . The growth of the reduced density of all uranium isotopes is linear, an indication that the source term is always dominant with respect to the fission rates.
After the two-day startup the uranium source was switched off and the value of the source of 113 Cd was dynamically assigned based on the ratio z of the current power density to the set power density, i.e., s x = 0.01 (z − 1). The ratios of the reduced densities of the uranium species to the corresponding pseudo-steady-state density are reported in Figure 2 for a period of 5 years, and the value of the average neutron flux in this time period is 2.42 × 10 18 m −2 s −1 . Figure 1. Plots of the ratio of densities at time t (in days) with respect to the corresponding stationary density: t 02 /t 02 (for 232 Th, lavender), u 23 /ū 23 (for 233 U, green), u 24 /ū 24 (for 234 U, red), and u 25 /ū 25 (for 235 U, blue) versus time. The period is the two-day startup with the uranium source operative for a reactor with radius R = 1 m.

Figure 2.
Plots of the ratio of densities at time t (in days) with respect to the corresponding stationary density: t 02 /t 02 (for 232 Th, lavender), u 23 /ū 23 (for 233 U, green), u 24 /ū 24 (for 234 U, red), and u 25 /ū 25 (for 235 U, blue) versus time. The period is five years following the two-day startup of a reactor with radius R = 1 m in the absence of a uranium sink.
The reduced density of 233 U at the time of switch-off of the uranium source is ≈20% above its stationary value, while the other isotopes, 234 U and 235 U, reach the equilibrium value just at the end of the 5-yr time period.
Operating the uranium sink for a period of 5 years the system stays critical by increasing the neutron flux to 4.18 × 10 18 m −2 s −1 and lowering the uranium state variables to the values given in Table 7. Plots of the ratios of the reduced densities of the various species to the corresponding pseudo-steady-state density are reported in Figure 3. The withdrawn uranium has the composition u 23 = 5.89 × 10 −3 , u 24 = 1.19 × 10 −3 , and u 25 = 1.29 × 10 −4 , corresponding to the molar fractions χ = 0.82, 0.17, and 0.018 for 233 U, 234 U, and 235 U, respectively. At these molar ratios and the number density of metallic uranium ρ U = 4.83 × 10 28 m −3 , the fission rate per fast neutron would be ηρ U χ = 1.32 × 10 8 , 1.92 × 10 7 , and 1.01 × 10 6 s −1 for 233 U, 234 U, and 235 U, respectively. All these values are able to sustain a fission chain reaction in a metallic uranium sphere with the above isotopic composition. The corresponding ratios between the reduced densities at the end of the five years and the pseudo-steady state reduced densities of uranium are u 23 /û 23 = 0.86, u 24 /û 24 = 0.34, and u 25 /û 25 = 0.32. The withdrawal of the uranium isotopes reduces the reduced density of 232 Th to 95.9% of its initial value.   Based on a cross section of 1.80 mb for the (n, 2n) reaction on 233 U by 6 MeV neutrons [3], the formation of 232 U in a period of five years through the (n, 2n) channel gives u 22 / (u 23 + u 24 + u 25 ) = 2.70 × 10 −4 , corresponding to a specific activity of 2.18 × 10 11 Bq kg −1 in the processed uranium with the composition of the withdrawn fuel. Kang and von Hippel estimated the exposure at a distance of one meter from a 5 kg uranium sphere with 0.4 % of 232 U in secular equilibrium with 208 Tl to be 0.076 Sv h −1 kg −1 . The other uranium isotopes would contribute a specific activity of 2.85 × 10 11 , 3.74 × 10 10 , and 1.41 × 10 6 Bq kg −1 for 233 U, 234 U, and 235 U respectively. For comparison, the specific activity of 238 U is 1.24 × 10 7 Bq kg −1 .
The ratio (s 23 + s 24 + s 25 ) startup t/ (s 23 + s 24 + s 25 ) sink = 2.02 × 10 −2 t/1.53 × 10 −4 of the total uranium input at the end of the two-day startup period to the uranium sink rate of a reactor with s f = 0.5 at the end of the five-year period following startup is 0.72 yr, that is, a reactor running with a uranium sink for 0.72 yr would breed enough uranium to start a new reactor of the same size in 2 days.
We conclude by comparing the yield of this simulated reactor with the more common reactors based on the uranium fuel cycle. The ratio of the fissioned mass to the energy produced for the homogeneous reactor described above as 1.97 × 10 −13 kg J −1 , while the corresponding value for the fission of 235 U enriched at 3% is 4.20 × 10 −13 kg J −1 .

Conclusions
This work explores the possibility of energy production by a fission reactor in the absence of a continuous supply of fissile isotopes, since the dynamics of a homogeneous light-water thorium reactor can be exploited for breeding the uranium isotopes necessary for fission. In fact, after a short startup period of two days requiring fissile uranium isotopes (for a total inventory of 12.26 kg m −3 ), the modeled reactor can operate at a nearly steady state with a constant source of 232 Th, both with and without a uranium sink, i.e., a steady withdrawal of all uranium isotopes. After the fissile 233 U or 235 U has been used for the startup of the first reactor, its operation with a continuous uranium sink for 0.72 yr can afford the necessary fissile uranium isotopes to start a new reactor, allowing the continuous production of energy from only 232 Th. However, the limitation to nuclear weapon proliferation does not appear to be ensured by the molar fraction of 234 U in the extracted uranium fuel.