A dynamic stochastic methodology for quantifying HAZMAT storage resilience

A disruption to hazardous (flammable, explosive, and toxic) material (HAZMAT) storage plants may trigger escalation effects, resulting in more severe storage performance losses and making the performance restoration more difficult. The disruption, such as an intentional attack, may be difficult to predict and prevent, thus developing a resilient HAZMAT storage plant may be a practical and effective way to deal with these disruptions. This study develops a dynamic stochastic methodology to quantify the resilience of HAZMAT storage plants. In this methodology, resilience evolution scenarios are modeled as a dynamic process that consists of four stages: disruption, escalation, adaption, and restoration stages. The resistant capability in the disruption stage, mitigation capability in the escalation stage, adaption capability in the adaption stage, and restoration capability in the restoration stage are quantified to obtain the HAZMAT storage resilience. The uncertainties in the disruption stage and the mitigation stage are considered, and the dynamic Monte Carlo method is used to simulate possible resilience scenarios and thus quantify the storage resilience. A case study is used to illustrate the developed methodology, and a discussion based on the case study is provided to find out the critical parameters and resilience measures.


Introduction
In the petroleum and chemical industry, hazardous material (HAZ-MAT) storage infrastructures play a critical role in industrial production by providing various chemical products such as petroleum, natural gas, and acrylonitrile. These chemicals are always stored in atmospheric tanks or pressurized tanks located nearby in an industrial area (e.g., oil depots, LNG terminals, and chemical storage facilities). Most of these chemical products are flammable, explosive, or toxic, making storage facilities vulnerable to disruptions, resulting in major accidents such as fire, explosion, and hazardous release [ [13], [27], [50], [62], [77], [80]]. Besides, these major accident scenarios in a storage installation may escalate to installations nearby, leading to a chain of accidents, resulting in the overall consequences more severe than the primary event, which is called "domino effects" or "escalation effects" [66].
A HAZMAT storage plant may encounter various disruptions. According to the nature of the disruption events, the disruptions may be divided into three categories: unintentional accidents, natural disasters, and intentional (cyber or physical) attacks [ [13], [16], [32], [71]].
Accidental disruptions may be caused by mechanical failure, corrosion, fatigue, and human errors, etc., such as the Intercontinental Terminal Company (ITC) Tank Fire in March 2019 at Deer Park, the U.S [15]. Compared with accidental disruptions, natural hazards may result in more severe consequences due to the damage of multiple chemical facilities, safety barriers, and other emergency response infrastructures [90,91]. The damage to industrial facilities caused by natural disasters is called the Natech event [ [17], [49]]. For instance, the hurricane of Harvey in 2017 led to the release of at least 18 hazardous storage tanks in Texas [ [60], [64]]. Both the accidental disruptions and the disruptions caused by natural hazards are unintentional, while intentional attacks may aim to cause damage to the attack objective by using external weapons besides the hazardous materials inside chemical facilities. For example, on June 26, 2015, two tanks in a France chemical plant were damaged due to an explosion attack [12].
Resilience engineering is becoming a more active and substantial research topic in the safety and security domain. Although no identical definition of resilience exists currently in the academic domain [37], the capabilities (metrics) of a resilient system for responding to unexpected disruptions can be summarized as follows [[14], [36], [82]]: (i) Absorptive capacity: the capability of a system to resist, absorb, or withstand the impact of disruptive events; (ii) Adaptive capacity: the capability of a system to adapt itself to maintain its operational performance without any recovery activity; (ii) Restoration capacity: the capability of a system to repair or restore damages from a disruption to recover the loss performance of the system, making the system reach a new stable state.
Safety management aims to take measures to reduce the likelihood and consequences caused by disruptions for avoiding and mitigating human loss, economic loss, environmental loss, etc. Different from safety management, resilience engineering is to enhance a system's capabilities to absorb, adapt, and recover from a disruption, reduce the impacts of the disruptions on the system's performance. Safety management may be used to enhance absorption capability while has no direct impacts on adaption and recovery capabilities. As a result, safety management is not as wide as that of resilience management/engineering. In light of the unpredictable or indefensible threats (e.g., intentional attacks and natural disasters), enhancing the resilience capability is an ideal approach to reduce the losses caused by disruptions and to quickly recover its storage performance [88].
The term "resilience" originated from the Latin word "resiliere" which refers to the capability to rebound [23]. Therefore, Resilience may be regarded as the ability to return to an equilibrium [26]. Holling [34] highlighted the role of resilience in absorbing changes. Walker et al. [89] defined resilience as the ability to absorb a disruption and reorganize while undergoing change while retaining the same state. Doorn et al. [23] suggested that resilience is a 'formal' concept and demonstrated that resilient performance is the ability to keep or enhance certain functions. In recent years, resilience has become an important concept in the safety and security domain [5]. The advancement of resilience engineering research will contribute significantly to chemical process safety [ [14], [22], [35], [62], [63]]. Past research attempts on resilience in the process industry identified the process resilience influence factors [ [22], [41]], resilience hazards [ [6], [42], [65]]. Besides, the Bayesian network [92] was used to quantify process resilience [ [1], [77], [88]]. However, Little attention has been paid to HAZMAT storage resilience in which escalation effects may play an essential role [ [14], [72]].
This study aims to develop a methodology for quantifying the resilience of hazardous storage systems, considering the dynamic stochastic evolution of disruptions due to escalation effects and the dynamic restoration process. This paper is organized as follows: Section 2 defines HAZMAT storage resilience and introduces the possible measures to enhance the resilience of hazardous storage resilience. a stochastic dynamic methodology for quantifying the storage resilience is elaborated in Section 3. A case study is provided in Section 4 and a discussion based on the case study is illustrated in Section 5. A discussion is provided in Section 6 and the conclusions drawn from this study are present in Section 7.

The definition of HAZMAT storage resilience
The resilience concept has been used in various industries and systems, and many resilience definitions for different domains are available in the academic domain [37]. For instance, the US National Academy of Sciences provides a general definition of resilience: "the ability to anticipate, prepare for, and adapt to changing conditions and withstand, respond to, and recover rapidly from disruptions" [58]. Carpenter et al. [9] stated that the operationalization of resilience depends on the answer to the question 'resilience of what to what?'. Thus, in this paper, we define the HAZMAT storage resilience as the capability of a storage plant to resist, mitigate, adapt and recover from disruptions, to maintain its storage performance. Safety and security measures aim to prevent undesired events and mitigate the consequences caused by the events. Resilience engineering measures intend to enhance a system's capability to anticipate and prepare for disruption and its ability to adapt and recover from the disruption. To improve the resilience of a HAZMAT storage plant, operators should apply measures in different stages to resist the impacts of an undesired event, mitigate the consequences by preventing possible domino effects, and adjust operation strategies to improve the storage performance before recovery and to rapidly recover the plants. Fig. 1 shows the storage performance changing over time. According to the storage resilience definition, a hazardous material storage plant may be in six stages when it comes to undesired events. Before the occurrence of undesired events, the storage plant is in the initial stage, i.e., storage performance is the maximum value S 0 . When a disruption (undesired event) occurs, the performance may decrease immediately and cause major accident scenarios due to the damage to one or more than one storage tank. The primary major accident scenarios may escalate to installations nearby, resulting in escalation effects. This will further reduce the storage performance. When the escalation is prevented (t 2 ), the residual storage performance reaches the minimum value S 2 . At that time, the storage plant may adapt its operation strategies to improve storage performance. For instance, the storage plant can utilize reserve tanks or speed inventory turnover ratio. The last resilience strategy is to recover the performance by repairing or rebuilding the damaged installations (t 3~t4 ). t 4 is the time that the storage performance of the plant is fully recovered rather than the end of the restoration since the final restored performance may exceed the original performance. In real cases, the storage performance of a recovered plant may be different from the plant in the initial stage. In this study, the entire resilience process is time-dependent, which is called a resilience evolution scenario. It is also a stochastic process due to the uncertainties in the vulnerability of tanks, hazardous scenario escalation, emergency response, etc.

Resilience metrics
As shown in Fig. 1, the storage performance varies over time. According to the resilience framework proposed by Bruneau et al. [7], resilience loss can be the expected degradation in performance over time. Based on the storage performance, we defined the storage resilience metrics as a dimensionless ratio, as follows: The numerator of formula (1) indicates the accumulation of storage performance S(t) between t 1 (disruption) and t 4 (fully recovered). The denominator represents the accumulation of initial storage performance S (t 0 ) between t 1 and t 4 . Although the resilience metrics is illustrated by the case of hazardous storage systems, it may be applied to other fields by substituting other performance functions for the storage performance function S(t).
There may be many possible resilience evolution scenarios in terms of the uncertainties in resilience (which can be seen as different performance curves). Considering N resilience evolution scenarios and the maximum value of t 4 is t max , then the resilience metrics can be adapted, as follows: S i is resilience evolution scenario i. It is worth noting that the resilience value R of different resilience evolution scenarios may change with the upper limit (t 4 ) of the integral in In Eq. (1) [ [73], [74]]. As a result, t 4 is substituted with t max in Eq. (2), to unify the time dimension and thus avoid overestimating the resilience with longer resilience evolution time t 4 . In the following calculation and analysis, t max should be consistent. The most resilient system (R = 1) is an ideal condition in which the disruption does not induce any performance degradation. In such case, the impact of the disruption on the system is fully absorbed. If the system is destructed and the recovery is impossible, R is equal to zero. The value of R is between 0 and 1. It should be marked that t 4 is the time the system is fully recovered rather than the end of the restoration. The performance at the end of the restoration may exceed the original storage performance while the maximum of restored storage performance at t 4 cannot exceed its original storage performance and R is no more than 1.

Capabilities of hazardous material storage resilience
According to the storage performance curve in Fig. 1 and the resilience metrics in Eq. (1), the capabilities of hazardous material storage resilience consist of resistance capability, mitigation capability, adaptation capability, and restoration capability, as shown in Fig. 2. Resistance capability is the capability to resist descriptive events to avoid failure and maintain operation. Various measures can enhance resistance capability, and different measures may be taken to tackling different disruptions. For instance, installing lightning masts around storage tanks and installing air terminals on the tank can prevent the damage of tanks caused by lightning strikes [ [61], [80]]; while security measures such as fence may be used to delay attacks and thus prevent intentional attacks [ [12], [84]]. Besides, early warning (detection) measures are essential for resisting intentional and unintentional disruptions since they may be the precondition of effective delayed measures and response actions. The commonly used warning measures in the chemical industry include smoke detectors, toxic gas detectors, closed-circuit television (CCTV) cameras, and intrusion detection  systems, etc. By applying these measures, resistance capability can be improved and thus increase the value of S 1 in Fig. 2, enhancing storage resilience. Fig. 3 Mitigation capability is the capability to prevent the escalation of possible major accident scenarios caused by disruptions. This capability is essential for hazardous storage infrastructures due to possible domino effects. Safety barriers are always used to prevent the escalation of domino effects in the process industries, such as active protection measures (e.g., pressure relief valve), passive protection measures (e.g., fireproof coating), procedural and emergency response (e.g., firefighting) [ [13], [45], [53]]. Since safety barriers can prevent the escalation of domino effects and avoid catastrophic events, they can increase S 2 and shorten the time between t 1 and t 2 , thus reducing the storage performance loss and improve resilience.
Adaptation capability is the capability to adapt a new operation state to fully or partly recover the performance before restoration. This capability can be achieved in hazardous material storage plants by adjusting operation strategies, such as utilizing reserve tanks, speeding inventory turnover ratio, or adjusting chemical inventory strategies. Enhancing the adaptation capability can increase the value of S 3 in Fig. 2 and improve the storage performance before the restoration of the damaged storage plant.
Restoration capability refers to the capability to quickly repair or rebuild the damaged installations to recover the performance of the storage plant. In this stage, the restoration capability mainly depends on the time to fully recover (TTR). Therefore, shortening the TTR (t 3 -t 4 ) can effectively reduce the lost performance and achieve a more resilient storage plant.
Modeling the resilience capabilities based on the performance curve is the key step to quantify the resilience of a hazardous material chemical plant. There, the next section is to develop a framework to quantify the resilience metrics by modeling these resilience capabilities.

Storage performance metrics
Hazardous material storage plants are industrial facilities for storing hazardous chemicals such as petroleum, benzene, and other chemical products. These products are delivered to end-users, process facilities, and other storage facilities. As a result, the average daily chemical flow rate at the initial stage (f in ) of a hazardous material storage plant is used to represent the system storage performance. It is the sum of the daily chemical flow rates of all hazardous materials, as follows: m is the number of hazardous materials in a storage plant and each hazardous material at least has one storage tank. If there are n storage tanks in the plant, the total storage volume (V) is the sum of the volume of all storage tanks, as follows:

Resistance modeling
The resistance capability is the ability to withstand and retain operation and avoid being damaged [82]. Resistance may be viewed as the antonym of vulnerability, representing the inability of an installation to withstand strains and the consequent failures [43]. The vulnerability of installations is often represented by the failure probability of installations exposed to a disruption. Therefore, the resistance capability of an installation can be obtained as follows: where P f is the failure probability of an installation exposed to a disruption, and C r is the resilience capability of the installation exposed to the disruption. In terms of a hazardous material storage plant, there are usually multiple storage installations. Duo to a disruption such as a terrorist attack [12], a natural disaster, a cyber-attack, multiple storage installations may be simultaneously damaged, resulting in the sudden decrease of storage performance (S 0 -S 1 ). S 1 is the storage capability of the undamaged tanks. The resistance capability of an installation mainly depends on different disruptions and the intensity of the disruptions. Therefore, the robustness aspect of resilience [75] is not considered in this study. Let us take an intentional attack using explosive devices as an example, the overpressure caused by the explosion is the main threat for installations, possible to be calculated by the TNT equivalency method [ [55], [78]]. In this method, the point-source TNT explosion model is adopted by transferring the net explosive mass to be an equivalent amount of TNT. Then the scaling distance can be determined as follows: where r s denotes the scaled distance; m TNT represents the equivalent mass of TNT, kg; r represents the distance between installations and the explosion. The overpressure is a function of the scaling distance, which can be read from the TNT blast chart or obtained from empirical formulas.
Since the blast chart needs to be read by humans and thus may not suitable for computer computation, many empirical formulas were developed in the past years. Eq. (6) provides an empirical formula of overpressure (P s ) based on the TNT blast chart [4].
To obtain the failure probability of storage vessels, the probit model [ [18], [19]] is adopted in this study, as follows: where Y p is the probit value. The failure probability (P f ) caused by overpressure is thus calculated using the cumulative distribution function of the standard normal distribution (ϕ), as follows: Fig. 3. Blast chart for calculating overpressure caused by TNT explosion [4].
In each resilience evolution scenario, the damaged tanks in the disruption stage can be determined by sampling random numbers according to P f (see the illustrations in Section 4). Therefore, the total storage volume (V t1 ) can be obtained according to the damaged storage volume in the disruption stage (V di ), as follows: It should be noted that the reflection and diffraction of shock waves due to the TNT explosion around the tanks have not been taken into account in this study.

Mitigation modeling
In hazardous material storage plants, hazardous storage installations are located nearby. Domino effects may occur due to possible major accident scenarios (e.g., fire and explosion) caused by disruptions. If domino effects occur, the consequences may be more severe than the primary disruptions. The mitigation capability in hazardous material storage plants thus refers to preventing or mitigating the escalation of domino effects. As shown in Fig. 2, enhancing the mitigation capability can raise S 2 and may decrease t 2 , to improve storage resilience.
When a storage installation is damaged and results in a loss of containment of hazardous materials, major hazards such as fire and explosion can occur, resulting in the nearby installations exposed to heat radiation or overpressure. Once the physical force damage the nearby installations, the major accident scenarios may propagate, resulting in a chain of accidents and decreased storage performance (t 1 -t 2 in Fig. 2). The damage likelihood caused by overpressure can be estimated by Eqs. (7) and (8). In terms of fire escalation, it is a dynamic process due to the damage of storage installations caused by heat radiation time to failure (TTF) [54] and the time to burn out (TTB) of flammable materials [44]. During the escalation process of domino effects, a storage installation may receive heat radiation from multiple fires, called a "synergistic effect" [11]. In that case, the propagation is triggered by the collaboration of multiple fires and the total received heat radiation can be simplified as the sum of the heat radiations caused by different fires: ([ [11], [44]]; Khakzad et al.) Q jk ,ti is the heat radiation from installation j to installation k at time t i . If installation j is not on fire, it is equal to zero. Q k,t i is the total heat radiation received by installation k at time t i . In the initial stage, the TTF can be calculated based on the total heat radiation, as follows [54]: where V is the volume of the storage equipment; c 1 , c 2 , c 3 , and c 4 are constants depending on equipment types, as shown in Table 1. A synergistic effect can lead to the increase of heat radiation (Q) and thus result in the reduction of TTF. Besides, the received heat radiation may change over time, and the heat radiation received in each period should be superimposed when calculating the TTF, which is called "superimposed effects" [12]. By the application of the dynamic graph approach [11], the TTF at time t i+1 (TTF k,t i+1 ) can be updated according to the TTF at time t i (TTF k,t i ), as follows: When the TTF equals zero, then the installation is considered "failed", and fire escalation occurs. To avoid the failure of installations, safety barriers such as passive barriers, active barriers, and emergency response barriers may be implemented [ [45], [53]]. The passive fire protection measures refer to these safety measures that do not need external activation to trigger the protection functions for containing fire or delaying fire escalation, such as fireproof coating, pressure safety valves, and fire-resistant walls. These protection measures are based on different mechanisms and thus have different performances for fire protection. The fireproof coating is a widely used passive protection measure for storage tanks in the process industry. Usually, the performance of the fireproof coating can be modeled by the fire-resistant time (FRT). In that case, fire escalation can be delayed by using the fireproof coating, providing time for firefighting. In terms of active protection measures, external activations are needed to trigger the protection function, such as the water spray system (WSS). The water spray system will be active by an actuation system (including fire detection equipment, power, pump, etc.) when a fire occurs. The reduced heat radiation due to WSS can be obtained using Eq. (11) [46]: (14) Where Q r is the reduced heat radiation; η is an effectiveness factor of the WSS, φ is the heat radiation reduction factor. The third safety barrier is emergency response. Emergency response actions such as firefighting are essential to prevent domino effects while a period is needed for the emergency response team to arrive. The emergency response system can be regarded as a socio-technical system with some uncertainties, and the performance can be represented by the time for emergency response (TER). As shown below, a lognormal distribution (with a mean value u and a variance σ 2 ) is usually used to model the TER considering its uncertainty.
If the emergency response teams arrive in time, domino effect escalation is assumed to be prevented [ [11], [83]] and t 2 is considered equal to TER. Emergency response capability can be improved and thus the TER may be decreased by having sufficient emergency drills, simulations, and training [ [57], [86]]. By applying one safety barrier or more, domino effects may be prevented and mitigated, reducing the number of possible damaged storage installations and/or the level of damage, and decreasing the performance reduction.
At the end of the escalation stage (t 2 ), the total storage volume (V t2 ) can be obtained according to the damaged storage volume in the escalation stage (V es ), as follows:

Adaptation modeling
The adaptation capability in this study refers to operation adjustments, which can lead to improved storage performance. Operation strategies include utilizing reserve tanks, speeding inventory turnover ratio and other chemical inventory strategies, etc. reserve tanks can directly increase the total storage volume and thus increase the daily chemical flow rate, as follows:

Table 1
The values of c 1 , c 2 , c 3 , and c 4 [54]. where V t3 is the total storage volume at t 3 ; f t3 is the daily chemical flow rate of hazardous materials; V re represents the increased storage volume due to reserve tanks. The daily chemical flow rate can be directly improved by speeding the inventory turnover ratio, as shown in Eq. (15).
where r t3 is the inventory turnover ratio at t 3 ; r t2 is the inventory turnover ratio at t 2 . According to these adaptation strategies. The loss of storage performance caused by the disruption and the sequential cascading effect can be partially recovered.

Restoration modeling
The loss of performance may be fully recovered by restoring the damaged installations. In this study, all the damaged tanks are considered to be reconstructed. In this stage, the time to full recovery (TTR) is a quantitative indicator. The reconstruction of a storage tank is a timeconsuming process. It includes several steps: installing the tank bottom, installing the hydraulic jacking system, installing the tank roof, assembling and lifting the first ring (top) of the tank wall, assembling and lifting the second ring of the tank wall, installing the accessories. The construction time of tanks depends on many factors, such as the tank volume, construction method, the number of people, and resources invested in the construction. Normally, the construction of a storage tank needs several months. If multiple storage tanks are damaged, the rebuilding sequence may also affect the TTR. The restoration capability is negatively correlated with the construction time. As a result, increasing the investments in construction can shorten the TTR and thus enhance storage resilience. Besides, a company may improve the level of preparedness to quickly recover from disruptions, such as the Fig. 4. Flow diagram of the algorithm for obtaining storage resilience.
availability of drawings, construction and maintenance teams, and financing, etc.

Simulation algorithm
This section provides a stochastic dynamic algorithm to obtain the resilience of the HAZMAT storage plant exposed to disruptions. Fig. 4 shows the flow diagram of the algorithm for obtaining storage resilience.
The algorithm based on the dynamic Monte Carlo method is an iterative procedure in which uncertainties in storage resilience are modeled. Iterative algorithms are widely used in dealing with uncertainty and decision-making issues. For example, Hausken [31] proposed an iterative procedure to illustrates how a threat impacts uncertainty in which multiple players have different thresholds for acceptable uncertainty. As shown in Fig. 4, firstly, we need to input the number of iterations N, the disruption time t 1 = 0, and the initial iteration n = 1. Given a disruption, vulnerability analysis will be conducted using the models in Section 3.2 to determine the failure probability of installations exposed to the disruption. Based on the failure probabilities, a set of random data (between 0 and 1) is sampled to determine the damaged storage installations. If a random number is less than the failure probability, the installation is considered damaged. Then, possible escalation is assessed using the escalation models in Section 3.3 to obtain the failure installations in the escalation stage. Again, random data will be generated and used to determine the end time of escalation t 2 . According to the results of vulnerability models and escalation models, the storage performance from t 1 to t 2 can be gained. Next, the improved performance due to adaptation measures needs to be determined. Furthermore, the restoration start time (t 3 ) and end time (t 4 ) should be determined based on the restoration strategy. Based on these calculations, the entire performance curve (t 1 -t 4 ) can be obtained. The above steps need to be repeated until n exceeds N. To calculate the resilience in each iteration, the maximum value of t 4 (t max ) can be found out. Setting the integral interval [0 t max ], the storage resilience in each iteration is calculated based on Eq. (1). Finally, the storage resilience can be obtained according to Eq. (2), considering the dynamic resilience evolution process and uncertainties in the disruption and escalation stages.

Case study descriptions
In this section, the resilience of a refined oil storage farm is examined using the resilience quantification methodology proposed in this study. The storage farm consists of 14 tanks (numbering T1-T14) and stores two hazardous materials: gasoline (T2-T5, T12-T14) and diesel (T1, T6, T7-T11). The characteristics of the tanks are listed in Table 3, and the layout of the storage farm is shown in Fig. 5. The total storage volume in the initial stage (initial) is 30500 m 3 . The flow rate in the initial stage of gasoline and diesel is 1088.6 m 3 /d (3 × 10 8 kg/y) and 658.6 m 3 /d (2 × 10 8 kg/y), respectively. The inventory turnover ratio of the storage farm is 20.9. Assuming that the average value of TER (μ) is 15 min and the variance (σ) is 5 min.
Assume that a disruption due to an intentional explosion occurs at the storage farm (represented by a red asterisk in Fig. 5) induced by a suitcase bomb with an improvised explosive device (IED), the explosion is assumed to be equivalent to 23 kg TNT [38]. The distances between different storage tanks and the explosion position are shown in Table 4.
The explosion may lead to fire on tanks, and the heat radiation generated by pool fire from one tank to other tanks can be obtained by using the ALOHA software [2]. The parameters for heat radiation calculation and the results are shown in Appendix. The mean time for emergency response (TER) of 15 min and the variance of the TER of 5 min are assumed for the storage plant. Once the attack results in tank damage, a suspended time (t 2 -t 3 ) of 30 days is assumed for incident investigation, preparation for restoration. No adaptation strategies are available in the current case. In the restoration stage, the damaged tanks are rebuilt according to the tank volume (descending order).

Results
According to the case study description, a stochastic dynamic resilience simulation can be conducted according to the methodology presented in Section 3 and the algorithm in Section 4. A desktop PC (CPU: Intel(R) Core(TM) i5, RAM: 8G) is used to carry out the simulation. If the number of iterations N is set to be 10 4 , the computation time is around 2s and the difference of storage resilience values between two calculations is less than 5/1000. While the computation time is around 105 s the difference is lower than 1/1000 when the N is equal to 10 5 . Since the accuracy difference between N = 10 5 and N = 10 4 is ignored, we select N = 10 5 for the computation in this paper. In this study, the uncertainties in resilience are considered. As a result, there are many (10 5 ) resilience evolution scenarios and each resilience scenario has a resilience value. Ranking these resilience values, we can obtain a minimum value, a mean value and a maximum value. Fig. 6 only shows three evolution scenarios (three curves). The black curve represents the resilience scenario with the maximum resilience value in which the disruption doesn't lead to any damage and performance reduction. The red curve represents a resilience evolution scenario with the mean resilience value while the blue curve denotes a resilience evolution scenario with the minimum resilience value. The storage resilience R (mean value) is equal to 0.822. The maximum value of R is 1 while the minimum value is 0.417. Fig. 6 shows the resilience distribution of the storage tank farm.
In light of the large difference between the minimum value and the maximum value, the stochastic characteristics of resilience cannot be ignored. Due to the overpressure caused by the explosion, two tanks (T3 and T12) close to the explosion position failed immediately in the two resilience scenarios. However, during the escalation stage, the residual 12 tanks are damaged in the minimum resilience scenario while only 5 tanks (T2, T5, T1, T4, T13) failed in the mean resilience scenario. As a result, the needed time to full recovery of the minimal resilience scenario (465 days) is much longer than that of the mean resilience scenario (945 days).
According to the escalation calculation results, the t 2 of the resilience evolution scenario with the minimum resilience is 38 min while that of the resilience evolution scenario with the mean resilience value is 15 min. Since the unit of the x-axis of Fig. 6 is the day, the difference is not obvious. Besides, in this case, the storage performance improved by adaptation is equal to zero (no available adaptation strategies). Consequently, the first significant performance improvement and the second significant improvement on both curves are caused by the restoration of two damaged tanks with the same volume. Suppose the escalation effects are not considered in the case study, the resilience increase from 0.821 to 0.905. Fig. 7 shows a typical resilience scenario without domino effects (red curve) and a typical resilience scenario considering escalation effects. Most of the red curve is lower than the blue curve, indicating that escalation effects have an ignored impact on the storage resilience of hazardous materials. Consequently, the resilience of hazardous material storage plants may be overestimated if escalation effects are neglected.
To further explain the role of escalation effects in storage resilience, Fig. 8 shows the failure probabilities of storage tanks exposed to only the attack and both attacks and possible escalation effects.
The blue bars represent the failure probability directly caused by the explosion, while the orange bars denote the failure probability considering escalation effects. It shows that T3 and T12 are prone to be directly damaged by the explosion attack since they are close to the explosion position. However, other tanks are more likely to be damaged by the escalation effects caused by the explosion attack. For example, T14 has a very low probability (0.00002) to be directly damaged by the attack while the failure probability is 0.17 due to the possible escalation effects caused by the damage of T12. If escalation effects are neglected, the failure probabilities of tanks would be underestimated, resulting in an overestimation of storage performance and resilience.

Discussion
Based on the case study in Section 5, this section will discuss the resilience model parameters and possible resilience measures to improve the resilience of hazardous material storage plants.

Resistance capability analyse
In light of explosion attacks, the resistance capability of storage tanks mainly depends on the TNT equivalent mass of the explosion, as shown in Fig. 9a. With increasing the TNT equivalent mass of the explosion, the

Mitigation capability analyse
In hazardous material storage plants, a disruption may lead to major accident scenarios and result in escalation effects. The mitigation capability refers to the safety barriers that can mitigate the consequences of disruptions by preventing or mitigating escalation effects. By applying a water spray system (WSS), the possible heat radiation can be partly reduced, thus increasing the resilience, as shown in Fig. 10a. The storage resilience increases with the reduced heat radiation increasing until escalation effects are almost prevented (60%). The minimal resilience value largely increases when the reduced heat radiation increases from 50 to 69% since escalation effects are almost impossible when the reduced heat radiation exceeds 60%. This can be seen in the failure probability curves in Fig. 10b. When the reduced heat radiation exceeds 60%, the failure probability of each tank remains unchanged (i.e., the tank failure can only be directly caused by the overpressure caused by  the attack rather than escalation effects). The results indicate that further decreasing the possible heat radiation may not be cost-effective when the reduced heat radiation is more than 60%. In Fig. 10b, the failure probabilities of T12 and T3 are not obviously decreased by the reduction of heat radiation since they are close to the explosion position. In that case, the dominant cause of the failure of the two tanks is explosion rather than domino effects caused by other tanks. Therefore, the reduction of heat radiation doesn't have obvious effects on the two tanks.
Similar to the water spray system (WSS), an emergency response such as firefighting can also be used to prevent escalation effects. The mean time for emergency response u is a critical parameter for mitigation capability, as shown in Fig. 11.
The storage resilience (i.e., the red curve in Fig. 11a) decreases with increasing μ because the failure probability of most of the tanks increases with increasing μ, as shown in Fig. 11b. In Fig. 11b, the failure probabilities of T3 and T12 are much higher than other tanks when μ is less than 10 min since they are more likely to be directly damaged by the blast overpressure caused by the explosion, as shown in Fig. 8. The failure probabilities of residual tanks increase with increasing μ since delayed emergency response can lead to more severe escalation effects. The failure probabilities of T12-T14 cannot exceed the failure probability of T12 since there are two domino islands where no domino effects can occur in between [70]. One island consists of T1-T11 and another consists ofT12-T14. As a result, T13 and T14 can only be damaged by the escalation effects caused by T12 when T12 is damaged by the attack. Therefore, the failure probabilities of T13 and T14 are almost equal to the damage probability of tank 12 (caused by the explosion) in which the probability of domino effects equals 1 when emergency response is not available on time.

Adaption capability analyse
Adaption measures such as utilizing reserve tanks and speeding inventory turnover ratio may partially compensate for the performance loss caused by a disruption before the storage plant is fully restored. The adaption capability is limited by the storage equipment (reserve tanks), loading, and unloading facilities. Fig. 12a shows the effects of speeding inventory turnover (represented by the increased inventory turnover rate (%) based on the storage capability at the end of the escalation stage) on storage resilience. As shown in Fig. 12a, the storage resilience (red curve) increases with an increasing inventory turnover rate in the adaption stage. The minimum resilience value is constant since all tanks are damaged in the minimal resilience scenario and speeding inventory turnover does not work. Besides the performance adaption, shortening the adaption time (t 2 -t 3 ) can also improve the adaption capability, as shown in Fig. 12b. Both the storage resilience (red curve) and the minimal storage resilience (blue curve) in Fig. 12b show a decreasing trend with increasing adaption time. As a result, reducing the adaption time and starting restoration early is also an adaption measure to enhance resilience.

Restoration capability analyse
Restoration is the final stage of a resilience process. The restoration time is a crucial parameter for the restoration capability, as shown in Fig. 13. As shown in the figure, the resilience is inversely proportional to the restoration time. Both the resilience value and the minimal resilience value increase with the decrease of restoration time. As a result, a quick restoration capability is essential for developing a resilient hazardous material storage plant. Different restoration strategies (e.g., restoration sequence) may lead to different resilience capabilities. In this study, the restoration sequence is based on the tank volume: ascending order (from small to large) and descending order (from large to small). There is no apparent difference (both are 0.822) between ascending order restoration and descending order restoration since the tank construction time is proportional to the tank volume. If the construction time is not proportional to the tank volume, the restoration sequence makes a difference. For instance, the resilience based on the descending sequence (0.861) is much higher than that based on ascending sequence (0.777).
According to the above analyse, there are many measures in different resilience stages to enhance the resilience of hazardous material storage, such as safety barriers in the escalation stage, speeding inventory turnover in the adaptation stage, and shortening restoration time in the restoration stage. Besides, these measures, inherent safety design [21] may also be used to improve storage resilience by preventing and mitigating the primary major accident scenarios and possible escalations. In this study, the costs of different resilience measures are not considered. In the future, resilience management approaches may be developed by combining the resilience quantification method developed in this study with economic tools such as cost-benefit analyse and cost-effectiveness analyse. Furthermore, with the digitalization of the process industry, storage plants may be vulnerable to cyber-threats due to a high level of automation and an increasing connection with external networks [40]. Therefore, storage resilience may be integrated with cybersecurity [30] in the future. Moreover, the robustness of resilience with respect to different types and intensities of threats may also be addressed in future research work.
In this study, we define the HAZMAT storage resilience as the capability of a storage plant to resist, mitigate, adapt and recover from disruptions, to maintain its storage performance. The resistance capability and mitigation capability concern damage aspects, while the adaptation and restoration capability consider the capability to maintain and recover from damages. All four capabilities are modeled in this study. This study pays attention to damage aspects of resilience since hazardous material storage plants are vulnerable to major accidents and domino effects, which is the main difference from other critical infrastructures. This study is the first quantitative study on the resilience of hazardous material storage plants; more work should be done in the future to improve the modeling on adaptation and restoration capabilities that depends more on human and organizational factors.
This study is the first work on HAZMAT storage resilience, so a widely used resilience quantification method ( [7]; Bruneau et al., 2009) is selected to show resilience mechanism and highlight the role of resilience in the protection of storage plants. In the future, we will conduct more research on restoration strategies. In that case, more advanced resilience metrics such as the "center of resilience" [ [73], [74]] may be selected to model the change of more restoration parameters such as restoration time and the final restored performance.
In terms of the mitigation capability assessment, the vulnerability models used in this study for fire escalation assessment may be further improved by experiments and numerical simulations, considering more influence factors, such as design factors, failure types, safety barriers, etc.

Conclusion
The storage resilience of HAZMAT is time-dependent and uncertain. This paper proposes a dynamic stochastic methodology to measure it, considering possible escalation effects and the recovery of damaged installations. In this methodology, the dynamic resilience process is divided into four stages: disruption, escalation, adaption, and restoration stages. This model considers the uncertainties related to the vulnerability of tanks exposed to disruptions, escalation effects, and emergency response time. The dynamic Monte Carlo method is used to simulate possible resilience evolution scenarios and thus obtain storage resilience. Compared with traditional safety and security risk assessment, the developed resilience methodology addressing the roles of adaptation and restoration, which is more suitable for tackling unpredictable disruptions. Finally, a case study is provided to demonstrate and test the proposed methodology and algorithm. The primary conclusions are: (1) the resilience values can range in a large interval (in the case study, they are between 0.4 and 1) due to the uncertainties in the dynamic resilience process. As a result, the uncertainties in the resilience process cannot be ignored in resilience modeling; (2) escalation effects play an essential role in hazardous material resilience; neglecting possible escalation effects may underestimate the storage resilience; (3). the storage resilience depends on the intensity of disruptions, the storage plant's resistance capability, mitigation capability, adaption capability, and restoration capability; (4) resilience measures such as safety barriers in the escalation stage, speeding inventory turnover in the adaptation stage, and shortening restoration time in restoration stage are effective for developing a more resilient storage plant; (5) economic tools such as cost-benefit analyse and cost-effectiveness analyse may be used in this study to develop a resilience management approach for hazardous material storage.