DOI: 10.3303/CET23104003 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Paper Received: 16 April 2023; Revised: 15 July 2023; Accepted: 2 August 2023 
Please cite this article as: Khakzad N., Amyotte P., 2023, Optimal Evacuation of Process Plants in Case of Tank Fires, Chemical Engineering 
Transactions, 104, 13-18  DOI:10.3303/CET23104003 
  

 CHEMICAL ENGINEERING TRANSACTIONS  
 

VOL. 104, 2023 

A publication of 

 

The Italian Association 

of Chemical Engineering 

Online at www.cetjournal.it 

Guest Editors: Sabrina Copelli, Ernesto Salzano 
Copyright © 2023, AIDIC Servizi S.r.l. 

ISBN 979-12-81206-06-9; ISSN 2283-9216 

Optimal Evacuation of Process Plants in Case of Tank Fires 

Nima Khakzada,*, Paul Amyotteb 

aToronto Metropolitan University, Toronto, Canada  
bDalhousie University, Halifax, Canada 

 nima.khakzad@torontomu.ca  

Effective firefighting and evacuation as integral parts of emergency response in petrochemical plants play a key 

role in protecting human lives in the event of major tank fires. Compared to firefighting, however, studies devoted 

to planning and optimizing evacuation plans in the event of tank fires, and particularly concurrent tank fires, 

have been very few. In the present study, considering the thermal dose as the main cause of casualties in 

outdoor fires, an innovative methodology is developed to identify proactive evacuation plans for credible fire 

scenarios. The methodology consists of three main parts: (1) For a given fire scenario (e.g., a single or multiple 

tank fires), the tank terminal is modelled as a thermal graph in which the weight of each node presents the 

corresponding heat flux, and the weight of each edge presents the thermal dose between the connected nodes; 

(2) Dijkstra’s algorithm is used to find the shortest paths (a series of connected edges with the least total thermal 

dose) to the safe spots (e.g., shelters); (3) Considering the limited capacities of safe spots, mathematical 

programming is used to identify the number of evacuees to be assigned to each safe spot so as to minimize the 

total risk of casualties during evacuation. Application of the methodology to an illustrative process plant resulted 

in intuitive evacuation plans, which is indicative of the methodology’s validity, particularly in the absence of 

similar studies for comparison and validation purposes.  

1. Introduction 

Among industrial plants, chemical and process plants are very susceptible to accidental and intentional fires 

due to the presence of large quantities of flammable chemicals and a variety of ignition sources. Fire at chemical 

and process plants may vary from small fires, which can safely be handled by portable fire extinguishers, to 

major tank fires, which demand emergency response measures, including plant evacuation and firefighting. 

Despite several studies on fire protection (Janssens et al., 2015; Landucci et al., 2015; Ovidi et al., 2021) and 

optimal firefighting of process plants (Khakzad, 2018, 2021), works devoted to evacuation of chemical and 

process plants in the event of major fires have relatively been very few.  

Previous research has usually come up with general recommendations for evacuation during major safety and 

security events (Hosseinnia, et al., 2018), has investigated the issue from an economic perspective – is it 

economically viable to shut down and evacuate the plant? (Reniers et al., 2007), or has been devoted to 

evacuation in the event of accidents other than fires (Zhang et al., 2017; Yoo and Choi, 2019). Also, available 

guidelines and standards such as CAAP (2010) and API RP 2001 (2011) are too general to be used to identify 

optimal evacuation plans in chemical and process plant fires.  

The present study is aimed at developing a methodology for optimal evacuation of process plants in the event 

of major tank fires. To this end, considering the thermal dose as the main cause of casualties, Dijkstra’s 

algorithm is used to identify the shortest evacuation path from an operating unit to a safe spot given single or 

concurrent tank fires, with the distance from node i to node i+1 being equal to the received thermal dose from i 

to i+1. This way, the shortest path would be a path with the least accumulated thermal dose, and would thus be 

the safest path. That being said, throughout this manuscript, shortest path and safest path are used 

interchangeably. In case of having multiple limited-capacity safe spots, mathematical programming can be used 

to optimally disperse the evacuees among the safe spots in order to minimize the risk of casualties. This way, 

pre-planned evacuation procedures can be identified and presented in the form of decision-making tables to 

facilitate evacuation of process plants in the event of fire scenarios. 

13

mailto:nima.khakzad@torontomu.ca


2. Methods and Materials 

2.1 Thermal grid 

Exposed to fire, the thermal dose (D) received by a human is dependent on the intensity of heat flux (q) and 

exposure time (te) (Assael and Kakosimos, 2010): 

 

𝐷 = 𝑡𝑒 . 𝑞4/3    (1) 

 

For a fire, q (W/m2) at a distance can be calculated using a variety of analytical and numerical models (e.g., the 

point source mode). The exposure time te(s) can be calculated using the relationship below (Assael and 

Kakosimos, 2010): 

 

𝑡𝑒 = 𝑡𝑟 +
𝐿

𝑢
   (2) 

 

where tr (s) is the reaction time; L (m) is the distance the exposed person should run to get the safe spot; u is 

the average speed of the person (~ 4 m/s). For trained personnel who expect fires and explosions at the 

workplace, tr ~ 3 s while for a lay person tr ~ 8 s.  

For a given fire scenario, such as a tank fire, a process plant can be modeled as a two-dimensional grid. The 

intensity of heat flux for each node of the grid can be calculated and assigned as the weight of that node. 

Subsequently, by combining Eq (1) and Eq (2), the thermal dose an evacuee may receive while escaping from 

their location (at node i = 0) to the safe spot (at node j) can be calculated as (Khakzad, 2023): 

 

𝐷 = 𝑡𝑟  𝑞0

  
4

3 + ∑ (
𝑞𝑖+𝑞𝑖−1

2
)

4

3
 
∆𝑥

𝑢

𝑗
𝑖=1     (3) 

 

where 𝑞𝑖 and 𝑞𝑖+1 are two adjacent nodes on the grid; ∆𝑥 is the distance between the adjacent nodes; q0 is the 

heat flux at the initial location of the evacuee. As an example, consider a thermal grid in Figure 1 where for a 

given fire, an evacuee should leave their initial location at node i = 0 and seek shelter at the safe spot at node i 

= 3. In doing so, the evacuee decides to traverse nodes i = 0, 1, 2, and 3 in sequence. The heat flux intensity at 

each node has been calculated and shown as q inside the node.  

 

 

Figure 1. Thermal grid given a fire in the process plant. Each node is indicated with two indices: i refers to the 

node number while q (kW/m2) shows the intensity of heat flux at that node. An evacuee is supposed to escape 

from their location at i = 0 to a safe spot at i = 3. The presented path here is for illustrative purposes and is not 

necessarily the shortest (safest) available path.  

Further, assume that the distances between the adjacent nodes are ∆01= 40 m, ∆12= 20 m, and ∆23= 40 m. 

Assuming tr = 3 s and u = 4 m/s, the total thermal dose the evacuee would receive while escaping from the ‘initial 

location’ to the ‘safe spot’ can be calculated as:  

D = tr q0

  
4

3 + ∑ (
qi+qi−1

2
)

4

3
 
∆x

u

j
i=1 = 3(12000)4/3 + (

12000+8000

2
)4/3 40

4
+ (

8000+8000

2
)4/3 20

4
+ (

8000+2000

2
)4/3 40

4
=

4,633,617 (W4/3m-8/3s).  

14



Given the thermal dose, dose-response relationships can be used to estimate the probability of injury (e.g., 2nd 

degree burns) or death. The following probit function can be used to estimate the death probability due to heat 

exposure (Assael and Kakosimos, 2010): 

 

Y =  −36.38 + 2.56 ln(D)   (4) 

𝑃 = 𝐹𝐶 . Φ(𝑌 − 5)   (5) 

Where Y is the probit value; FC is the clothing coefficient, accounting for the role of clothing in protecting an 

exposed human, with typical values of 0.85 and 0.18 for summer and winter outfits, respectively; Φ(.) is the 

cumulative density function for standard normal distribution. 

2.2 Dijkstra’s algorithm  

Dijkstra’s shortest path algorithm (Dijkstra, 1959) is a greedy algorithm to identify the shortest path between two 

nodes on a graph (or mesh). In Dijkstra’s algorithm, the state of a node is identified by two variables: the distance 

value of the node dj, which is its distance from the starting node, and the status label of the node, which can be 

permanent (p) or temporary (t). The status label is permanent only if the distance value of the node is the shortest 

distance from the starting node. Dijkstra’s algorithm measures and updates the distance values of the nodes 

until all the status labels become permanent. The algorithm for finding the shortest paths from a starting node 

to other nodes of the graph includes three stages: Initialization, Updating, and Termination.  

• Initialization  

Assign zero to the distance value and ‘permanent’ to the status label of the starting node i = 0: The state of this 

node is now (0, p). Assign to every other node a distance value ‘ꚙ’ and a status label ‘temporary’. The states 

of the other nodes are now (ꚙ, t). Designate the starting node i = 0 as the ‘current’ node. 

• Updating distance values  

Find the set J of nodes with temporary labels that can be reached from the current node via a single link. For 

each j ∈ J, the distance value of node j is updated as: dj = min {dj, dj + wij}, where wij is the weight of the edge 

linking node i to node j. Among the nodes, the one with the smallest updated distance value (node j*) is added 

to the set of current nodes and its status is changed to ‘permanent’. The state of this node is thus updated to 

(dj*, p). 

• Termination 

The previous step should be repeated by finding the set of nodes with temporary labels that can be reached 

from any of the current nodes i or j∗ via a single edge. The algorithm terminates when the status labels of all the 

nodes (or the destination node of interest) that can be reached from the starting node i = 0 turns to ‘permanent’.  

An application of Dijkstra’s algorithms has been exemplified in Figure 2 to find the shortest path between nodes 

A and D. 

 

A

B C

D

(0, p)

(ꚙ, t)(ꚙ, t)

(ꚙ, t)

2

2

5

4 1

(a)
  

A

B C

D

(0, p)

(ꚙ, t)(2, p)

(ꚙ, t)

2

2

5

4 1

(b)
  

(c)

A

B C

D

(0, p)

(4, p)(2, p)

(ꚙ, t)

2

2

5

4 1

  

A

B C

D

(0, p)

(4, p)(2, p)

(5, p)

2

2

5

4 1

(d)
 

Figure 2. Application of Dijkstra’s algorithm to find the shortest path and distance from node A to node D. Panels 

(a)-(d) present the sequential steps via the algorithm (Khakzad, 2023). 

2.3 Mathematical programming 

Mathematical programming includes a number of techniques to maximize or minimize the value of a function 

(known as objective function) by finding the best values of the function parameters. Maximum or minimum value 

of the function and thus the optimal values of its variables would be subject to some constraints. In some 

mathematical programming techniques, a feasible solution (i.e., the set of optimal values of variables) can be 

obtained only if all the constraints are satisfied. However, some other techniques such as goal programming 

15



can render a feasible solution even if some less important constraints are violated in favor of more important 

constraints. Mathematical programming in the form of minimization, for instance, can be presented as: 

 

𝑚𝑖𝑛𝑖𝑚𝑖𝑧𝑒 𝑍 = 𝑓(𝑋)   (6) 

𝑆𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜: {
𝑔(𝑋) ≤ 0

ℎ(𝑋) = 0
     

 

where 𝑋 = {𝑋1, … , 𝑋𝑛} is the set of variables, or function parameters; Z is the objective function as a function of 

X; and g(X) and h(X) are the constraint functions.  

3. Methodology  

3.1 Case study  

Consider a process area with two operating units, two gasoline storage tanks and two fire shelters in Figure 2. 

Given concurrent tank fires at Tank #1 and #2, the process area can be modelled as a thermal grid showing the 

intensity of heat flux at each node of the grid. The heat flux intensity at each node can be calculated for each 

tank fire via a number of models and techniques (e.g., the point source model) knowing the tank size and 

weather conditions, and then in case of multiple simultaneous fires these heat flux intensities can be 

superimposed to find the total heat flux intensity at each node. For the sake of exemplification, assume that the 

heat flux intensities for each tank fire were already calculated, superimposed, and presented inside each node 

with q. This way, each node of the grid can be indexed with a number ‘I’ and its corresponding heat flux intensity 

‘q’. 

 

         

(a)  (b) 

Figure 2. (a) A process area, consisting of two operating units, two gasoline storage tanks, and two fire shelters. 

(b) Modelling the process area as a thermal grid given tank fires at the storage tanks. Each node has been 

identified with a node number i and the intensity of heat flux q at that node (in kW/m2).   

Knowing the heat flux intensities and the distance between adjacent nodes in Figure 2, the thermal dose one 

would receive while traversing between two adjacent nodes, either horizontally or diagonally, can be calculated 

using Eq (3). The results can then be summarized in the form of a thermal-dose graph in which the weight 

assigned to the edge connecting node ‘I’ and ‘i+1’ shows the corresponding thermal dose when one moves from 

the former to the latter or vice versa. Regarding the process area in Figure 2, since the aim is to find the shortest 

paths from the operating units (nodes i = 18, 25) to the shelters (nodes i = 1, 10), 𝑡𝑟  𝑞0

  
4

3 on the right-hand side 

of Eq (3) can be ignored as this term contributes equally to all the paths and thus its omission would be fine for 

comparison purposes (Khakzad, 2023). Figure 3 shows the thermal-dose graph for the process area where the 

weight of each edge shows the respective thermal dose in 105 W/m2, presuming that an evacuee would not 

run toward the tank fires while escaping to the shelters. Due to this presumption, the edges running from tank 

fires’ adjacent nodes toward the tank fires have been deleted from the graph for simplification.  

 

16



   

Figure 3. Modelling the process area as a thermal-dose graph. The weights assigned to each edge show the 

corresponding thermal dose while an evacuee traverses between the nodes connected by that edge. Presuming 

that an evacuee would not run directly toward a tank fire, such edges have been deleted from the graph. The 

node numbers are presented with bold numbers.   

3.2 Identification of shortest paths 

Since a smaller thermal dose is attributed to a smaller probability of injury or death, an optimal evacuation path 

between the operating units (nodes i = 18 & 25) and the shelters can be considered as a path (a sequence of 

connected edges) with the least accumulated thermal dose. 

To this end, the shortest paths from Unit #1 (node i = 18) to the shelters (nodes i = 1 & 10) can be identified by 

applying Dijkstra’s algorithm as:  

• From Unit #1 to Shelter #1: i = 18 → i = 22 → i = 21 → i = 16 → i = 11 → i = 6 → i = 1, with the total thermal 

dose 77 × 105 kW/m2, and 

• From Unit #1 to Shelter #2: i = 18 → i = 14 → i = 9 → i = 10, with the total thermal dose 74 × 105 kW/m2.  

Similarly, the shortest paths from Unit #2 (node i = 25) to the shelters (nodes i = 1 & 10) can be identified as: 

• From Unit #2 to Shelter #1: i = 25 → i = 20 → i = 15 → i = 10 → i = 4 → i = 3 → i = 2 → i = 1, with the total 

thermal dose 43.5 × 105 kW/m2, and 

• From Unit #2 to Shelter #2: i = 25 → i = 20 → i = 15 → i = 10, with the total thermal dose 41 × 105 kW/m2.  

Using Eq (4) and Eq (5), and considering a cloth effect of FC = 0.5, the probability of death for a single evacuee 

who may escape from unit i to shelter j (Pij) is calculated as P11 = 10.5%, P12 = 9.5%, P21 = 1.5%, and P22 = 

0.5%. 

3.3 Optimal allocation of evacuees  

Having the shortest paths identified by applying Dijkstra’s algorithm, mathematical programming can be 

employed for optimal allocation of evacuees to the shelters. To demonstrate this, assume that the number of 

evacuees at Units #1 and #2 are, respectively, 5 and 10, while the capacity of Shelters #1 and #2 are, 

respectively, 15 and 8. The goal is thus to determine the number of evacuees from each unit that should be sent 

to each shelter in order to minimize the number of casualties. Considering Xij as the number of evacuees that 

should be sent from Unit i to Shelter j, the optimization problem can be formulated as: 

 

𝑀𝑖𝑛𝑖𝑚𝑖𝑧𝑒 𝑍 = 𝑃11𝑋11 + 𝑃12𝑋12 + 𝑃21𝑋21 + 𝑃22𝑋22    (7) 

 

𝑆𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜: {

𝑋11 + 𝑋21 ≤ 15
𝑋12 + 𝑋22 ≤ 8  
𝑋11 + 𝑋12 = 5   
𝑋21 + 𝑋22 = 10

    (8) 

 

17



The first two constraints in Eq (8) are to ensure that the number of evacuees at each shelter will not exceed the 

shelter capacity while the last two constraints are to ensure that no one would be left out. Solving Eq (7) and Eq 

(8) in the Microsoft Excel Solver Tool Pak, the optimal values of the evacuation variables Xij are identified as: 

X11 = 5, X12 = 0, X21 = 2, X22 = 8, with the total number of casualties as Z = 0.6. The result is completely intuitive 

as evacuees from each unit are first sent to the nearest (safest) shelter (Unit #1 is near Shelter #1 and Unit #2 

is near Shelter #2) and then, anyone left is sent to the second nearest shelter.  

4. Conclusions 

Evacuation is an integral part of fire safety and protection, which plays a key role in protecting human lives in 

the event of major accidents in chemical and process plants. Among the major accidents, however, the majority 

of the previous studies have developed evacuation methodologies for release and dispersion of toxic chemicals 

with too little attention, if any, to major fires and particularly concurrent multiple fires. 

In the present study, considering the thermal dose as the main cause of casualties during major industrial fires, 

a methodology was developed by integrating Dijkstra’s algorithm and mathematical programming for optimal 

evacuation. For this purpose, Dijkstra’s algorithm was employed to identify the safest path – a path with the 

least accumulated thermal dose – from a unit to a safe spot while mathematical programming was used to 

identify the number of evacuees that should be assigned to each safe spot.  

The applicability of the developed methodology can further be improved by considering a wider range of process 

plants and fire scenarios and the combined effect of heat and smoke in the form of a heat-smoke-graph. 

Acknowledgments 

Financial support provided by the Natural Sciences and Engineering Research Council of Canada (NSERC) via 

a Discovery Grant (RGPIN/03051-2021) and a Launch Supplement (DGECR/00220-2021) is appreciated. 

Publication Grant, Faculty of Community Services, Toronto Metropolitan University is gratefully acknowledged.  

References 

API RP 2001, 2011, Fire Protection in Refineries. Candidate Ballot Draft for API 2001 9th Edition, API Publishing 

Services, Washington, DC. 

Assael M.J., Kakosimos K.E. (Eds), 2010. Fires, Explosions, and Toxic Gas Dispersions: Effects Calculation 

and Risk Analysis, CRC Press, Boca Raton, US.  

Dijkstra E.W., 1959, A note on two problems in connexion with graphs, Numerische Mathematik, 1, 269-271. 

Hosseinnia B., Khakzad N.,Reniers G., 2018, Multi-plant emergency response for tackling major accidents in 

chemical industrial areas, Safety Science, 102, 275-289.  

Janssens J., Talarico L., Reniers G., Sörensen K., 2015, A decision model to allocate protective safety barriers 

and mitigate domino effects, Reliability Engineering & System Safety, 143, 44-52. 

Khakzad N., 2018, Which fire to extinguish first? A risk-informed approach to emergency response in oil 

terminals, Risk Analysis, 38(7), 1444-1454. 

Khakzad N., 2021, Optimal firefighting to prevent domino effects: Methodologies based on dynamic influence 

diagram and mathematical programming, Reliability Engineering & System Safety, 212,107577. 

Khakzad N., 2023, A methodology based on Dijkstra’s algorithm and mathematical programming for optimal 

evacuation in process plants in the event of major tank fires, Reliability Engineering & System Safety, 236, 

109291. 

Landucci G., Argenti F., Tugnoli A., Cozzani V., 2015, Quantitative assessment of safety barrier performance in 

the prevention of domino scenarios triggered by fire, Reliability Engineering & System Safety, 143, 30-43.  

Ovidi F., Zhang L., Landucci G., Reniers G., 2021, Agent-based model and simulation of mitigated domino 

scenarios in chemical tank farms, Reliability Engineering & System Safety, 209, 107476. 

Reniers G., Pauwels N., Audenaert A., Ale B., Soudan K.,2007, Management of evacuation in case of fire 

accidents in chemical industrial areas, Journal of Hazardous Materials, 147, 478-487. 

The Canadian Association of Petroleum Producers (CAPP), 2010, Atlantic Canada Offshore Petroleum Industry 

Escape, Evacuation and Rescue. St. John’s, NL, Canada, June 2010, 1-30. 

Yoo B., Choi S.D., 2019, Emergency evacuation plan for hazardous chemicals leakage accidents using GIS-

based risk analysis techniques in South Korea, International Journal of Environmental Research and Public 

Health, 16, 1948. 

Zhang J.H., Liu H.Y., Zhu R., Liu Y., 2017, Emergency evacuation of hazardous chemical accidents based on 

diffusion simulation, Complexity, 4927649.  

18


	60khakzad.pdf
	Optimal Evacuation of Process Plants in Case of Tank Fires