The Pearson correlation coefficient

The Pearson correlation coefficient [30] measures the linear correlation between input X and output Y and is represented by "r." It ranges between -1 and 1, with -1 indicating a perfect negative correlation, 0 indicating no correlation, and 1 indicating a perfect positive correlation. The Pearson correlation coefficient is used to determine the strength and direction of the relationship between two variables and is widely employed in statistical analysis to assess the association between two numerical variables. However, it is sensitive to outliers and does not account for the possibility of one parameter being related to another. The formula for calculating the Pearson correlation coefficient is as follows: $r_{\text{Pearson}}=\frac{\mathrm{cov}\left(X,Y\right)}{{\sigma }_X{\sigma }_Y}=\frac{E\left[\left(X-{\mu }_X\right)\left(Y-{\mu }_Y\right)\right]}{{\sigma }_X{\sigma }_Y}=\frac{{\displaystyle \sum _{k=1}^N\left(x_{i,k}-{\overline{x}}_i\right)}\left(y_{j,k}-{\overline{y}}_j\right)}{\sqrt{{\displaystyle \sum _{k=1}^N{\left(x_{i,k}-{\overline{x}}_i\right)}^2}}\sqrt{{\displaystyle \sum _{k=1}^N{\left(y_{j,k}-{\overline{y}}_j\right)}^2}}}$ (2) where: $\begin{array}{l}{\overline{x}}_i={\displaystyle \sum _{k=1}^N{x}_{i,\text{k}}/N}\\ {\overline{y}}_j={\displaystyle \sum _{k=1}^N{y}_{j,k}/N}\end{array}$ (3) where N is the sample size; ${\mu }_X$, ${\mu }_Y$, ${\sigma }_X$ and ${\sigma }_Y$ represent the mean and standard deviation of variables X and Y, respectively; $x_{i,k}$ and $y_{i,k }$ represent the corresponding values of input and output parameters for the i-th parameter in the k-th sampling; $\overline{x}$ and $\overline{y}$ represent the sample means.

The Spearman rank correlation coefficient

The Spearman rank correlation coefficient [31] is based on the ranks of the inputs and outputs, rather than their actual values, to calculate the correlation between them. Therefore, even if there are large differences in the magnitudes between the inputs and outputs, the Spearman rank correlation coefficient is still effective, thus overcoming the limitations of the Pearson correlation coefficient. ${\rho }_{\text{Spearman}}=\frac{{\displaystyle \sum _{i=1}^{N}R}\left(X_i\right)R\left(Y_i\right)-N{\left(\frac{N+1}{2}\right)}^2}{{\left({\displaystyle \sum _{i=1}^{N}R}{\left(X_i\right)}^2-N{\left(\frac{N+1}{2}\right)}^2\right)}^{1/2}{\left({\displaystyle \sum _{i=1}^{N}R}{\left(Y_i\right)}^2-N{\left(\frac{N+1}{2}\right)}^2\right)}^{1/2}}$ (4) where R(Xi) and R(Yi) represent the rankings of the i-th variable among all samples, and N represents the total sample size.

The partial rank correlation coefficient (PRCC)

The Spearman rank correlation coefficient [31] calculates the correlation between two variables based on the ranks of their inputs and outputs, rather than their actual values. This approach allows it to remain effective even when there are significant differences in magnitude between the inputs and outputs, overcoming the limitations of the Pearson correlation coefficient. The Partial Rank correlation coefficient [9] measures the degree of the linear relationship between inputs and outputs while accounting for the linear effects of other parameters. In situations where we want to assess the correlation between two variables while considering the potential influence of other variables on their relationship, the partial correlation coefficient helps us determine the independent relationship between them by removing the influence of other variables.

Because the PRCC isolates relationships between variables from the influence of other variables, it provides a more accurate analysis and prediction of each variable's behavior and helps in better understanding the causal relationships between them [32]. In cases where multicollinearity exists among the parameters, the PRCC outperforms the Pearson correlation coefficient. ${\rho }_{{\text{PCC}}_{Y{X}_1{X}_2}}=\frac{r_{Y{X}_1}-r_{Y{X}_2}{r}_{X_2{X}_1}}{\sqrt{1-r_{X_2{X}_1}^2}\sqrt{1-r_{Y{X}_2}^2}}$ (5) where, ${\rho }_{{\text{PCC}}_{{YX}_1{X}_2}}$ represents the partial correlation coefficient betweern X₁ with Y after removing the influence of X₂ variable, and $r_{Y{X}_1}$, $r_{Y{X}_2}$ and $r_{X_2{X}_1}$ represent the linear correlation coefficient between two variables.

The standardized regression coefficient (SRC)

The Standardized regression coefficient [33] is a commonly used statistic in multiple linear regressions. It quantifies the average change in the dependent variable when the independent variable changes by one unit of standard deviation. Using X_j and Y to represent the inputs and outputs, respectively, we can establish a multiple linear regression model: $\widehat{y}=b_0+{\displaystyle \sum _{j=1}^J{b}_j{X}_j}.$ (6)

The absolute range of the standardized regression coefficients falls between 0.0 and 1.0, and their sign signifies whether there is a positive or negative correlation between X_j and Y. The formula for calculating SRC is as: ${\rho }_{{\text{SRC}}_j}=b_j\sqrt{\frac{\text{Var}\left(X_j\right)}{\text{Var}\left(y\right)},}$ (7) where b_j represents the regression coefficient of X_j on Y and Var(X_j) and Var(Y) represent the variances of X_j and Y, respectively.

Another crucial coefficient for assessing a regression model's capacity to fit the data is the model determination coefficient R². It assumes values between 0 and 1, where a value approaching 0 indicates a suboptimal fit to the data, while a value nearing 1 signifies a strong fit.

Hydrogen production

Figure 6a shows the transient change curve of hydrogen mass under various operating conditions, while Fig. 6b shows a scatter plot of total hydrogen production. As apparent from the figures, hydrogen was primarily generated during the initial stages of large-break loss of coolant accidents. This is because, immediately after the accident was triggered at 0.0 s, a substantial amount of coolant was expelled, leading to a rapid decrease in the core water level, coupled with an increase in fuel temperature. Simultaneously, the CMT and LHSI systems promptly initiated water injection into the high-temperature-exposed core, resulting in the generation of a significant amount of steam. Zirconium alloys undergo a vigorous oxidation reaction with water vapor when temperatures exceed 1500 K, leading to the early production of a substantial amount of hydrogen. Thanks to the timely activation of safety systems and core cooling, further hydrogen production through oxidation was averted. By the time 1000 s had elapsed, hydrogen production had reached a steady state.

Fig. 6Full-screen View

(Color online) Uncertainty analysis of hydrogen production. a Hydrogen production in pressure vessel; b hydrogen mass scatter point; c probability density distribution of hydrogen mass

As the accident developed, the uncertainty band for hydrogen production gradually widened. Through calculations, it is evident that the total mass of hydrogen generated falls within the range of 182.784 kg to 330.664 kg. Darnowski [13] calculated that hydrogen production from zirconium oxidation accounted for approximately 97% of the total hydrogen mass in a double-ended break LOCA accident in a PWR reactor. According to the design parameters of a million-kilowatt nuclear power plant, hydrogen production from 100% zirconium-water reaction is approximately 1010.1 kg. Based on the calculations in this article, it can be inferred that approximately 18.10% to 32.736% of zirconium reacts with water to produce hydrogen within the reactor core.

Examining the various percentile curves in Fig. 6a, it is evident that the calculated result of 278.348 kg for the reference operating conditions closely aligns with the median value of 253.242 kg. Figure 6c illustrates the probability density distributions (PDF) of hydrogen production under all operating conditions. After data processing, the overall sample calculation results fit a normal distribution with an R² value of 0.8968, and the one-sided upper limit of the hydrogen mass associated with a 95% cumulative probability density is 306.621 kg.

Thermal power is directly linked to the size of the reactor core and, to some extent, is proportionate to the mass of zirconium. Due to variations in reactor design parameters and calculation programs, a metric (the ratio of hydrogen production to thermal power within the reactor core) was introduced for comparative analysis, as outlined in Table 3. According to the table, the MELCOR simulation results for LBLOCA in different reactor types project hydrogen-to-thermal power ratios ranging from 0.04 to 0.11 kg H₂/MWth. Zhao [34] and colleagues [34] conducted a sensitivity analysis of the COR node in a Westinghouse three-loop nuclear power plant during an LBLOCA, predicting similar values within the range of 0.07–0.1 kg H₂/MWth. In this study, the best estimates and uncertainty analysis of the hydrogen-to-thermal power ratio obtained through ISAA were 0.0881 and 0.0578-0.105, respectively, which align with other simulation results.

Fission product release fraction

Figure 7a displays the curve representing the release fraction of fission products during the accident. Cs and I symbolize the release behaviors of highly volatile fission products, with their behaviors being relatively similar. The primary factor influencing their release from the fuel is the fuel temperature [22]. Shortly after the PCS pipeline experienced shear rupture, the rapid loss of coolant resulted in a swift decline in the core water level and a rapid increase in fuel temperature. During this phase, fission products were rapidly released at a high rate as the fuel temperature increased. Subsequently, the release rate of fission products decreased as the safety injection system consistently injected water to cool the core materials. Ultimately, further progression of the accident was arrested by the continuous cooling provided by the Low-Head Safety Injection (LHSI) system, resulting in the cessation of fission product release, which stabilized at a certain mass level.

Fig. 7Full-screen View

(Color online) Uncertainty analysis of FP release fraction. a Release fraction of Cs/I from fuel; b scatter plots for Cs/I fraction; c probability density distribution of Cs/I release fraction

Analyzing the data from Figs. 7a and 7b, it becomes evident that the range of the proportion of nuclide release to the initial reactor inventory is quite extensive, spanning from 15.6% to 84.3%. In comparison to the release fraction (less than 35%) obtained by Ahn et al. [9] in their uncertainty calculation for the Station Blackout (SBO) accident in the OPR1000 advanced reactor, the LBLOCA results in considerably more severe radioactive consequences. This paper primarily focuses on the fraction of fission products released after the core has been cooled. The instantaneous data distribution function at 15,000 s is shown in Fig. 7c. The fraction of fission products released from the fuel follows a normal distribution (R²=0.94), with the calculated value (43.956%) closely aligning with the median value (48.1%), indicating a high degree of consistency with the reference case.

The influence of Willks order

Non-parametric statistical methods, employing Wilks' formula, were used to calculate one-sided tolerance upper limits for hydrogen production and fission product release fractions at a 95% confidence level and 95% probability levels. The upper limits obtained through resampling are listed in Table 4.

The results reveal that non-parametric statistical methods inherently exhibit conservatism when calculating one-sided tolerance limits at lower orders. Calculations using the higher-order Wilks' formula yield more stable results. Studies have shown that the one-sided tolerance upper limits calculated with Wilks' formula gradually converge after the third order [24, 35]. Fig. 8 underscores this consistent trend. Mean values calculated using Wilks' formula with third-order or higher sampling have already converged, and increasing the order has relatively minor effects on the results.

Fig. 8Full-screen View

Wilks’ method analysis. a The 95/95 upper limit of hydrogen production; b the 95/95 upper limit of the FP release fraction

Sensitivity analysis of hydrogen production

Various sensitivity analysis methods were employed to assess the degree of influence of each uncertain input parameter on hydrogen production. The calculation results are presented in Table 5 and Fig. 9a. According to the results, the rupture temperature of the zirconia oxide layer (RT-ZrO₂) exhibits a strong positive correlation with hydrogen production, while the maximum flow rate of molten material after the oxide layer rupture (MaxFR-Mat) and the equivalent diameter of the particulate debris in the core (EDPD-CORE) show weak negative correlations with hydrogen production. Other parameters did not display significant correlations.

Fig. 9Full-screen View

(Color online) Sensitivity analysis of hydrogen production. a Importance of input parameters; b scatter plot of RT-ZrO₂ and hydrogen production; c scatter plot of EDPD-CORE and hydrogen production; d scatter plot of MaxFR-Mat and hydrogen production

Zirconium alloy forms an oxide shell layer on the exterior of the cladding due to the oxidation reaction, which prevents the molten material from undergoing relocation. It ruptures and relocates only when the oxide shell layer reaches its failure temperature. The failure temperature of the oxide shell layer determines when the molten material is released and impacts the reaction between water vapor and the zirconium alloy, subsequently affecting hydrogen production. A higher failure temperature results in a higher material temperature when the fuel rod collapses into particle debris, leading to a faster zirconium-steam oxidation rate and a greater amount of hydrogen production. To visually represent the linear relationship between this parameter and hydrogen production, Fig. 9b shows a scatter plot of temperature values and hydrogen production for each operating condition and conducts a linear fitting (R² = 0.593). The linear fitting outcome is consistent with the correlation coefficient calculations, indicating a significant positive correlation between the two. Wang et al. [15] conducted a sensitivity analysis of a Nordic boiling water reactor and obtained results that align with the conclusions of this study.

As shown in Fig. 9a, there is a negative correlation between EDPD-CORE and hydrogen production, consistent with the overall trend observed in the linear fitting shown in Fig. 9c. Smaller particulate debris likely have a larger contact area with water vapor, making oxidation reactions more likely to occur and produce more hydrogen, resulting in a negative correlation.

There is a weak negative correlation (ρ≈0.20) between the maximum flow rate of molten material after oxide layer rupture (MaxFR-Mat) and hydrogen production, as shown in Fig. 9d. This parameter partially determines the amount of molten zirconium exposed to oxidation outside the core. The analysis results in this study are consistent with the sensitivity analysis results of Choi et al. [10] for the OPR1000 reactor mitigation strategies. This could be attributed to the safety injection system rapidly injecting water into the reactor, with higher flow rates enhancing the heat exchange between the molten material and water, rapidly cooling it and resulting in less hydrogen generation.

Sensitivity analysis of FP release fraction

Table 6 and Fig. 10a show the sensitivity analysis results for the 17 model parameters concerning the total fission product release from the fuel. Similar to hydrogen production, a correlation was observed between RT-ZrO₂, MaxFR-Mat, and EDPD-CORE and the fission product release fraction. However, the porosity of the debris (POR-DB) exhibited a significant negative correlation (-0.853, -0.632) with the release of fission products, distinct from the sensitivity analysis of hydrogen production. Other model parameters did not exhibit clear relationships in this study.

Fig. 10Full-screen View

(Color online) Sensitivity analysis of FP release fraction. a Importance of input parameters; b scatter plot of RT-ZrO2 and release fraction; c scatter plot of MaxFR-Mat and release fraction; d scatter plot of POR-DB and release fraction; e scatter plot of EDPD-CORE and release fraction

As previously discussed, RT-ZrO₂ and MaxFR-Mat can influence the release of molten materials from the rod and the heat exchange process with water. A higher failure temperature of the oxide layer results in higher debris temperatures, delaying their cooling, and thus a more evident positive correlation between RT-ZrO₂ and the fission product release. After reaching the failure value of RT-ZrO₂, molten material is ejected from the oxide layer. At this point, due to the activation of the safety system and external water injection, a larger MaxFR-Mat enhances the heat exchange and cooling of the fragments, affecting the release of fission products.

A higher porosity of the fragments increases the heat transfer surface area within the core debris bed, enabling more effective cooling by the injected water. This, in turn, is closely related to the temperature, with lower fragment temperatures leading to reduced fission product release. Fig. 10d illustrates a strong negative correlation. The PRCC calculation results, corrected for the effects of other parameters, exceed 0.8, further indicating a strong correlation between porosity and fission product release. Additionally, there is a weak negative correlation between the equivalent diameter of the PD in the core (EDPD-CORE) and the release of fission products, as shown in Fig. 10e. In general, a smaller equivalent diameter of the PD allows for better contact with the surrounding water or steam, enabling more effective heat exchange and cooling, ultimately reducing fission product release.

Important parameters for FOMs

Through uncertainty and sensitivity analysis of LBLOCA with severe accident mitigation measures in a 1000 MWe APWR, several parameters significantly impacting the target output have been identified and organized in Table 7. Notably, RT-ZrO₂, MaxFR-Mat, EDPD-CORE, and POR-DB play more prominent roles among the selected 17 parameters. Importantly, these parameters are all user-defined empirical values. To minimize uncertainties in numerical simulations, further improvements are necessary in future research.