Numerical Modeling of the Macroscopic Behavior of a Crowd of People under Emergency Conditions Triggered by an Incidental Release of a Heavy Gas: an Integrated Approach

In the present study, we focus on modeling and simulations of macroscopic behavior of a crowd of people under emergency conditions. We present a newly developed integrated numerical approach, which combines the time-dependent RANS method in predicting the temporal and spatial evolution of a heavy gas turbulent dispersion with the macroscopic model of a crowd movement. The coupling is imposed through locally determined concentration of the heavy gas which directly impacts the movement of the crowd of people. The potential of the developed approach is demonstrated in a series of cases that include the crowd behavior on an open railway platform and inside a train station, as well as predictions of field studies of turbulent dispersion of a heavy gas under windy conditions. Finally, various scenarios have been analyzed of the coupling between the local concentration of a potentially harmful heavy gas on crowd behavior through different forms of the cost functions. In conclusion, the version of the integrated model is recommended for future optimization studies of crowd dynamics following an incidental release of heavy gas.


Introduction
The constant increase in population density and higher population mobility in urban areas leads to an increased risk of various situations that can trigger a sudden movement of large crowds of people. In the present literature, there are various approaches dealing with predictions of the crowd behavior. Among these approaches, the most common techniques are based on: 1. a cellular automata model [1,2], 2. a lattice gas model [3], 3. a social force model [4], 4. a fluid-dynamics similarity [5][6][7], 5. an agent-based model [8], 6. the game theories [9], and finally, 7. the studies and experiments with animals [10,11].
In the present work, we focus on a continuum description of a crowd of people movement, i.e. we consider modeling of the flow of human crowds. Our choice is primarily based on our expertise in classical fluid mechanics and because of the fact that there are numerous similarities between movement of large crowds and a fluid flow, i.e. a large crowd can be considered as a kind of 'thinking fluid', as discussed in great details in [7]. Same author also postulated three main mechanisms that determine the motion of crowd: (1) the density of crowd, (2) a common sense of crowd to reach exit (so called potential), (3) a tendency of crowd to minimize travel time (evacuation time) and to avoid extreme densities, [7]. All these parameters are then mathematically defined and included in the time-dependent crowd continuity equation, as shown in [6,7,12].
Here, we present details of development of an integrated in-house computation code for predictions of a crowd movement under influence of an incidental release of a heavy gas. The first part of the code is based on a macroscopic description of a crowd dynamics (fluiddynamics similarity) and discretised form of the crowd continuity equation, [6,12]. This code is then integrated with an existing in-house code based on solving a full set of the time-dependent Navier-Stokes equations (T-RANS) [13,14] with additional terms due to concentration buoyancy in order to properly mimic the behavior of a heavy gas. The novel contribution is in a fully coupled approach, where we propose how to redefine some of the model parameters (so called cost functions) used in the crowd dynamics model, in terms of the local concentration of a heavy gas.

Crowd evacuation model
The conservation of a total crowd mass is used for macroscopic behavior of the flow of a crowd of people [6,7,12]: where ρ crowd ( x, t) (in people/m 2 ) is the local crowd density, and v(x, t) represents the local mean velocity of the crowd, which will be calculated from v = V n, n = − ∇φ |∇φ| (2) where n is the unit vector that describes the direction in which the crowd will move, V is the local travel speed, and φ is the potential that is a measure of the travel costs which pedestrians want to minimize, given by the following Eikonal equation: where c(ρ crowd ) is the cost function, is the domain interior, 0 is the exit, w is the wall, |∇φ| is the magnitude of the gradient of the potential (in s/m). The equation indicates that the potential (φ) (defined as φ = 0 at the exit) increases with distance from the exit by the local travel costs per unit distance. The potential can be simply interpreted as a time needed by a moving crowd to reach the exit. Since we are dealing with the dynamically adjusting local conditions (the spatial and temporal variation of the crowd density), the gradient of the potential will be expressed in terms of the cost function (c (ρ crowd )). The system of equations (1)-(3) will be closed if this cost function is known, which is calculated as: where V max and ρ max are maximum velocity and crowd density, respectively, defined as the model inputs. It can be seen from (4) that the local velocity of crowd approaches zero when the local crowd density is equal to the maximum crowd density, which generates infinitely large cost function. As a consequence, the rest of crowd will try to avoid these regions. The situation is opposite for the low crowd density regions where the local travel speed will be close to the maximum velocity.

A heavy gas turbulent dispersion
The mass and extended momentum equation takes into account the concentration buoyancy term are written as: with where density of the gas mixture is calculated as Here C is a heavy gas concentration fraction (C = 1 for a heavy gas and C = 0 for air) calculated from the following transport equation: The concentration buoyancy terms are also included in modeled transport equations of the turbulent kinetic energy (k) and its dissipation rate (ε) as: where The turbulent stress and turbulent concentration flux are expressed through the Boussinesq assumption and generalized gradient diffusion hypothesis (GGDH), respectively: where the model coefficients are given in Table 1.

Coupling of the crowd evacuation model and a heavy gas dispersion: the redefined cost function
The coupling between the behavior of the crowd and the local concentration of a heavy gas is implemented through two main mechanisms: redefinition of the perceived cost function (i.e. the crowd tries to move away from the dangerous cloud) and movement reduction (i.e. to take into account possible incapacitating effects). The first contribution is evaluated as c ρ crowd , C gas = 1 V (ρ crowd ) + αC gas (14) where the parameter α is calculated as Table 1 The turbulence model coefficients with C ref gas as the referent value of the concentration when the crowd starts to avoid a direct movement through the defined region. Finally, the second mechanism (movement reduction) is modeled as where f is the fraction of the crowd which will slow down immediately when it encounters a gas cloud defined with a local referent concentration C ref gas . This reduction of the crowd movement can be caused by various effects, e.g. reduced visibility and/or breathing difficulties due to presence of the heavy gas cloud. Finally, as the very last step in the coupling mechanism, we introduce the critical levels of the heavy gas concentration (C crit gas ) for the point at which the people in the crowd start to lose consciousness and collapse, i.e. the movement speed is automatically set to zero when the critical level of heavy gas concentration reached.

Numerical Method
Here we provide the most important details of numerical methods developed and used in the present research for predictions of crowd movement and the turbulent dispersion of heavy gas. In essence, we have two separate algorithms (crowd dynamics and CFD of heavy gas turbulent dispersion), which are coupled through a/the two-dimensional projection of the local distribution of the heavy gas concentration at the particular distance from the ground. These local distributions of the heavy gas concentration are then used as the input parameters to redefine the particular cost functions of the crowd movement.

Algorithm used for evolution of crowd behavior
The implementation of the numerical algorithm for a crowd movement can be split into three major parts: (i) solving of the Eikonal equations in the domain of interest to obtain the potential; (ii) calculation of fluxes along the cell-faces of two-dimensional control volumes; (iii) calculations of the crowd density. The entire algorithm can be split into the following operations: 1. A numerical mesh of the domain is created and an initial crowd density distribution is set. 2. Using the crowd density distribution, the speed and cost distribution are calculated. 3. The Eikonal equation is solved using the speed and cost distributions, resulting in the potential φ. 4. The gradient of φ is calculated, and, combined with the speed distribution, the velocity field is obtained. 5. The fluxes in each cell are calculated from the density distribution and the velocity field in the center of control volumes.
6. Using the calculated fluxes, the density distribution at the next time (t + t) is calculated. 7. Steps 2 to 6 are repeated until the desired time is reached.
A graphical representation of the algorithm is shown in Fig. 1 (enclosed in red rectangle). Next, we provide some numerical details of the algorithm. The Eikonal equation is solved by applying a fast sweeping method and the Godunov upwind difference scheme, as proposed in [15]. The discretised form of the Eikonal equation (3) can be written as: where φ (i,j ) indicates the potential in the control volume cell (i, j ). Here, the plus sign indicates the positive part: where φ x min and φ y min are defined as the smallest neighboring values in the xand ydirection, respectively, and are calculated as: Boundary conditions at walls are implemented by assigning numerically very high values of the potential. At exits, zero values of the potential are specified and surrounding cells are extrapolated with a radius of 2h into the domain. The discretised form of Eq. 17 has a unique solution: Fig. 1 Graphical representation of the integrated algorithm used for crowd movement (evacuation) (red) and a heavy gas turbulent dispersion (blue), with corresponding remeshing which translates the three-dimensional heavy gas results to a usable two-dimensional form for the updated cost functions of crowd movement where D = 2c 2 (i,j ) h 2 − φ x min − φ y min 2 (21) and the potential can be updated from: By sweeping the domain, alternating in four different orderings, the potential φ is updated until convergence is reached. The alternating ordering used is: With the potential found, the velocity field along which the crowd moves can be solved. The direction is provided by the negative gradient of the potential, whereas its magnitude in each cell is given by the speed distribution provided from: The gradients of the potential are calculated from the central-differencing scheme as With a known velocity field, the corresponding crowd density flux can be calculated as However, we have to evaluate the flux across the cell faces, which can be done by applying the Lax-Friedrich's scheme as: where Finally, the discretised form of crowd conservation can be written as: where x and y are the sizes of the control volume (i, j ) which are taken to be equal to the grid spacing h. The time rate of change is calculated using the implicit 3rd order Runge-Kutta method (here we applied so called Radau IA scheme, [16]) as follows: with where all constants are given in Table 2. With this, we complete the algorithm of the crowd people movement. The particular adjustments for the coupling between a heavy gas dispersion and the crowd evacuation is done through adaptation of discretised forms of particular terms, which follow from Eqs. 14-16.

Numerical details of the finite-volume code for a heavy gas turbulent dispersion
For the CFD part of the algorithm, we have used in-house finite-volume computer code based on discretization of the mass, momentum, energy and concentration equations within the time-dependent Reynolds-Averaged Navier-Stokes (RANS) approach, for the non-orthogonal geometries with a structured non-uniform numerical mesh based on the hexagonal elements, [13]. This computer code is used as a platform for numerous additional extensions to mimic some essential features of atmospheric and environmental flows. Some of the most important features include: (i) possibility to simulate a complex terrain orography [17], (ii) an efficient mesh generation and visualization tools to take into account the presence of buildings within the complex urban areas [13,18,19] (iii) the presence of vegetation and surface roughness effects [20], (iv) ability to use a highly skewed mesh for air flows around hilly or mountain terrains [21], (v) inclusion of complex atmospheric chemistry for proper simulations of the ozone distribution [22], (vi) as well as the radiation and ground/soil conductivity effects as a part of the global energy balance for simulation of the urban heat island phenomenon [23,24].
Here, we will provide just a short description of the most essential numerical details of the CFD algorithm. For the numerical algorithm, the non-staggered grid was employed for all variables with the Rhie-Chow interpolation to prevent decoupling between velocity and pressure fields. The standard iterative predictor-corrector SIMPLE algorithm was used to couple the instantaneous velocity and pressure to satisfy both continuity and momentum equations. The time-rate of change of transport variables was discretised by the fullyimplicit three-consecutive time-step method. The diffusive term of discretised transport equations was approximated with the second-order central-differencing scheme (CDS). The convective term in momentum equation was represented by the second-order quadratic upwind scheme (QUICK), while the second-order linear upwind scheme (LUDS) was applied for the scalar equation.

Numerical details of the integrated approach
To couple the CFD heavy gas turbulent dispersion model and crowd evacuation model, the pre-specified horizontal cross-sections from a three-dimensional concentration field are extracted and projected onto the underlying two-dimensional mesh by using standard bilinear interpolation, as a part of the 'Remesh' procedure shown in Fig. 1. Similarly, because of Table 2 Constants used in the Runge-Kutta discretization, the implicit Radau IA scheme, [16]  different characteristic time-steps used in CFD and crowd evacuation algorithms, a simple linear interpolation is applied to get values of the concentration at particular time-instants that are required to calculate characteristic redefined cost functions for a fully coupled approach.

Validation of the crowd evacuation model
To validate the implementation of the crowd movement model, first we select a welldocumented case from literature, which was previously analyzed in [12]. The geometry used Fig. 2 Snapshots of contours of calculated crowd density (in people/m 2 ) at different times for the railway platform test case of [12] in this case is displayed in Fig. 2a. The case study consists of a railway platform with an obstacle at its center. Pedestrians enter the platform and make their way across towards two exits. The inflow of pedestrians through the entrance is uniform and given by the following equation: where 'f ' is the flux of pedestrians (in pedestrians/ms) and 't' is time (in s). This flux gives a total of 15000 pedestrians that enter this railway platform over a period of 120 s. It is assumed that the maximum walking speed of pedestrians is V max = 2 m/s and that they come to a complete standstill when the pedestrian density reaches ρ max = 10 people/m 2 . The computational domain is represented by a numerical mesh of 200 × 100 regular rectilinear cells (in the xand y-direction, respectively) which gives a grid spacing of h = 0.5 m. The characteristic time-step of t = 0.01 s is used. The spatial and time resolutions are selected to match values of [12]. Contours of the crowd density distribution at different time instants are shown in Fig. 2b-e. At the initial time instant, t = 0 s, the  [12] railway platform is empty and the maximum of the incoming flux is reached at t = 60 s, and some initial signs of a local increase of the crowd density at the corners of the central obstacle can be observed, Fig. 2b. At t = 90 s, the jamming (i.e. the local concentration of crowd density reaches its maximum value of ρ max = 10 people/m 2 ) is observed around the front corners of the central obstacle. The process of jamming is further intensified at t = 120 s, where the jamming front extends over the entire upper and lower passages around the central obstacle, as well as around the edges of two outlets. The railway platform region is empty at t = 240 s. To compare the current simulations with results of [12], we extracted two characteristic profiles (fixed at x = 36 m and y = 33 m, respectively) at the particular time instant of t = 120 s, Fig. 3. it can be seen that a good agreement between our results and simulations of [12] are obtained at both locations. We also performed both temporal and spatial grid independency studies. The time step variation in the 0.0025 ≤ t ≤ 0.02 s range, produced negligible differences for profile distributions shown in Fig. 3b and c. Similarly, the further grid refinement of 400 × 200 (h = 0.25 m) or coarsening of 100 × 50 (h = 1 m), proved that results obtained with 200 × 100 (i.e. with h = 0.5 m) are grid independent, which is also in good agreement with the grid-dependency studies of [12].
To further illustrate the capabilities of developed algorithm, we analyze the more complex situations characterized by multi-directional outlets, the non-orthogonal boundaries of the simulation domain, as well as the presence of multiple obstacles -all these features are embedded in the analysis of the crowd movement predictions for the Rotterdam central railway station, Fig. 4a. The simulated geometry is marked by the dashed rectangle in Fig. 4a, and includes the gate system used to check entry tickets to the railway platforms. Here, we apply the V max = 2 m/s and ρ max = 6 people/m 2 . We analyze a few different scenarios: (i) the first case is with simplified straight walls and with 15 fully open gates with a typical width of 1 m, (ii) the second case includes the slanted walls and with 21 gates that fit across the corridor and integrated spacing of 1 m, with all gates open, and finally, (iii) the third case which is geometrically identical to the second case but where only each 2nd gate is open. All geometries have an elevator shaft with dimensions of 3 × 3 m 2 located 10 m behind the gates, which is also treated as an obstacle. The people crowd is initially located at 3 m distance from the elevator shaft and occupies an area of 20 × 23 m 2 with a uniform crowd density of 3 people/m 2 , containing a total of 1380 people. The time-dependency of people crowd size for three analyzed scenarios is shown in Fig. 4b. It is interesting to see that the inclination of the surrounding walls compared to the original straight configuration ('simple') results in almost identical evacuation time despite the fact that each second gate is closed. The results also indicate that for the case of slanted walls and all gates open, the evacuation time is approximately 62 sec, whereas for other scenarios at this particular moment in time, more than 275 persons are still not evacuated. The movement of the crowd is shown in Fig. 4c-h where the contours of crowd density are plotted at different time instants. Here we show only results for the first ('simple' geometry and all gates open) and third (slanted walls with each 2nd gate open) case. The slanted walls speed up the initial evacuation progress, as the simplified geometry is slower at the start. The closed gates appear to quickly become a limiting factor, and it can be seen that at t = 30 s, significant jamming is generated, Fig. 4e and f. However, as the rate of evacuation for the 'closed' geometry case is lower, its crowd size is surpassed by the 'simple' geometry just shortly before the evacuation is completed at t = 65 s, as seen in Fig. 4b.
We conclude that the crowd evacuation algorithm applied here showed a good agreement with similar studies in literature and the potentials of its application are demonstrated for real-life cases.

Validation of heavy gas turbulent dispersion model
For the validation of the heavy gas turbulent dispersion model, we have selected field experiments based on the Thorney Island Trials studies, [25][26][27]. The Thorney Island Trials experiments addressed the temporal evolution of heavy gas cloud under various wind conditions. The geometrical setup of these experiments contains the cylindrically shaped tent with a diameter of D = 7 m and height of H = 13 m, inside of which a heavy gas (31.6% Freon-12 and 68.4% Nitrogen mixture) is contained, Fig. 5. The ratio of this heavy gas and air densities is ρ gas /ρ air = 2. The second object is a solid cubical element with dimensions L × L × L = 9 3 m 3 , where a few measuring concentration sensors are placed. In the present investigation, we focus on the so-called 'Number #26 Trial case', which is characterized by the upwind location of the heavy gas release at the distance of 50 m from the front side of the solid cubicle. The experiments were performed in such a way that the quasi-steady conditions were reached initially for a given direction and intensity, after which the tent was removed, causing the instant release of the heavy gas. Note that due to a relatively high ratio of densities of the gas mixture and air, locally strong concentration buoyancy effects will take place, which will also affect the local velocity distributions, resulting in a twoway coupling between air flow and heavy gas concentration (concentration of the heavy gas is an active scalar). The entire simulated domain covers a volume of approximately 300×250×55 m 3 , in the x-, yand z-direction, respectively. At the inlet, approaching wind conditions were specified such that the streamwise velocity component follows the function U x (z) = U 0 (z/z 0 ) λ , where U 0 is the wind speed at the height of z 0 = 10 m, and the 'λ' is the non-dimensional parameter, both specified in Table 3. The inlet turbulence parameters  [25]. The red cylinder represents the tent where the heavy gas was initially confined, and the blue block indicates the solid obstacle where the measuring concentration sensors were located were specified as k(z) = (I t U x ) 2 and ε(z) = C μ k(z) 2 /(ν t /ν)i, where I t = 1% is the intensity of turbulence, and the ratio between turbulent and molecular viscosity (ν t /ν) = 10 3 . At the outlet boundary, the zero gradient boundary condition was applied. Along the lateral as well as the top boundary, the symmetry boundary conditions were applied. The specification of such boundary conditions proved to be sufficient since during the total simulated  time of the turbulent dispersion, the local concentration field was still away from the outlet and lateral boundaries. Similar types of the boundary conditions were also applied in [28][29][30]. The ground surface with the typical roughness of z 0 = 0.005 m was imposed throughout the extended wall-function for velocity as U + = 1/κ ln min (E, E R /e) z + , with e = U τ z 0 /ν, U τ = C 1/4 μ k 1/2 p , and coefficients E = 8.432, E R = 30, [14]. To mimic the transition between the fully developed flow conditions before and after the tent collapse, we applied an approach similar to the treatment of flow around the vegetation. In the initial stage, when the tent is the non-permeable object, the cylinder was represented as 40 × 40 prismatic elements with very low permeability (both the drag coefficient and porosity of the elements were specified to an order of magnitude higher values than normally used for very dense vegetation canopy) and the quasi-steady flow conditions were obtained. 1 An illustration of the fully-developed flow and turbulence conditions is shown in Fig. 6. To mimic the moment of the tent collapse and instantaneous release of heavy gas, at the initial time step 1 Additional term in the momentum equation to account for the presence of a porous structure is F D i = −C d a|U |U i , where C d is the drag coefficient and a is the leaf area density. Similarly, additional term in k and ε equations are S D k = C d a β p |u| 3 − β d |U |k and S D ε = C d a C ε4 β p |U | 3 ε/k − c ε5 β d |u|ε , respectively. The original model coefficients are given in [20]. of simulation the drag coefficient in the porous vegetation terms was changed to C d = 0, providing the entrance of the air flow to the initially blocked tent region and producing the turbulent dispersion of the heavy gas. To represent the solid cubicle, we applied the passive element approach, where the walls of the cubicle are treated with the extended wall functions, [14,18,20]. The numerical domain was represented with a structured locally refined (in the proximity of the ground and walls of the tent and solid obstacle) numerical mesh with approximately 1.65×10 6 CVs, which details are provided in Table 4. To check the grid dependency, a finer mesh with a total of approximately 6.2 × 10 6 CVs, was generated and results were compared to the coarser mesh, as shown in Fig. 7. It can be seen that relatively small differences were obtained (since the QUICK scheme was used to represent the convective terms of the momentum equation), confirming that the coarser numerical mesh can also be used in order to reduce the total computational costs. The time evolution of the heavy gas turbulent dispersion process is shown in Fig. 8. Here, we plotted the isosurface of the characteristic non-dimensional concentration C/C 0 = 0.02 at different times, including the initial state and time instants after tent removal. These distributions indicate a rapid collapse Fig. 9 a A sketch of the approaching wind condition and the front location of the measuring sensor at the cubical obstacle for the 'Thorney Island Trial Case Number #26'. b The details of the front measuring sensor which is placed at the height of 6.4 m and 4.5 m distance from the obstacle edge. c Time evolution of the non-dimensional concentration (C/C 0 ) at the front measuring location shown above: comparison with experiments of [25,26] and RANS simulations of [30] of the heavy gas column and its spreading in all directions. Note that concentration cloud also propagates against the wind, which is especially visible in the initial phase of the heavy gas column collapse, Fig. 8b and c. At later time instants, the local concentration of the heavy gas is reduced due to a rapid turbulent dispersion, making the concentration buoyancy mechanism weaker, and convective and turbulent diffusion mechanisms of concentration transport take over, generating strong lateral and downwind mixing, Fig. 8d.
The time-response of the concentration sensor located at the front side of the cubicle is shown in Fig. 9. Here we compare the present results with experiments of [25,26] and the earlier RANS study of [30]. The exact location of the measuring sensor is indicated in Fig. 9a and b, i.e. with its height of 6.4 m from the ground and 4.5 m distance from the side edge of the cubicle. It can be seen that a overall qualitatively good agreement is obtained between the present RANS and experiments. This is in contrast with previous RANS results of [30], which were not able to properly capture the general behavior of measurements. The present model correctly predicts the characteristic behavior of a rapid increase in the value of the non-dimensional concentration (C/C 0 ) by characteristic maximum peak, followed with a rapid decay and appearance of the secondary peak, Fig. 9c. We also analyze the influence of the concentration-buoyancy terms in momentum and equations of turbulence parameters (i.e. with G C k terms, Eqs. 10 and 11). It can be seen that the inclusion of the concentration buoyancy terms in turbulence parameters improved agreement with experiments, i.e. the first peak value is triggered earlier and the secondary peak shows a closer agreement with the measured value. The quality of agreement is also reasonable for the second measurement location, placed at the back side of the cubicle, now much closer to the ground at a height of 0.4 m, as illustrated in Fig. 10a and b. For this location, the inclusion of G C k terms also improved agreement with measurements, especially for the time instants above 120 s, Fig. 10c.
The physics of flow responsible for the appearance of this sudden concentration peak is easy to understand from characteristic spatial and temporal evolution of the heavy gas cloud illustrated in Fig. 11. The first peak is generated by the impact of heavy gas cloud on the cubicle (splash like mechanism) and a sudden increase of the vertical momentum, as shown in Fig. 11c. This is followed with a rapid decay in the local concentration due to the action of gravity, Fig. 11d.
Finally, from the observed behavior of the CFD results for a heavy gas turbulent dispersion, we conclude that the present model (with additional concentration buoyancy terms and   . 13 The contours of crowd density distribution (ρ crowd ) at different times without gas coupling for t = 0, 36, 125 and 180 s. The crowd density is displayed in color contours matching the scale to the side. The concentration of the heavy gas is displayed as a see-through gray-scale with black contour lines at normalized concentrations C/C 0 = 0.01, 0.02, 0.03 and 0.05, respectively. Simulations are performed with V max = 2 m/s GGDH representation of the turbulent concentration fluxes), was able to reasonable predict available field scale concentration measurements at various locations. This is in contrast to some earlier RANS studies reported in the literature (for the identical test case), which failed to show a good agreement with experiments.

Integrated approach
As the last step, we are going to integrate the crowd behavior dynamics with an incidental release of heavy gas. A sketch illustrating different stages of the crowd movement/gas release interaction is shown in Fig. 12. The initial crowd distribution (indicated by block rectangle) is directed toward the gas cloud, Fig. 12a. Upon the first contact with the critical gas threshold (C crit ), a part of the crowd is incapacitated (indicated by red block), Fig. 12b, causing the rest of the group to split into two parts trying to avoid regions with critical gas levels, Fig. 12c. The entire process of the crowd/heavy gas interaction is described by Eqs. 14-16, mentioned in Section 2.3. Here, the main input information for the integrated approach is the local concentration of a heavy gas at a particular height. We select the previously analyzed 'Thorney Island Trial #26' case as a basis for testing various coupling mechanisms. The results from the CFD simulation at the height of 1 m above the ground are projected onto the uniformly distributed two-dimensional mesh with a typical cell size of 0.25 × 0.25 m 2 . For the crowd behavior analysis, we defined a scenario where outlets are placed at the inlet plane (entire length) and at the front side of the cubicle with an entrance We start our analysis by presenting results of the crowd density distribution for a case without coupling, Fig. 13, i.e. the presence of the heavy gas does not affect the behavior of the crowd, similarly to the cases presented in Section 4.1. Despite the fact that there is no coupling, we also show a background concentration of the heavy gas (indicated in semitransparent gray-scale levels with black isolines). It can be seen that already at t = 36 s, jamming occurs around the exit at the front side of the cubicle, Fig. 13b. At t = 125 s, a crowd is split into two major groups: the first group still trying to escape through the front side of the cubicle, whereas the second group is directed toward the large outlet at the west side of the domain, Fig. 13c. It can be seen that the second part almost entirely left the domain at the t = 180 s, whereas the first part of the crowd is still jammed around the first exit, Fig. 13c. We also monitor the time-evolution of the total number of people still present in the simulation domain, Fig. 15a. It can be seen that the total evacuation time for the case without coupling is approximately t = 250 s.
Next, we introduce various scenarios of coupling between the background concentration of the heavy gas and corresponding dynamics of the crowd behavior. The first scenario includes testing of different ratios of the perceived and regular cost function, i.e. parameter α in Eq. 14 is in 0 ≤ α ≤ 250 range, whereas corresponding ratio values are shown in Table 5. Note that the selected referent concentration values (C ref ) are arbitrarily selected. The results of changing parameter α are shown in Fig. 15a. The change of parameter α can be interpreted as a ratio between a tendency of the crowd to avoid the gas cloud (perceived The last column shows the ratio between the two terms in the cost function (14) for typical gas concentration (C gas = 0.02) with no crowd (ρ = 0), i.e. c 2 /c 1 = αC gas / (1/V (ρ crowd )) cost function) and to reach the exit (regular cost function). It takes between 30 and 50 s for the very first part of the crowd starts to leave the domain, as indicated by the first deviation from the straight initial level of 15000 persons. Note that a crowd leaves the domain initially at an identical rate despite the different values of parameter α. This is due to the tendency of the crowd to try to escape through the nearest exit, i.e. the first small exit through the front side of the cubicle, as shown in Fig. 13b. Because of jamming, after some time a large part of the crowd starts to move through the now mostly dispersed gas (consequently reduced perceived costs) cloud to reach the large secondary exit at the west side of the domain. Upon reaching this exit, a large drop in the crowd size is observed between 180 and 240 s, depending on the value of parameter α, Fig. 15a. For a relatively small value of α = 1, crowd behavior is still dominated by a mechanism of finding the shortest escape route and the effect of the heavy gas concentration is negligible. Larger values of the parameter α ≥ 25 show significant changes in evacuation behavior, including the part of the crowd taking shelter behind the cubicle and moving around the heavy gas cloud, both in an effort to avoid a high local concentration of gas -through the significantly increased perceived cost functions, resulting in 40% longer evacuation time. In the second scenario of coupling, we test the effects of model parameters that directly influence the crowd velocity, i.e. the parameter f , defined in Eq. 16. The tested combinations are given in Table 6. Here, we assume that the reference value of the non-dimensional concentration is C ref = 0.02. Results are shown in Fig. 15b. It can be seen that the overall tendency of increasing parameter f is, as expected, in extending evacuation time. The time of reaching the first exit at the front side of the cubicle is longer than 100 s, and is a direct consequence of the reduced speed of movement. By reducing the reference concentration threshold (C ref = 0.01), the sensitivity of the crowd is significantly enhanced, resulting in a further increase of evacuation time. In all the above-presented coupling scenarios, we assumed a 100% survival ratio. In particular cases with a toxic heavy gas leakage, we also define a critical gas threshold (C crit ) at which the individuals in the crowd will start to lose consciousness and collapse, making them unable to leave the simulated domain, later to be treated as potential casualties. We will gradually increase the level of complexity involved in coupling with the critical concentration threshold, with previously analyzed increased perceived costs and movement reduction, as listed in Table 7. In the present study, we assume that the critical non-dimensional concentration threshold is C crit = 0.05. Results of the coupling involving a critical gas threshold are shown in Fig. 15c. It can be seen that active coupling with critical threshold (C-CP-CM-CPM) increases the evacuation time. In the case of critical threshold coupling (C) only, due to casualties collapsing at the front of the crowd when the critical threshold is reached, the rest of the crowd needs to go around them increasing the walking/running path and, consequently, total evacuation time. Similar behavior is also observed when both critical threshold and perceived cost coupling (CP) is performed. In the case of critical threshold and movement reduction coupling (CM), the crowd is initially slowed down and people take shelter behind the cubicle, which prevents some early casualties. For the final case, when the entire set of interactions is taken into account, i.e. (CPM) coupling, snapshots of the crowd density contours are shown in Fig. 14. It can be seen that different distributions are obtained at particular time-instants in comparison with the case without coupling, as previously already discussed in Fig. 13. Already at t = 36 s, the crowd starts to split into three major parts, because part of the crowd already starts to feel the presence of the critical concentration level. Next, in order to avoid contact with the gas cloud, the crowd moves farther away along the upper and lower boundaries of the simulated domain, Fig. 14c. With time, the local concentration of gas decreases, making the part of crowd previously jammed around the front cubicle entrance also to move straight to the large outlet at the west side, as indicated by contours between 50 and 100 m in the x-direction, Fig. 14d. The inclusion of the concept of the critical gas concentration also provides information regarding the number of incapacitated persons as shown in red lines in Fig. 15c. The simplest model, where the only critical concentration threshold is used (C), predicts the highest number of casualties (about 960). The models which include separate effects of the perceived (CP) and movement (CM) cost function, predict approximately 650 and 705 casualties, respectively. The full model (CPM), reduces this number significantly -to approximately 130 casualties. This is a result of combined effects of the reduced movement speed, which gives more time to respond to the critical concentration. This is caused by a reduction of the speed of movement and an enhanced tendency of the crowd to avoid contact with regions characterized with even a slight increase in background concentration. Consequently, the evacuation time is longest for the fully coupled model (CPM). It is interesting to observe that the majority of casualties occurred in a rather short interval (approximately 30 s), after which the crowd starts to avoid the critical regions. Lastly, we also analyze the effects of the initial crowd size on the evacuation time, Fig. 15d (with the parameters provided in Table 8). Based on different values of initial crowd density (0.1 ≤ ρ 0 ≤ 8 people/m 2 ), crowd size varies between 300 and  24000 persons, respectively. Note that in contrast to the social force based model of crowd dynamics, the macroscopic evacuation model presented here is numerically equally efficient for small and large crowd sizes. It can be seen that the total evacuation time increases with crowd size, although the differences obtained for larger initial crowd densities are rather small. Simulations reveal a change from monotonic (for small crowd size) to non-monotonic behavior (for large crowd size) of the evacuation times. For smaller crowd sizes (ρ 0 ≤ 0.5 people/m 2 ), the entire crowd still evacuates through the small exit at the front size of the cubicle, whereas for the larger crowd sizes (ρ 0 ≥ 2.5 people/m 2 ), due to jamming, a larger part of the crowd will move towards the wider exit. A sudden drop in distribution indicates the moment that this part of the crowd starts to cross the large exit.

Conclusions
We have presented a detailed mathematical framework behind the newly developed integral numerical approach that mimics the macroscopic response behavior of a crowd of people when exposed to an incidental release of heavy gas. We have shown that a relatively simple macroscopic evacuation model can produce qualitatively sound results which were in good agreement with similar studies in the literature. The crowd evacuation model we have developed has demonstrated great flexibility and numerical robustness in simulating complex floor configurations which include multiple walls (both straight and slanted) and outlets. We have demonstrated that CFD based on a time dependent two-equation eddy viscosity RANS approach can serve as a simple and reasonable accurate approach in simulations of a heavy gas turbulent dispersion. Here we also stressed the importance of the additional concentration buoyancy terms in both momentum and turbulence parameters equations, as well as the use of the GGDH in modeling of turbulent concentration fluxes. Finally, we presented some originally developed coupling schemes between the evacuation model and CFD, in which we applied various forms of the redefined cost functions and crowd movement reduction, which were based on the local referent (or critical) values of the heavy gas concentration. The presented results revealed the importance of various contributions to the redefined costs functions. The inclusion of the critical heavy gas threshold (which marks a toxic concentration level of the released gas, and consequently, an estimate of potential casualties) additionally improved the predictive capabilities of the fully integrated approach. Integrated simulations of the presented scenario of a hypothetic incidental release of a heavy gas in environmental conditions and its impact on the dynamics of the crowd behavior showed many features which are expected in realistic situations. We conclude that the approach presented here can serve as a solid foundation for further developments of fully coupled crowd dynamics/heavy gas release interactions.

Compliance with Ethical Standards
Conflict of interests The authors declare that they have no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.