Using Machine Learning to Predict the Fuel Peak Cladding Temperature for a Large Break Loss of Coolant Accident

In this paper the use of machine learning (ML) is explored as an efficient tool for uncertainty quantification. A machine learning algorithm is developed to predict the peak cladding temperature (PCT) under the conditions of a large break loss of coolant accident given the various underlying uncertainties. The best estimate approach is used to simulate the thermal-hydraulic system of APR1400 large break loss of coolant accident (LBLOCA) scenario using the multidimensional reactor safety analysis code (MARS-KS) lumped parameter system code developed by Korea Atomic Energy Research Institute (KAERI). To generate the database necessary to train the ML model, a set of uncertainty parameters derived from the phenomena identification and ranking table (PIRT) is propagated through the thermal hydraulic model using the Dakota-MARS uncertainty quantification framework. The developed ML model uses the database created by the uncertainty quantification framework along with Keras library and Talos optimization to construct the artificial neural network (ANN). After learning and validation, the ML model can predict the peak cladding temperature (PCT) reasonably well with a mean squared error (MSE) of ∼0.002 and R2 of ∼0.9 with 9 to 11 key uncertain parameters. As a bounding accident scenario analysis of the LBLOCA case paves the way to using machine learning as a decision making tool for design extension conditions as well as severe accidents.


INTRODUCTION
Deterministic safety analysis has traditionally been utilized to demonstrate the robustness of nuclear power plants, usually adopting a conservative approach. However, the conservative approach relies on a number of assumptions that do not necessarily reflect the real plant performance (Queral et al., 2015). On the other hand, the best estimate (BE) approach provides a more realistic system response based on detailed thermal-hydraulic mechanistic models provided it is accompanied with uncertainty quantification (UQ). The integration of BE and UQ is known as best estimate plus uncertainty (BEPU) and is built upon a statistical foundation to provide a more realistic estimation of the safety margin and hence ensure that the safety limit is met.
Utilities were given an ample opportunity to apply the best estimate plus uncertainty (BEPU) methodology following the United States Regulatory Commission (USNRC) amendment of 10CFR50.46 Appendix-K in 1988. Accordingly, the USNRC assisted in the steady transition from the conservative to BEPU methodology by introducing the USNRC Regulatory Guide 1.157, "Best Estimate Calculations of Emergency Core Cooling System Performance" and the demonstration of code scaling, applicability and uncertainty (CSAU) methodology in 1989 aiming to quantify the uncertainty parameters (USNRC, 1989). In addition to the CSAU methodology, various international collaboration projects had been undertaken to propose and validate other uncertainty quantification methodologies such as the uncertainty method study, UMS, (OECD, 1998), the BEMUSE project (OECD, 2007a) and the SM2A study within the SMAP framework (OECD, 2007b).
The BEPU methodology has been used to predict key safety parameters such as the peak cladding temperature (PCT), departure from nucleate boiling ratio (DNBR), etc. for critical accident scenarios. In BEPU analysis, a BE code is used to simulate the plant response given the variations in a multitude of uncertain parameters (UPs) that can be propagated within the thermal-hydraulic system code. The process of uncertainty propagation is however lengthy and hence BEPU analysis has so far been limited to the analysis of bounding design basis accident (DBA) scenarios, e.g. large break loss of coolant accident (LBLOCA) (Chang et al., 2020) and only recently to the analysis of a station blackout (SBO) (Musoiu et al., 2019) and to the main steam line break (MSLB) (Petruzzi et al., 2016).
This concern can be addressed by using data-driven approaches that provide a prediction based only on the database previously obtained from experimental, or simulation results. A data-driven model tries to learn the salient characteristics embedded within the system by developing a mathematical relationship between the system parameters rather than solving the physics-based models to describe the system performance. This process is known as machine learning (ML).
ML is one of the branches of artificial intelligence (AI). Currently, there are many machine learning tools that can be used for prediction or classification such as artificial neural network (ANN), support vector machine (SVM), Naïve-Bayes algorithm, random forest, decision tree, logistic regression (LR), K-nearest neighbors (KNN), etc. Any of these tools may be used to develop a machine learning algorithm. Each algorithm is based on its unique strategy in making predictions. Generally speaking, ML algorithms learn from the datasets and try to decipher the salient characteristics within the data that reflect the relationship between the inputs and outputs. Based on the datasets, a mathematical relationship can be generated between the input vector variables and the scalar output variables. The learning process helps improve the relationship by constantly changing the learning parameters to tune the model until the objective function is optimized. The objective functions for each machine learning algorithm is different and needed to be specified accordingly.
Recently, the International Atomic Energy Agency (IAEA) has urged the nuclear community to integrate ML in the industry within the framework of emerging technologies, given its superior capability in handling big-data (IAEA, 2020). In fact, the potential of using ML technology has been explored to estimate some key figures of merit such as the power pin peaking factor (Bae et al., 2008), the wall temperature at critical heat flux (Park et al., 2020), the flow pattern identification (Lin, 2020), to detect anomalies and warn of equipment failure (Ahsan and Hassan, 2013;Chen and Jahanshahi, 2018;Devereux et al., 2019); to determine core configuration and core loading pattern optimization (Siegelmann et al., 1997;Faria and Pereira, 2003;Erdogan and Gekinli, 2003;Zamer et al., 2014;Nissan, 2019), to identify initiating events and categorize accidents (Santosh et al., 2003;Na et al., 2004;Lee and Lee, 2006;Ma and Jiang, 2011;Pinheiro et al., 2020;Farber and Cole, 2020) and to determine of key performance metrics and safety parameters (Ridlluan et al., 2009;Montes et al., 2009;Farshad Faghihi and Seyed, 2011;Patra et al., 2012;Young, 2019;Park et al., 2020;Alketbi and Diab, 2021), and in radiation protection for isotope identification and classification (Keller and Kouzes, 1994;Abdel-Aal and Al-Haddad, 1997;Chen, 2009;Kamuda and Sullivan, 2019), etc. However, it is worth noting that the application of ML in nuclear safety is still limited despite its potential to enhance performance, safety, as well as economics of plant operation (Chai et al., 2003) which warrants further research (Gomez Fernandez et al., 2017). For a more comprehensive review of the status and development efforts utilizing data-driven approaches in nuclear industry, the reader may consult (Gomez Fernandez et al., 2017;Gomez Fernandez et al., 2020).
In this study, an artificial neural network (ANN) is developed to predict the PCT under LBLOCA conditions as a bounding accident scenario. The goal is to develop a fast and cost-effective tool for uncertainty quantification of PCT under LBLOCA conditions using ML. This is achieved by using a database to train the ML algorithm, and once trained and tested, the metamodel can be used as a predictive tool. The database required to train and test the model is generated via the thermal hydraulic system code MARS-KS (KAERI, 2004) within an uncertainty quantification framework using Dakota (Adams et al., 2020). Once proven, the ML technology may be used to help the nuclear designers and/or operators to expedite the decision making process particularly in those situations that involve complex interconnected phenomena during design optimization or in the event of a nuclear accident.

ARTIFICIAL NEURAL NETWORK
ANN is a machine learning model inspired by the biological network of the nerve cells that make up the brain. Fundamentally, the ANN behaves in a way similar to the nerve cells. However, each biological structure is replaced with layers of neurons with a pre-defined architecture that communicates data between the input signals and output signals via weights, biases, and activation functions to find the best weight matrix that best describes the relationship between the inputs and outputs. The ANN structure can be split into three different classes; artificial neural network (ANN), convolutional neural network (CNN) and recurrent neural network (RNN). This research focuses only on ANN.
The ANN generally refers to the modelling of the data through a stack of computational layers. The ANN utilizes the back propagation based on the stochastic gradient descent (SGD) technique that approximates the loss function optimal points which guarantees convergence and terminate at the optimal solution (Dawani, 2020). Currently, the improvised gradient descent techniques such as the adaptive moment estimate (Adam) is being widely used for many ANN applications. The SGD teaches the ANN how to tweak the connection weights and biases in order to converge to the closest mathematical representation of the data at hand. A number of activation functions, such as the hyperbolic tangent function and the rectifier linear unit function (ReLU), may be used to provide signal transformations for each input layer and hence provide better representation of the underlying non-linearity of complex systems.
The ANN is based on a multi-layer perceptron model and can be used for both regression and classification problems. The training process of an ANN is a two-step process. The first is a forward propagation step and involves evaluating the error or loss function. In the second step, the resulting error is propagated backwards through the network to adjust the weights and biases. This process is repeated until no further improvement in the error between the predicted outputs and desired values is achieved. The complexity of the network is determined by the number of hidden layers, the number of nodes in each layer, the type of activation function. The output is predicted by summing the functions within the hidden layers to produce a net input function.
The goal of the training process is to tune the model hyperparameters for better prediction. Those hyper-parameters include: the number of neurons, number of hidden layers, network structure, activation function, type of optimizer, loss function, etc. Common hyper-parameters associated with the regression problems are described in Table 1.
During the optimization process, provisions should be made to ensure a global minimum is achieved rather than a local minimum or saddle point, for which small neural networks are prone. A balance between generalization and fitness to the training data should be achieved to ensure that over-fitting or retaining redundant features in the neural network is avoided.

METHODOLOGY
To achieve the goal of this paper, three main objectives can be identified which are, 1) the thermal hydraulic model development, 2) the uncertainty quantification and database creation and 3) machine learning model development. Each will be delineated in the next subsections.

Thermal-Hydraulic Model Development
This section focuses on the details of the thermal hydraulic model development. In this investigation, the best estimate system code, MARS-KS version 1.4, is used to simulate the nuclear power plant response under LBLOCA conditions. MARS-KS is a multidimensional two-phase thermal hydraulic system code developed by KAERI (2009). The model representation in MARS-KS, including nodalization, boundary and initial conditions as well as the main assumptions will be presented next.

APR1400 Nodalization
First, the details of the APR1400 reactor are described using the system nodalization shown in Figure 1 to represent the key systems and components which includes a reactor pressure vessel, a pressurizer, two loops with four cold legs and two hot legs. The nodalization also includes detailed description of the two steam generators (SG), the main steam lines and associated valves (MSSV, MSIV, ADV, etc). The turbine however, is modeled as a boundary condition. The safety injection system (SIS) is modeled to represent the emergency core cooling system (ECCS) of the APR1400. The SIS modelling is necessary to understand the reactor dynamic behavior in dealing with the LBLOCA scenario and to ensure the safety parameters stay within the safety limits during the different phases of the accident.
The containment building is modelled as a boundary condition and the containment spray system (CSS) are excluded from the nodalization process as it is not required in the analysis of the reactor response under the LBLOCA scenario. The core is divided into a hot channel and an average channel. The hot channel represents the hottest fuel assembly; while the average channel represents the remaining 241 fuel assemblies in the APR1400 nuclear power plant core configuration.
The core channels are modeled in MARS using a pipe hydrodynamic component with a vertical orientation using 20 axial nodes. A solid structural element is attached to the hydrodynamic component to represent the fuel assembly with 20 axial nodes and nine radial nodes. For the downcomer, an annulus hydrodynamic component is used. The annulus is divided circumferentially into six channels at different angle (0°, 60°, 120°, 180°, 240°, 300°, 360°), each having five axial nodes with a vertical orientation.
The same nodalization scheme is used for the upper annulus region, which hosts the entry location of the emergency coolant from the SIS that should be directed to the core. The general pathway of emergency coolant should start from the upper annulus to the downcomer, to the lower plenum, to the bottom of the core and finally to the core itself.
The inlet and the outlet of the core are nodalized using branch components. The core inlet hosts the lower support structure (LSS) and in-core instrumentation while, the core outlet hosts the

Hyper parameters
Typical values fuel alignment plate (FAP), the upper guide structure, the core shroud and the fuel barrel assembly. The crossflow is permitted between the downcomer regions, the upper annulus, the downcomer upstream, between the average channel and the hot channel and between the loop elevation reference regions. The cross flow allows the inventory to move sideways between the hydrodynamic components to represent the secondary cross flow.
At the middle part of the reactor pressure vessel (RPV), a bypass is used to connect the bottom region of the core towards the fuel alignment plate without passing through the core site. The bypass allows liquid to move vertically from the lower plenum to the fuel alignment plate (FAP). To simulate the safety injection system (SIS), both the safety injection tanks (SITs) and the safety injection pumps (SIPs) are modelled in the nodalization.
There are four units of SITs that are modelled using the accumulator hydrodynamic component. The SIT are connected to the upper annulus and is controlled using a combination of logical and variable trips. A set point that refers to the low pressurizer set point pressure is used to navigate the turning on of the SITs. The SIT valve is divided into two part for each train to represent the low-flow and highflow conditions of the actual APR1400 fluidic device. The fluidic device act as a flow regulator in the actual APR1400 plant to optimize the usage of the emergency inventory during the loss of coolant scenario. Meanwhile, the SIPs are comprised of 4 units however, only two units of SIP are available in accordance to the conservative assumption adopted by the APR1400 Design Control Document (DCD) for LBLOCA evaluation. The SIPs will be available only for the break side and the opposite direction across from the break, i.e. located at 60°and 240°, respectively. The SIPs are modeled using a time dependent volume and a time dependent junction instead of pump hydrodynamic components. This configuration allows the user to impose flow boundary conditions and control the velocity of the coolant. Each SIP and SIT is connected to the upper annulus model at different circumferential angles based on APR1400 description. Consistent with APR1400 DCD the decay heat model is based on the ANS-1973 model and a reactivity table is also provided in the code to account for the negative reactivity insertion due to control rod insertion. However, based on the APR1400 DCD, the negative reactivity contribution from the control rod is discredited for conservatism when conducting the LBLOCA analysis. This will allow the SIS capability in managing the accident and maintaining the core integrity to be fully tested during the LBLOCA accident. Regarding the reactor internals, the heat structure components are used and attached to the related hydrodynamic volumes to reflect the heat transfer boundary conditions and architecture of the APR1400 design. Meanwhile to describe LBLOCA scenario, a double ended guillotine break is placed on the cold leg after the pump discharge. This is achieved by incorporating two trip valves to divert the coolant from the vessel to the time dependent volume attached to each of the trip valves when the accident is initiated.

LBLOCA Model Assumptions
Following the DCD recommendation, a double ended guillotine break (DEGB) equivalent to double the area of the pipe with the largest cross section of the RCS, i.e. the cold leg piping is used in this work. The standard internal diameter of the connecting pipes between the pump discharge and the reactor pressure vessel inlet nozzle is 762 mm (30 inches) which corresponds to the break area of 0.456 m 2 . The thermal hydraulic model development is based on several assumptions similar to those reported in the DCD document: 1) LBLOCA occurs at loop B1 near the pump discharge site. 2) Break type is double ended guillotine break (DEGB). 3) Loss of offsite power (LOOP) for the RCPs. 4) No negative reactivity insertion from the control rods. 5) Single emergency diesel generator (EDG) is not functioning causing two out of four safety injection pumps (SIPs) to be non-operable. 6) All safety injection tanks (SIT) are in operation.

Uncertainty Quantification Framework Development
The statistical tool, Dakota (Adams et al., 2020), is used in this work to propagate the uncertainty parameters into the thermal hydraulic model. Dakota is an open source statistical software tool developed by Sandia National Laboratory. It can be used for optimization, sensitivity analysis and uncertainty quantification. The uncertainty propagation process is achieved by developing the uncertainty quantification framework by loosely coupling the best estimate system code, MARS-KS, and the statistical tool, Dakota, via a python script to manage the data exchange process. Several important files such as, the Dakota input file, the python interface script, the MARS steady state file and the MARS transient file are necessary for the uncertainty quantification framework to run smoothly and propagate the uncertainty parameters.

Uncertain Parameters Identification
As indicated in the introduction, the current work explores the possibility of using ML to predict the PCT under the conditions of a LBLOCA, being an important bounding accident scenario. LBLOCA was used in nuclear safety as a design basis for the emergency core cooling system, ECCS, to provide assurance that the ECCS would not violate any of the safety limits and hence preserve the fuel integrity during a loss of coolant accident (LOCA). For LBLOCA, the key performance measure of the ECCS, as defined by the 10CFR50.46 Appendix-K guideline is that the PCT does not exceed the safety limit of 1477 K (2,200℉) to ensure the integrity of the fuel under the accident conditions (Martin and O'Dell). Now for the ML algorithm to be developed and trained, it is necessary to use a database of the most important system parameters (features) that impact the safety parameter of interest, in this case the PCT. Generally speaking, the model can use a database originating from the plant historic data, which is not possible under DBA conditions. Alternatively, it can be generated using simulation results produced by system codes. In this work, the database is created, using the latter approach.
A BEPU analysis is undertaken to generate a database of the system response under LBLOCA. In general, uncertainty quantification can be achieved using either the input uncertainty propagation approach or the output uncertainty propagation approach (Martin and O'Dell, 2008). The former approach will be followed in this work.
To conduct the BEPU analysis, the uncertainty quantification process requires the identification of the uncertain parameters that can impact the PCT. Those uncertainty parameters can be derived from the phenomena identification and ranking table (PIRT). The PIRT describes the key phenomena and processes relevant to the plant's thermal-hydraulic response for a specific accident condition. Most of the PIRT developed throughout the years centered on the LBLOCA cases due to its importance to nuclear safety as a bounding DBA scenario. Several PIRTs have been developed for LBLOCA scenario such as: Westinghouse PIRT (USNRC, 1988) ). For the current project, the investigation will focus on the APR1400 PIRT which is based on the KNGR PIRT, which in turn is derived from the Westinghouse PIRT.
Based on the work of (Lee et al., 2014;Kang, 2016), eight key phenomena were considered in this study. The key phenomena underlying the LBLOCA scenario are gap conductance, energy stored in the fuel, decay heat, rewetting process, reflooding heat transfer, critical flow, pump performance and core reflooding as shown in Table 2. A total of 19 uncertainty parameters have been derived from these key phenomena. These key uncertain parameters and the statistical information associated with each (range and distribution) are listed in Table 2.

Data Pre-Processing
Before propagating the uncertainty parameters into the best estimate thermal-hydraulic system model, it is essential for the uncertainty parameters to undergo a normalization process. The normalization is done with respect to the statistical information available for each uncertain parameter derived from the PIRT. First, the mean value is calculated using the following expression: Next, the upper and lower limit can be scaled using this mean value as follows: And finally for the standard deviation can be calculated as follows: where x i represents value of uncertain parameter in the sample, and x high and x low are the upper and lower values of the uncertain parameter, respectively. For a normal distribution function, the mean, standard deviation and upper and lower limits are required; while, for a uniform distribution, only the upper and lower limits are required. Table 3 shows the uncertain parameters after scaling.

Uncertainty Propagation and Database Generation
With the key uncertain parameters identified and scaled appropriately, they are randomly propagated into the thermal hydraulic system code, MARS-KS using DAKOTA. The goal is to generate a large enough sample that can be representative of the realistic system performance in accordance with the USNRC requirement specified in 10CFR50.46 Appendix-K, i.e. the safety criteria should be satisfied with a probability of 95% and a confidence level of 95%. The 95/95 rule has been recognized by the USNRC to have sufficient conservatism for LBLOCA analyses. Usually a large number of samples are required which can be achieved using the Monte-Carlo random sampling technique. To determine the minimum number of samples required for Monte Carlo method to

ANN Model Development
In this investigation, six ANN model development steps are applied. The six steps are 1) input selection, 2) data splitting, 3) architecture selection, 4) structure selection, 5) model optimization and 6) model validation.
The input data for the ANN model will be a selected set of the uncertain parameters identified earlier and derived using the PIRT for LBLOCA. The feature selection process is important for ANN model development since too many variables will slow down the optimization process and may prevent the model from finding a good solution (Geron, 2019); whereas a few features may not be sufficient for the model to properly learn the system characteristics embedded in the data. A correlation matrix based on Spearman's method is therefore used to identify the key features from the 19 uncertain parameters that impact the PCT the most.
Next, the random sampling technique is applied to split the database into three main categories: one for training (3,022 samples), validating and testing the model (202 and 332 samples, respectively). In order to improve the ANN performance metrics during training, the input and output parameters should have the same scale. Before the propagation of uncertainty, all input parameters have been scaled; hence only the output parameter (peak cladding temperature, PCT) needs to be normalized using the minmax scaling function: where y max and y min are the maximum and minimum values of PCT in the dataset, respectively; while y i represents the temperature to be scaled and y scaled represents the scaled temperature.
Architecture selection refers to the choice of ANN hyperparameters. Since there are a lot of hyper-parameters that can be tuned, choosing the suitable set is a delicate task given the large number of degrees of freedom that the user can manipulate during the ML tuning process. The tuning process of the ML model is therefore more an art than a science and depends on the problem at hand as well as the characteristics of the data. Hence, finding a set of optimal hyper-parameters that provide the best model performance without compromising either its predictive accuracy or generalization capability can be computationally challenging. This is particularly true for hyper-parameters as opposed to other model parameters since the former are not learnt by the model during the training process but must be set manually. Various techniques can be employed to search for the most appropriate hyper-parameters: grid search, random search and Bayesian optimization. In this work, the random search method is used to expedite the convergence. Table 4 shows all the hyper-parameters that are tuned in this study.
ANN tuning is an important step to enhance the model predictability by converging on the most optimum combination of hyper-parameters. In this study, an automatized optimization tool, the Talos (Autonomio Talos, 2019) software, is used. Talos is an open sources software written in Python language. It is compatible with Keras (Chollet, 2015) application programming interface (API) that is suited for the development of artificial neural network (ANN) models. Currently, Talos does not support any other machine learning model other than the ANN architecture and it only supports Keras backend machine learning algorithms.
Initially, the user needs to define the Keras for the ANN algorithm development. Then the user needs to define the search space boundary in the format of key-value pair python dictionary. Afterwards, the scan function is used to run the Talos experiment. The arguments of the scan function include the type of search method (grid or probabilistic), model's name, number of epochs, batch size and search constraints. Talos will generate a list of possible hyperparameters combinations along with their corresponding values that can be analyzed using the built in command such as report and predict functions. The size of the results list depends on the number of parameters defined in the search space boundary dictionary. The analysis process can be done by analyzing the whole or the specific model combinations.
If the user is satisfied with the value of the performance metric generated from the results list, then the deploy function is used to save and call the model from the defined python dictionary path and hence the Talos experiment is complete. As mentioned earlier, the random quantum search method is used to optimize the model. To further reduce the computing time to find the optimized model, both early stopping and window reduction strategies are used. Early stopping prevents the Talos tool from evaluating models that shows unproductive permutation when its metrics are no longer improving; whereas, the window reduction strategy allows the Talos tool to compare the upcoming model with the previously evaluated model based on the specific metrics. Once any of the two criteria is satisfied first, the computation will stop and the results list is generated. Table 5 shows all the parameters used to reduce computation time.
It should be noted that the Talos ability to find the optimized architecture based on hyper-parameters combinations relies heavily on the defined search space boundary. As such, if the final results do not provide reasonable performance metrics, it is essential to redefine the search space dictionary by adding new hyperparameters or retuning their corresponding values in order to improve the model accuracy. Even though the Talos is an automated optimization tool, it is still based on a trial and error process that required extensive knowledge in regards to the behavior of each hyper-parameter towards the ANN model. However, the Talos can expedite the optimization process.
It is worth noting that neither under-fitting nor overfitting is desirable in machine learning. To prevent under-fitting, a large database has been generated to train the model, as many input parameters are used to develop the model and the training time was increased until the cost function is minimized to an acceptable value. To mitigate over-fitting, number of techniques can be used for example: cross validation, regularization, dropout, and early stopping. An overfitting model tends to have good learning metrics during training but performs poorly during the validation process. To avoid this, the data is split into "training", "testing" and "validation" data subsets. The model uses the "training" subset during learning and used the unseen subsets for validation and prediction, respectively. Validation metrics were therefore generated using a subset of the database unseen during the training process to ensure that the ANN model is not overfitting the data. In addition to cross-validation, regularization and early stopping were also used to make sure the model does not memorize the data. Further a dropout layer was placed between the input and the hidden layer, whereby the drop rate is determined solely based on random search algorithm.
The regression type ANN was evaluated using the mean squared error (MSE) that represents the squared difference between the predicted and actual value as shown in Eq. 6: where y pred and y act are the predicted and actual or known value of the dependent variables, respectively; while, n is the number of samples in the dataset. Another important performance metric is the determination coefficient or the R 2 value that measure how well the model predict under sporadic unseen data. The highest value is 1.0 and indicates that the model has strong generalization capability. The R 2 metric is expressed as Eqs 7, 8:

RESULTS AND DISCUSSION
This section is dedicated to the obtained results. The first subsection focuses on the results of the thermal hydraulic model. This is followed by the results of the uncertainty quantification and post-processing of the generated database. Finally, the results of the ML model are presented.

Thermal Hydraulic Model Results
The thermal hydraulic system response is validated against values reported in the APR1400 DCD (KHNP, 2014) for both steady state and transient simulations. The comparison between the MARS and the APR1400 steady state response are shown in Table 6. Based on Table 6, the steady state simulation agrees reasonably well with the plant reference data and the calculated variables are considered to be within the acceptable error limit of less than 5% when compared to the corresponding values reported in the DCD for APR1400. For the transient simulation, key events for the LBLOCA scenario are listed and compared in Table 7.
The progression of the LBLOCA is characterized by three different phases; that is, the blowdown, refill, and reflood. At time t 0 where blowdown phase is taking place, the RCP discharge piping starts to break with double ended guillotine break (DEGB) condition. Instantaneously, the loss of offsite power occurred causing all RCPs to coast down. For this scenario the loss of a single emergency diesel generator (EDG) unit is assumed. This leads to the loss of a single safety system train and hence the loss of two SIPs.  During the blowdown phase, the core uncovers and as a result the heat transfer coefficient drops significantly which causes the fuel cladding temperature to increase reaching a maximum of 1165.1 K (891.95°C) at ∼10 s which is still below the acceptance criterion of 1447.6 K (1204.44°C) as illustrated in Figure 2. Later on, as the decay heat drops the cladding temperature starts to decrease. Further reduction in the cladding temperature is observed due to the large condensation that occurs at the upper guide structure (UGS) and reactor vessel upper head. The condensate passes through the reactor core in reverse direction (from top to bottom) into the lower plenum and headed towards the downcomer region.
The cool down effect continues until the coolant inventory from the upper head of the reactor vessel is depleted. Once the top quenching is over, the fuel-cladding starts to reheat again approximately ∼17 s after the break due to the accumulation of decay heat. After this time, there is no longer any cooling mechanism for the core except the one initiated by the safety injection system (SIS) during the reflooding phase.
The refill period starts when the emergency coolant reaches the lower plenum of the vessel and stabilizes till is completely filled and ends when the water level in the lower plenum vessel reaches the core inlet. The SIT high-flow injection starts when the core pressure reaches the set point pressure, approximately 4.2 MPa (∼43 kg/cm 2 ) at ∼ 17 s after the pipe break.
The emergency coolant will flow from the downcomer, towards the lower plenum and up into the core. However, even though the SIP is initiated earlier compared to the SIT because of the higher set point pressure approximately, 12.5 MPa (128 kg/cm 2 ), the SIP starts to inject the emergency coolant at 48.0 s with a 42.0 s delay time. As such after the blowdown phase, the core is cooled initially by the UGS inventory before being assisted by the injection from the SIT followed by the SIP. This delay accounts for the time required for signal actuation time as well as the time needed to start the SIP.
The reflood phase is further subdivided into two phases: the early reflood and late reflood phases. During early reflood phase, sufficient injection from the SIS helps the downcomer to be filled up relatively quickly. However, due to the limited space available in the downcomer combined with the excess amount of emergency coolant, some of the inventory is bypassed through the pipe break regions causing inventory loss. Nonetheless, the amount of coolant available is still sufficient to maintain the core integrity.
The downcomer is nearly filled at around 50 s after the break. At this time, the water level in the downcomer region starts to stabilize as the fluidic device shifts from the high-flow injection to the low-flow injection. The low flow injection stabilizes the downcomer level as the water rises up into the core. The quench front moves vertically upwards to quench the whole core during this times. The maximum flow rate is reached approximately 30 s after the blowdown phase; while the lowflow injection will continue until 200 s. The SIP assist the low flow injection to cover and quench the core. This phase ends when the entire core is quenched from the bottom up gradually and the fuel rod temperature is slightly above the coolant saturation temperature.
During the early refill phase, the steam binding phenomenon may occur which may slow down the process. However, this effect diminishes once the vapor from the upper section of the core no longer received the de-entrainment liquid from the bottom part of the vessel at the surface of the quench front. Afterwards, the steam binding effect starts to disappear after some time which allows the reflooding phase to resume and the reactor core to be filled with water again.
The late reflood phase is marked by the SIT depletion as the SIT low-flow injection comes to an end. During this time, the task of replenishing and providing the emergency inventory for core cooling and core coverage process is achieved solely by the SIPs. The downcomer water level is maintained at a relatively constant value. Both the hot core and average core are finally quenched around 150 s without violating the fuel acceptance criterion of 1477.0 K.

Uncertainty Quantification Results and Database Post Processing
Using the uncertainty quantification framework, the databases for the machine learning model was generated. The Monte Carlo random sampling technique is used to generate 5,000 runs in order to acquire a large data base for the machine learning model. The simulation is conducted using a single PC platform with 3 GHz Intel ® Xeon ® Gold CPU processor, with 64.0 GB random access memory (RAM), 24 parallel processors and a Windows 10 platform. The time taken to complete the simulation is approximately 3 days.
It is worth noting that, out of the 5,000 samples, only 3,556 samples were successful and used to train the machine learning model. To ensure the number of datasets is enough to represent the 95% probability and 95% confidence criteria, the mean value for the PCT is averaged over the number of samples. As seen in Figure 3, the Monte Carlo simulation stabilizes after approximately 2000 runs, hence, a sample size of the 3556 is adequate for to meet the criteria. Figure 3 shows the spread of PCT which follows a normal distribution with a mean value of 1169 K. The majority of the data are well below the PCT safety criterion of 1477 K. However, there are two data points that lie very close to the safety limit with values of 1462.3 K and 1451.8 K.
As discussed earlier, 19 uncertain parameters are propagated into the thermal hydraulic model using the Dakota uncertainty framework in order to generate the PCT response under LBLOCA scenario. The independent variables are the uncertain parameters (UPs); while the dependent variable is the peak cladding temperature (PCT). However, each UP has different degree of influence towards the PCT. Hence, a sensitivity study is conducted to assess the correlation between the uncertain parameters and the PCT. Spearman's rank correlation, which is a non-parametric measure of the statistical dependency between the variables, is used for the sensitivity analysis. Using the Spearman's correlation coefficients, the strength and the direction of the relationship between the independent and the dependent variables can be evaluated using the following expression: where x i refers to the input variable, i, x is the input variable's mean, y i is the output variable,i, and y is the output variable's mean. Figure 4 shows the Spearman's correlation matrix generated from the database. A positive correlation coefficient means when the independent variable is increased, the dependent variable also increases. While, a negative correlation coefficient means that when the independent variable is increased, the dependent variable decreases. The sensitivity study proves that the uncertainty parameters are independent from each other.
From the sensitivity analysis, a threshold value needs to be defined in order to choose the most significant parameters to reduce the number of inputs for the machine learning model. By selecting a threshold of ±20%, any correlation coefficient higher than the threshold value is deemed to be strongly correlated with the PCT either proportionally or inversely. However, selecting the threshold value is subjective as such it should be tested to find the best value possible.

ANN Model Results
The ANN algorithm has been successfully developed and trained using the database created via the uncertainty quantification framework. After tuning, the model was deployed using the Talos optimization tool. 20 models that differ in architectures and hyper-parameters are generated. Among those 20 models, the best model is selected based on the lowest validation metric, i.e. MSE. The final ANN structure recommended by the Talos optimization tool is composed of an input layer, three hidden layers and an output layer Table 8 shows the hyper-parameters for the selected ANN model achieving the best performance.
The model is trained using the training subset of the available PCT database. Next, the validation dataset is used to measure the model ability to learn the salient characteristic of the data. Once the meta-model has been trained, its performance is tested using an unseen subset of data to make the predictions. Figure 5 shows a comparison between the model predictions of the PCT and the known values produced using the uncertainty quantification  framework for all 19 key uncertain parameters. The ML metamodel predicts the PCT with reasonable accuracy (MSE 0.0039); however, it tends to underestimate the high temperatures which is problematic from a safety point of view. This may be due to the fact that unnecessary data from other uncertainty parameters may confuse the ML algorithm and hence impact the model accuracy. However, the obtained MSE depends on the chosen input parameters. To assess the model sensitivity to the number of input parameters, the model is tested with multiple sets of inputs (5, 7, 9, 11, 17, and 19 UPs) that correspond to different threshold values for the Spearman's correlation coefficient (i.e. level of importance to PCT). Table 9, summarizes the performance metrics of the model with different number of uncertain parameters used as inputs to the meta-model. After dimension reduction, the various cases investigated are compared to each other using a number of performance metrics: the determination coefficient (R 2 ), the mean square error (MSE), the mean arithmetic error (MAE) and the mean logarithmic squared error (MLSE).
Judging by both R 2 and MSE, the model with 9 to 11 parameters achieves reasonable performance. Considering the results presented in Table 9, the lowest possible MSE is approximately 0.00185 which is obtained using 11 inputs with an R 2 value of ∼0.89. When, the machine learning model is tested with nine inputs, approximately similar results are obtained with a loss function, MSE, of ∼0.00186 at an R 2 value of ∼0.93. Outside this range (9 to 11 inputs) the ML model performance deteriorates. Given, the aleatory nature of the ANN model optimization which is based on random optimization, the optimum number of input parameters is expected to be within the range 9 to 11 variables with slight variation in performance metrics results (MSE ∼0.002, R 2 ∼ 0.9). It is worth noting that from a safety perspective, it is better to tune the model for high temperatures. Figure 6 shows the prediction results in comparison to the actual known PCT values with different number of input parameters. Clearly, the lower number of input parameters does not capture fully the relationship between inputs and outputs embedded in the database. On the other hand, the higher number of input parameters may include unnecessary details that may confuse the model. One would suspect an optimum number of input parameters may exist for better prediction capability as evidenced by the results shown in Table 9.

CONCLUSION
The aim of this work is to develop a machine learning (ML) algorithm that is capable of accurately predicting the key safety parameter, PCT, under LBLOCA scenario. The algorithm was trained using a database created using the best estimate code MARS-KS V1.4, with uncertainty quantification using the statistical tool, Dakota to propagate the uncertainty parameters through the thermal hydraulic model. A Monte Carlo sampling method was used to generate 3,556 successful samples to train the ML model using 19 key uncertain parameters. The Monte Carlo simulation converged after 2,000 samples yielding the required average PCT in consistency with the USNRC requirement of 95% probability and 95% confidence interval. An ANN model was developed, trained and optimized using the provided database. The ANN model was successfully tuned using the Talos optimization tool to predict the PCT with high accuracy. The optimum model is chosen based on the desirable objective function and the validations metric, MSE. A model with 9 to 11 inputs best represents the data and can be used to predict PCT accurately with a MSE of ∼0.002 with R 2 value of ∼0.9.
This study successfully shows that ANN can be used as a surrogate to the thermal hydraulics MARS-KS model to predict the PCT value for the LBLOCA scenario using only the key uncertain input parameters with reasonable accuracy. For future work, the framework developed for this project can be used for uncertainty quantification of other key safety parameters such as Departure from Nucleate Boiling Ratio (DNBR) under LBLOCA or other critical scenarios. This is a preliminary step towards developing an expert support system that can be used to guide the operator actions under the stressful accident conditions. As a bounding accident scenario, the analysis of the LBLOCA case paves the way to using machine learning as a decision making tool for design extension conditions as well as severe accidents.

DATA AVAILABILITY STATEMENT
The datasets presented in this article may be available based on personal communication with the corresponding author.
Requests to access the datasets should be directed to aya.diab@kings.ac.kr.