Study on the Calculation Method of Molten Pool Decay Heat Distribution under IVR Condition

As an important strategy to mitigate severe accidents, the in-vessel retention (IVR) technique has been applied to the new generation of pressurized water reactor (PWR). However, under IVR conditions, the decay heat distribution in the molten pool is very uncertain because of the complexity of the molten pool and the calculation method limitations. To explore the calculation method and distribution of the decay heat of lower head molten pool under IVR conditions, the decay heat calculation method is developed based on Reactor Monte Carlo Code (RMC). The verification results show that the relative error of calculation result is generally within ± 0.25%. In addition, geometric modeling for lower head molten pools has been carried out, and distribution of the decay heat in two-layer and three-layer structures has also been accurately calculated. The calculation results indicate that the decay heat power spatial distribution is relatively uniform in the two-layer molten pool structure. The decay heat power at the center of the lower head decreases from 0.71° W/cm 3 to 0.023° W/cm 3 within 1d-5d. In the three-layer molten pool structure, the spatial distribution of the decay heat power is severely uneven due to the precipitation of heavy metal uranium. Besides, in actual engineering calculations, it should lay emphasis on the heat transfer characteristics and design margin of the upper part of the heavy metal layer and the lower part of the oxide layer because the maximum decay heat power appears at these two positions.


INTRODUCTION
In-Vessel Retention (IVR) technology is an important strategy for serious accident mitigation. It arrests relocated molten core materials in the vessel during the severe accident and has been singled out as an appealing accident management approach to many reactors (Theofanous, 1989;Su, 2016). The specific process of IVR technology is as follows. During the severe accident in PWR, the core melts if it cannot be effectively cooled. The melting core relocates to the lower head of the pressure vessel and then forms a liquid molten pool with a variety of pool configurations. Since the IVR strategy confines all the fission products in the vessel, complex ex-vessel severe accident phenomena can be avoided.
At present, a commonly used method to achieve IVR is by carrying out External-Reactor-Vessel-Cooling (ERVC) strategy to remove the in-vessel decay heat and ensure the vessel integrity by fully depressurizing Reactor Coolant System (RCS) and flooding the cavity to the designed elevation. When ERVC occurs, the singlephase convective heat transfer and boiling heat transfer outside the pressure vessel will carry away most of the decay heat of the molten pool and ensure the integrity of the pressure vessel. The schematic diagram of the IVR + ERVC combination strategy is shown in Figure 1.
In Figure 1, the heat transfer imposed by the in-vessel molten pool is a vital part of IVR success, considering the difficulty of significantly altering ex-vessel CHF. However, the boiling crisis will occur when the current heat flux is higher than the CHF, which will lead to the rapid rise of the temperature on the lower head and even cause pressure vessel failure. Therefore, the IVR safety strategy requires that the decay heat of the molten pool must be taken away via a cooled vessel wall. To learn a modest safety margin of IVR strategy, the physical phenomena like heat transfer characteristics relevant to IVR must be studied systematically and methodically to ensure the heat flux imposed by the molten pool is below the CHF of a watercooled vessel wall, which is known as thermal criteria.
From Table 1, the previous study of the IVR has several basic deficiencies. First, the experiments often used water or liquid mixture as a substitute for the molten pool and used electric heated wires or rods to carry out uniform heating. The simulations cannot reflect the actual physical process in the lower head. Secondly, there is no calculation or estimation of the decay heat generation and its temporal and spatial distribution in the lower head. Therefore, the above experiments have not comprehensively proposed an analysis methodology that is sufficient enough to quantify molten pool characteristics and support subsequent heat transfer analysis.
In terms of calculation methods, due to the complexity of the physical and chemical characteristics of the molten pool, it is almost impossible to calculate the decay heat through the theoretical analysis method, making the estimation of the inner heat load of the lower head greatly uncertain. Besides, the existing methods for calculating the decay heat are also limited, and the existing programs for calculating the decay heat are often unable to deal with the problem of calculating the large-scale burnup in engineering. For that reason, there are few studies on the heat distribution of decay heat generation in the lower head after core-melting. Therefore, the evaluation of the IVR strategy lacks unified results.
At present, the commonly used decay heat calculation standard is ANSI / ANS-5.1 (ANSI/ANS, 2005) issued by the United States, and the relevant domestic standards (e.g., NB / T 200056-2001(Nuclear Industry Standardization Institute China Institute of Atomic Energy, 2011 ) are also based on their calculation methods. In the latest ANSI / ANS-5.1 2005 edition, the decay heat calculation formula is as follows.
In the formula, The reactor shuts down after running for T seconds, and the total decay heat generated hereafter in t seconds is P d (t, T), P dFP (t, T) is the decay power of fission products, P dAP (t, T) is the decay power after activation of fission products, and P dHE (t, T) is the decay power of transuranic elements produced by capturing neutron. Although the decay heat calculation standard provides an available calculation method for calculating the decay heat in the molten pool, the spatial distribution and decay power of each element are usually unknown in the calculation of the decay heat in the IVR. And in this calculation standard, a large number of approximate and simplified methods (e.g., the transuranic element only considers the U239 and Np239) are used, which cannot realize the accurate calculation of the decay heat distribution in the molten pool.
To solve the above problems, based on the Depth program coupled with Monte Carlo code RMC, this paper has developed the decay heat calculation function. The program has realized the function of large-scale grid element decay heat calculation and can output a variety of decay heat information. It can be used for lots of scenarios like reactor shutdown, molten pool calculation. In this paper, the BEAVRS benchmark (Horelik et al., 2017) is taken to validate the decay heat generation of molten pool in the lower head with time and space by using the above calculation program.
The reactor shutdowns after running for T seconds at the power P , and the decay heat power generated hereafter in t seconds P d (t, T) can be calculated using the following empirical formula: However, this formula can only be used to roughly estimate the decay heat power at certain times after the reactor shutdowns, and it is completely impossible to accurately calculate the decay heat distribution of the molten pool in the lower head This work developes the new feasible method to analyze the decay heat generation of the molten pool in the circumstance of severe accident with Monte Carlo code RMC based on the BEVARS benchmark. RMC is used for molten pool configurations modeling and criticality calculation and the decay heat generation is also obtained using the DEPTH program. Furthermore, a comparison has been made between this work and previous results to validates the method.

DEPTH PROGRAM DECAY HEAT CALCULATION MODEL
At present, the commonly used decay heat calculation program is represented by the ORIGEN-S program developed by ORNL Laboratory (Hermann and Westfall, 2011) . The ORIGEN-S program couples an external third-party interface to the Monte Carlo program (such as MCNP) and completes data transmission by reading the output file while generating input files. However, as this coupling method cannot go deep into the bottom of the Monte Carlo transport module for method improvement and efficiency optimization, there is a general lack of computing power for large-scale burnup problems. When there are more burnup zones, it is impossible to accurately calculate the burnup information of the whole molten pool. The DEPTH program (She, 2013) was developed by the REAL team of Tsinghua University. Based on the Monte Carlo program RMC (She et al., 2010;Wang et al., 2011), the internally coupled burnup calculation and decay heat calculation functions were realized. The large-scale burnup calculation ability (such as the automatic expansion of the repetitive structure burnup zone, the large-scale parallel efficiency optimization of the Monte Carlo, and the parallel calculation of the ignition burnup) could be used to calculate the whole PWR burnup and other examples.
Because of the common decay heat calculation requirements in reactor shutdown and molten pool scenarios, this paper proposes the corresponding decay heat calculation method based on the DEPTH program coupled with the Monte Carlo code RMC, which enriches the functional modules of the RMC program. Besides, as the RMC program being the completely independent intellectual property rights of the domestic Monte Carlo program, the development of the decay heat calculation module is of great significance to break the monopoly of foreign computing software.
The flow chart of decay heat calculation in the DEPTH program is as Figure 2.  (Theofanous et al., 1997(Theofanous et al., ) 1997 The melt is heated directly to a higher temperature and poured into a pressure vessel RASPLAV (Theofanous et al., 1997;Asmolov, 2000;Asmolov, 1999) Korchatov Nuclear Energy Research Institute of Russia

2015-2019
water/ 350°C/ the mole fraction of NaNO 3 is 20% and KNO 3 is 80% Electric heating rod heating PWR LIVE-L4 (Miassoedov et al., 2007;Fluhrer et al., 2008;  In the DEPTH program, first of all, the transport calculation is needed to obtain the neutron flux in the lower head, and the power distribution is carried out according to the total power to obtain the power of each burnup zone. Then, in each time step, the burnup is calculated for the decay reaction to obtain the nuclide density in the burnup zones. Finally, the decay heat in the burnup zones is calculated with the formula below: where Q is the total amount of decay heat generated in the burnup zone, R is the reaction rate, and q is the decay heat released of each nuclide. Compared with the existing methods, the method proposed in this paper has the following innovative points: 1. Compared with the existing ORIGEN-S program, it realized the parallel decay heat calculation of the large-scale burn-up area, and achieved the precise distribution result of the decay heat release inside of the lower head by calculation; 2. The decay heat release was calculated by Monte Carlo method, which is more accurate than the empirical formula methods that are commonly used.

CALCULATION PROCESS OF MOLTEN POOL DECAY HEAT
When a severe accident happens, the core is completely exposed and melted, and then the melt penetrates the core bottom plate. Finally, various elements are redistributed in the lower head of the pressure vessel, and a molten pool is formed. Figure 3 shows the calculation flow chart of the decay heat of the molten pool in this paper. The specific process is shown as follows: (1) Calculation of total density of core nuclides. The total nuclide density of all fission products in the core is obtained through the calculation of the whole core burnup. In this paper, the total burnup of commercial nuclear power reactors described in the BEAVRS benchmark is calculated by using the Monte Carlo program RMC and the TianHe II supercomputing platform. The total nuclide density of fission products in the core is obtained. The calculation time of the reactor is 85°days, and the thermal power of the reactor core is 3411MW, which is divided into 520,000 burnup zones. In addition, calculating the total nuclide density of control rods and other core structural parts refer to the BEAVRS benchmark manual.
(2) Extract important nuclide density. The DEPTH program can output the density information of 1487 nuclides in the burnup chain. However, previous studies have shown that about 50 nuclides can envelope more than 95% of the decay heat contribution after shutdown (Li et al., 2014) . Therefore, in this paper, only the 50 nuclides that contribute greatly to the decay heat, as well as the main heavy isotopes and fission products in the main burnup chain are considered. And the zirconium-water reaction was not considered in this paper.
(3) Calculation of the specific volumes of three layers is based on several assumptions. In this paper, it is assumed that the oxidation degree of the zirconium claddings is 0.5, the molten pool is well-mixed the moment it falls into the lower head, and the densities of each component in the molten pool are consistent. Based on the phase diagram and researches on previous studies, the mass of the main components in each layer is obtained, and herein the volumes and thickness are also determined (Liu et al., 2019) (Liu et al., 2019) . (4) Calculation of nuclide density in each layer of the molten pool. In this paper, it is assumed that the atomic number of each nuclide in the lower head at the initial time is equal to the value of the core burnup calculation of the end time.
Because the decay heat generated will decrease with time, so this assumption can simulate the maximum decay heat power released by the molten pool of the lower head, Then, dividing the molten pool into an approximate number of burnup zones along the radius direction, height direction and circumference direction, in order to accurately calculate the decay heat generated in the molten pool. (6) Calculation of decay heat of the molten pool. The DEPTH program was used to calculate the decay heat power distribution in the molten pool.

Hierarchical Structure of Molten Pool Model
In the previous IVR safety strategy analysis, it is generally believed that the structure of the molten pool in the lower head is a two-layer model (Theofanous et al., 1996) . That is, the light metal layer is located above the oxide layer. However, the MASCA project carried out by Russian research center of Kurchatov institute has revealed that the physical-chemical interactions in the molten may cause participation of metallic uranium and zirconium initially located in the oxide phase.
The change of density of stratification will lead to sinking of heavy metal layer. Previous experiments have found that, a great amount of uranium element is extracted by unoxidized zirconium form oxidic phase to metallic phase. With the increment in density of metallic phase, it sinks and be located below the oxide phase, thus forming a three-layer corium pool (Kim and Olander, 1988;Hayward and George, 1996). So the above-mentioned two-layer molten pool model cannot completely envelope the possibility of occurrence (Parker et al., 1990) and INEEL proposed a more enveloping three-layer melt configuration based on the twolayer one (Gu et al., 2017) (Tsurikov et al., 2007) . In this configuration, the heat flux distribution is different from that of the two-layer molten pool structure, resulting in changes in the safety margin for the IVR scheme. So, the two configurations need to be analyzed separately. Figure 4 shows two kinds of molten pool structures.  As shown in Figure 4A, in the three-layer molten pool model, the components of the oxide layer and the light metal layer are the same as those in the two-layer molten pool model. The difference is that the elements in the partial oxide layer and the metal layer will enter the bottom of the oxide layer after the reaction, forming a heavy metal mixture layer.
As shown in Figure 4B, in the two-layer molten pool model, the pool can be divided into the bottom oxide layer and the top light metal layer. The light metal layer consists of stainless steel and unoxidized zirconium, and the oxide layer consists of uranium and zirconium oxides.
For the two-layer and three-layer structures of the molten pool, this paper uses the geometric module of the RMC program to carry out the three-dimensional modeling. It divides the burnup zones of the light metal layer, oxide layer, and heavy metal layer in the two-layer and three-layer structures of the molten pool to accurately calculate the decay heat power distribution inside the molten pool of the lower head.

Modeling of Two / Three-Layer Molten Structures
According to the above hierarchical model of the molten pool, the components and heights of each layer need to be obtained first. In the molten pool, the main components are the core fuel and the melts of structural components in the core, which form the U-Zr-O multivariate system after the core melts. Assuming that the density of each component is constant in the melting process and the molten pool is instantaneously uniform, the components of each layer can be obtained by multivariate phase diagram analysis (Chevalier et al., 2004), and then the volume of each layer can be calculated by the density of each component.
The heights of each layer can be obtained by a multivariate phase diagram as shown in Table 1. And the material components and typical nuclides of each layer are shown in Table 2. This stratification is based on the following assumptions: 1. The main components of the light metal layer are Zr and stainless-steel SS; 2. The main components of the oxide layer are UO 2 and ZrO 2 ; 3. The main components of the heavy metal layer are U and Zr; 4. The Ag-In-Cd elements and fission products in the control rod are uniformly distributed in the oxide layer and the heavy metal layer.
In this paper, the three-dimensional modeling of the twolayer and three-layer molten pool structures (Figure 4) is carried out by using the RMC program, and the light metal layer, heavy metal layer, and oxide layer in the molten pool structure are divided into burnup zones, to accurately calculate the decay heat power distribution inside the molten pool. The components of each layer and typical nuclides are shown in Table 3, 4:

Three-Layer Molten Pool Model Burnup Zones Division
To accurately calculate the decay heat generated in the molten pool, it is necessary to finely divide the molten pool into an approximate number of burnup zones. In this paper, the molten pool of the lower head is first stratified with each 10°cm along the Z-axis height direction. So, for the three-layer pool, the light metal layer is divided into five layers; The oxide layer is divided into 12 layers; The heavy metal layer is divided into four layers. Each layer is divided into 11 layers along the radius direction and 18 layers along the circumference direction. Finally, the light metal layer is divided into 990 burnup zones; the oxide layer is divided into 2016 burnup zones; The heavy metal layer is divided into 288 burnup zones, and the complete partition and filling of the light metal layer, oxide layer, and heavy metal layer in the three-layer molten pool model are then realized. For the two-layer molten pool, the light metal layer is divided into six layers; The oxide layer is divided into 16 layers. Each layer is divided into 11 layers along the radius direction and 18 layers along the circumference direction. Finally, the light metal layer is divided into 1288 burnup zones; The oxide layer is divided into 2304 burnup zones. The complete partition and filling of the oxide layer and heavy metal layer in the two-layer molten pool model are then realized.
By calculating the decay heat generated in each burnup zones in the molten pool, the complete spatial distribution of the decay heat release in the lower head can be obtained. Figure 5 is the side views and overlooking view of the molten pool model.

CALCULATION RESULTS AND ANALYSIS OF DECAY HEAT Verification of the Decay Heat Calculation Function
To validate the decay heat calculation function, the researchers used DEPTH and ORIGEN-2 to calculate the nuclear density and decay heat of all the main products of 237 Np after 106°years' decay (She, 2013) . To facilitate the comparison of calculation results, both programs applied the decay database of ORIGEN-2. Table 5 shows the comparison of the decay calculation results of 237 Np by DEPTH and ORIGEN-2. The relative error of the calculated decay heat of DEPTH and ORIGEN-2 fluctuates within ± 0.25%. The maximum relative error in the entire decay chain is -1.14%, which was produced by the last nuclide 209 Bi in the decay chain.

Decay Heat Distribution Calculations of the Melt in Lower Heads
The decay heat power of the melt within 5°days is calculated. The decay process is divided into 30 burnup steps. About 1,320 core hours are used for each calculation. The fission source iteration method is adopted for the critical calculation, specifically, 200,000 for neutrons in each generation, 20 for inactive generations, 120 for the total generations, and the initial fission source is the point source at the center of the lower head. The total power in the molten pools is calculated by the decay heat power empirical formula (2) after the reactor shutdowns. As the power change shown in Figure 6, the accuracy of the total power assumption is ensured.
For the lower head with a three-layer molten pool structure, the light metal layer is divided into 990 burnup zones; the oxide layer is divided into 2,016 burnup zones; the heavy metal layer is divided into 288 burnup zones; for the lower head with twolayer molten pool structure, the light metal layer is divided into 1,288 burnup zones; the oxide layer is divided into 2,304 burnup zones.
Previous researchers usually choose analysis programs such as IVRSA and MAAP to analyze and calculate the IVR characteristics of reactors under severe conditions to obtain the corresponding parameters. For example, Zavisca et al. used MELCOR, MAAP, etc. to analyze and calculate the serious accidents caused successfully by a pressure vessel's external submergence Zavisca et al., 2003); Esmaili et al. used the same method to calculate the initial IVR parameters of the AP1000 nuclear power plant in the severe accident (Esmaili and Khatib-Rahbar, 2004) . The decay heat parameters calculated by the predecessors and obtained herein were listed below. It can be observed that the calculation results of the decay heat are roughly the same when the redistribution time is similar, which verifies the validity of decay heat calculation results of the molten pools in this research.

Analysis of Decay Heat Distribution Calculation Results of the Melt in Lower Heads
For time distribution of the melt in lower heads, Figure 7 exhibits the change of the decay heat power generated by each layer in the two-layer molten pool structure and the three-layer molten pool structure after 30 burnup steps within 0d-5d.
It can be seen from Figure 7: (1) In the two-layer molten pool structure, the decay heat generated in the light metal layer, the oxide layer and the heavy metal layer all decrease with time. It decreases the most on the first day, and then the trend becomes more and more gentle, which meets the decline trend of the decay heat power after shutdown. (2) The decay heat power of the oxide layer is about 10-15 times that of the heavy metal layer, while the decay heat power of the light metal layer is only about onethousandth of that of the heavy metal layer and the oxide layer. Therefore, in the evaluation of the IVR safety strategy, more attention should be paid to the heat flux density at the position of the oxide layer and the heavy metal layer.
Regarding the spatial distribution of the melt decay heat in lower heads, due to the symmetry of lower heads, it can visually observe the distribution of the decay heat power density of the entire lower heads with the cross-sectional side view. The decay heat power density shown in the Figure 8 is obtained by the following formula: there, E cell is the power density in the burnup zone, that is, the decay heat power released per unit volume in the burnup zone; W cell is the decay heat power generated in the burnup zone; and V cell is the burnup zone's volume. Figure shows a half cross-sectional side view of the decay heat power density of the melt in lower heads of the oxide layer and heavy metal layer at θ 0°, and it is colored according to the magnitude of the decay heat power in two single molten pools. The left is the decay heat power generated by the three-layer molten pool structure in 1d, 2d, 3d, 4d and 5d; and the right is the decay heat power generated by the two-layer molten pool structure in 1d, 2d, 3d, 4d and 5d. It can be seen from Figure 8 that the spatial distribution of the decay heat power of the oxide layer is relatively uniform and at the same order of magnitude in the two-layer molten pool structure. It is because oxides such as uranium dioxide and zirconium dioxide are uniformly distributed in the oxide layer in the two-layer structure, forming a stable configuration. Notably, because the initial fission source is the point source at the center of the lower head, the maximum value of the decay heat power density in the lower head appears at the center of the lower head, and the value gradually decreases as it approaches the wall of the lower head. In terms of time distribution, the decay heat power at the center drops rapidly from 0.71 to 0.033 within 1d; it slowly decreases to 0.023 and gradually reaches a steady state within 2-5d; while the decay heat power at the inner wall surface is relatively stable with little fluctuation, from 0.051 to 0.022 within 1d-5d.
In the three-layer molten pool structure, the spatial distribution of the decay heat power of the oxide layer and the heavy metal layer is very uneven, and the maximum value of the decay heat power is increased by about 4 orders of magnitude compared with the minimum value over the same time frame. This is because a complex chemical reaction occurred between molten uranium dioxide and molten zirconium under severe accident conditions, resulting in the precipitation of heavy metal uranium that caused a very uneven distribution of elements inside the lower head. In addition, since the initial point source is located in the middle of the molten pool. Therefore, in the three-layer molten pool structure, the decay heat is mainly concentrated at the bottom of the oxide layer and the upper part of the heavy metal layer. At the same position of the heavy metal layer, the decay heat power of the threelayer structure is 2-3 times that of the two-layer structure; on the upper part of the oxide layer, the decay heat power of the threelayer structure is 2-3 orders of magnitude lower than that of the two-layer structure. It is noted that the maximum decay heat power density in the lower head appears in the middle and lower parts of the oxide layer. In terms of time distribution, the maximum decay heat power density decreases from 2.9°W/cm 3 to 1.5°W/cm 3 within 1d-5d; while the decay heat power at the inner wall surface is relatively stable with unobvious declination, from 0.094°W/cm 3 to 0.085°W/cm 3 within 1d-5d. It is noteworthy that the maximum decay heat power of the inner wall surface appears at the top of the heavy metal layer and the lower part of the oxide layer or corresponding positions in both two-layer and threelayer structures. Therefore, the design margin of these two locations should be laid emphasis on in actual engineering calculations.

CONCLUSION
All along, the calculation and estimation of the decay heat in molten pools under IVR conditions are controversial due to the limitations of calculation tools and the complexity of physical processes in molten pools. To solve the above problems, the decay heat calculation module is developed based on DEPTH coupled with RMC, and its functional verification is carried out. Then, the program developed in this paper is taken to calculate the distribution of the decay heat in the two molten pool structures, and the following conclusions can be drawn: 1. The decay heat calculation results meet the physically decline change trend after shutdown. By comparing with the results calculated by previous programs such as MELCOR and MAAP, the validity of calculation results of is verified. 2. In the two-layer molten pool structure, due to the uniform distribution of nuclides, the spatial distribution of the decay heat power of the oxide layer is relatively uniform. At the center of the lower head, the decay heat power slowly falls from 0.71°W/cm 3 to 0.023°W/cm 3 . Besides, in the three-layer molten pool structure, the spatial distribution of the decay heat power is severely uneven due to the precipitation of heavy metal uranium. 3. The maximum decay heat power appears at the top of the heavy metal layer and the lower part of the oxide layer or corresponding positions in two structures. Therefore, it should lay emphasis on the heat transfer characteristics and design margin of these two locations in actual engineering calculations.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.