Dynamic vulnerability assessment of process plants with respect to vapor cloud explosions

Abstract Vapor cloud explosion (VCE) accidents in recent years such as the Buncefield accident in 2005 indicate that VCEs in process plants may lead to unpredicted overpressures, resulting in catastrophic disasters. Although a lot of attempts have been done to assess VCEs in process plants, little attention has been paid to the spatial-temporal evolution of VCEs. This study, therefore, aims to develop a dynamic methodology based on discrete dynamic event tree to assess the likelihood of VCEs and the vulnerability of installations. The developed methodology consists of six steps: (i) identification of hazardous installations and potential loss of containment (LOC), (ii) analysis of vapor cloud dispersion, (iii) identification and characterization of ignition sources, (iv) explosion frequency and delayed time assessment using the dynamic event tree, (v) overpressure calculation by the Multi-Energy method and (vi) damage assessment based on probit models. This methodology considers the time dependencies in vapor cloud dispersion and in the uncertainty of delayed ignitions. Application of the methodology to a case study shows that the methodology can reflect the characteristics of large VCEs and avoid underestimating the consequences. Besides, this study indicates that ignition control may be regarded as a delay measure, effective emergency actions are needed for preventing VCEs.


Introduction
In petroleum and chemical industrial plants, fire, explosion and toxic release arising from loss of containment (LOC) are concerned as major hazards [1][2][3]. Fire is most common while explosion may impact a wider area and cause severe consequences, leading to multiple fatalities and extensive damage to property [4]. With respect to the amount and rate of vaporization, large releases often result in vapor cloud explosion (VCE) rather than fires [5]. Abdolhamidzadeh et al. [6] investigated 224 domino accidents that occurred in the process industries and indicated that explosion is the most frequent cause of domino effects (57%). VCE has been responsible for 84% of the domino effects induced by explosions. Several catastrophic accidents occurred in recent years due to VCEs, such as the Puerto Rico explosion (2009, USA), the Sitapura explosion (2009, India), and the Amuay explosion (Venezuela, 2012). The Amuay disaster caused by a large VCE at the Amuay refinery, situated in northwestern Venezuela, led to over 50 fatalities and more than 100 injures, damaging 1600 houses and resulted in financial losses up to $1 billion [7,8].
Although VCEs have frequently occurred in the petroleum and chemical industries, the mechanism of the blast is not well understood [9]. For example, the VCE in the Buncefield depot, in the UK, brought about an unexpected overpressure with the maximum value of more than 2000 kPa [10]. A release of hazardous substances can induce a fire if the released substance is immediately ignited while it is more likely to result in a VCE when the ignition is delayed. The subsequent fast expansion of flames produces the overpressure or so-called shock wave, resulting in damaging effects [11]. Many factors influence the evolution and the intensity of a VCE, including the type and quantity of the released flammable substance, the delayed time to ignition (DTI), the space figuration of the release position, the position and the number of ignition sources in the affected area, etc. [12].
Many attempts have been conducted to model the vapor cloud dispersion or estimate the overpressure created by VCEs. The TNT equivalent method [13] is the most widely used method in risk analysis [14,15]. This method provides a simple method for estimating a farfield blast effect, neglecting the space configuration where the explosion takes place, ignition sources and the dispersion of the vapor cloud, thus usually underestimating the overpressure [16,17]. The Multi-Energy method based on the gas explosion mechanism, which considers the VCE as a number of sub-explosions inside special obstructed areas, is recommended as an alternative method for the TNT equivalent method [11]. This method is more suitable for estimating the significant overpressure produced by a large VCE in fuel storage plants [7,10,18]. Other widely used overpressure evaluation methods for VCEs include the Baker-Strehlow method [19] and CFD simulation [20,21].
Compared with fire scenarios, the VCE phenomenon is more difficult to assess due to the uncertainty of ignition position, the uncertainty of delayed ignition time (DIT) and the complexity of overpressure intensity calculation. In fact, the VCE induced by the LOC of hazardous substances in chemical plants is a dynamic process along with the vapor cloud dispersion. However, previous risk analysis methods [22,23] for VCE always assume that the explosion takes place immediately at the release place, which is inconsistent with the observations from large VCEs in the recent years. The position of a vapor cloud explosion depends on the vapor cloud dispersion and ignition sources that can be inside or outside chemical plants. On November 28, 2018, a vapor cloud explosion outside a chemical plant at Zhangjiakou (China) was for instance caused by a Vinyl chloride release inside the plant, leading to 23 fatalities and 22 injuries [24].
The present study aims to establish a dynamic risk assessment methodology based on a discrete dynamic event tree (DDET) to integrate plant physical models and ignition sources into a stochastic simulation engine so as to model the timing dependencies and ignition uncertainty in the evolution of VCEs. The overpressure induced by VCE is calculated by the Multi-Energy method while the damage probability of installations is calculated using probit models. VCE, its characteristics, and previous studies are represented in Section 2. Section 3 illustrates the dynamic accident evolution methodology. A case study is presented in Section 4 and a discussion based on the results is present in Section 5. Finally, conclusions are drawn in Section 6.

Explosion mechanism
A flammable vapor cloud (FVC) is formed by mixing released flammable gases or evaporated flammable liquids and air during the leakage of flammable substances. Flash fire (FF) and VCE are the possible consequences of vapor cloud ignition. The ignition may take place if the concentration of flammable gases lies within the flammability limits (between the lower flammability limit and the upper flammability limit) and an ignition source is present for supplying the required energy (usually of the order of 10 J). FFs can result from the sudden ignition of a FVC, where the flame is not accelerated due to insufficient obstacles or the influence of turbulent dispersion. Alternatively, the flame speed may accelerate to sufficiently high velocities and produce significant blast overpressure [12].
The expansion mechanism of a VCE can be analyzed according to the flame speed which is proportional to the developed overpressure. Following an ignition, the flame starts to propagate away from the point of ignition, with a speed of 5-30 m/s, producing very low overpressure. Next, a wrinkled-flame front appears due to the unstable nature of the flame and large turbulent eddies, resulting in an increase of the flame surface and thus an acceleration of flame speed (30-500 m/s.) and forming an overpressure of up to 2-3 mbar. The presence of obstacles in the flow induces a further increase in the flame speed (500-1,000 m/s), leading to an overpressure of up to 1 bar. This physical process of flame speed acceleration is a deflagration. If the flame speed continues to increase, and the reactive mixture in the front zone of turbulent combustion is compressed and heated due to mixing with combustion products, a shock wave can be created when the reactive mixture's temperature is higher than the self-ignition. This physical effect is called detonation, resulting in a flame speed up to 2,200 m/s and overpressures up to 20 bar. Johnson and Tam [25]  Conditional probability of the ignition of ignition 1 at t 2 given no ignition before time t 1 P I I ( | ) IS t 1 1 Conditional probability of no ignition before time t 1 given no immediate ignition at time t 0 P I LOC Conditional probability of no immediate ignition at time t 0 given a LOC event P IS1 (t 2 ) Ignition probability of source 1 before time t 2 P IS1 (t 1 ) Ignition probability of source 1 before time t 1 P IS2 (t 2 ) Ignition probability of source 2 before time t 2 P IS2 (t 1 ) Ignition probability of source 2 before time t 1

Impact assessment of vapor cloud explosions
The available models to simulate or predict the effects of vapor cloud explosions can be categorized as empirical analytic models and numerical models. The empirical analytic models include Congestion Assessment method, TNT Equivalent method, Multi-Energy method, Baker-Strehlow method, etc. [14] while numerical models are mainly based on CFD codes, such as Flacs and Fluent [21,27].
(1) Multi-Energy method The Multi-Energy method is based on gas explosion mechanism that regards the VCE as a number of sub-explosions inside special obstructed areas [11]. The layout of the space where the cloud is spreading is characterized as a strength coefficient in this method. The value of the coefficient is proportional to the blast overpressure which increases with augmenting the obstacle density in the area. In general, the TNT equivalent method can be used to quickly calculate the overpressure as a function of the distance while the results of Multi-Energy method is usually more accurate and closer to actual conditions [12].
(1) Baker-Strehlow method The Baker-Strehlow method [19], based on the Multi-Energy method, takes into account the flame propagation speed. The flame propagation speed depends on the way the flame front propagates, the reactivity of the fuel, and the density of the obstacles [14].
(1) Numerical models Numerical methods based on CFD codes have received much attention in recent years [20,21,28]. These methods are able to model the effects of terrain shape and the presence of obstacles on the dispersion of a vapor cloud [29]. The CFD simulation of the whole industrial facility can lead to more accurate results, but it is very complex, timeconsuming and expensive, thus may be unsuitable for risk analysis of large and complex process plants [12]. The accuracy of CFD codes especially in the case of VCE simulation in congested environments depends upon the adopted combustion models, turbulence closure models and the constants for the computation of turbulence and for the description of the complex interaction between fame front and turbulent flow field [28].

Frequency assessment of vapor cloud explosions
Event tree has been widely employed to analyze scenarios and frequencies of accidents induced by a LOC event in the process and chemical industries [23,[30][31][32][33]]. Fig. 1 shows a general event tree analysis for a release of hazardous liquid substances from an atmospheric storage tank.
Following a LOC event, a pool fire scenario can occur if the released substance is ignited immediately. Otherwise, the released substance would vaporize and form a vapor cloud. In case of a delayed ignition, A FVC can induce a VCE or FF during the dispersion process. If there is no immediate ignition and delayed ignition, the release event may form a large hazardous vapor cloud that may be harmful to surrounding people or damage the environment.
Assessment of fire accidents triggered by immediate ignition based on the general event tree analysis is reasonable since there is no delay, and the fire can be regarded to occur at the release position. However, it ignores the uncertainty of delayed ignition time (DIT) and the uncertainty of ignition and the uncertainty of delayed ignition position that is essential for the assessment of VCEs. First, the size of the vapor cloud increases over time, which has a great impact on the ignition likelihood and the explosion intensity. Conversely, the ignition position also influences the explosion intensity and the damage effects on other installations. Table 1 lists the DIT values of several large VCE accidents occurred in the process and chemical industry [7,[34][35][36][37][38][39].
As shown in Table 1, the DIT values range from 20 s to 4500 s and ignition source areas can be inside (e.g., pump house, wastewater treatment areas) or outside (e.g., vehicles) chemical plants. Therefore, neglecting the uncertainties caused by vapor cloud dispersion may result in significant errors. Besides, the most of these accidents occurred in no-wind or low-wind conditions, which indicates that stable and large vapor clouds are more likely to form in these weather conditions [26,34].

Methodology
Although a lot of work has been done on the vulnerability of installations subject to vapor cloud explosion caused by LOC, the spatialtemporal evolution of such accidents has been overlooked. However, the vulnerability of installations depends on the intensity of overpressure caused by the VCE which is relevant to the spatial-temporal evolution of vapor clouds before ignition. This section thus aims to establish a dynamic VCE evolution assessment (DVEA) methodology,  integrating the dispersion process of vapor cloud and ignition uncertainty into a stochastic simulation engine in order to assess the vapor cloud explosion risk in process industrial areas. The flow chart of the DVEA methodology is shown in Fig. 2. The subsequent steps of the methodology are more thoroughly explained in the following subsections.

Step 1: identification of hazardous installations and characterization of LOC scenarios
In process industrial areas, large quantities of hazardous (flammable/ explosive/ toxic) substances are handled, transported and stored via all kinds of installations, such as process vessels, pipelines, valves, flanges, heat exchangers, pumps, storage tanks, etc. The inherent hazard of an installation depends on the quantity of substance present, the hazardous properties of the substance as well as the specific operation conditions (e.g., temperature, pressure) [40]. Since this study mainly focuses on fire and explosion accidents, the toxic effects of hazardous substances are ignored. Therefore, only the hazardous installations that may become a release source of flammable or explosive substances are identified in the first step of the developed methodology.
Following the identification of hazardous installations, the loss of containment (LOC) events should be characterized to obtain the LOC scenarios and the corresponding frequency of each scenario. Both LOCs caused by unintentional events and intentional events should be identified in this step. The former should include generic LOCs, externalimpact LOCs, loading and unloading LOCs and others [11]. Generic LOCs cover all failure causes not considered explicitly, such as corrosion, construction errors, welding failures and blocking of tank vents. External-impact LOCs are tailored for transport units. Loading and unloading LOCs are those that occur during loading and unloading operations [41], such as overfilling. To estimate the frequency of LOCs, some specific information may be employed. In terms of intentional LOCs (e.g., deliberately opening valves [42]), security risk analysis according to available information should be conducted, including threat analysis, attractiveness analysis, vulnerability analysis, and consequence analysis, etc. [43][44][45]. Based on the analysis of LOC scenarios, the parameters used for vapor cloud analysis such as initial pressure and temperature of hazardous substances in facilities, the mass of hazardous substances, and the possible leak sizes should be characterized.

Step 2: analysis of vapor cloud dispersion
The results of hazardous installation identification and LOC scenario characterization such as release position, maximum release time and mass flow rate, are the prerequisite for the formation and dispersion analysis of a vapor cloud. This step serves, therefore, to model the vapor cloud dispersion process over time, achieving the vapor cloud volume and position over time. The total release mass (M t ) at time t (the initial release time is zero) for a time-varying release scenario can be expressed as the integral of the mass flow rate (m t ) with respect to time: where m t is the mass flow rate which can be represented as a function of time t. In order to simplify the calculation, the time period t can be divided into n discrete segments, then m t can approximately be expressed as the sum of the masses in each segment, as shown in Eq (2). If the mass flow rate can be regarded as a constant (m) independent of release time, Eq. (2) can also be simplified as Eq. (3) The mass flow rate m of a leakage from a hole can be obtained by using Eq. (4) as follows [11]: where C d represents the discharge coefficient; Α h (m 2 ) denotes the crosssectional area of the leakage hole; P o (Pa) denotes the initial gas pressure in the vessel (for each time step); W g (kg/mol) represents the molecular weight of the gas; γ denotes the Poisson ratio; R represents the universal gas constant (8.314 Jmol −1 K −1 ); T (K) denotes the temperature of the gas; Pa (Pa) represents the ambient pressure,; ρ (kg/ m 3 ) represents the density of the liquid,. Other methods for the calculation of leak rate have also been developed [12,17]. The mass of the flammable substance in the vapor cloud M f,t is represented as the released mass multiplied by the ratio of the evaporation rate to flow rate α which is equal to 1 if the released substance is a gas. In that case, the volume of the flammable gas (V f, t ) is represented as: where ρ f is the density of the flammable gas. A vapor cloud is deemed as a mixture of air and flammable gas. As a result, the total volume of the vapor cloud (V t ) can be obtained as: where V a,t is the volume of air mixed in the vapor cloud. V a,t can be determined by considering that the flammable gas is fully mixed with oxygen [12].
In the Multi-Energy method, the vapor cloud shape is modeled as a hemisphere [13]. Fig. 3 shows a sketch of the vapor cloud hemisphere model and the possible hazardous installations, ignition sources, and obstacles covered by the vapor cloud. In that case, the boundary of the vapor cloud can be characterized by the radius of the hemisphere (R t ), as shown in Eq. (6).
In terms of dense gas, a cylindrical shape [11,46,47] can be used to model the dispersion process, as shown in Fig. 3b. The dispersion is characterized by a consistent height h and a radius of R, as follows: However, the height of vapor cloud varies with the cloud location and evolves over time. Thus a simplified gravity-driven model developed by Atkinson [48] is adopted to model the dispersion of dense gas, as shown in Eq. (8).
C E is an empirical constant and varies from 0.91 to 1.15. Q is the vapor cloud flow rate, g' is the relative density of the vapor. This method is developed based on experiments and simulations results in low-wind conditions in which most large VCE accidents occurred; it is therefore more conservative to model the vapor cloud dispersion when wind velocity is very low. (From a risk assessment perspective, it is more conservative as it would result in more devastating shock waves). Besides the analytical methods based on the hemisphere or cylinder assumption, the available software for dispersion modeling includes ALOHA [49], PHAST [50], EFFECTS [51], etc. CFD software such as FLUENT [21], CFX [52], and FLACS [27] may obtain more accurate results by addressing meteorological parameters, salient relief, and plant layout. In this study, analytical methods are adopted to model vapor cloud dispersion since a large number of release scenarios may involve risk assessment, overcoming the time-consuming aspect of CFD software.

Step 3: identification and characterization of ignition sources
Identification and characterization of ignition sources is a critical step for the dynamic accident evolution assessment given a LOC event. The first task of this step is to identify ignition sources that may contribute to immediate ignition or delayed ignition, such as flare, boiler, and vehicles. In the chemical and process industries, measures for eliminating possible ignition sources are regarded as a significant and practical way to reduce the risk of fire and explosion; such measures may include decreasing the flow rate during loading and unloading operations for preventing static electricity and ground rods for preventing lightning. However, it is impossible to eliminate all ignition sources in an industrial environment. Besides, the ignition sources outside the industrial area have been responsible for some large vapor cloud explosion accidents that occurred in the chemical and petrochemical industries [7,24,37]. Therefore, this step should identify as many as possible ignition sources within the chemical plant as well as outside the chemical plant. In other words, the possible ignition sources within the maximum area of the vapor cloud should be identified, no matter whether they are (in the chemical plant or not). The maximum vapor cloud can be determined by the dispersion model recommended in Step 2, under the premise that the hazardous substance inside the installation is completely released.
As shown in the event tree in Fig. 1, ignition can be divided into immediate ignition and delayed ignition. The immediate ignition is defined as ignition at or near the release source and occurring quickly enough to preclude the formation of an appreciable vapor cloud [32]. Thus, the immediate ignition depends upon both the likelihood of autoignition and the likelihood of static discharge, which can be deemed as irrelevant to release time and vapor cloud dispersion. The probability of autoignition (P Aut ) and the probability of static discharge (P Sta ) can be determined as follows [32]: . This method assumes that there is always a likelihood of non-immediate ignition if T is no more than 200 °C and higher than the AIT. These correlations were developed based on a combination of ignition data and expert judgments. Therefore users should use these correlations with discretion, select conservative values of input parameters, and should not read more accuracy into their predictions than is warranted [32]. Certainly, other correlations and data can be easily included in the developed methodology according to different users and applications.
A delayed ignition can be defined as any ignition other than immediate ignition, where there is a delayed time that allows the formation and dispersion of a vapor cloud. More ignition sources may involve in the accident evolution due to the vapor cloud dispersion. The cumulative probability of ignition caused by an identified source (s) can be modeled as a function of time when the ignition source is present in the vapor cloud (t IS ) and the ignition effectiveness (ω), as shown in Eq. (9) [11].
The ignition effectiveness ω (s −1 ) depends on a lot of factors, such as ignition source types, ignition energy, ignition control measures, etc.  [53,54]. The estimation of ω is a key step in the assessment of delayed ignition probability of a flammable vapor cloud. In this study, we adopt the ignition probability estimation method developed by HSE [53], considering the type of hazardous area, properties of on-site ignition sources (strength, frequency and duration of activity, and density), and ignition control measures in place. Since the ignition probability is expressed by an exponential time-independent function, the probability of ignition is equal to the probability of ignition in one minute which in turn can be used to calculate the ignition effectiveness. It should be marked that this equation can only be used when the ignition source is active and covered by the vapor cloud. Therefore, the cumulative probability should be equal to zero before the ignition source is active or before the vapor cloud arrives. In terms of the ignition caused by vehicles on a road or railway near the plant, the ignition probability can be determined by the average traffic density d. The average traffic density d is defined as: N is the number of vehicles per hour, L is the length of a road or railway section, v is the average velocity of the vehicle. Therefore, the ignition probability caused by vehicles on a road or a railway can be calculated using Eq. (15)

Step 4: explosion frequency and delayed time assessment
This step aims to assess the explosion frequency and the temporal dependencies caused by a LOC event in process plants using a dynamic probabilistic tool. Different safety analysis tools for different research domain are available in the literature [76][77][78][79][80][81]. The widely used dynamic probability tools include dynamic event tree [55], dynamic fault tree [56], dynamic bow-tie [57], dynamic Bayesian network [58] and Monte Carlo simulation [59]. Siu [60] classified dynamic risk assessment methods into three categories: digraph-based methods (e.g., dynamic event tree), explicit state-transition methods (e.g., explicit Markov chain models) and implicit state-transition approaches (e.g. discrete event simulation). Dynamic event tree is recommended as a typical digraph-based tool for modeling system evolution while considering its stochastic behavior and possible dependencies among failure events [60,61]. In order to directly present the accident evolution process and the possible accident scenarios, a discrete dynamic event tree (DDET) is employed in the present study. The DDET is used as a framework to simulate and analyze the dynamic interactions among the vapor cloud dispersion and ignition sources, as shown in Fig. 4. Fig. 4 shows a VCE evolution process with four ignition sources: autoignition, static discharge, ignition source 1, and ignition source 2. Taking the scenarios marked in bold as an example, the probability of VCE/FF caused by the ignition of source 1 at t 2 , can be obtained:

P P P I LOC P I I P I I P VCE I
= ( ) ( )

P P P I LOC P I I P I I P FF I
where P t IS VCE , 1, 2 is the probability of VCE caused by ignition source 1 at t 2 , P t IS FF , 1, 2 is the probability of FF caused by ignition source 1 at t 2 , is the conditional probability of VCE given delayed ignition caused by ignition source 1 at t 2 , is the conditional probability of FF given delayed ignition caused by ignition source 1 at is the conditional probability of the delayed ignition caused by ignition source 1 at t 2 given no ignition before time t 1 , P I I ( | ) t t 1 0 is the conditional probability of no ignition before time t 1 given no immediate ignition at time t 0 , P I LOC ( | ) t0 is the probability of no immediate ignition at time t 0 given a LOC event, and P LOC is the probability of the LOC event. According to Eq. (11), P I I ( | ) t IS t , 1 2 1 can be calculated as: where P IS1 (t 2 ) is the ignition probability of source 1 before time t 2 , P IS1 (t 1 ) is the ignition probability of source 1 before time t 1 , P IS2 (t 2 ) is the ignition probability of source 1 before time t 2 , P IS2 (t 1 ) is the ignition probability of source 2 before time t 1 .
To simplify the calculation, a constant time step (Δt = t i+1 -t i ) is recommended in dynamic event tree analysis. The value of Δt should be determined based on the required calculation accuracy and the needed calculation time. The event tree will end when the probability of the Fig. 4. A discrete dynamic event tree for accident evolution assessment of a LOC event. vapor cloud is less than a threshold, which means that the accident scenarios with a probability lower than the threshold can be ignored.

Step 5: overpressure calculation
In this section, the Multi-Energy method is introduced to calculate the overpressure of VCE scenarios identified in Step 4. To apply the method, two parameters should be determined: (i) strength coefficient and (ii) scaled distance. The coefficient of the strength which characterizes the strength of the explosion blast depends on the obstacle density of the explosion area. The coefficient ranges from 1 to 10 and increases with the increase of obstacle density. The obstacle density is used to characterize the congestion level of the area covered by a vapor cloud. Low obstacle density is defined for areas in which there are few obstacles in the flame path, or the obstacles are widely spaced and there are only one or two layers of obstacles. High obstacle density areas have three or more closely spaced obstacle layers with a blockage ratio of 40% or more [17]. Since obstacle density is the most difficult to quantify in the Multi-Energy method, uncertainty exists in the determination of the strength coefficient. As a result, the coefficient value may be determined by using expert judgment. Appendix A describes a method adopted in this study for estimating the strength coefficient (SC). Since it may be difficult for users to determine the value of strength coefficient, a conservative value is recommended by TNO (The Netherlands Organization for Applied Scientific Research) to avoid underestimating the blast strength [17]. The scaled distance (r sc ) is calculated by Eqs. (18) and (19).
where E (J) is the total combustion energy; P a (Pa) is the ambient pressure; r (m) is the distance from the center of the explosion; ΔH (J/ kg) is the combustion heat of the flammable gas. The scaled overpressure (P sc ), as a function of the scaled distance and the strength coefficient of the explosion blast, can be read from a blast chart [62], as shown in Appendix B. As a result, the overpressure can be obtained as: It should be marked that the uncertainty exists in each commonly used calculation method for vapor cloud explosion. For example, the main uncertainty parameter in the Equivalent TNT Mass method is the fraction of energy released as shock wave (coefficient f E ), while in the Multi-Energy method the unknown parameter is the coefficient of strength of the explosion blast.

Step 6: damage assessment
In order to address the uncertainty of domino escalation and support for vulnerability assessment of installations subject to domino effects, probability models were used to assess the vulnerability of installations [82][83][84]. Bagster and Pitblado [63] proposed a probability approach defining a damage probability function based on the distance from the center of primary scenarios and the safety distance. Khan and Abbasi [64] adopted a probit function to model the damage probability caused by overpressure, considering peak overpressure (static pressure) and dynamic pressure. The probit function was firstly developed by Eisenberg et al. [65] and only peak overpressure was considered in the literature, as shown in Eq. (21).
where △P (Pa) is the peak overpressure; Y is the probit value; a and b are constants. Then the damage probability P r can be calculated using the cumulative standard normal distribution (Φ), as shown in Eq. (22).
Cozzani and Salzano [15] developed probit models for each category of equipment (atmospheric, pressurized, elongated and small) rather than using a general model for all equipment. The equipmentspecific models significantly reduced the error caused by the general probit model, presenting the important difference between the damage probabilities and the damage threshold of different categories of equipment. Therefore, this study adopts these special probit models to estimate the damage probability of installations.

Case study
The large VCE accident in the Buncefield oil storage and transfer depot, 4.8 km from the town center of Hemel Hempstead, Hertfordshire, on 11 December 2005 [66], is used as a case study to illustrate the developed methodology.

Description of the plant and the VCE accident
The layout of the Buncefield oil depot is shown in Fig. 5. It typically stores 150,000 tons of fuel (gasoline, fuel oil, kerosene) with a total capacity of 273,000 m 3 . The main 35 storage tanks are numbered and shown in Fig. 5.
Before the explosion, overfilling occurred during a delivery operation of unleaded petrol via a pipeline to Tank 2, forming a large vapor cloud that covered part of the plant. The overfilling lasted about 23 minutes before the vapor cloud was ignited, resulting in a powerful VCE and the following fires. The accident damaged 23 storage tanks and injured 43 people [29,66,67]. Since the primary VCE event escalated and resulted in overall consequences more severe than the primary event, a domino effect was involved in the Buncefield accident. This study mainly focuses on the assessment of the VCE accident and thus ignores any second-level or higher-order escalation of domino effects [68,69].

Methodology application
To apply the methodology illustrated in Section 3 to the Buncefield plant, we firstly should identify hazardous installations and characterize the possible LOC scenarios. The main hazardous installations include 35 ground storage tanks, pipelines linking with these installations, loading and unloading facilities, and other components such as valves and pumps. The possible LOC scenarios may be releases from storage tanks, tank overfilling during loading and unloading operations, and leakages from pipelines. To illustrate and validate the proposed methodology, only the overfilling scenario is considered, i.e., the excessive liquid flowed down from the vents in the fixed tank roof. The mass flow rate is estimated as a constant of 115 kg/s [47]. The probability of valve failure (i.e., the leakage) is considered to be 5 × 10 −2 per year referred to the failure frequency estimation during loading and unloading operations in the process industry [70].
According to the characteristics of the LOC scenario, a vapor cloud analysis in step 2 can be conducted to obtain the vapor cloud dispersion over time. The ratio of evaporation rate to flow rate, α, is approximately equal to 0.17 given an evaporation rate of 19.5 kg/s. Consequently, the vapor addition rate to the cloud is estimated as 199 m 3 /s based on empirical formulas developed by Atkinson and Coldrick [47]. Since the gasoline vapor is denser than air, the gravity-driven model is used to analyze vapor cloud evolution, considering C E = 1 and g' = 0.5 [48]. As a result, the vapor cloud radius R can be calculated according to Eq. (7) given a delayed time t. Fig. 6 shows the vapor cloud contour at each time slice. The vapor cloud contour is idealized neglecting the effects of site topology and obstacles, etc. It should be noted that by considering these parameters via more advanced CFD methods one may obtain more accurate contours which may be irregular in shape, offset in one direction, or of different thickness [29].
Since the ambient temperature was very low during the accident, the autoignition probability P ia is considered to be zero. The ignition probability of static discharge P is is 0.0156 given a minimum ignition energy of 0.23 mJ. According to the accident investigation [71,72], there are two possible ignition sources in the oil storage depot, i.e., the pump house (IS1) and the Northgate emergency generator (IS2). Fig. 7 shows the dynamic event tree up until 25 min with 66 possible VCE scenarios.
For illustration, other possible ignition sources such as the vehicles parked nearby are not considered in this case study. Both the primary ignition probabilities of the two sources are considered to be 0.1/min (with 'good' ignition controls [53]) since the two equipment items were not in operating condition during the accident. As a result, the parameters of ω in Eq. 10 of the two ignition sources is equal to 0.0018. IS1 is active after t = 1.5 min while IS2 is active at t = 5 min. The next step is to estimate the explosion probability and possible delayed time using a dynamic event tree for a time step Δt = 1 min.
Based on the event tree, the cumulative probability of ignition over time can be obtained, as shown in Fig. 8. The cumulative probability increases over time while the increase rate decreases after the vapor cloud reaches the second ignition source. At t = 1 min, the delayed ignition probability is zero since the vapor cloud hasn't covered any ignition source. The cumulative probability reaches 0.98 at 23 min and still increases over time, finally approaching 1.
There are four possible ignition causes: immediate ignition caused by static discharge, delayed ignition at the pump house (SI1), delayed ignition at the Northgate emergency generator (SI2) or simultaneous ignitions at the pump house (SI1) and the Northgate emergency generator (SI2). Fig. 9 shows the ignition probabilities of different ignition causes over time. The maximum ignition probability is in the 6th min when the vapor cloud reaches Northgate emergency generator and the two ignition sources are active. The ignition probability in 2 nd min is lower than that in the 3rd min because the vapor cloud reaches the pump house at t = 1.5 min (i.e., the ignition source active from 1.5 min). The ignition probability decreases rapidly over time as the cumulative ignition probability increases, so a VCE with a long DIT may be regarded as a low-probability event. Based on the spatial-temporal dispersion results, Step 5 can obtain the overpressure caused by VCE at each discrete time using the Multi-Energy method. Since the industrial area was blocked by various buildings, tanks, and plants, a conservative strength coefficient of 10 (the maximum value) is considered for all the VCE scenarios. Fig. 10 shows the overpressure caused by VCE scenarios at different distances.
The overpressure increases over the delayed ignition time (DIT) due to the increase of total explosion energy. The maximum overpressure can reach 2000 kPa which is consistent with the previous study on Buncefield accident investigations [10,18]. But the overpressure rapidly decreases with increasing the distance from the center of the explosion. So, equipment with a large distance from the center of the explosion may survive from the VCEs, such as T34 and T35. The final step assesses the vulnerability of hazardous installations, obtaining the damage probability of installations and the likelihood of domino effects caused by possible VCE scenarios. The parameter values of a and b are considered as −9.36 and 1.43, respectively, for atmospheric tanks [73,82]. Consequently, the damage probability of installations subject to these VCEs can be calculated using Eqs. (21) and (22) . Fig. 11 shows the conditional probability of damage to T20 and T35. The damage probability of T35 is lower than that of T20 for a VCE since the distance from the explosion to T35 is larger than that to T20. The conditional probability of damage for each tank increases with the increase of delayed ignition time (DIT). The explosion caused by simultaneous ignitions at SI1 and SI2 leads to more severe consequences than the explosion caused by single ignition at SI1 or SI2. Fig. 12 shows the damage probability of tanks caused by the VCE scenario at t = 23 min (the Buncefield explosion accident in 2005). The results indicate that Tanks 1-21 are very likely to be destroyed by the explosion due to the high damage probabilities (> 0.9). Fig. 13 shows the layout of the Buncefield plant after the accident in 2005. The real damaged tanks can easily be identified and these 21 tanks (marked by yellow circles) are very likely to have been destroyed by the explosion due to the high damage probabilities (> 0.9). Among the 21 tanks, only T5 and T9 were not really damaged in the accident, which indicates that the results obtained by the developed methodology are almost in agreement with the real Buncefield accident.
Finally, we can obtain the conditional probability of damage of installations subject to possible VCEs given an overfilling scenario at Tank 2, as shown in Fig. 14.
Tanks 25-35 have a lower damage probability than other tanks since they are situated at a substantial distance from the release tanks. The number of expected damaged tanks (the sum of damage probability of each tank) is 16. In that case, domino effects may be inevitable. The results indeed indicate that a large vapor cloud explosion can lead to the damage of multiple tanks or knock-on effects, resulting in severe consequences. This was the case for instance with the Buncefield VCE accident in 2005, the San Juan VCE accident in 2009, and the Jaipur VCE accident in 2009.

Discussion
The case study indicates that the developed dynamic assessment methodology can model the influence of the spatial-temporal evolution of a vapor cloud and the uncertainty of delayed ignition on the vulnerability of installations subject to the major accident scenario "vapor cloud explosion" (VCE). The obtained results are consistent with the accident investigations for the past large VCEs. This section analyzes the critical parameters in the methodology and the possible future   research issues to further improve this study. In a real chemical plant, multiple ignition sources may be present within or outside the plant. As a result, collecting and characterizing the main ignition sources is a critical step in the developed methodology. The ignition probability of a single source depends on a lot of factors, such as the type of the source, the ignition energy of the source, and also the control measures for the source. Fig. 15 shows the cumulative probability of ignition over time with different values of ignition effectiveness (ω).
As shown in Fig. 15a, the ignition probability increases with augmenting the ignition effectiveness of single sources. Therefore, ignition control is widely used to decrease the ignition effectiveness so as to prevent major accidents in the process and chemical industries. However, ignition cannot be completely eliminated by ignition control measures. Besides, due to ignition control measures, a VCE may be delayed, resulting in a larger vapor cloud and thus a larger VCE of more severe consequences. As shown in Fig. 15b, the conditional damage probability of the tanks decreases with increasing the distance between the tanks and the VCE center (the dips in the figure are due to situations where the higher numbered tanks are actually further away from the ignition source). Moreover, the conditional damage probability of the tanks increases with decreasing the ignition effectiveness. Therefore, ignition control can be considered as a delay measure, which may not be adequate to prevent VCEs. To prevent VCEs, ignition control measures may be integrated with emergency response actions such as diluting oil vapor by water vapor. In other words, ignition control may be used to provide enough time for emergency response actions to prevent VCEs.
The Multi-Energy method was adopted to calculate the overpressure caused by VCEs in this study. The key issue in the application of this method is to determine the strength coefficient based on the congestion (obstacle density) of process plants. Obstacle density is the most difficult parameter to quantify in the application of the Multi-Energy method. Although TNO has already published the yellow book to guide the application of the Multi-Energy method, it is still difficult for users to determine the value of the strength coefficient due to the uncertainty of obstacle density. Since the actual overpressure can easily be underestimated according to experimental results, it is recommended to be conservative in the determination of the strength coefficient [17]. Taking the case in Section 4 as an example, the maximum conditional damage probability subject to VCEs caused by overfilling at T2 decreases to 0.35 if the strength coefficient of 10 is substituted by 3. Thus, the strength coefficient should be determined by meticulously analyzing the layout of a process plant in the application of the developed methodology. Otherwise, the worst strength coefficient should be adopted in order to obtain conservative results in vulnerability assessment.
To make the developed methodology user-friendly, we adopted an analytic method to predict the dispersion of vapor cloud, neglecting the VCE dilution with distance. As a result, the application of this methodology would lead to more conservative results in risk assessment. To account for VCE dilution with distance as well as upper and lower  explosion limits, CFD methods may be integrated into this methodology to obtain more accurate results in future studies. With the rapid improvement of computational resources, applying CFD methods in risk assessment may become easier and acceptable for engineers in the future. Monte Carlo simulation can also be integrated into this methodology to evaluate the frequency of each scenario when the number of possible accident scenarios becomes too large due to the increase in the number of ignition sources.
The developed methodology in this study was illustrated and verified by the VCE at the Buncefield oil storage facility as a case study. The results agree with the observations that more than 20 tanks were damaged by the VCE. Besides the real accident scenario, other possible scenarios were also obtained by the application of the developed methodology which shows the effectiveness of the methodology in considering the uncertainties (more than one accident scenarios could have occurred).
In this study, second or higher-level escalation was neglected, which could be considered for future work. Besides the application in risk assessment, the developed method combined with Bayesian theory may be used in accident investigations to identify the most likely ignition source based on evidence-based reasoning.

Conclusions
This study introduced a new methodology based on dynamic event tree (DET) to model the vulnerability of process plants to VCEs, considering both the spatial-temporal dispersion of vapor cloud and the uncertainty of delayed ignition time (DIT). This work demonstrated how DET can effectively be used to assess the damage probability of critical installations exposed to VCEs caused by loss of containment. The dynamic methodology can address severe consequences caused by a large VCE due to a long DIT, such as the Buncefield accident in 2005.
Different from previous work on vulnerability assessment for vapor cloud explosions, the key outcomes of the present study can be summarized. Firstly, the time dependencies in vapor cloud dispersion and the uncertainty of delayed ignition should be considered to assess the VCE; this is crucial for reflecting the characteristics of possible large VCEs and for avoiding the underestimation of their consequences. Secondly, the vulnerability of installations to VCEs depends on the congestion of the plant layout and DIT. A long-delayed explosion may result in multiple-failure of installations, resulting in catastrophic disasters. Thirdly, the DIT is related to the distance between the release position and the ignition sources, the type of ignition sources, and the ignition control measures in place. The ignition control measures can decrease the ignition probability of single sources and may delay (if not   completely eliminate) the VCE. However, a delayed ignition (which might be considered a good thing) could actually lead to a larger VCE and more severe consequences. Lastly, combining ignition control measures with emergency response actions (e.g., diluting oil vapor by water vapor) may be a way to prevent VCEs in process plants since ignition control might provide enough time for emergency response actions to prevent VCEs.

Declaration of Competing Interest
The authors declare that they have no known competing for financial interests or personal relationships that could have appeared to influence the work reported in this paper.

A. strength coefficient
There are several methods based on qualitative factors available in the literature [85], such as the index method developed by Kinsella [74]. The method based on three factors: (i) degree of obstruction by obstacles inside the vapor cloud, (ii) ignition energy, (iii) degree of confinement. The first factor is divided into three levels: high (obstacles in the gas cloud with a volume blockage fraction no less than 30% and with spacing between obstacles no more than 3 m), low (obstacles in gas cloud with a blockage fraction less than 30% and/or spacing between obstacles in excess of 3 m) and none (no obstacles within gas cloud). The factor of parallel plane confinement is divided into two levels: confined (gas clouds, or parts of it, are confined by walls/barriers on two or three sides), and unconfined (gas cloud is not confined, other than by the ground). The factor of ignition strength is divided into two levels: high (high energy source), and low (low energy source). The strength coefficient then can be estimated according to