Quantitative resilience assessment of chemical process systems using functional resonance analysis method and Dynamic Bayesian network

The emergent hazards of chemical process systems cannot be wholly identified and are highly uncertain due to the complicated technical-human-organizational interactions. Under uncertain and unpredictable circumstances, resilience becomes an essential property of a chemical process system that helps it better adapt to disruptions and restore from surprising damages. The resilience assessment needs to be enhanced to identify the accident's root causes on the level of technical-human-organizational interactions, and development of the specific resilience attributes to withstand or recover from the disruptions. The outcomes of resilience assessment are valuable to identify potential design or operational improvements to ensure complex process system functionality and safety. The current study integrates the Functional Resonance Analysis Method and dynamic Bayesian Network for quantitative resilience assessment. The method is demonstrated through a two-phase separator of an acid gas sweetening unit. Aspen Hysys simulator is applied to estimate the failure probabilities needed in the resilience assessment model. The study provides a useful tool for rigorous quantitative resilience analysis of complex process systems on the level of technical-human-organizational interactions.


Resilience
The technology is moving forward, and the process safety measures are increasingly becoming automated to facilitate the safe environment in the chemical process plants. However, automation has also led to complexities and accident occurrence due to technical-human-organizational interactions. Nowadays, the safety goal of industrial plants is not only preventing accidents but also dealing with the daily disturbances and variabilities to sustain the system in the normal operational state. To achieve this goal, the industrial system management should be more flexible and resilient [ [18], [43]]. Resilience engineering application in the process plants is particularly important because of the high level of complexity that arises from the non-linear interconnections and high level of variability and uncertainty which arises from this complexity [9].

Dynamic Bayesian network (DBN) for resilience assessment
Yodo et al. [24] developed a DBN model for resilience assessment. They viewed resilience as a constant value though a probabilistic resilience metric was used. Cai et al. [7] considered the performance and time-related properties in resilience assessment, in which high availability and a short time to recovery characterized the system's high resilience. They applied DBN to predict the future values of availability and their temporal variation by employing the data available at the current state. However, they evaluated resilience as a parameter dependent on the failure and repair rates of the system components (i.e., the technical factors). Tong et al. [44] proposed to measure resilience using system functionality. A DBN model was developed to quantify the probabilities of system functionality states based on the four attributes of resilience, namely absorption, adaptation, restoration, and learn. Zinetullina et al. [48] used a bow-tie (BT) diagram to support the development of the DBN model and applied it to analyze the resilience of a two-phase separator operating in a harsh environment. In both approaches, resilience is defined as the ability of the system to sustain high functionality under disruptive effects or restore the low functionality state to a high functionality state. Kammouh et al. [25] also used DBN for resilience assessment, where resilience was defined in terms of fast recoverability, reduced vulnerability, reduced failure likelihood, robustness, or reduced failure consequences. Kammouh et al. [25] measured resilience at pre-, intra-, and post-disruption stages.
The majority of the previous studies considered resilience as a static property of the system. It is an oversimplified assumption for complicated process systems when interacting with human operators, organizational factors, and various disruptions.

Functional resonance analysis method (FRAM) for resilience assessment
FRAM, firstly introduced by Hollnagel [17], can qualitatively model how a complex system functions considering socio-technological interactions [35]. Sawaragi et al. [41] conducted a FRAM analysis for the aircraft accident in Columbia in 1995 and investigated the main reasons and conditions behind that accident. De Carvalho [10] used FRAM to analyze the reasons behind the mid-air collision between two aircraft in Amazonian Sky in 2006. Later, FRAM was employed in other domains. Shirali [42] exploited FRAM for the identification of the emergent risk in a process plant. Rosa et al. [39] used FRAM for identifying the hazards associated with the performance variability of the couplings in the construction industry. The couplings are the interconnections between the main functions (or activities) of the system that contribute to the accomplishment of the main goal of the system. Patriarca et al. [35] presented a semi-quantitative approach for FRAM analysis with the employment of Monte Carlo Simulation (MCS) for air traffic management systems. Having a quantitative FRAM model facilitates the identification of the most critical coupling that needs prominent attention.
Furthermore, it helps to maintain the safe state of complex systems effectively and cost-efficiently. Before Patriarca et al. [35], Rosa et al. [39] proposed to quantify FRAM and identify critical couplings with the employment of AHP and Subject Matter Expert's (SME) judgment. However, the limitation still exists in the subjectivity of the judgements for the variability of the function parameters. This limitation could be mitigated to some extent by the application of MCS to the variables.
The existing resilience assessment methods are ineffective in modeling emergent disruptions and complicated technical-human-organizational interactions. The nonlinear interactions among the system components, human and organizational factors, and their interdependencies with function variabilities are not well treated in the resilience assessment models (Andersson et al., 2002;Levenson, 2004;[43]). Conducting FRAM analysis at the initial stage of the quantitative resilience assessment of a complex system could facilitate a more rigorous analysis of the system's operational state and identify concrete solutions to enhance system safety [39].
FRAM could be considered as a viable method for resilience assessment due to two main reasons: (1) FRAM identifies the root causes on the level of non-linear technicalhuman-organizational interactions [43]. (2) A FRAM model consists of activities performed both on and by the system. The resilience parameters can better be characterized by actions of the objects rather than by their presence ([35]b).

Objecitve of the study
The present study aims to develop a quantitative resilience assessment method for chemical process systems, in which technical-humanorganizational interactions are considered. To the authors' knowledge, this paper will, for the first time, integrate FRAM and DBN to facilitate the assessment of the ability of a chemical process system, its operators, and associated organization to respond to unprecedented disruption, recover from unexpected damages, and learn from those events.
In this paper, the FRAM is mainly applied to describe system functions and identify potential variabilities and their critical coupling. This study helps to advance the knowledge of the approaches where the FRAM has been applied with other methods for the resilience engineering of chemical process systems. The DBN is used to quantify the dynamic resilience of process systems. The DBN includes nodes to represent the absorption, restoration, adaptation, and learning parameters. Nodes are structured according to the outcomes of the FRAM analysis. The main steps of the methodology will be described, and its application will be demonstrated on a two-phase separator." Fig. 1 describes the methodology employed in this study. First, the FRAM analysis on the process system will be accomplished in conjunction with the MCS. Resultantly, the critical coupling will be identified. For the critical coupling, the attributes of the resilience will be developed, and the DBN model will be built. The conditional probability tables of the nodes in the DBN will be filled with probabilities of failures (PoF) collected from the literature and Subject Matter Experts (SMEs), and calculated with the Aspen Hysys simulation. DBN model will be updated with additional safety measures, and the generated resilience curve will be compared with the previous resilience profile. After that, the first developed DBN model will undergo sensitivity analysis to identify the most influential safety measures. The following sections describe the main steps of the proposed method.

Step I. FRAM modeling
At the beginning of the FRAM analysis, it is necessary to identify the goal of the FRAM analysis, specifically whether it is accident investigation analysis or the risk assessment [35].

Step I.1: develop FRAM model for the system
Firstly, the main activities or functions performed by or on the system to reach a particular outcome are defined. Functions are interconnected via the aspects placed in the corners of the hexagon, according to the Structured Analysis and Design Technique (SADT) [35]. Each function consists of six aspects: input, output, precondition, resources, time, and control. The output of one function (upstream function) serves as one of the five aspects of the other function (downstream function) in a coupling. The other function is the downstream function. The six aspects of a function are as follows [35]: ■ Input (I)-the aspect that initiates the function ■ Output (O)-the outcome of the function, which serves as the input for the downstream function ■ Precondition (P)-the condition that should be performed before the initiation of the downstream function ■ Resource (R)-that the function needs to produce the output ■ Time (T)-temporal requirements or constraints obliged to the function, such as start time or finish time. ■ Control (C)-what performs the control and monitoring of the function to achieve the set outcome.

.1 Identify the variability
Each function has its variability. Variability could be classified as positive if it reduces the risk of possible accident development or negative, otherwise. According to Hollangel [19], variability could be classified based on multiple phenotypes, such as time, precision, flow rate, speed, duration, direction, object, and force. In this study, variability manifestation based on time and precision will be used. The reasons are user-friendliness and readability of the analysis and the universal application of those phenotypes for any function.
The variabilities of "time" are classified as on time, too late, too early, and not at all. Not at all is used for the cases when the function is performed that late that its outcomes are no longer attractive. For "precision", the variabilities are classified as precise, acceptable, imprecise, and wrong. Based on the SME's judgment, those classifications are ranked with numbers and further used in MCS [35].

.2 Aggregate the variability
This step deals with the variability of the couplings. Variability of the couplings may arise as a result of each function's internal variabilities or the impact of the upstream function. When the coupling produces large negative variability, this results in the resonance of the connected functions and identifies the critical path. These couplings are called critical. The critical coupling defines the leading cause of the accident in the accident investigation scenario, or the main factors contributing to the hazard development in risk assessment.
If the upstream function generates positive variability, this results in the dampening of the variability of the downstream function [35]. In step I.2.4, the MCS is integrated to identify the exact critical coupling, which results in the malfunction of the two-phase separator of the acid gas sweetening unit.

(3) Step I.2.3 Manage the variability
This step aims to suggest measures to increase positive variabilities and decrease negative variabilities. Improvement recommendations need to be made for the identified critical coupling to avoid the accident from happening, or in case of an accident, to restore the system to the high functional state. Based on the critical coupling identified in the next step, we will develop the DBN model consisted of the absorption, adaptation, and restoration parameters for resilience assessment. (4) Step I. 2

.4 MCS analysis
MCS is the tool for the mathematical modeling of multiple systems using random sampling to obtain outcomes as close as possible to reality (Zio, 2013; [31]). MCS method is mostly employed for solving complicated and multidimensional problems and equations. For instance, it is used for estimating differential and integral equations, for addressing the deterministic and stochastic problems, for simulating the random case development scenarios. One of the frequent applications of the MCS is in reliability and risk analysis of the engineering systems. The outcome of the MCS analysis facilitates the estimation of the probability of failures (PoFs) [ [12], [29], [33], [47]].
In this research work, the MCS is used to estimate cumulative variabilities in the FRAM model. The idea of MCS for estimating the cumulative variability in the FRAM Model was taken from the work of Patriarca et al. [35]. In their work, they identified the critical couplings as ones with a high tendency to cause an accident. They also provided the resonant path that emerges as a result of the variability accumulation for the neighboring functions to the critical couplings. The critical couplings were determined by setting the critical cumulative variability value and the confidence level of 95%. The couplings for which 5% of the cumulative variabilities were greater than the setpoint value, equal to 24, was classified as critical. Hence, the corresponding recommendations were provided by Patriarca et al. [35] to avoid accident occurrence.
In this paper, a similar approach of MCS for cumulative variability estimation is employed to identify the critical couplings in FRAM. However, in this paper, MCS also serves as a bridge between the FRAM model and the DBN model for the quantitative resilience assessment of the process unit.
From the SMEs' judgment in the work of Patriarca et al. [35], the discrete probability for each of the four variability states based on the timing and precision were defined (Tables 1-3).
In this study, the variability of timing was assumed to be too late, and variability of precision was assumed to be of acceptable precision. For illustrative purposes, the random samples from normal distribution were assigned to the values of the variabilities based on timing and precision.
The amplification factors are distributed based on the timing ( ) ij T and precision ( ) ij P relations of the upstream and downstream couplings. The output of a function may vary with respect to timing (e.g., being too early, on time, too late, or omitted). It may also vary for precision (e.g., being precise, imprecise, or acceptable). Based on the effect of the variability of the upstream function conveyed to the downstream functions, the amplification for both timing and precision are assumed in the ranges in Table 4. Amplification factors describe the impact that the upstream function has on the downstream function connected to it. For example, if the amplification factor on "timing" is equal to 0.5 (i.e., the dampening effect applies), the time aspect of the upstream function being "too late" may result in the downstream function to be on time. If the amplification factor is equal to 10, then the upstream function that is too late will result in the downstream function being too late either. If the amplification factor is equal to 1, the upstream function will not affect the variability of the downstream function. For illustrative purposes, the normal distributions were assigned to the values of amplification factors.
The overall variability is defined as a product of the variabilities of upstream functions based on timing (V j T ) and precision (V j P ) (Eq (1)) To calculate the overall variability, each value of too late timing ( Table 2) and acceptable precision (Table 3) were multiplied correspondingly.
The cumulative variability (CV ij ) is the joint variability of the coupling. It is estimated by multiplication of the overall variability and amplification factors for timing and precision assigned for the corresponding coupling.
In this study, we conducted the MCS in Excel to estimate the OV and CV. For the variables, including timing, precision, and associated amplification factors, we assumed they follow normal distributions. We derived the mean and standard deviation values for the scores of "too late timing" and "acceptable precision" based on Tables 2 and 3. First, using Eq. (1), the OVs were calculated based on randomly generated samples of variability values of timing and precision. Then, using Eq.
(2), the CVs were obtained based on a random sampling of the amplification factors. 1000 iterations were applied. It is worth noting that the number of iterations may affect the results. This aspect should be investigated in future work.
At the next stage, the set-point and the confidence level are assigned in Excel for 1000 generated random numbers with normal distribution using Eq. (2). The confidence level in this work is set to be equal to 95%. After the generation of 1000 random numbers, five percentiles value of those random numbers is compared with the predefined set-point (in the case study, a set-point of CV = 24 was assigned). If the five-percentile number is higher than the set-point, the corresponding coupling is set as the critical. The critical value of 24 can be customized by the model user, and the confidence level of 95%. It depends on the user's requirement for the criticality level. It is also worth investigating how the choice of CV and confidence level may affect the results of the entire resilience assessment model.
The further recommendations and work on risk prevention or accident recovery are prioritized, starting from the coupling with the highest five percentile number.

Step II Hysys simulation
The PoFs are obtained from Aspen Hysys Dynamics simulation for the process unit using the strip charts for the parameters characterizing the performance of the process system of interest [21]. From the strip chart, the time range of the system malfunction or being in the out of normal operating range is divided by the total time range of operations. The PoF characterizing the functioning of the process system could also be extracted from Aspen Hysys Sensitivity analysis. In this case, simulation is run at the specified range of the several input parameters. The output is monitored for the number of scenarios out of the normal range scope. Then, the PoF is estimated as the ratio of the amount of the abnormal outcomes to the total number of outcomes.

= PoF
No of abnormal outcomes Total no of outcomes The estimated PoFs are then inputted into the corresponding nodes of the DBN model developed in the next step.

Step III generate DBN parameters
DBN is a probabilistic graphical model. It has the same structure and principles as the Bayesian network, with an exception that it allows estimating the joint probabilities of the variables with time. Thus, using DBN allows estimating the joint probabilities at time t, t + 1, t + 2, and so on ( [23]; Neapolitan, 2004). It consists of the nodes (variables) and arcs that link the nodes and allow deriving the interdependencies of the nodes based on the conditional probabilities. The nodes are classified as the parent nodes, child nodes, and the root nodes. The root nodes are assigned for marginal probabilities. The rest are assigned for the conditional probabilities. The joint probability of conditional nodes is calculated at the child node dynamically (Eq. (1)) [23].
In DBN and BN, the joint probability could be updated with the addition of new nodes.
Having these characteristics, DBN enables predicting the spatial and temporal evolutions of systems probabilistically [27]. In resilience   Table 4 The ranges for amplification factors [37].
The variability of the downstream function is amplified by the output of the upstream function =1 The variability of the downstream function is not affected by the output of the upstream function <1 The variability of the downstream function is dampened by the output of the upstream function engineering, it was first applied by Yodo et al. (2016) for the assessment of the reliability and restorative properties of the system. However, their study was limited to an estimation of restoration only; however, resilience consists of four parameters, which are absorption, adaptation, restoration, and learning. Furthermore, in their research, the resilience was considered as a constant term. In contrast, the resilience is a characteristic that changes following the changes of variabilities of components and either internal or external disturbances [44].

Step III.1 establish linkages between the identified critical couplings and the attributes of resilience
The Absorption, Adaptation, and Restoration nodes are complemented with the nodes corresponding to the most critical coupling obtained as a result of the FRAM and MCS. For example, if the identified most critical coupling are the inaccurate performance of the level controllers and interlock system and late or incorrect response of the operator in the control room, the nodes for Absorption, Adaptation and Restoration can be selected corresponding to only the performance of the level controllers and interlock system and the performance of the operator in the control room.

Step III. 2 develop conditional probability tables (CPTs)
The CPTs of the nodes can be filled based on the outcomes of the simulations, historical data, and data from the surveys.
For the nodes in DBN, two states were assigned in the CPTs, "High Functionality" and "Low Functionality". As for assigning the marginal probabilities, the "Low Functionality" section in CPT was filled with the PoF values, estimated from Aspen Hysys simulation, historical data, and data from the surveys. "High Functionality" section was filled by subtracting the estimated probabilities (in the corresponding "Low Functionality" section of the CPT) from 1.
As for the parent nodes of "Absorption", "Adaptation" and "Restoration", consisting of the different combinations of the functionality states of the corresponding nodes, the conditional probabilities were assigned based on the expert judgment for each separate scenario. They are starting from the highest conditional probability when all constituting nodes have a high functionality and finishing with the lowest conditional probability for the case when each constituting node has a low functionality state.

Steps IV and V resilience assessment before and after improvement
The DBN model will be developed based on the Absorption, Adaptation, Restoration, and Learning parameters for the most hazardous case identified by the FRAM model.
For Absorption, Adaptation, and Restoration nodes, additional safety nodes are added. Then the resilience profile is compared with the previous one. The updated resilience model should present the enhanced resilience properties for a longer time until the functionality drops. A higher value of the functionality should be attained after the restoration stage.

Step VI sensitivity analysis for identification of the most influential nodes in DBN model
Sensitivity analysis facilitates identifying the nodes that have a high impact on the state of the system's functionality. This will assist in the application of the specific safety measures for the system. First, the sensitivity analysis was conducted in GeNie (https://www.bayesfusion. com/genie). As a result, the most influential nodes for the "State of Functionality" target node are identified. Next, for each influential node, the evidence was changed to the low functionality state for the ten-time steps sequentially. Based on the obtained results, the parameters having the highest impact on the system's resilience reduction of the system are identified.

Case description
This case study will focus on the quantitative resilience assessment of a two-phase vertical separator of the acid gas sweetening unit. Sour gas treating unit is a part of the oil and gas preliminary treating plants. It accomplishes the absorption of the sour gasses (mainly carbon dioxide (CO 2 ) and hydrogen disulfide (H 2 S)) with diethanolamine (DEA) from the gas coming after the crude oil stabilization. Crude oil stabilization is the process during which oil coming from the well is separated into three phases (mainly gas, oil, and water). The pressure is reduced in the stabilizer from the high upstream pressure to the normal operating pressures not to harm gas, oil, and water treating process units. Fig. 2 presents the process flow diagram of the acid gas sweetening. Aspen Hysys, a chemical process simulator, was used to simulate this process. First, sour gas is depleted from process hydrocarbons in the two-phase separator FWKO TK. Then, dehydrated sour gas goes to the absorption column, where it is scrubbed using the DEA solution and leaves the column as a sweet gas for further production of the sales gas. The rich DEA stream, with CO 2 and H 2 S, then is separated from the hydrocarbons in the separator Flash TK and passes to the regenerator column, where the products, lean amine and acid gas are produced. Lean amine stream is refluxed to the absorption column and acid gas serves as the feed for the Sulfur Recovery Unit [21]. In this study, the simulations were conducted with a focus on the analysis of the valve failure scenarios. Based on the simulation results, the PoFs of valve actuator failure (as a node contributing to the internal disruption of the system) were estimated.

FRAM-DBN model development
The resilience assessment was conducted on a two-phase oil and gas separator located before the coalescing filter and acid gas contactor column. The two-phase separator has the level controllers and interlock system for the liquid outlet stream, pressure controller and interlock system for the gas outlet, pressure safety valves activating in case of the overpressure, alarms interconnected with 2-out-of-3 alarm counting mechanism. Methanol is injected at the inlet of the feed stream to avoid the hydrates forming in the pipes. The performance indicators are monitored and controlled distantly by an operator in the control room and by an operator at the unit area.

FRAM analysis
The FRAM analysis was developed according to the steps described in the methodology part (Section 2). First, the functions connected with the operational functions of the two-phase separator was listed. Then, the functions listed in Table 5 were interconnected with each other with the assistance of the six aspects. The values of the amplification factors in Table 5 were assigned based on expert judgment. For example, coupling "Pressure indicators" being an input parameter to the "Pressure Safety Valves" has an amplification factor based on the timing of 60 ± 10, where 60 is the mean value and 10 is the standard deviation. The values of mean and standard deviation were applied in MCS to define the normal distribution formula in the Excel file. The value of the amplification factor based on timing, in this example, means that the function of "pressure indicators" acts "too late". This results in the output of the "pressure safety valve" being "too late".
Then the MCS method was employed following the process described in the work of Patriarca et al. [35] with the application of myFRAM (https://functionalresonance.com). The scenario of having too late response of the acceptable precision was assumed. Next, the variabilities for each of the function couplings and amplification factors were assigned for each coupling. The cumulative variabilities have been estimated for each coupling according to the Eq. (2).
For each of the variables in Eqs. (1) and (2), including V , and , ij P 1000 random positive numbers were generated by sampling from the normal distributions with the parameters defined in Table 5 to compute the cumulative variability. Setting the 95% confidence level and set-point value of 24 for the overall variability, the outcome of each of 31 couplings was compared with the set-point value. If the 5-percentile value of the given coupling was greater than 24, it was considered as critical (Fig. 3 developed using the FMV, https:// functionalresonance.com). In the following case, the coupling of the highest criticality was considered, which is the coupling of the operator of the control room controlling the level and interlock system, resulting in the cumulative variability value at 5 percentiles of 38. Thus, imprecise and/or too late a response of the operator at the control room to the level alarms and indicators has the highest potential to cause an accident. Furthermore, the imprecise performance or too late response of the level alarm and interlock systems has the highest potential to result in the accident or two-phase separator failure.

DBN model
Based on the outcome of the FRAM analysis, the DBN model is developed for the coupling with the highest criticality. The node "Disruption" lists the potential malfunctions that may cause accident development. The nodes "Absorption" and "Adaptation" contain the current safety absorption and adaptation measures applied for the twophase separator. "Restoration node" includes the recommendations for the enhancement of the 2-phase separator performance.

Application of simulation for extraction of PoFs
The PoF for the "2-phase separator performance" node was estimated using the sensitivity analysis module in the Aspen Hysys simulator. The sensitivity analysis is a one-at-a-time (OAT) method. The input parameters were the ranges of the pressure, temperature, and flow rate of the feed. The "Nested" state-input type was specified for the input parameters. Thus the sensitivity analysis was automatically run throughout the specified ranged first for pressure, then for temperature and finally for flow rate. The monitored parameter for each simulation run was the dihydrogen sulfide concentration in the "Sweet Gas" stream, leaving the amine contactor. As a failure, the case of having the H 2 S concentration higher than 4 ppm was assumed. Consequently, the PoF equivalent to 0.122 was estimated using Eq. (4). A total of 700 simulation runs were carried out, which was the maximum number of runs that the program allowed to perform for the specified ranges. It was a pragmatic decision to stop at 700 iterations due to the software limitation and may affect the accuracy of the results. It is worth noting that the users may consider conducting a sensitivity analysis of the number of simulation runs to identify an appropriate number. The PoFs of the other nodes were identified using Aspen Hysys Dynamic simulation scenarios.

Scenario 1. Valve actuator Failed Open.
The Aspen Hysys Dynamic simulation for estimation of PoF was completed for the two-phase separator located before the absorption column (Fig. 5). The separator products were the sour gas routing to the absorption column and processed hydrocarbons routing to the stabilization unit. The PoF evaluation was accomplished with the strip charts generated as a result of the dynamic simulation (Fig. 6). For filling the disruption node "Valve actuator Failed Open", the valve VLV-102 failopen scenario was simulated in Aspen Hysys Dynamic (Fig. 5). The normal liquid percent level was set to the value of 20-50%. However, the level of liquid in the tank dropped to 4.13% when the valve VLV-102 failed open (Fig. 6). Employing Eq.

Scenario 4. Heat introduction and fire mitigation scenario
In this scenario, the PoF is estimated for node "Adding valve and interlock at the inlet stream of the separator" given fire in the twophase separator (Fig. 7). The heat was introduced to the vessel. With the employment of the event scheduler and cause and effect matrix, the closing of all valves during the fire extinguishing was simulated with subsequent fire relief via the pressure safety valves before the pressure is returned to the normal operating pressure. After that, the closed valves were opened, and the operation continued at the recovered conditions.
The normal operating pressure range for this scenario was the pressure range of 490-520 kPa. The valves VLV-100, VLV-101, and VLV-102 were set to close at the pressure of 655 kPa. Based on the Given the PoFs obtained from the simulation, the resilience curve for the 2-phase separator is obtained. Fig. 8 presents the resilience of the separator given the disruption in terms of the malfunctions that may occur due to inaccurate performance or too late a response of the level alarm and interlock systems or the operator working in the control room. It takes 9 h for disruption to bring the separator to the lowest level of resilience, i.e., 0.48. Furthermore, it takes 56 h for the system to restore to 0.94 -the highest resilience after the adaptation and restoration measures applied to the separator.
Additional nodes were considered for "Absorption"," Adaptation" and "Restoration" (Table 6). A new DBN model ( Fig. 9) for the twophase separator was built. With the additional nodes, the resilience profile is expected to enhance. The resilience profile displays stronger "Absorption", "Adaptation" and "Restoration" characteristics, the structure of the DBN model can be validated.
From a comparison of the results presented by Table 7, it is observed that the lowest value of resilience has increased from 48.34% to 50.00% with the addition of the new parameter for the "Absorption" node. This designated a slightly more robust absorption properties of the system.  The final resilience state is slightly increased at the updated DBN model, and the time for 90% recovery is shorter than in the initial DBN model. Both characteristics are evidence of the strengthening of the "Adaptation" and "Restoration" parameters of the system. Thus, the DBN model is partially validated. It is worth noting that the obtained results (regarding the safety measures recommendation) are for demonstration purposes because no independent validation was performed.

Sensitivity analysis for identification of the most influential nodes in DBN model
Sensitivity analysis was carried out to identify the nodes having the highest impact on the state of system resilience. Identification of the critical nodes assists in the allocation of specific safety measures to the system. Table 8 lists 11 scenarios that may contribute to the variation of the resilience of the system and the results of the sensitivity analysis. Table 8 indicates that the highest impact on the resilience drop occurs at low absorption conditions. However, with the strong adaptation and restoration characteristics, the system restores to the normal operational state (90% of recovery) in the shortest time. The most significant impact on the resilience drop is caused by the effect of the poor operator skills and not functioning alarm and interlock system. Scenarios 4-11 show the resilience reduction to the value of around 40% in approximately 9 h, representing the strong absorption properties or well-thought inherent safety design of the developed model.

Methods comparison
In this section, the previous work on the quantitative resilience assessment [48] is compared with the current study. In the previous work, the elements of BT were developed based on their association with the resilience attributes, namely, absorption, adaptation, and restoration. The developed BT was then converted into the DBN model with the nodes of absorption, adaption, and restoration as the immediate parent nodes linked to the node of "state of functionality", which was used to measure the system resilience. In the present study, instead of subjectively identifying the factors, FRAM was applied to identify the critical factors affecting the system resilience. It is a more objective approach.
Initially, it was assumed that disruption occurs due to the electricity outage and subsequent malfunction of the electric heat tracing of the three-phase separator operating under the harsh cold conditions. Comparatively, in this study, the possible reason for the separator malfunction was not assumed but found through the FRAM analysis with the MCS approach. Thus, the current research enables identifying the causes of the possible accidents, inventing and applying the particular measures for either avoiding the accident from happening or in case the accident occurred to adapt to it or recover from it in a fast manner. With this approach, the safety measure currently applied and newly proposed could be analyzed for their effectiveness for the specific failures with the employment of this method. Thus, the approach offers both a practical and cost-efficient approach to ensure safe operations in the plant.
Furthermore, some studies (e.g., Kammouh et al. [25]) did not develop DBNs for disruption and recovery attributes. They did not consider the system's learning capability as an attribute of its resilience. The present work identifies the most vulnerable technical-human-organizational interactions within the system to support DBN development, which involves the sub-networks for modeling disruptions and recovery and considers learning to be one of the attributes of resilience. Moreover, the present work enables investigating the effectiveness of the supplementary safety measures to enhance the resilience attributes.
Bellini et al. [4] have developed a Q-FRAM approach that quantifies the system's resilience at a specified time by using the System Resilience Index (SRI). SRI is estimated by applying the Choquet integral for the capacitive variabilities assigned for the function couplings, referring to the corresponding resilience "cornerstones" (anticipate, respond, monitor, learn). The value of SRI > 1 designates the risk of accident   development as a result of the resonance in the poor performance of the functions, and thus characterizes the system as non-resilient.
With the integration of FRAM and DBN, the present study provides the opportunity to obtain the resilience profile as a probabilistic and time-dependent evolution of a system's functionality state, affected by the critical function couplings.

Conclusions
Dynamic quantitative resilience assessment is a research problem of great interest in the sphere of resilience engineering for process systems. This requires the rigorous analysis of the root causes of the accidents and mishaps, and the resilience attributes of the process systems on the level of technical-human-organizational interactions. The current study demonstrated the applicability and effectiveness of the integration of the FRAM and DBN for the quantitative resilience assessment of complex process systems consisting of the technical-humanorganizational aspects. The proposed approach can identify the root causes of an accident with the employment of FRAM. The essential resilience attributes strongly associated with these causes are considered in the development of the DBN to assess the temporal changes of the resilience of a complex process system. The most important safety measures can also be identified with sensitivity analysis. This study supports that resilience should be assessed in both probabilistic and temporal terms. It investigates the usefulness of process simulators as a process data source for the estimation of failure probabilities. For future work, we are planning to investigate the opportunity of learning the structure of the Bayesian network from various data (e.g., process data, management and organizational data, operational environmental data) for dynamic resilience assessment.

Declaration of Competing Interest
The research being reported in this paper titled "Quantitative resilience assessment of chemical process systems using Functional Resonance Analysis Method and Dynamic Bayesian Network" was supported by Nazarbayev University. The authors of this paper have the Fig. 9. Updated DBN model with additional parameters (orange nodes presented in Table 6). IP ownership related to the research being reported. The terms of this arrangement have been reviewed and approved by the Nazarbayev University in accordance with its policy on objectivity in research.