Simulation of Neptunium Extraction in an Advanced PUREX Process—Model Improvement

ABSTRACT Routing neptunium to a single product in spent nuclear fuel reprocessing is a significant challenge. In this work, we have further improved the simulation of neptunium extraction in an advanced PUREX flowsheet by applying a revised model of the Np(V)–Np(VI) redox reaction kinetics, a new nitric acid radiolysis model, and by evaluating various models for the nitrous acid distribution coefficient. The Np disproportionation reaction is shown to have a negligible effect. The models are validated against published ‘cold test’ experimental results; the ‘hot test’ simulation suggests that high neptunium radiolysis could help to achieve high recoveries using this flowsheet.


Introduction
Nuclear energy plays an important role in the supply of sustainable and secure electricity. However, reducing the impact on the environment is a critical challenge in the development of next generation nuclear reactors and the associated fuel cycles. The spent nuclear fuel, which is highly radioactive and toxic, must be treated carefully before disposal. Due to the complexity of the nuclear reactions, many fission products, as well as transuranic actinides, are present in spent nuclear fuel. Among them, neptunium has been highlighted to have a potential environmental impact, at least under certain conditions, due to its long half-life, high mobility, and radiotoxicity. [1] Therefore, there are potential benefits if neptunium (as well as other minor actinides) is removed from spent nuclear fuel before storage in a geological disposal facility. [2][3][4][5] The PUREX process has been successfully applied in spent nuclear fuel processing since the 1950s. [6] However, this process generally focused on the separate recovery of plutonium and uranium. In the solvent extraction process, neptunium is commonly distributed between aqueous and organic streams in the first cycle of PUREX flowsheets and thus requires specific stages to purify products from neptunium contamination. Routing neptunium from spent nuclear fuel to a single product will, therefore, simplify the process and this can lead to the overall reductions in the radioactive waste volumes and plant size. The preferred method is to adjust the operational parameters of the primary separation stage of the PUREX process to route neptunium with uranium and plutonium as part of an 'advanced PUREX' process designed for future closed fuel cycles. [7,8] However, the neptunium reduction-oxidation reactions are complicated due to the dual role of HNO 2 , which reduces Np (VI) to Np(V), but at low concentrations catalyzes the oxidation of Np(V) and the extractabilities of Np(IV), Np(V), and Np(VI) are significantly different. [8][9][10] Therefore, it is not easy to fully recover neptunium in this advanced PUREX process. [11,12] Thus, a mathematical model of the process, implemented as a simulation code, is needed to predict the neptunium extraction behavior and to guide experimental design and flowsheet operation.
In a previous work, [13] a process model for the flowsheet simulation of an advanced PUREX process [11] was developed and implemented in gPROMS which is a software environment allowing users to build, validate, and execute steady-state and dynamic process models. [14] In this simulation code, the choice of rate equations for neptunium redox reactions has an important effect on the outcome of the model. The main equilibrium reaction of neptunium in the nitric acid solution is expressed as: [15] NpO þ 2 þ 1:5H þ þ 0:5NO À 3 Ð NpO 2þ 2 þ 0:5HNO 2 þ 0: This reaction is one of the most complex reactions in spent nuclear fuel processing. Many studies of the reaction kinetics of this reaction are reported in literature studies. [15][16][17][18][19][20] However, these studies do not provide clear and consistent descriptions of the reaction rate. [21] Our previous simulation work adopted the kinetics of Koltunov [18] for this reaction. Based on the research of neptunium redox in aqueous-only phase and two-phase extractions, [11,22] a revised description of neptunium kinetics is applied in this work to improve the simulation model. During this reaction, the role of nitrous acid is important: it catalyzes oxidation at low concentrations, but acts as a reductant at high concentrations, reducing Np(VI) to Np(V). Hence, predicting the distribution of nitrous acid in the advanced PUREX process simulation is a key requirement. In this work, a more accurate model of nitrous acid distribution coefficients (also known as distribution ratios) replaces the simple method of Uchiyama [23] used in our previous work. [13] Due to the influence of nitrous acid on the redox reaction, nitrous acid generated by radiolysis of nitric acid [24] also needs to be considered in flowsheet simulation of spent nuclear fuel re-processing. Therefore, a new model for calculating radiolytic nitrous acid production was integrated into the simulation of neptunium extraction. Finally, the effects of the disproportionation reaction of Np(V) were evaluated and are also reported in this article.

System overview
The Advanced PUREX flowsheet for neptunium co-extraction with uranium and plutonium reported by Taylor et al. [11] and modeled by ourselves [13] is illustrated in Fig. 1. In this flowsheet, the neptunium in the aqueous feed F1 is assumed to be Np(V) because this is the most stable valence of neptunium in aqueous solution for up to 5M HNO 3 without nitrous acid present. [11,25] As the distribution ratio of Np(VI) in the extraction agent, tributyl phosphate (TBP), is much larger than that of Np(V), this flowsheet was designed to oxidize most of the feed Np(V) to Np(VI) by nitric acid in the presence of nitrous acid and then to extract Np(VI) by TBP. So, the aqueous feed A2 was designed to simulate radiolytically generated nitrous acid in the active feed to the cold test (using a surrogate solution without significant radiation power emission); this feed is not required in the subsequent 'hot test'(using real spent nuclear fuel) simulation. (Note that the nitrous acid in the aqueous feed F2 to the flowsheet is retained in the hot test simulation.) 14 Figure 1. Neptunium extraction flowsheet being simulated (primary extract-scrub section of an advanced PUREX process) [11] .

Neptunium redox reaction model
In the previous work [13] , we applied the reaction kinetic model of Koltunov [18] for the redox reaction equation (1). The current work replaces this kinetics model with expressions for the forward and reverse neptunium redox reaction kinetics as shown in Eqs. (2) and (3) by Edwards et al. [26] which is based on the work of Moulin: [16] v F ¼ À dC NpðVÞ;aq dt ¼ 1:24 Â 10 8 e À 8540 Tþ273 ð Þ C H þ ;aq C NO À 3 ;aq C HNO 2 ;aq C NpðVÞ;aq C HNO 2 ;aq þ C NpðVÞ;aq Tþ273 ð Þ C NpðVIÞ;aq ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi C HNO 2 ;aq p where T is the temperature in°C, C is the concentration in mol·L -1 , v F is the forward reaction rate, and v B is the reverse reaction rate. v F and v B are both in mol·L -1 ·s -1 . The presence of the organic phase (TBP-diluent) can change the rate of the neptunium redox reaction. [16,17,27] Considering this effect, the rate constants and the activation energies in Eqs. (2) and (3) have been regressed against our single-stage two-phase neptunium extraction experimental data. [11] The organic phase redox reaction kinetic model is the same as that used previously. [13] 2.3 Disproportionation of neptunium (V) The forward and reverse disproportionation reaction of Np(V) are expressed as Eqs. (4) and (5), respectively: [28] 2NpO Np(IV) can be oxidized by nitric acid according to the reaction [16] : The kinetics of these reactions have been investigated by several researchers. [16,[28][29][30][31][32] Koltunov et al. [29] found that, in the aqueous phase, the rate of disproportionation of Np(V) is dependent on the concentration of Np(V) and the acidity: In the organic phase, the reaction kinetic model for this reaction given by Sarsfield et al. is: [31] À dC NpðVÞ;or dt ¼ 0:74e À 93000 8:314 In the aqueous phase, Tachimori [33] (citing Rykov et al. [34] ) proposed a reaction kinetic equation for the reproportionation reaction (Eq. (5)): where C N; aq ¼ C H þ ;aq þ C NpðVÞ;aq þ 2 C UðVIÞ;aq þ C NpðVIÞ;aq À Á þ 4 C PuðIVÞ; aq þ C NpðIVÞ;aq À Á is the total concentration of nitrate in the aqueous phase.
In the organic phase, Wehrey et al. [32] showed that the reaction rate of the reproportionation reaction is given by: Together with Eqs. (2) and (3), Eqs. (7)(8)(9)(10)(11) fully define the kinetics of redox reactions in the aqueous and organic phases in our experiments [11] and are applied in simulation in this work. Kinetic expressions similar to Eqs. (7)(8)(9)(10)(11) were used by Tachimori [33,35] but with different Np(V)-Np(VI) redox reaction kinetics which was also derived from the work of Moulin. The main modification to the reaction models in this work is, for the first time, to use the Np(V)-Np(VI) redox reaction kinetics of Eqs. (2) and (3). Furthermore, the parameters of Eqs. (2) and (3) were regressed against single-stage two-phase experimental data.

Distribution coefficients
This work extends the previous work, [13] improves the model for determining nitrous acid distribution coefficients, and adds a distribution coefficient model for neptunium (IV). As before, the distribution coefficients of nitric acid and uranium (VI) are calculated using the SEPHIS model, [36] and distribution coefficients of neptunium (VI) are calculated using the model of Kolarik. [37] The distribution coefficient of neptunium (V) is set to 0.01, consistent with the approach of Tachimori in the EXTRA.M code. [35] The distribution coefficient of neptunium (IV) is calculated based on the relationship of Rozen for distribution coefficients [38] : The distribution coefficients of plutonium (IV) are calculated by the SEPHIS method. [36] The mean absolute percentage error (MAPE) of this method when applied to the 266 experimental data points of Kolarik and Dressler [39] is about 40%. This agreement is better than that with the model of Kolarik [37] (MAPE: 50%) or of Tachimori [35] (MAPE > 100%).
The βand γ-radiation powers are usually treated together, therefore: In spent nuclear fuel reprocessing, radiative power originates from radioactive isotopes in both aqueous and organic solutions. We can calculate the total emitted radiation power from the sum of the isotopes' specific output power P i (per unit mass) and their amount in each phase. This calculation is given in Eq. (23), where m i is the atomic weight of isotope i, N i is the molar amount of i. Values for P i for some isotopes are listed in Table 3. (Note that Table 3 does not distinguish between the specific output powers of β and γ) We assume the same absorbed power per unit volume in aqueous and organic phases. The absorbed radiation power of aqueous and organic solutions then can be calculated as in Eqs. (24) and (25), where V aq and V or are the holdup volumes in the centrifugal contactor. (The volume of connections between two contactors is treated as negligible in the model).

Single-stage neptunium extraction simulation
The actual neptunium extraction experiments carried out in a single-stage centrifugal contactor [11] are depicted in Fig. 3. The single-stage neptunium extraction simulations emulate these real experiments. In these single-stage experiments, the organic phase is passed through the contactor once while the aqueous phase is cycled via a reservoir.
(i) Neptunium redox reaction kinetics Figure 4 shows the results of simulations applying the new redox reaction kinetics described above for the single-stage contactor experiments. In these simulations, the radiolysis reaction is included in the simulation model; nitrous acid distribution coefficients are calculated by Eqs. (13)(14)(15)(16); all the mass transfer coefficients are set to 2×10 -5 m·s -1 . Other parameters and models are the same as in our previous work. [13] The results in Fig. 4 show reasonable agreement with the experimental results [11] . Figure 5 compares simulation predictions applying revised (Moulin) kinetics and Koltunov kinetics, together with the above revisions to the model, to experimental results. The mean absolute error (MAE) is taken to be the average value of the absolute error over all sample points of an experiment.

(ii) Neptunium disproportionation reactions
In nitric acid solutions, the neptunium (V) disproportionation reaction rate increases as the acidity increases. [25,28] The rate of this reaction is slow, due to the requirement to break Np-O bonds when reducing neptunium (V) (NpO 2 + ) to neptunium (IV) (Np 4+ ). Therefore, this reaction usually makes a negligible contribution to neptunium speciation in nitric acid. [21] However, the nitric acid concentration is around 5 mol·L -1 in the advanced PUREX process; therefore the simulation of the single-stage experiments [11] was carried out with and without including the disproportionation reaction to check the effects of this reaction. Figure 6 presents the deviation between simulation results and experimental results with and without considering disproportionation reactions. It may be seen that the disproportionation of neptunium (V) does not significantly change the simulation results. The concentration of neptunium (IV) is significantly lower than that of other neptunium valences in both phasesan example is shown in Fig. 7. This low concentration means that the disproportionation of neptunium (V) occurs, but its influence on neptunium extraction is negligible at the experimental conditions.
the simulated experiments were cold test experiments so the radiation power was negligible. Hence, the radiolysis of nitric acid in both phases is negligible in these simulations and integrating the radiolysis model does not change the cold test simulation results. The following section explores the potential impact of radiation on the process by considering radiolysis.

Flowsheet (multi-stage) simulation
Flowsheet simulation with various redox kinetics. The model of Koltunov for neptunium oxidation kinetics [18] was used in our previous work. [13] The oxidation reaction rate can be calculated as given in Eq. (26): [51] À  The value of the index z was 0.5 originally; it is found that, for nitrous acid concentrations around 10 −3 mol·L -1 , a value of z = 2 gives a better fit to experimental data. [51] This change in the value of z might be due to the accelerating effect of the organic phase on the neptunium redox reaction [16,17,27] or because the original z value was obtained at different experimental conditions to those studied here.
The results are presented in Fig. 9; these show that simulation with the revised kinetic model based on parameters obtained from single-stage experiments gives the best simulation of experimental results from the neptunium extraction flowsheet test. Simulation with Koltunov kinetics with z = 2 also gives good agreement with experimental data in HA banks (High Active banks, stages 5 to 14 in Fig. 1), but deviates from the experimental results in the HS bank (Hot Scrub bank, stages 1 to 4 in Fig. 1). Simulation with Koltunov kinetics with z = 0.5 are quite different to experimental results, especially in the HA banks. As can be seen in Fig. 9, Eqs. (2) and (3) provide a more realistic prediction than the kinetic model used by Tachimori. [33,35]   Effect of nitrous acid distribution coefficients. In our previous work, [13] the method of Uchiyama [23] was applied to predict nitrous acid distribution coefficients for the nitric acid/TBP-OK (odorless kerosene) system. Figure 10 presents two sets of simulation results for a multi-stage flowsheet, using nitrous acid distribution coefficient models based on the work of Uchiyama [23] and of Tachimori; [35] these results are compared with the experimental data. Figure 10 shows that the HNO 2 concentration profile simulated using the distribution coefficient model of Tachimori is in better agreement with the experimental profile than the simulated profile generated using the method of Uchiyama. Note that concentrations on stages 9-14 are below detectable limits of the relevant instruments. Figure 11 shows the simulated nitrous acid concentration in the organic phase. Figure 10 also indicates that, with the method of Tachimori, the HNO 2 concentrations predicted in the aqueous phase are closer to experimental data than the HNO 2 . Calculated concentrations of neptunium species present in aqueous and organic phases, based on simulation of singlestage experiment 12 in ref. [11] .  concentration predicted using the method of Uchiyama. However, Fig. 10 shows that in the HS bank, discrepancies between the simulation results and the experimental data are still significant. These discrepancies may also cause deviation between predicted and measured neptunium concentration profiles in the HS bank shown in logarithmic scale in Figs. 12 and 13. Figures 12 and 13 also suggest that although the two approaches predict different nitrous acid concentration profiles in HA banks (particularly in stages [6][7][8][9][10][11][12][13][14], predicted neptunium concentration profiles are relatively similar; differences in predictions are relatively small and only occur in the HA2 and HA3 banks (stages [9][10][11][12][13][14]. That different nitrous acid concentration profiles could correspond to rather similar Np profiles is unexpected based on an intuitive analysis of kinetic equations such as Eqs. (5) and (28); this result implies that the neptunium concentration profiles are affected by phenomena other than the reaction kinetics, such as extraction equilibrium and mass transfer effects.
Hot test simulation. Using the model presented above and validated against the results of the singlestage contactor experiments and cold test flowsheet, a flowsheet for reprocessing spent fuel was    designed by simulation (here referred to as the 'hot test' flowsheet). In the cold test flowsheet experiment, to replicate the radiolytic generation of nitrous acid in the extraction process, an extra flow of NaNO 2 (stream A2 in Fig. 1) was added into centrifugal contactor stage 7. In the hot test flowsheet, this A2 feed was removed; the nitrous acid content of the feed (F2) was maintained at the same level as in the cold test.
For the purpose of developing and demonstrating this radiation model within the overall flowsheet simulation, spent nuclear fuel of the composition given in Table 4 is assumed as the active feed. The specific output powers of the species based on a reference spent fuel (40 GWd per tonne, 5 year cooled) are listed in Table 3. Presently, during the simulation of the extraction process, the plutonium (IV) distribution coefficient is calculated by SEPHIS model, [36] while the fission products (FP) and all other isotopes including trivalent minor actinides (MA) are treated as nonextracting nitrates [52] . To simplify the calculations, it is also assumed that the average molecular weight of these non-extracting nitrates is 100 g·mol -1 and that the nitrate salt formed with any of them is FPNO 3 or MANO 3 . These assumptions would be refined in future work. Further work would also be needed to consider redox reactions of plutonium for hot test simulation.
Simulation results for the hot test are shown in Figs. 14 and 15. Figure 14 shows that the γradiation output power is stronger in the HA banks of contactors than in the HS bank. This effect is due to the strong γ-radiation from fission products which are modeled as being inextractable into TBP and are hence routed to the aqueous raffinate. The γ-radiation output power in the HS bank is consequently only due to the extracted (organic phase) plutonium. The α-radiation output power changes less significantly than that of γ-radiation. This may be explained based on the understanding  of the two main contributions to α-radiation: about half the contribution is from plutonium extracted from the aqueous to the organic phase by TBP and then enters the HS bank, and about half is related to inextractable minor actinides which remain in the aqueous phase and flow to the HA bank. Figure 15a shows the simulation results for nitrous acid in both cold and hot tests. It may be seen that nitrous acid concentrations are higher in the HS and HA1 banks (stages [1][2][3][4][5][6][7][8] in the cold test than in the hot test; this result implies that the amount of NaNO 2 added in A2 in the cold test is greater than the radiolysis yield calculated for the hot test. In the HA2 and HA3 banks, the modeling results show the nitrous acid concentration in the hot test is higher than in the cold test, and that the nitrous acid concentration reduces slowly from stage 9 to stage 14. This concentration profile is the result of radiolytic generation of nitrous acid in those stages. Although the concentrations of nitrous acid in the HA2 and HA3 banks in the hot test are higher than those in the cold test, the concentration of nitrous acid in both phases is still less than 0.1 mmol·L -1 . At these low nitrous acid concentrations, the reduction of Np(VI) to Np(V) by nitrous acid is negligible; instead, the nitrous acid mainly acts as a catalyst to accelerate the   oxidation reaction of Np(V) to Np(VI) by nitric acid. As a result, the low concentration of nitrous acid in HA2 and HA3 banks can improve the oxidation of Np(V) and the extraction of neptunium in these banks, thus reducing the leak of neptunium to the aqueous raffinate. Figure 15b shows that in stages 9-14, the aqueous phase neptunium concentration in the hot test is less than that in the cold test, while in the hot test the organic phase neptunium concentration is slightly higher than that in the cold test. The simulation of the hot test shows that the nitrous acid yield through radiolysis can improve the neptunium extraction in the advanced PUREX process and that 99% of neptunium would still be routed to the organic product when radiolysis occurs.

Conclusions
This article reports an extension to a previously published (open source) flowsheet simulation code designed to model neptunium extraction in an Advanced PUREX flowsheet. [11] The purpose of the modifications above are to improve simulation accuracy by applying a new neptunium redox reaction kinetic model, by modifying the nitrous acid distribution coefficient model and by integrating a model for predicting radiolytic nitrous acid generation into the flowsheet model. The effects of disproportionation of neptunium are also investigated and these were proved to be negligible.
Simulations of flowsheet performance with and without radiolytic nitrous acid generation (known as hot and cold tests, respectively) are presented. The results of the cold flowsheet simulation are in reasonable agreement with experimental test results [11] and the accuracy of simulation is shown to have improved. The hot test simulation results indicate the likely influences of radiolysis on neptunium extraction under radiation conditions that will be present when reprocessing spent nuclear fuel. The hot test simulation shows that, without extra NaNO 2 added to feeds, that the flowsheet used in the cold test should still reach about 99% recovery for neptunium (for the reference fuel and activity used). This result is supported in part by results from the CEA (France) who obtained an improved percentage recovery of neptunium from 90% in a cold test to >99% in a hot test when testing a similar flowsheet using pulsed columns as the contacting equipment in their Atalante facility. [53] CEA also used a specific nitrous acid feed in the cold test but omitted this feed and relied on radiolytic generation of nitrous acid in the hot test. Therefore, it appears that this flowsheet design can be used as the basis for a future hot test of the primary extract-scrub section of an advanced PUREX process that uses centrifugal contactors for separations.