A dynamic multi-agent approach for modeling the evolution of multi-hazard accident scenarios in chemical plants

In the chemical industry, multi-hazard (toxic, flammable, and explosive) materials such as acrylonitrile are stored, transported, and processed in large quantities. A release of multi-hazard materials can simultaneously or sequentially lead to acute toxicity, fire and explosion. The spatial-temporal evolution of hazards may also result in cascading effects. In this study, a dynamic methodology called “ Dynamic Graph Monte Carlo ” (DGMC) is developed to model the evolution of multi-hazard accident scenarios and assess the vulnerability of humans and installations exposed to such hazards. In the DGMC model, chemical plants are modeled as a multi-agent system with three kinds of agents: hazardous installations, ignition sources, and humans while considering the un- certainties and interdependencies among the agents and their impacts on the evolution of hazards and possible escalation effects. A case study is analyzed using the DGMC methodology, demonstrating that the risk can be underestimated if the spatial-temporal evolution of multi-hazard scenarios is neglected. Vapor cloud explosion (VCEs) may lead to more severe damage than fire, and the safety distances which are implemented only based on fire hazards are not sufficient to prevent from the damage of VCEs.


Introduction
The past decades have witnessed an increase in the number, size, and diversity of chemical plants due to the increasing population and the increasing requirement for products (energy, chemicals, commodities, and food, etc.) [1,2]. The rapid expansion of the process plants and infrastructures brings huge economic benefits while unavoidably increasing the exposure to major hazards caused by hazardous materials in chemical industrial areas, resulting in human losses, environmental damage and economic losses [3][4][5][6][7]. Major hazards such as fire, explosion, and toxic release arising from loss of containments may occur due to intentional or unintentional causes [8][9][10][11]. Intentional hazards are security-related threats, including terrorist attacks, sabotage, thief, etc. Unintentional hazards consist of accidental hazards (e.g., corrosion, fatigue, mechanical damage) and natural hazards (e.g., earthquake, flood, lightning). In hazardous chemical areas, fire is the most frequent hazard (44%), followed by explosion (36%). Toxic release without fire and explosion accounts for 20% of all major accidents and toxic substances are involved in almost 30% of these accidents [12]. Besides, chemical industrial areas are usually congested with hazardous storage tanks, complex piping, high-pressure compressors, and separators in which a loss of containment (LOC) event may lead to cascading effects and multiple hazards.
All the major hazards of fire, explosion and toxic release can be simultaneously or sequentially present in one disaster due to the evolution of hazards. Many catastrophic disasters in the past two decades originated from the hazardous release of process vessels and evolved to VCEs and finally fires. On October 23, 2009, a large VCE happened at the Caribbean Petroleum Corporation (CAPECO) terminal in Bayamón, Puerto Rico, during the offloading of gasoline from a tanker [13]. The subsequent fires triggered by the explosion lasted about 60 h and resulted in significant damage to 17 of the 48 petroleum storage tanks and other equipment [13]. On November 28, 2018, a Vinyl chloride release in a chemical plant at Zhangjiakou (China) caused a VCE outside the chemical plant, triggering fires on tank trucks and leading to 23 fatalities and 22 injuries [14].
In light of these past disasters and due to the severe consequences of unpredicted hazards, modeling the spatial-temporal evolution of hazards originating from release of hazardous materials in industrial areas is essential for protecting staff, nearby residents and emergency rescuers [15][16][17]. For example, in the Tianjin port disaster in 2015, which was caused by a spontaneous ignition of nitrocellulose, many of the emergency rescuers were killed in the disaster due to an unpredicted evolution of the fire to an explosion. Besides, the disasters caused by natural hazards (Na-tech) can make emergency response more difficult due to the damage to safety barriers and other infrastructures, resulting in more severe consequences [18][19][20][21][22]. To avoid such catastrophic disasters, many post-accident analyses have been conducted to predict the overpressure induced by explosions [23][24][25][26][27] and vapor cloud dispersion [28][29][30]. Besides, a lot of work has been done on vulnerability assessment of installations to VCEs, risk assessment of domino effects caused by VCE [31][32][33][34][35][36] and domino effects triggered by fire [37][38][39]. Regarding the evolution of fire, the time to failure of equipment exposed to fire is critical for assessing the vulnerability of installations. As a result, dynamic methods were used to assess the vulnerability of installations exposed to fire and fire-induced domino effects [40][41][42][43][44]. Among these dynamic tools, the dynamic graph approach and the dynamic Bayesian network approaches are able to model the spatial-temporal evolution of domino effects caused by fire and visualize the escalation paths of fire [8,11,40]. Monte Carlo simulation has also been applied to address the evolution uncertainty of domino effects [45].
Compared to the research devoted to VCE or fire evolution, little attention has been paid to the evolution of possible toxic release, VCE and fire in a catastrophic disaster and the assessment of human exposure to multiple hazards. Dai et al. [46] analyzed the vulnerability of tanks exposed to fire and explosion in different accidents, but overlooked the possible evolution between different hazards. He and Weng [47] studied the synergic effects of multi-hazard on vulnerability assessment but ignored the dynamic evolution process. The evolution of toxic release to VCE and vice versa is a dynamic process along with a vapor cloud dispersion [48]. The present study, therefore, aims to establish a dynamic methodology for human and facility vulnerability assessment considering the spatial-temporal evolution of multiple hazards: toxic release, VCE, and fire. In our study, chemical plants are modeled as a multi-agent (component) system [49,50] through the application of dynamic graphs. Graph-based methods have been used for domino effect analysis, vulnerability and reliability analysis, and resilience assessment. Besides, Monte Carlo simulation [51] is used in this study to solve the dynamic multi-agent model. Consequently, both the uncertainty of ignition and the uncertainty of evolution of different hazards are taken into consideration in the present study. The model and algorithm are developed in Section 2. A case study is provided in Section 3 to show the application of the developed methodology. A discussion based on the results of the case study is presented in Section 4. The conclusions of this study are summarized in Section 5.

Methodology
Graph-based methods are commonly used and are effective tools to analyze multiple interacting agents in a system [52,53]. In a graph, the agents are modeled by nodes and their dependencies are represented by edges [38,44,54]. As a result, graph-based methods provide a visible structure (graph or network) to represent the complex agent interactions while agent-based modeling focuses on agent behaviors (e.g., attributes and interactions), making it very flexible to model socio-technical systems [50,55,56]. Graph metrics such as betweenness and closeness have been used to assess domino effects and the vulnerability of installations subject to fire and explosion hazards [38,57]. The time-dependent escalation of fire can also be modeled by dynamic graphs [41]. The dynamic graph approach is able to model the dynamic evolution of fire hazard while static graph methods provide merely a snapshot of the whole process at once. When it comes to modeling the evolution between different hazards, it is difficult to address the uncertainty in hazard evolution by merely using the dynamic graph approach.
Compared with analytical methods, Monte Carlo simulation is widely used to model uncertainties that cannot be easily accounted for due to the intervention of random variables, avoiding complex mathematical calculations [58][59][60]. In this study, a new methodology is developed based on dynamic graphs and Monte Carlo simulation to model the complexity and uncertainty of hazard evolution. The methodology is called Dynamic Graph Monte Carlo (DGMC). The DGMC is defined as a dynamic graph with time-dependent parameters and random parameters in which Monte Carlo simulation is used to solve the model.

Modeling
In order to model the hazard evolution process and thereby dynamically assess human vulnerability exposed to the possible toxic cloud, heat radiation, and overpressure, we define a Hazard Evolution Graph (HEG) based on the developed DGMC method. The HEG can be defined as a dynamic graph with nine-tuple, as shown in Eq. (1).
represents the evolution time of hazards starting from a hazardous release (t 1 =0). The dynamic graph HEG is sliced into G static graphs by these time nodes. The HEG parameters are updated at each time node t g due to the change of hazards, human states, and installation states. If an ignition occurs, t 2 is equal to the ignition time (IT). IT is a random variable that depends on the number of ignition sources, ignition effectiveness, and vapor cloud dispersion, etc. The IT is equal to zero if the released materials are immediately ignited. In the chemical industry, the likelihood of immediate ignition is always determined by the autoignition of flammable substances and the static discharge caused by the release [48]. If the ignition is delayed, the possible ignitions caused by different ignition sources are considered as independent events, and the ignition probability of a single ignition source depends on the ignition effectiveness and the period that the ignition source is covered by the flammable vapor, as follows [61]: f I,k represents the cumulative ignition probability caused by the ignition source k. ω denotes the ignition effectiveness of the ignition source. t I,k represents the time that ignition source k (within the flammability limit) is covered by the flammable vapor. To determine t IS , the vapor cloud dispersion model developed by [62] can be adopted, as follows: where R t is the radius of the area in which cloud might be ignited at time t; c E is an empirical constant approximately equal to 1; V is the volume flow rate of the flammable gas; ρ is the vapor density relative to air. This dispersion model is suitable for low-wind conditions and neglects the effects of obstacles on dispersion.

Numbering hazardous installations
is a set of nodes representing the hazardous installations that may be involved in the evolution of hazards.

Numbering human positions
N = [m + 1, m + 2, m + 3, …, m + n] is a set of nodes denoting the human position that may be affected by toxic release, fire, or overpressure hazards.

Numbering ignition sources
K = [m + n + 1, m + n + 2, m + n + 3, …, m + n + k] is a set of nodes denoting the ignition sources that may cause the ignition of a flammable vapor cloud.

Node states
S is a node parameter indicating the state of installations, humans and ignition sources at the evolution time T. According to possible major hazards in industrial areas and the vulnerability characteristics of hazardous installations and humans, five states of hazardous installations, three human states, and three states of ignition sources are defined, as shown in Tables 1-3. As shown in Table 1, the state of "operational" is an initial state while the state of "extinguished" is a terminal state. The states of "release", "fire" and "VCE" are harmful to other installations and humans. If a release occurs at an installation, the state of the installation changes from "operational" to "release". The state of "fire" is caused by an immediate ignition while the state of "VCE" results from delayed ignition. Table 2 shows human states including one initial state "safe" and two terminal states "injured" and "dead". Table 3 lists the four states of ignition sources, including one initial state (inactive), one terminal state (ignited) and one intermediate state (ignited). It should be noted that all the foregoing states are time-dependent and may be updated with the spatial-temporal evolution of the hazards.

Physical effects
E is a set of directed edges denoting the physical effects that may cause damage to hazardous installations or be harmful to humans. In this study, the heat radiation induced by fire, the overpressure caused by VCEs and the toxicity induced by toxic vapor are considered. There are six kinds of directed edges: the heat radiation from installation nodes to installation nodes or human nodes, the overpressure from installation nodes to installation nodes or human nodes, and the toxic effects from installation nodes to human nodes or ignition nodes.

Acute intoxication
C is a set of edge parameters from release source to human denoting the concentration of toxic vapor at human positions. The acute intoxication of exposed humans caused by a toxic cloud depends on the toxic concentration (C t ) and exposure time (t e ). The probit function for acute intoxication is used to quantify the death probability due to human exposure to toxic vapor, as follows: where c 1 , c 2 and c 3 are constants that vary with different toxic substances. These constants for different toxic substances can be adopted from the Green Book [63]. Y h is the probit value of human vulnerability exposure to toxic gas. Aa a result, the death probability (f t ) caused by acute toxicity can be obtained using Eq. (5): ϕ is the cumulative distribution function of the standard normal distribution. It should be remarked that besides the toxic concentration and exposure time, other factors such as demographics (e.g. ages) and Personal Protection Equipment (PPE) are not considered in this formula.

Damage induced by VCEs
P is an edge parameter denoting the overpressure generated by VCEs when the flammable vapor is ignited by an ignition source. The commonly-used overpressure estimation methods include the TNT equivalent method, the Baker-Strehlow method and CFD simulation, etc. The TNT equivalent method is a simple approach based on TNT explosion mechanism to calculate overpressure, which neglects the effects of space configuration, ignition sources and flammable gas distribution and thus may underestimate the overpressure. The Multi-Energy method is developed for gas explosions, dividing the explosion as a number of sub-explosions and addressing the effects of congestion levels, ignition and gas distribution in obstructed areas. Netherlands Organization for Applied Scientific Research (TNO) recommended the Multi-Energy method for overpressure calculation in quantitative risk analysis [61]. Consequently, the Multi-Energy method [62] is adopted to calculate the overpressure P obtained by different installations and humans. According to this method, the overpressure induced by VCEs depends on the strength coefficient (c s ) and the scaled distance (r s ). The c s is a constant parameter (1-10) determined by the congestion of the chemical plant. In view of the severe consequences caused by past VCE accidents in chemical industrial areas, the most conservative value of 10 is usually assigned to c s in risk assessment [48]. The scaled distance r s is characterized by four parameters: the mass of the flammable vapor (M f ), the combustion heat of the vapor (ΔH) and the distance between the calculation point (installations or humans) and the explosion center (r), as shown in Eq. (6): Then the overpressure obtained by different installations and humans can be derived from a blast chart [64]. The damage probability of hazardous installations and death probability caused by overpressure (f p ) can also be calculated by the application of probit functions, as follows: where c 4 and c 5 are constants, as shown in Table 4.

Table 1
States of hazardous installations.

State Description
Operational The hazardous installation is not physically damaged and is operational.

Release
The hazardous installation is physically damaged, resulting in the loss of containment of hazardous materials and/or poisoning humans nearby. Fire The installation is on fire due to immediate ignition, causing heat radiation on humans and/or other installations. VCE The installation's loss of containment induces a vapor cloud explosion due to delayed ignition. Extinguished The installation is physically damaged but does not generate any hazardous effects. The human is injured due to exposure to toxic gas, heat radiation, or overpressure. Dead The human is decreased due to exposure to toxic gas, heat radiation, or overpressure.

Table 3
States of ignition sources.

State Description
Inactive Flammable vapor is not present at the ignition source, or the concentration of the vapor is out of the flammability limit. Active Flammable materials are present at the ignition source, and the concentration of the vapor is between the lower and upper flammability limits. Ignited The ignition source has ignited the flammable vapor. The damage probability of installations or the death probability of humans induced by overpressure is calculated by the application of Eq. (7) while substituting Y p for Y t .

Damage induced by fires
H is a m × (m +n) matrix representing the heat radiations generated by nodes in "fire" states. h i,j is an element of the matrix denoting the heat radiation induced by an installation i in a "fire" state to installation or human j, as follows: H is not a square matrix because people can only receive but not generate heat radiation. Considering possible synergistic effects [8] induced by multiple installations in the "fire" states, the heat radiation received by a node j (Q j ) can be calculated as: Considering that installation i starts receiving effective heat radiation at evolution time t g , the residual failure time T f tg i of the installation can be calculated as [65]: where T f tg i is the residual failure time of installation i at t g (min); The values of c 6 , c 7 , c 8 , and c 9 are shown in Table 5.
If the received heat radiation Q i,t g evolves to Q i,Tg+1 at t = t g+1 and the T f tg i > t g+1 − t g , the installation i will not be physically damaged at t g+1 , and the residual time to failure of installation i in the "heat up" state at the time t g+1 will be updated according to the superimposed effect model [41].
The vulnerability of humans exposed to heat radiation is estimated by using the exposure time and the received heat radiation (Q). Consequently, the probit value of human vulnerability exposure to multiple fires can be estimated as: The heat radiation received by people (Q) varies with the number of hazardous installations in the "fire" state during the spatial-temporal evolution of hazards. Therefore, the hazardous effects caused by heat radiation on humans at different time periods should be superimposed (e.g., the superimposed effects of heat radiation on human vulnerability). At evolution time t G , the probit value (Y tG f ) can be estimated as: Finally, the death probability (P f ) can be calculated by the cumulative distribution function ϕ in Eq. (5).
During fire escalation, emergency response measures such as firefighting may effectively extinguish the fires and thus prevent the evolution of hazards. Thus a cumulative log-normal distribution function is adopted to model the uncertainty of the time needed for an effective emergency (t e ) [8]: u is the mean value, and σ 2 is the variance of log t e . Fig. 1 shows the state transitions and physical effects due to different states in HEGs. As shown in Fig. 1, Dotted lines represent state transition of nodes and solid line denote physical effects caused by a node to other nodes. Humans may be injured due to acute toxicity caused by installations in the "release" state, heat radiation induced by installations in the "fire" state, as well as overpressure caused by installations in the "VCE" state. Since human vulnerability to heat radiation and toxic vapor depends on the intensity of the physical effects and the exposure time, humans may die after a period time of exposure. But the overpressure induced by VCEs may induce an immediate death since the death likelihood caused by explosions is determined by the intensity of explosion regardless of the exposure time. In a hazard evolution caused by overfilling, a human may suffer from different hazards at different subsequent times.

Graph update rules
To further illustrate the graph update rules, Fig. 2 shows an example of a HEG with 9 static graphs. As shown in Fig. 2, a hazardous installation's state changes from "operational" to "toxic release" when a toxic release occurs at the installation (e.g., T1 in Fig. 2a) due to accidental events, natural events, or intentional attacks. If the released material is ignited immediately, the installation's state will change to "fire" and induce heat radiation on humans and other hazardous installations. Otherwise, a vapor cloud may form and disperse along with the vaporization of the release material, resulting in acute toxicity, as shown in Fig. 2b.
As the vapor cloud continues to spread, the ignition source may change from an "inactive" state to an "active" state (Fig. 2c). As a result, a VCE may occur when the vapor cloud is ignited, resulting in the damage of hazardous installations and casualties. As shown in Fig. 2d, T2 and T4 are damaged by the VCE while T3 is not. At time t 4 , both H1 and H2 become exposed to the overpressure caused by the VCE but the injury of H1 may be more severe than that of H2 since the former suffers from acute toxicity as well. Simultaneously, T2 and T4 may be on fire because the explosion can release a lot of heat and energy which increases the likelihood of immediate ignition at the two damaged tanks. As shown in Fig. 2e, the two tanks in "fire" states induce synergistic effects of heat radiations on H1, H2, and T3. H1 may die at time t 5 (Fig. 2f) while H2 may die at time t 6 (Fig. 2g). The state of T4 changes from "fire" to "extinguished" at t 5 due to the burn out of flammable substances. Finally, the fire at T2 is extinguished and T3 survives since the escalation is blocked by firefighting actions. The evolution ends since there is no escalation. Fig. 2 uses a dynamic graph to represent an evolution process of hazards originating from a toxic release. Evolution uncertainties such as the ignition time and the death probability of humans cannot be fully considered by listing all the possible hazard evolution paths. Besides, the evolution may be more complex when it comes to real chemical storage areas with multiple ignition sources and many hazardous installations. As a result, Monte Carlo simulation is employed to generate DHEGs, addressing the time-dependencies and uncertainties in the hazard evolution process. Fig. 3 shows the developed algorithm based on the HEG model and the Monte Carlo simulation.

Simulation algorithm
According to the simulation algorithm, at first, basic data is inputted, including the number of iterations (NI max ), the industrial layout, release Table 5 The values of c 6 , c 7 , c 8 , and c 9 [65].
Atmospheric tank − 2.67 × 10 − 5 1 − 1.13 9.9 Pressurized tank 8.845 0.032 − 0.95 0 scenarios, possible ignition sources, etc. Next the ignition type is determined by random sampling. If it is a delayed ignition, the ignition time (IT) and the ignition source are determined by random sampling based on Eq. (2). Death caused by toxic vapors can also be determined using Eq. (4). At the IT, the HEG updates and the curves in the graph represent overpressure. The overpressures suffered by humans and hazardous installations can be calculated based on Eq. (6). As a result, the fatalities and the subsequent fires caused by the overpressure can be obtained by random sampling based on Eq. (7). The graph will be updated again at the IT, and the curves in the graph will then represent heat radiation. The heat radiation can be calculated by the application of ALOHA, and then the time to failure of hazardous installations can be determined using Eqs. (10) and (11). If there is an immediate ignition, the calculation procedures for explosion and toxic vapor are neglected, and the heat radiation is immediately calculated. During the fire escalation period, the graph is updated when a new fire occurs, or an existing fire is extinguished when the evolution is over. The above calculations are repeated NI max times. Finally, the death probability and failure probability of installations during the dynamic hazardous evolution are obtained. Besides, the possible evolution paths, evolution time nodes, expected DIT can also be determined using the simulation.

Application of the methodology
The DGMC methodology consisting of the HEG model and the simulation algorithm was developed in Section 2. To demonstrate its application to a dynamic hazard evolution, an illustrative case study is used in this section.

Case study
A chemical storage facility including 37 chemical storage tanks (T1-T37) in three tank areas (I, II, III), 5 possible human positions (H1-H5), and two possible ignition sources (S1-S2) is considered in this study. The layout of the chemical storage facility is shown in Fig. 4. Table 6 summarizes the main characteristics of the storage tanks considered in the dynamic vulnerability analysis.
An overfilling of acrylonitrile at T1 with a filling rate of 100 kg/s was considered as the primary scenario. The release of acrylonitrile can result in acute intoxication and the subsequent explosions and fires may lead to human's exposure to overpressure and heat radiation. The released acrylonitrile vaporizes and disperses around. The ambient temperature is 0 • C and the wind speed is equal to zero. According to the vapor cloud dispersion model in Eq. (8), the ignition source S1 is active after 5.1 min while S2 is active after 2.8 min. The autoignition probability P ia is zero and the ignition probability due to static discharge P is is estimated as 0.02 given the minimum ignition energy of 0.16 mJ [66] and the autoignition temperature of 481 • C [67]. The possible heat radiations induced by tanks and the burning rates of fires are calculated through the ALOHA software. The number of iterations (NImax) is set as 10 5 in which the computation time is 3.9 min, and the average deviation of two computations is lower than 0.001. In terms of large chemical plants with many installations, the computation time will increase though being still acceptable. Taking the case with 150 tanks [8] as an example, the computation time would be 21 min. This computation is conducted by an ordinary personal computer (Intel (R) Core (TM) i7 CPU, 8GB RAM). The computation time can be deviously reduced by using a computer with better computation performance.

Results
Due to the spatial-temporal evolution of hazards, humans may get exposed to different hazards. The death probabilities caused by different hazards at H1-H5 are shown in Fig. 5.
As shown in Fig. 5, both the total death probabilities at H1 and H2 are around 0.99, indicating that people at H1 and H2 will die due to the toxicity and fire. The total death probabilities at H1 and H2 are approximately equal to their death probabilities caused by toxicity or fire. It demonstrates that people at H1 and H2 may simultaneously or sequentially receive multiple hazards. Humans at H1, H4 and H5 are mainly threatened by toxicity and fire while the main hazards at H2 include toxicity, VCE, and fire. The hazards at H3 are dominated by toxicity since it is far from T1, and it is not in the storage tank area. Although there is a long distance between T1 and H5, fire can escalate from tanks nearby T1 to the tanks close to H5. The results can be explained by analyzing the spatial-temporal evolution between the hazards. In this case, explosion is not the main cause of fatalities, but it can cause injures such as eardrum rupture. The probability of eardrum rupture at H1 and H2 is 0.34 and 0.68, respectively. According to these results, different kinds of PPEs can be assigned to humans in different positions. For example, humans at H3 only need to take a gas mask while protective clothing to protect against potential heat is also needed for humans at H1, H2, H4 and H5 due to possible multi-hazard. More details of PPE are presented in the Discussion. Fig. 6 shows the failure probabilities caused by fires and explosions of the 27 hazardous tanks due to the overfilling of acrylonitrile at T1. The failure probabilities of T1-T6 in tank area I is around 0.98 since they are close to the release source. The failure probabilities of tanks in storage area II obviously decrease with increasing the distance between the release sources and the tanks (e.g., T7-T11, T12-T15). The failure probabilities of tanks in tank area III are much lower than those of the tanks in storage areas I and II since they are located farther from the release source T1. Besides, the fires cannot escalate from area II to area III due to the safety distance between these areas. However, the explosion at tank areas I and II can damage the tanks in area III, possibly resulting in fires. It indicates that the safety distances provided for preventing fire escalation is not sufficient to prevent the damage caused by VCEs.
The tanks around T1 (i.e., T1-T8 and T12) are more likely to be directly damaged by explosions (VCEs). These damaged tanks may result in multiple fires, and the fires may escalate spatially as well as     temporally, resulting in the damage of other tanks such as T10 and T11. Fig. 7 exemplifies one of the possible hazard evolutions: acrylonitrile starts releasing from T1 at t = 0; the released acrylonitrile forms a flammable vapor cloud, and the vapor cloud disperses around. People at H1 and H2 get exposed to the toxic gas at t = 1.3 min and t = 3.0 min, respectively. The ignition source S2 is active at t = 2.8 min while S2 is active at T = 5.1 min. Then the vapor cloud is ignited by S2 at t = 7.6 min, resulting in a VCE. During the vapor dispersion process, people at H1 and H2 die due to acute toxicity. At t = 7.6 min, 11 tanks (T1-T9, T11, T12, and T15) are damaged and catch fire due to the overpressure caused by the VCE. The fires rapidly escalate to T10, T13 and T14. Finally, people at H4 die due to heat radiation. As shown in Fig. 7, T10, T13 and T14 are damaged by multiple fires at t = 10.4 min due to synergistic effects. The result demonstrates that fire escalation after a delayed VCE is usually inevitable. Despite the fact that the fire cannot escalate from tank storage area II to tank storage area III due to the physical distance between them, there is a low probability for the tanks in storage area III to get damaged (0.0-0.1) and for people at H5 to die (0.09) in a fire: This case may happen if a VCE causes damage to the tanks in the storage area III and subsequently triggers fire in the area. It indicates that the tanks are more vulnerable to VCEs and the safety distances solely based on fire risk assessment would be ineffective for VCEs.
Due to the hazard evolution, the death probabilities of humans at H1, H2 and H4 change over time. Fig. 8 shows the cumulative probabilities of death at H1, H2, and H4. People at H1 start inhaling toxic vapor at t = 1.3 min when the toxic vapor spreads to H1 and the death probabilities increase over time due to the amount of inhaled toxic vapor. The cumulative probability of death at H2 increases from t = 3.0 min when the people at H2 when people at H2 begin to be exposed to toxic gas. At t = 7.6 min, multiple tank fires are induced by the VCE, resulting in hazardous effects at H1, H2 and H4. As a result, the cumulative probability of death at H4 starts to increase due to the induced heat radiation and that at H1 and H4 further increase due to toxicity and heat radiation.

Discussion
The case study in Section 3 demonstrates that multi-hazard chemicals (e.g., acrylonitrile) in the process industry can simultaneously or sequentially lead to multiple hazards to humans and installations due to the cascading effects. The results are consistent with the characteristics of the disaster in Zhangjiakou, China (2019). Based on the case study, this section discusses parameters that may have considerable effects on human vulnerability and the multi-hazard evolution.

Atmosphere parameters
Atmosphere parameters such as ambient temperature and wind have an impact on the evaporation of liquid hazardous materials and the dispersion of vapor clouds. Both wind speed and direction have been shown to affect the vulnerability of humans and facilities during cascading effects [41]. Since large disasters caused by overfilling usually occur at low wind conditions [68], only the effects of temperature on the death probability of humans and the failure of hazardous tanks are discussed in this study. As shown in Fig. 9, both the death probability of humans and the failure probability of tanks increase with the rise of ambient temperature. As shown in Fig. 9a, the death probability of humans at H3 only slightly increases with rising temperature, indicating that acute toxicity is not sensitive to temperature. The reason is that the rise of ambient temperature increases the toxic gas concentration while decreases the delayed ignition time (DIT). The failure probabilities of tanks in tank area III (T16-T27) are much lower than those in area II since the physical distance between tank area II and III becomes a barrier for fire escalation between the two areas. When the ambient temperature rises, the likelihood of tank failure in tank area III increases rapidly since the overpressure caused by VCEs increases so does the damage likelihood of tanks in area III. In general, humans with lower death probabilities and tanks with lower failure probabilities are more sensitive to ambient temperature. It can be demonstrated that the damage effects increases with increasing ambient temperature.

Flow rate
The flow rate is a key parameter for characterizing a loss of containment. Since the amount of hazardous material released is proportional to the flow rate, the death probability increases with increasing flow rate, as shown in Fig. 10a. The death probability at the five positions slightly rises with the increase of flow rate. Fig. 10b shows the effects of flow rate on the failure probabilities of hazardous tanks. All the failure probabilities of the 27 tanks display an increasing trend when increasing the flow rate. Therefore, it can be demonstrated that a larger release is more likely to result in more severe consequences, no matter what hazards the release causes. By comparing Fig. 10 with Fig. 9, it can be demonstrated that the effect of flow rate is much smaller than that of ambient temperature since the evaporation rate of toxic vapor is more sensitive to ambient temperature than the flow rate. For example, the concentration of toxic vapor increase by 34% if the flow rate increases from 80 kg/s to 160 kg/s while that increases by 170% when the ambient temperature increases from 20 • C to 40 • C.

Probability of immediate ignition
An immediate ignition can lead to fire but the increase of the probability of immediate ignition (PII) decreases the likelihood of VCEs and toxicity. Fig. 11 illustrates the effects of FII on the total death probability of humans and the total failure probability of the tanks. Both the death probability and failure probability decrease with an increase of the PII. It indicates that the damage caused by the fire is lower than the damage caused by VCEs and toxicity. The total risk caused by hazardous materials towards individuals and facilities may be underestimated if only Fig. 7. One of the hazard evolution processes including toxic release, VCE and fire. fire hazard is considered in the hazard evolution. The tanks close to the release source are more sensitive to the PII because an increase of PII sharply decreases the damage probability caused by VCEs, as shown in Fig. 11b.

Emergency response
The time needed for effective emergency response (t e ) is essential for mitigating the consequences of a disaster by cutting or delaying the spatial-temporal escalation of hazards. In this study, a log-normal distribution is used to model the uncertainty of t e . The effects of the average values of t e (μ) on the total death probability and failure probability are shown in Fig. 12. As shown in Fig. 12a, The death probability of humans far from the release source increases with increasing the time needed for emergency response while the death probability at H1 and H2 close to the release source is not sensitive to the t e . It is more difficult to rescue people near the release source via emergency response. The failure probability of tanks also increases with the increase of tee. Besides, tanks are more vulnerable to fire are more likely to survive since the emergency response can largely decrease the possibility of fire escalation.

Personal protection equipment
PPE is used to minimize human exposure to hazards, such as respirators, protective clothing, helmets, goggles, glasses and other garments that are designed to prevent the wearer from hazards [69]. In the case study, the main hazards for humans are toxic gas and heat radiation. Respirators are commonly used for protecting humans from toxicity by filtering out toxic gases in the air [70]. In case PPE is available in chemical plants, a human response time of 5 s for humans to respond to hazardous scenarios is considered [63]. Fig. 13 shows the death probabilities of humans with respirators in different positions. Compared to Fig. 5 (without respirators), the death probability caused by acute intoxication is largely reduced. For example, the death probability at H1 decreases from 0.99 to 0.05, decreasing by 95%. However, the total death probabilities at H1, H2 and H4 don't decrease since people at these positions are also threatened by heat radiation. As a result, thermal protective clothing is needed to provided enough time for humans to escape from fire scenarios [63,71]. Fig. 14 shows the death probabilities under the protection of respirators and thermal protective clothing.
As shown in Fig. 14, the total death probabilities at H1-H5 are decreased due to the application of respirators and thermal protective clothing. In the case scenario, both respirators and thermal protective clothing are needed for people at H1, H2, H4 and H5 while only respirators are needed in H3 since people may only be exposed to toxic hazards. Consequently, this study can facilitate the allocation of PPE for    Besides the discussed issues, in the future, CFD [72] may be used to model the dispersion of toxic gas, considering the influences of wind velocity and obstacles on dispersion, thus improving the accuracy of the proposed method. In terms of Monte Carlo simulation, the application of more advanced computers may be needed for the case with hundreds of installations and multiple ignition sources.

Conclusions
There are many installations for storage, transport, or process of hazardous materials in chemical process plants. Once a release occurs at an installation, hazards as toxic release, VCE and fire can simultaneously or sequentially occur, and the generated hazards can evolve spatially and result in a cascading disaster.
In this study, a dynamic methodology based on dynamic graphs and Monte Carlo simulation is developed to assess the vulnerability of humans and facilities in exposure to multiple hazards while considering the spatial-temporal evolution of the hazards. A case study was used to illustrate the application of the methodology and its capabilities in modeling the occurrence and evolution of time-dependent multi-hazard under uncertainty.
The main achievements of the present study can be summarized as follows: (i) The methodology can effectively model simultaneous and sequential multiple hazards caused by the release of hazardous materials; (ii) only considering one type of hazard in vulnerability assessment may largely underestimate the risk, possibly resulting in ineffective allocation of personal protection equipment (PPE); (iii) humans in different locations may be threatened by different hazards, thus different protection strategies may be formulated for people within and around the chemical plants; (iv) VCE and toxic release may result in more severe consequences than fire as long-delayed ignition can result   in the damage of multiple installations and acute toxicity of people around the release source; (v) the concurrent fires resulting from a VCE may be inevitable due to a rapid escalation rate and limited emergency resources; (vi) hazardous installations are more vulnerable to VCEs, and the safety distances based on fire hazards are not sufficient for VCEs; (vii) people close to the release source are prone to multi-hazards while the deaths outside the hazardous storage areas are mainly caused by acute toxicity and VCEs.