An improved version of the consequence analysis model for chemical emergencies, ESCAPE

Abstract We present a refined version of a mathematical model called ESCAPE, “Expert System for Consequence Analysis and Preparing for Emergencies”. The model has been designed for evaluating the releases of toxic and flammable gases into the atmosphere, their atmospheric dispersion and the effects on humans and the environment. We describe (i) the mathematical treatments of this model, (ii) a verification and evaluation of the model against selected experimental field data, and (iii) a new operational implementation of the model. The new mathematical treatments include state-of-the-art atmospheric vertical profiles and new submodels for dense gas and passive atmospheric dispersion. The model performance was first successfully verified using the data of the Thorney Island campaign, and then evaluated against the Desert Tortoise campaign. For the latter campaign, the geometric mean bias was 1.72 (this corresponds to an underprediction of approximately 70%) and 0.71 (overprediction of approximately 30%) for the concentration and the plume half-width, respectively. The geometric variance was


Introduction
Considerable quantities of toxic industrial chemicals (TIC's) are produced and stored in the chemical and process industries and transported by railway, road or sea (e.g., Cunha et al., 2015). These can be stored either in liquid form under pressure at near ambient temperatures, or liquefied by cooling the chemical to or below its boiling point. Accidental releases of TIC's may cause significant hazard to both human health and the environment. Reliable methods are therefore needed for estimating the release rates and the atmospheric dispersion of TIC's following, for example, a process failure, accident or terrorist attack (e.g., Britter et al., 2011).
The initial cloud of a TIC release is often observed to behave as a dense gas (i.e., it is heavier than the ambient air). A denser-than-air cloud can disperse substantially differently to those observed for trace contaminants or neutrally buoyant gases. The dispersion of so-called heavy gas clouds has been studied fairly extensively over the last four decades. Reviews of the heavy gas dispersion models and experiments have been written by, for example, Koopman and Ermak (2007) and Hanna et al. (2008).
Heavy gas dispersion (HGD) models can be distinguished into four classes based on the details of physical processes modelled: (1) workbook, or phenomenological, models, (2) integral models, (3) shallow layer models and (4) computational fluid dynamics (CFD) models. Koopman and Ermak (2007) divided the models into Navier-Stokes models, modified Lagrangian puff or particle models, similarity profile models, and modified Gaussian or box models. Majority of the HGD models can be classified as either integrals or CFD models. In particular, integral models have proved to be successful hazard analysis tools, performing well within their design scope. Integral models describe the evolution of one-dimensional overall integral (bulk) properties of a cloud with time or distance. However, the effects of non-flat or obstructed terrain cannot be readily modelled, although attempts have been made to account for simple topographic effects in integral dense gas dispersion models (e.g., Kunsch and Webber, 2000).
The treatments of various processes of source term evolution and HGD have been presented in the literature fairly extensively. However, the literature is sparse regarding the development of modelling systems that are based on the use of web browser technologies. In particular, the extraction of reliable and representative meteorological input data for such modelling systems can be challenging, especially in case of non-expert users. There are very few modelling systems for accidental releases that can operationally use real-time and forecasted meteorological data. However, the meteorological input data for the WebCHARM model (Complex Hazardous Air Release Model;CHARM, 2016;Eltgroth, 2016) can be obtained from instruments on the site, from a meteorological station or from estimates prepared by a weather center.
First, this study aims to describe how an integral consequence analysis model, ESCAPE (Expert System for Consequence Analysis and Preparing for Emergencies), has recently been refined to include improved mathematical descriptions of the relevant processes. New mathematical treatments have been included for atmospheric boundary layer scaling, and new submodels have been introduced for dense gas and passive atmospheric dispersion. Second, we address the verification of the model (i.e., its fidelity to model equations and the accuracy of the numerical solutions) against the Thorney Island field trials, and its evaluation against independent experimental data measured at the Desert Tortoise field campaigns. Third, we describe a new computer implementation of the model that can readily be used over the internet on laptop computers, tablets and mobile phones. In particular, the relevant meteorological input variables can be automatically provided operationally in real time for the ESCAPE model, as predicted or forecasted by the HIRLAM numerical weather prediction model.

Overview of the ESCAPE modelling system
The ESCAPE model can be used for evaluating releases, source terms, atmospheric dispersion and consequences of hazardous materials. The model is applicable to evaluate both steady continuous and instantaneous ground-level releases of toxic and flammable gases into the atmosphere. The model can be utilized in release cases, which are a consequence of a sudden rupture of a container or the rupture of a pipe or container wall, either from the liquefied or gaseous state of a refrigerated or pressurized container. For previous descriptions of the modelling system, the reader is referred to Kukkonen (1990) and Riikonen et al. (1999Riikonen et al. ( , 2002. The model can also be used for estimating the effects of blast wave overpressures resulting from a boiling liquid expanding vapour explosion (BLEVE; Riikonen et al., 2002). The model allows for an aerosol or gaseous plume or puff to be formed (e.g., Kukkonen et al., 1994). In certain conditions, if the release scenario results in the formation of a liquid pool, its evaporation may subsequently produce either a gas puff or a plume. The current version of ESCAPE includes a library of property data of 36 commonly used and transported hazardous chemicals; the user of the program may add new chemicals by specifying their required physico-chemical properties .
The treatments for source terms and atmospheric dispersion can be classified as integral models. Integral models are characteristically formulated in terms of differential equations for the bulk properties and processes within the dispersing puff or plume. However, the effects of non-flat or obstructed terrain cannot be readily modelled. The released material is also assumed to be chemically inert and consisting of only one component.
An overview of the modelling system has been presented in Fig. 1. The system contains also semi-empirical criteria for the automatic selection of various source term and dispersion models, depending on the type of the release and its further predicted evolution. These criteria determine whether the system will use, for instance, a two-phase jet model or a liquid pool model for a particular release scenario. The modelling principles applied at all the transitions from a sub-model to another are: (a) the contaminant mass (flux) is conserved, (b) the maximum concentrations of the contaminant are continuous and (c) the physical dimensions of the puff or plume do not exhibit any apparent discontinuities. However, it is not always possible that all of the above mentioned conditions would be simultaneously valid. For such cases, the conservation of contaminant mass is considered to be the most important criterion, and the requirement of continuous maximum concentrations is the second most important one.
Two-phase or liquid-only releases through a pipe or a breach in a container wall is estimated with the two-phase pipe flow model TPDIS (Two-Phase DIScharge of liquefied gases through a pipe; Kukkonen, 1990). It is originally based on the model of Nyren and Winter (1987). Outflow through pipes and orifices from the vapour phase of a vessel is estimated with the GASFLOW model . Source terms for catastrophic failures of vessels (instantaneous releases) are evaluated by the model MIST (Model for an Instantaneously released Source Term). The model for liquid pool vaporisation POOL (Riikonen et al., 1999) is based in part on the GASP model (Webber, 1990). The two-phase jet model FLJET (FLashing JET) has been described by Kukkonen (1990). Atmospheric heavy gas and passive dispersion is estimated with the model CRISTA (Contaminant Releases Involving Spreading and Transport within Atmosphere). It is a modified version of the model DRIFT (Webber et al., 1992a;1992b).
The model described above can automatically use the real-time predictions and forecasts of the numerical weather prediction model HIRLAM, after their pre-processing. In the modelling of dispersion, the atmospheric conditions, such as wind velocity and air temperature, are assumed to remain spatially and temporally constant in the course of each hour. The dispersion modelling is currently not applicable for time-dependent or transient releases. However, updated meteorological parameters can be derived from the HIRLAM model also for historical reanalyses or forecasting in time.
The modelling system that uses the real-time meteorological data of HIRLAM has recently (in 2016) been taken to operational use, and it is described for the first time in this article. The predicted results can be post-processed using geographic information systems.

Refinement of the mathematical structure of the ESCAPE model
New integral models for the atmospheric dispersion of dense and passive gas plumes and puffs have recently been developed and included to the ESCAPE modelling system. These new submodels include the latest state-of-the-art knowledge. All of these treatments, for plumes and puffs, and for dense gas and passive dispersion regimes, have been derived using harmonised scientific principles. The new dispersion models are variations of the DRIFT model (Dense Releases Involving Flammables or Toxics) that was originally developed by Webber et al. (1992a and1992b, see also Tickle and Carlisle, 2008).
However, the treatments of the vertical profiles of wind speed, temperature, pressure and density in the atmosphere used in the model are different from those in the DRIFT model; these have been presented in detail by Kukkonen et al. (2014). The most recent knowledge of the atmospheric boundary layer scaling has been taken into account in the treatments of the profiles used in this study. The thermodynamic description in the model CRISTA is also different from that in the DRIFT model. There is also a difference in the treatment of the gravitational spreading between the two models. Finally, the computer implementation and the numerical procedures for solving the set of differential equations are different in these two models.
The source term models used within the ESCAPE system have been presented in Appendix A in the supplementary materials; however, these models are not discussed in detail in this article. A summary of the main treatments of the atmospheric dispersion are presented in the following.

Atmospheric dispersion of continuous plumes and instantaneously released puffs
For a dense gas cloud, vertical turbulence is damped due to stable density stratification within the cloud, reducing the vertical mixing with ambient air. Another effect is the horizontal gravity spreading that is generated due to density gradients in the horizontal direction. Thermodynamic phenomena are also often more pronounced for TIC releases, compared to releases and dispersion of nonaccidentally released contaminants.
Within the updated ESCAPE modelling system, the treatments for estimating the dispersion of instantaneously released puffs and steady continuous plumes near the ground level are mostly based on, but not identical compared with those in the model DRIFT. In addition to the different boundary layer scaling scheme described above, the present model is different compared with the DRIFT model regarding the assumptions of the contaminant liquid evaporation. In the treatments of the ESCAPE model, the contaminant is treated as a mixture of vapour and liquid in the source term dispersion regimes. However, it has been assumed that all the contaminant liquid has been evaporated after the source term evolution. This assumption simplifies substantially the model structure; however, it has commonly a negligible effect on the concentration predictions.
The model comprises a calculation of the evolution of integral scales for the size, density, temperature, composition and advection of the cloud with downstream distance. The cloud evolves smoothly in all respects from a sharp-edged heavier-than-air cloud to a diffuse one which ultimately behaves passively. The model equations describing the integral scales of the cloud can be derived from a fairly general profile, which comprehends the sharp edged (top hat) and Gaussian profiles. This method allows the writing of the equations for integral quantities all the way into the passive regime, without any discontinuities. Variations of the above approach have been applied in many current HGD models, for example, in the model HEGADAS (Witlox, 1994a;1994b), which is based on Colenbrander (1980) and te Riele (1977), and in the SLAB model (Ermak, 1990).
A definition is therefore needed for quantities such as the amount of air considered to be in the cloud, when the cloud has not sharply defined edges. The required equations that allow for the modelling of the cloud profiles are presented in Appendix B in the supplementary material.
The present sub-model for atmospheric dispersion is called 'Contaminant Releases Involving Spreading and Transport within Atmosphere' (CRISTA). The effects of the humidity of ambient air and the deposition of contaminant vapour have not been taken into account. In the following, we address first the atmospheric dispersion of continuous plumes, and thereafter instantaneously released puffs.
3.1.1. Continuous plumes 3.1.1.1. Concentration profiles. The time-averaged concentration field c for a continuous, steady plume is assumed to be of the form: (1) where x, y and z are the Cartesian spatial coordinates with x in the direction of the mean wind and z in the vertical direction, C m (x) is the mean ground-level, centre line concentration, F v (x,z) and F h (x,y) are the vertical and horizontal profiles, respectively, and a, b, s, and w are functions of downwind distance x. The vertical and horizontal profiles satisfy F v ðx; 0Þ ¼ F h ðx; 0Þ ¼ 1.
The functions a and b are measures of the plume dimensions, determining the vertical and crosswind extents, respectively. The functions s and w define the sharpness of the plume edge. The distribution is Heaviside ("top-hat"), as s and w tend to infinity in the dense gas limit, and approaches Gaussian forms with the value 2, or even smoother for values less than 2. The parameter s is determined from atmospheric profiles of vertical eddy diffusivity, while w depends on density effects.
3.1.1.2. Conservation of contaminant mass flux. Assuming no deposition or chemical reactions, the conservation of contaminant mass flux through a cross section of the plume on a plane perpendicular to wind direction q g for a steady release is simply where the flux is defined by q g ¼ R CU dA, where C is the concentration of the plume, U is the velocity of the plume and A is the cross-sectional area of the plume.
3.1.1.3. Lateral spreading. The lateral spread is defined by an equation reflecting gravitational slumping, while the plume is denser than surrounding air. For a dilute, neutrally buoyant plume spreading is determined by passive (atmospheric) diffusion. In order to ensure continuity between gravity dominated and passive dispersion, the lateral spreading is modelled in terms of both processes. The change of the effective integral half-width of the plume W can be written in the form where b gravity and b passive are the gravitational and passive lateral spreading rates. The plume half-width is defined in terms of the lateral concentration profile of the plume, where G is the gamma function. The gravitational and passive lateral spreading rates are given by where U f is the gravity spreading velocity, U is the effective (integral) plume velocity and u E,passive is the edge entrainment velocity in the passive regime. The form of gravity spreading (b gravity ) is different from that in the model of Webber et al. (1992a), as the spreading is assumed to be orthogonal with respect to the centre line of the plume. Webber et al. (1992a) assumed that the gravity spreading is orthogonal with respect to the edge of the plume. The difference of these two formulations is significant only for low wind speeds and close to the source. The gravity spreading velocity (U f ) is customarily modelled as where K is an experimentally determined constant, g is the acceleration due to gravity, D ¼ ðr À r a Þr À1 and D0 ¼ ðr À r a Þr À1 a are the relative integral (bulk) density differences, r is the effective (integral) density of the plume, r a is the density of ambient air, and H is the effective (integral) height of the plume (e.g., Deaves, 1987). The values evaluated by these two equations of the gravity spreding do not differ substantially, unless the plume is very dense.
The approach chosen is the latter equation (gD0H), which represents an approximate balance of the forces between gravity (density difference) and air resistance, whereas the former relation (gDH) does not allow such an interpretation (Webber et al., 1992a). However, U f may become arbitrarily large, as the density difference D0 increases. An estimate of the upper bound of U f can be found by noting that the spreading velocity is bounded by the potential energy of the plume, i.e., the potential energy is converted only to kinetic energy of plume spread. Therefore, the spreading velocity is approximated with (Webber et al., 1992a) where K f ¼ 1.05 is the densimetric Froude number, obtained by fitting the modelled plume spread with the Thorney Island experiments. The second term on the right-hand side of Eq. (7) represents a theoretical limit value, due to the conservation of energy, based on the analysis of Webber and Brighton (1986) on shallow water equations for liquid pool spread. However, the second term in U f will rarely be significant, only for a short time in the initial stages.
3.1.1.4. Air entrainment. The air entrainment rate q a per unit plume length into a ground based plume is given in terms of the volumetric top (Q ET ) and edge (Q EE ) entrainment rates, The volumetric entrainment rates per unit length of the plume can be written in terms of the top (u T ) and edge (u E ) air entrainment velocities The top entrainment velocity in both the gravity dominated and passive regions is modelled by where u T,passive is the top entrainment velocity in the passive regime, u *c is the friction velocity within the plume, u * is the atmospheric friction velocity and fðRi * Þ is a function of the bulk Richardson number accounting for the inhibition of entrainment due the density difference between the cloud and the air. u *c is estimated by rescaling the atmospheric friction velocity to account for the perturbation of a dense cloud on u * . The bulk Richardson number is defined by Ri * ¼ gDHu À2 * . The modelling of u T,passive closely follows the approach of the passive puff model of Wheatley (1988) (Webber et al., 1992a). Wheatley's model is an approximate solution of the atmospheric diffusion equation (gradient-transfer or K-theory) describing the diffusion of a passive puff released at the ground level into a turbulent stratified shear flow. As the density difference tends to zero (i.e. in the passive limit where fðRi * Þ/1 and u *c /u * ), the top entrainment velocity approaches the passive limit. The form of the top entrainment velocity is based on the analysis of the Thorney Island experiments by Wheatley et al. (1986). Forms of fðRi * Þ have also been proposed by, for example, Havens and Spicer (1985) and Puttock (1987), based on laboratory experiments. The relation applied here is very close to that of Havens and Spicer (1985).
Future development of the model should include an investigation of the performance of the top entrainment formulation adopted. For example, Briggs et al. (2001) suggests that u T ¼ 0:65u * ð1 þ 0:2Ri * Þ with Ri * up to 20, fitted to laboratory measurements of vertical diffusion of continuous, dense gas plumes over rough surfaces in neutral boundary layers. The entrainment equation presented here reduces to u T ¼ ku * (with k ¼ 0.4) for a passive plume in neutral atmospheric conditions, implying about 38% smaller entrainment velocity compared to the equation by Briggs et al. (2001) in the passive limit, as Ri * ¼ 0. A comprehensive analysis of the theory and experimental datasets related to the vertical entrainment of passive clouds dispersing in a neutral boundary layer is presented by Britter et al. (2003).
Again, in order to ensure continuity as the heavy plume becomes passive, the edge entrainment velocity is modelled as so that u E tends to the standard edge entrainment model u E ¼ a E U f in the gravity-dominated regime (w/∞) which has been discussed in detail by, for example, Wheatley et al. (1986). The edge entrainment model reduces to the passive model as w/2. The empirical value chosen for the entrainment constant (a E ¼ 0.7) optimises fits to the Thorney Island experimental data on instantaneous releases (Wheatley et al., 1986).
3.1.1.5. Plume temperature. We assume the plume to be gaseous and disperse in dry air, with no deposition. Therefore, the temperature of the plume is influenced by the enthalpy flux into the plume due to entrained dry air and the heat transfer between the plume and the ground. Heat transfer is most significant for initially cold plumes moving over a warmer substrate (ground or water). The temperature of the plume T can be solved from an enthalpy balance equation, where C q is the total heat capacity flux of the plume, c pa is the specific heat capacity of air, T a is the temperature of ambient air, Q a ¼ Q ET þ Q EE is the volumetric air entrainment rate per unit length of the plume into the plume and S H is the mean vertical turbulent heat flux per unit area. The first part of the right-hand-side of this equation is the enthalpy flux into the plume due to entrained air. The second part (S H ) describes the amount of heat transfer between the plume and the ground. The heat transfer is determined from the maximum of the heat fluxes due to forced and free convection. Free convection is assumed to occur only for clouds colder than the ground temperature.
3.1.1.6. Evolution of plume momentum. The advection model consists of specifying the momentum balance of the plume. Including the effect of momentum transfer due to entrainment, the change of the momentum flux of the plume is where q p ¼ rU 2 A≡rU 2 2WH, U AT and U AE are two velocities characteristic of the ambient flow. The advection formula may be interpreted as saying that top entrainment is bringing in material with velocity U AT while edge entrainment is bringing in material with velocity U AE . The velocities are set to where l E and l T are empirical constants and U w is the wind speed at the height of plume centroid. We have adopted the values as given in Appendix C in supplementary material. This reproduces the empirical model, which Puttock (1987) has found successfully to reproduce the cloud advection speeds in the Thorney Island instantaneous release data. Briggs et al. (2001) defines the mean transport speed of a plume (u p ) to be a concentration-weighted mean of the in-plume velocity (u) as u p ¼ R uCdz= R Cdz, where C is the concentration. However, the definition requires knowledge (measured or modelled) of the velocity profile within the plume and is therefore not readily suitable for modelling purposes. Some dense gas models describe the transport speed with the concentration-weighted average and assume u to be equal to the out-of-plume, ambient wind speed (e.g. HEGADAS; Witlox, 1994a). Obviously, the latter assumption does not allow for e.g. the initial acceleration of an initially stand still puff (for instance, the initial situation in the Thorney Island experiments). Briggs et al. (2001) presents also an alternate method for evaluating the transport speed (if there is no information inside the plume on C and u). It is an empirical relation derived from (i) a mass balance equation assuming no deposition or chemical transformation, and (ii) laboratory measurements of in-plume and out-of-plume vertical velocity profiles in neutral conditions.

Instantaneously released puffs
3.1.2.1. Concentration profiles. For the purpose of defining profiles of a puff, it is useful to define a relative coordinate x' ¼ x -X(t), where X(t) is the centroid position of the puff at time t. The steady, time-averaged concentration field C can be written as Cðx; y; z; tÞ ¼ C m ðtÞF v ðx; zÞF h ðx0; yÞ≡C m ðtÞexp À z a s exp where C m (t) is the centre, ground-level concentration at time t, F v (x,z) and F h (x',y) are the vertical and horizontal profiles, respectively, and a, b 1 , b 2 , s and w are functions of travel time t. The concentration contours are thus elliptical in any horizontal plane.
The parameter a is a length scale determining the vertical extent of the puff and b 1 and b 2 are length scales determining the puff size in x and y directions. The parameters s and w define the sharpness of the puff edge, giving a sharp edged puff and in the dense gas limit and approaches Gaussian forms in the passive limit. As for the plume model, the parameter s is determined from atmospheric profiles of vertical eddy diffusivity, while w depends on density effects.

Conservation of contaminant mass.
Assuming no deposition or chemical reactions the conservation of contaminant mass m g is simply dm g dt ¼ 0: 3.1.2.3. Lateral spreading. The lateral spread of a puff is modelled using the same principle as for a plume, ensuring continuity between gravity dominated and passive dispersion regimes. The rate of change of the integral (effective) radii of the puff in downwind (semi-major-axis) and crosswind (semi-minor axis) directions R i can be written as where U i,gravity and U i,passive are the lateral spreading rates in the gravity dominated and passive regimes, respectively. We shall refer with i ¼ 1 to the longitudinal spreading and with i ¼ 2 to the crosswind spreading. Spreading in the gravity dominated regime is assumed to be equal in crosswind and downwind directions.
Allowing also a delay time for the onset of radial acceleration, the gravity spreading velocity is modelled with where q is the Heaviside function (q(x) ¼ 0 for x < 0; q(x) ¼ 1 for x > 1), l g is an experimental constant and t G is a gravity spreading time scale given in terms of the initial state of the puff. The delay by l g t G is a very simple model of the initial radial acceleration, and is introduced to aid comparisons with experimental data. However, its influence is minor in most circumstances (Webber et al., 1992b). The passive spreading is based on the puff dispersion model of Wheatley (1988) allowing for different degrees of wind shear. Incorporating the Wheatley's model with the concentration profiles given above, leads to (Webber et al., 1992b) where G is the gamma function, the root mean square turbulence velocities in the longitudinal and crosswind directions are taken to be s u ¼ 2:5u * and s v ¼ 2:0u * , respectively, and s y is the standard deviation of the crosswind concentration distribution. The passive puff dispersion model of Wheatley (1988) allows the puff to lean over (skew) towards downwind rather than assuming a (vertically) right cylinder. The current model does not allow the effect of passive puff skewness.
3.1.2.4. Air entrainment. The air entrainment rate into the puff is given in terms of top (Q ET ) and edge (Q EE ) entrainment. The rate of change of the amount of air in the puff m a is written as where the volumetric top and edge entrainment rates are modelled with where A is the integral (effective) puff area (the horizontal ellipse section) and p A is the perimeter of the puff ellipse. The top entrainment velocity (u T ) is modelled as for the continuous plume model. For both the gravity and passive regime, the edge entrainment is set to The edge entrainment velocity reduces to u e ¼ a E U f in the dense gas limit (w large) and to the passive model of Wheatley (1988) as w/2. The edge entrainment model for gravity dominated region has been discussed by Wheatley et al. (1986), who found that a E ¼ 0.7 optimises fits to Thorney Island data.
3.1.2.5. Puff temperature. We assume the puff to be purely gaseous and disperse in dry air, with no sources or deposition. The rate of change of the puff temperature T with time is then expressed in the form (similarly to the evolution of plume temperature) where C c is the total heat capacity of the puff, Q a ¼ Q ET þ Q EE is the volumetric air entrainment rate into the puff and S H is the heat flux per unit area modelled as for the continuous plume.
3.1.2.6. Evolution of puff momentum. The downwind travel model consists of a specification of the force acting on the puff as a whole. Including the effect of momentum transfer due to entrainment, the rate of change of the momentum of the puff p c is dp c dt where p c ¼ mU and m is the total mass of the puff. The characteristic velocities (U AT , U AE ) are scales modelled as for the plume model and the volumetric entrainment rates (Q ET , Q EE ) are given in Eqs. (24) and (25).

Evaluation and verification of the model
We have compared the predictions of the refined ESCAPE model with the measured data of six well-known field campaigns: Landskrona, Sweden, 1982, 1984, Skelleftehamn, Prairie Grass, USA, 1956, Landskrona, Sweden, 1993e1994, Thorney Island, UK, 1982e1984 and Desert Tortoise, USA, 1983. The comparisons with the Thorney Island and Desert Tortoise measurements are addressed in this article. We have selected these two campaigns to be considered here, as (i) these experiments are amongst the most important datasets in use for the development and evaluation of dense gas dispersion models, and (ii) the organization of the campaigns and the results have been described with sufficient detail to facilitate an accurate model evaluation.
However, some aspects (entrainment velocities and cloud spread) of the original DRIFT model in the dense gas regime have been tuned against the Thorney Island trials (Webber et al., 1992b;Tickle and Carlisle, 2008). Comparison of the predictions of the present model with the Thorney Island data should therefore be considered primarily as a verification (of the writing of the computer program and its numerical accuracy), instead of an evaluation. Generally, there are so few available full scale hazardous gas dispersion experiments that it is in many cases not possible to find a truly independent dataset (e.g., Hanna et al., 1991a).
For readability, we first present a brief overview of these campaigns, and then the results of the model evaluation and verification are discussed.

An overview of the Thorney Island trials
The Heavy Gas Dispersion Trials at Thorney Island (U.K.) were carried out between June 1982 and June 1984. For a detailed description of the trials, the reader is referred to McQuaid and Roebuck (1985) and McQuaid (1985McQuaid ( , 1987. The main aim of the trials was to study the dispersion of fixed-volume (z2000 m 3 ), isothermal, instantaneously released dense vapour clouds. As a sequel to the main series, three continuous release trials were performed over unobstructed, flat terrain.
The trials of the main series were divided into two phases. Phase I comprised of 16 trials of instantaneous release, where cloud dispersion took place over flat and unobstructed terrain. In Phase II, 10 instantaneous release trials were organised in the presence of isolated buildings, and solid and porous fences. The release gas was a mixture of dichlorodifluoromethane (CCl 2 F 2 ; freon-12) and nitrogen (N 2 ) to simulate gases over a range of cloud densities. A 30 m tall meteorological tower with wind and temperature sensors installed at five elevations was located 150 m upwind of the release point. Measurements of gas concentration were taken at 38 masts at the elevations from 0.4 up to 14.5 m. The masts were located on a rectangular grid with distances up to approximately 800 m from the release point.
During phase I the release and meteorological conditions had the following characteristics. The initial cloud was cylindrical in shape with 14 m in diameter and at maximum approximately 13 m high. Initial relative densities compared to ambient air varied from 1.0 to 4.2. Estimated atmospheric stability ranged from moderately unstable (Pasquill B) to stable (Pasquill F). Measured wind speeds, averaged over the duration of each experiment, at the height of 10 m varied from 1.7 to 7.5 m s À1 .

An overview of the Desert Tortoise experiments
The series of Desert Tortoise field experiments were conducted in a desert at the Frenchman Flat area of the Nevada Test Site, US in 1983; a more detailed description has been presented by Goldwire et al. (1985) and Koopman et al. (1986). The measurement series was designed to document the transport and diffusion of ammonia vapour, resulting from four large scale cryogenic releases of liquid ammonia.
Meteorological instrumentation included eleven cup-and-vane anemometers at a height of 2 m at various positions within the site. Additionally, a 20 m tall meteorological tower was located upwind of the spill area, with temperature measured at four levels and wind speed and turbulence at three levels. Pressure, flow and temperature measurements were made on the release system. NH 3 concentrations and temperatures were measured on arcs of seven and five measurement towers, distributed at the distances of 100 m and 800 m downwind of the source. The towers were equipped with concentration and temperature sensors at three elevation levels, ranging from 1.0 to 8.5 m.
The wind speeds were fairly strong, with 5e8 m s À1 at the height of 3.4 m. Estimated atmospheric stability ranged from neutral (Pasquill D) to moderately stable (Pasquill E). Ambient air was hot and dry, with the ambient temperatures ranging from 29 to 34 C and the relative humidity from 13 to 21%. However, a heavy rainfall prior to the start of the experiments covered the site by a shallow layer of water that prevailed during most of the experiments.
Pressurized liquefied ammonia was released during the four trials from a spill pipe pointing horizontally downwind at a height of 0.8 m above the ground. The experiments were carried out with emission rates of 80e130 kg s À1 and the durations of 2e6 min. The liquid NH 3 flashed as it exited the pipe, resulting in about 18% of the liquid changing phase to gas (Hanna et al., 1993). The rest of the liquid (82%) broke up into an aerosol by the turbulence inside the jet. Very little, if any, of the unflashed liquid was observed to form a pool on the ground. Dispersion close to the release point was dominated by the dynamics of a turbulent two-phase jet. At later stages downwind of the jet, slumping and spreading of the cloud indicated dense-gas dynamics. In most cases, most of the plume was below the height of 6 m and within the lateral domain covered by the measurement towers.

Statistical assessment measures for model performance
For simplicity, we address here only two performance measures: the geometric mean bias (MG) and the geometric variance (VG), defined as where c o and c p are the observed and predicted concentrations, respectively; and the overbar denotes an average over the dataset. A perfect model compared against perfect observations would be placed at MG ¼ VG ¼ 1.0. Geometric mean indicates the central tendency or typical value of a set of numbers. Values of geometric mean MG > 1 imply model underprediction and those of MG < 1 overprediction. Typical magnitudes of performance measures (MG, VG) have been summarized by Chang and Hanna (2004) and Hanna et al. (2004). They concluded that well performing models appear to have the following typical characteristics: the geometric mean bias is within ±30% (0.7 < MG < 1.3) and the random scatter is smaller than from a factor of two (VG < 1.6) to three (VG < 3.3) of the predicted mean.
The logarithmic form of MG ensures equal weight for over-and underestimation, and is less sensitive to outliers. However, the value MG ¼ 1 obviously does not mean that predictions coincide with measurements. Geometric variance (VG) is a variance measure related to, or being the variance counterpart of MG. For predictions having only mean bias without any random scatter present, ln (VG) ¼ ln 2 (MG); this condition defines a minimum possible value of VG for a given MG.

Model verification for the Thorney Island experiments
Model performance measures MG and VG for the instantaneous dense gas releases at Thorney Island trials are presented in Fig. 2 for numerous models. Performance measures for the maximum ground level concentration are shown for 10 models that were evaluated previously by Hanna et al. (1991a), for the IITHG-I model by Mohan et al. (1995), for the LDM model by Lee et al. (2007), and for the latest version of the ESCAPE model, the latter evaluated in this study. The Thorney Island dataset includes a total of 61 concentration samplings from nine trials. The experimental concentration data represent averages of 0.6 s.
The results for all the models, except for the DEGADIS model, show fairly little random scatter. Most of the points in the diagram are located near the parabola that marks the minimum values of the mean geometric variance.
The results for the ESCAPE model can be interpreted only as model verification, i.e., in terms of the testing of the accuracy of the computer implementation and numerical solutions. The entrainment and cloud spread parameters in the model were determined from studies of the Thorney Island data (Webber et al., 1992a;1992b). The agreement of the predictions of ESCAPE and the measured data is excellent: there was a very slight overprediction of the observed concentration values (MG ¼ 0.92), and a small random scatter (VG ¼ 1.19). The above mentioned value of MG represents an overprediction of approximately 8%. The random scatter is clearly lower than a factor of two; a scatter of a factor of two would correspond to a value of VG of 1.6. A fraction of 87% of the concentration predictions of ESCAPE were within a factor of two of the observations. These results add more confidence that the model equations are probably correctly implemented and sufficiently accurately solved in the computer code, at least for cases that are similar to the Thorney Island trials.

Model evaluation for the Desert Tortoise experiments
The Desert Tortoise releases involved strong two-phase momentum jets. The evaluation therefore addressed both the jet and the heavy gas dispersion regimes. The model performance measures VG and MG for the Desert Tortoise experiments are presented in Fig. 3aeb. The Desert Tortoise dataset includes a total of eight concentration and plume half-width samplings from four experiments. The results are presented for 14 and eight models in Fig. 3aeb. The results for the other models except for ESCAPE have been presented previously by Hanna et al. (1991b). The half-width of the plume was defined as the standard deviation (s y ) of the horizontal crosswind concentration distribution.
The results for maximum concentrations (Fig. 3a) are different from the corresponding results for the Thorney Island experiments (Fig. 2) in one respect. The geometric mean variances were in general clearly lower for the Thorney Island experiments (except for one model). This implies that there was relatively more random scatter in either or both of the measured data and the model predictions in the Desert Tortoise experiments.
The ESCAPE model tends to somewhat underpredict the maximum centre line concentrations (Fig 3a; MG ¼ 1.72, corresponding to 72%), and slightly overpredicts the plume half-widths (Fig 3b; MG ¼ 0.71, corresponding to 29%). Predictions of the ESCAPE model for both the concentration and the plume half-width have a small geometric variance (VG < 1.5). The random scatter is within a factor of two, if VG < 1.6. The ratios of ESCAPE-predicted to observed maximum centerline concentrations (c p /c o ) range from 0.39 to 0.99. The ratios of predicted to observed plume half-width (s yp /s yo ) range from 1.2 to 1.6. In relation to the model evaluation measures of other participating models, the statistics for the ESCAPE model can be judged as fairly good.
The values of the model evaluation measures can be dependent on the averaging times of the measured concentrations. The concentration observations applied for the performance measures range for the longest available averaging time from 80 to 300 s (Hanna et al., 1991a;1991b). However, in the treatments of the ESCAPE model, the averaging time is related to the entrainment velocities adopted by the model, and it is therefore not possible to adjust these without changing the model structure. The mean concentration predictions in the heavy gas region are considered to represent averages of about 6 s, decreasing smoothly to averages of approximately 10e15 min in the passive cloud region (Webber et al., 1992a;1992b;Tickle, 2001;Webber, 2002). However, during the gravity dominated phase the fluctuations of concentrations tend to be relatively smaller compared with corresponding trace dispersion, and the results are therefore in that stage less dependent on the precise averaging time.

The internet-based version of ESCAPE
An operational, web browser based version of ESCAPE was introduced in January 2016. The program can use either real-time or forecasted meteorological data produced by the numerical weather prediction (NWP) models, the results of which are available at the Finnish Meteorological Institute. User-defined meteorological data can also be used. The results of the following NWP models are continuously available: HIRLAM (HIgh Resolution Limited Area Model) (Und en et al., 2002), HARMONIE (HIRLAM ALADIN Research on Mesoscale Operational NWP In Euro-eMediterranean Partnership) (Nielsen et al., 2014), ECMWF (European Centre for Medium-range Weather Forecasts) and LAPS (Local Analysis and Prediction System) (Hiemstra et al., 2006).
A ready-made operational system has been developed for obtaining the relevant meteorological data from the output of the NWP model HIRLAM. The HIRLAM model was selected for two main reasons. First, this model is known to provide fairly accurate weather forecasts for the whole of Europe, and its treatments have specifically been adopted for the conditions in the Northern European region. Second, the NWP computations are performed operationally in-house, which simplifies the transfer of data between the ESCAPE and HIRLAM models.
In the following, we first present a short introduction of the structure and application of the HIRLAM version used in this study, and then summarize the structure of the operational implementation of ESCAPE.

The numerical weather prediction model HIRLAM
HIRLAM is a hydrostatic numerical weather prediction model. The modelling domain includes Europe, Northern Atlantic, and part of western Asian and northern polar regions. In this study, an operational version was used, called HIRLAM RCR (Regular Cycle with the Reference). Meteorological forecasts are continuously produced four times a day, with a temporal resolution of 1 h and a forecast length of two days ahead in time. The horizontal grid spacing of the model is 0.068 (approximately 7.5 km). Vertically, the model spacing ranges from tens of metres near the ground surface to approximately 300 m at the height of 5 km (Lauri et al., 2012). The vertical grid consists of 65 levels; the lowest level is approximately at the height of 12 m.

The implementation of the operational version of ESCAPE
A new version of the ESCAPE model has been developed, based on web browser technology. This program is currently applicable on laptop computers, tablets and mobile phones. The program has been tested to be compatible with the latest versions of the most commonly used internet browsers. The main aim of developing this version was to provide a user-friendly tool of assessment for the emergency response personnel. This operational model version can be used for emergency contingency planning, for the safety analysis of industrial installations and for the training of emergency personnel.
This application uses a graphical user interface both for input data entries and for the presentation of computational results. Navigation through the model input screens has been designed to be as user-friendly and straightforward as possible. The user only needs to specify input data regarding the characteristics of the release and the environment (the type of release and its associated details, the location and time of the accident, and the contaminant). The program also automatically checks that the input data values and their combinations are not unphysical, and that these are within the applicability range of the model. The meteorological input data can be extracted and preprocessed automatically for the spatially and temporally most relevant coordinates (e.g., for the geographical location and time of the accident). Generally, the predicted meteorological values are also more representative for the accident location and time, compared with, e.g., local measurements at a station in an industrial complex. The on-line use of the products of a NWP model also simplifies and speeds up the input of data for any specific accident scenario; this also minimizes potential human errors. In case of non-expert users, the determination of meteorological variables, such as, e.g., atmospheric stability, could otherwise be a challenging task. In case of long-term accidents, the user can also forecast the hazard areas for two days ahead in time.
Currently, a set of detailed digital geographical maps provided by the National Land Survey are contained within the modelling system. These maps cover the whole of Finland. As output, the model provides, e.g., spatial concentration distributions near the ground level, presented on these maps. The user can compile dedicated accident data bases consisting of input data for a range of potential accident scenarios for a selected industrial installation. Such cases can easily be retrieved and edited as necessary.

Conclusions
We have discussed the physical and chemical refinement of the consequence analysis model ESCAPE, the verification and evaluation of the model against experimental field data, and a new internet-based implementation of the model. We have recently refined the mathematical structure of the model, especially regarding the atmospheric dispersion modelling. The models for treating dense gas and passive dispersion have been replaced by new modules, both for plumes and puffs. The atmospheric dispersion modelling has been mainly based on the mathematical treatments by Webber et al. (1992a and1992b) and a boundary layer scheme presented by Kukkonen et al. (2014).
We have compared the predictions of the refined ESCAPE model with the measured data of six field campaigns; the comparisons with the Thorney Island and Desert Tortoise measurements were addressed in this article. As some of the parameters of the model have been adjusted based on the results of the Thorney Island trials, these comparisons need to be considered only as a verification of the writing of the computer program and the numerical solutions, instead of a model evaluation. In case of the Thorney Island trials, the agreement of the predictions and the measured data was good, with a very slight model overprediction (the geometric mean bias ¼ 0.92), and a small random scatter (the geometric variance ¼ 1.2). For the Desert Tortoise campaign, the geometric mean biases were 1.72 and 0.71 for the concentration and the plume half-width, respectively, and the geometric variance was smaller than 1.5. These values for the Desert Tortoise campaign indicate a good agreement of the model predictions and the measured data.
However, the actual evaluation of the model that has been presented in this article is based only on one measurement campaign. Clearly, the two campaigns addressed in this study represent only a limited range of spill scenarios and atmospheric conditions. The atmospheric stabilities for the Thorney Island experiments were reported to range from moderately unstable (Pasquill B) to stable (Pasquill F), and for the Desert Tortoise campaign from neutral (Pasquill D) to moderately stable (Pasquill E). The measured wind speeds at the considered Thorney Island experiments at the height of 10 m varied substantially (from 1.7 to 7.5 m s À1 ), and those for the Desert Tortoise campaign were commonly strong, 5e8 m s À1 at the height of 3.4 m.
These two campaigns therefore include only little information on light wind speed cases, and no information on calm conditions (defined as u < 1.0 m s À1 at the height of 10 m). In particular, the Desert Tortoise experiments include solely strong wind speed cases. The information is also scarce regarding both stable and extremely stable atmospheric conditions. In case of hazardous materials, it is commonly necessary to exclude calm and nocturnal conditions from the campaign for safety reasons. In general, the stable and low wind speed cases are the most demanding conditions for atmospheric boundary layer scaling and dispersion models (e.g., Kukkonen et al., 2012). A more extensive model evaluation needs to be performed in the future. The field experiments of Kit Fox, for carbon dioxide (Hanna and Steinberg, 2001) and those of Jack Rabbit, for ammonia and chlorine (Hanna et al., 2012(Hanna et al., , 2016 provide well-controlled data regarding dense gas dispersion. Evaluations of some models (HEGADAS, Phast, SLAB) against the Kit Fox and Jack Rabbit experiments have been presented by Hanna and Chang (2001), Hanna et al. (2004Hanna et al. ( , 2012 and Witlox et al. (2014).
A web browser version of the ESCAPE model has been developed, which is applicable on laptop computers, tablets and mobile phones. The model can automatically use the real time predictions and forecasts of the numerical weather prediction model HIRLAM. This commonly results in more reliable and accurate predictions, for several reasons. First, the accurate evaluation of the required meteorological variables (e.g., atmospheric stability) can be a challenging task for a non-expert user. Second, the automatic provision of the best available predictions of a NWP model simplifies and speeds up the evaluation of an accident scenario. Third, the predicted meteorological values tend to be more reliable and better representative for the accident location and time, compared with, for example, a local measurement at some specific point in an industrial site. Fourth, in case of long-term accidents, the user can also forecast in time the hazard areas.
The internet-based program also includes other characteristics in order to make the evaluation of accident scenarios as fluent and flexible as possible. For instance, a set of detailed digital geographical maps are contained within the modelling system, and the results can be presented and processed directly on a userselected geographical domain. The user can also compile and archive ready-made accident data bases, evaluated specifically for a selected target installation or industrial area.
The ESCAPE model is currently widely used nationally in Finland by the regional rescue services and by other expert users (in addition to its use in-house). The new model version has already proved to be a useful tool of assessment for the needs of emergency response authorities. The model could also be used throughout Europe or globally after some adjustments, such as an extension of the digital map archive, and the use of the forecasts of a global NWP model, instead of the European-wide HIRLAM model. and Antti Hellsten (FMI). Fire Protection Fund is thanked for the financial assistance (grant number: SMDno/2012/1187) of this study.