An emergency response model for evaluating the formation and dispersion of plumes originating from major fires (BUOYANT v4.20)

A mathematical model called BUOYANT has previously been developed for the evaluation of the dispersion of positively buoyant plumes originating from major warehouse fires. The model addresses the variations of the cross-plume integrated properties of a rising plume in a vertically varying atmosphere and the atmospheric dispersion after the plume rise regime. We have described in this article 10 an extension of the BUOYANT model to include a detailed treatment of the early evolution of the fire plumes, before the plume rise and atmospheric dispersion regimes. The model input and output consist of selected characteristics of forest or pool fires, and the properties of a source term for the plume rise module, respectively. The main model structure of this source term model is based on the differential equations for low-momentum releases of buoyant material, which govern the evolution of the plume 15 radius, velocity and density differences. The model is also partially based on various experimental results on fire plumes. We have evaluated the refined BUOYANT model by comparing the model predictions against the experimental field-scale data of the Prescribed Fire Combustion and Atmospheric Dynamics Research Experiment, RxCADRE. The predicted concentrations of CO2 agreed fairly well with the aircraft measurements conducted in the RxCADRE campaign. We have also compiled an operational 20 version of the model. The operational model can be used for emergency contingency planning and for the training of emergency personnel, in case of major forest and pool fires.


Introduction
The dispersing fire plumes can represent a substantial hazard on the health of people and the state of the 25 environment, in addition to the direct effects of major fires at the accident sites. Major fires include, e.g., fires in warehouses and industrial sites, and wildland fires. The latter category includes, e.g., heath, moorland and forest fires. Major wildland fires can result in substantially more extensive and intensive https://doi.org/10.5194/gmd-2021-88 Preprint. Discussion started: 25 May 2021 c Author(s) 2021. CC BY 4.0 License. evaluation of the required input data is a challenging task. Clearly, this has also up to date been the case for all other models for the dispersion of buoyant plumes from major fires. The application of such models has therefore been possible only for expert users.
The overarching aim of this study has been to develop a comprehensive model for evaluating the 95 dispersion of plumes originating from major fires, including all the relevant dispersion regimes. An objective has been to develop an operational and user-friendly model for this purpose. The objective has been that such an operational model could be used also by well-trained meteorologists and by trained emergency rescue personnel. We have therefore developed a novel semi-empirical model for evaluating the initial conditions (i.e., source term) for models that treat the buoyant plume rise and the subsequent 100 atmospheric dispersion of plumes from major fires. The input data required by the source term model is substantially simpler than the corresponding input required by the original model version.
The specific objectives of this study are the following. (i) The first objective is to present a new model for evaluating the initial properties of fire plumes in terms of fairly simple characterizations of major 105 fires. We have also included the developed source term model to the BUOYANT model. (ii) The second objective is to compare the predictions of the latest version of the BUOYANT model against available field-scale experimental data. (iii) The third objective is to present the structure and functioning of the developed operational emergency response model. Both the original research model and its operational application could in the future be used worldwide, the latter after some slight modifications, for a better 110 preparedness and rescue operations in case of major fires.

Model
For readability, we present first a brief overview of the whole modelling system. Second, we address the extension of the system to include source term evolution.

Overview of the modelling system 115
The BUOYANT model is applicable for evaluating the initial plume rise and atmospheric dispersion of pollutants originated from buoyant sources. An overview of the modelling system is presented in the following. For a more detailed description of the mathematical model, th e reader is referred to Kukkonen et al. (2014) and in case of the atmospheric dispersion after the plume rise regime, to Nikmo et al. (1999).

120
The relevant flow regimes and an overview of the currently applicable modelling system have been presented in Figs. 1a-b. The model includes treatments (i) for the initial plume properties immediately above the fire (source term), (ii) for the dispersion of the buoyant plume and (iii) for the dispersion after the plume rise regime. All of these sub-models constitute a model called BUOYANT. The model can be used to predict the spatial concentration distributions of pollutants originated from fires. 125 The BUOYANT modelling system also includes an atmospheric dispersion module, which applies the gradient transfer theory and Gaussian equations in the vertical and horizontal directions, respectively (Nikmo et al., 1999). The source strength and the atmospheric conditions are assumed to remain constant 135 in time during the atmospheric dispersion. However, after the plume rise regime, the plume properties can also be used as input information for any other dispersion model.
The plume rise predictions of the BUOYANT model have previously been evaluated against two experimental field data sets regarding prescribed wild-land fires (Kukkonen et al., 2014). These were the 140 "Smoke, Clouds and Radiation -California" experiment (SCAR-C) in Quinault in the US in 1994 (e.g. Kaufman et al., 1996;Gassó and Hegg, 1998) and an experiment in Hyytiälä in Finland in 2009(e.g., Virkkula et al., 2014a. The plume rise treatment of the BUOYANT model has also been evaluated against wind tunnel experiments (Martin et al., 1997;Kukkonen et al., 2000). The dispersion module (after the plume rise regime) has also been tested against the Kincaid field trials (Olesen, 1995;145 Nikmo et al., 1999). Sofiev et al. (2012) has compared the BUOYANT plume rise predictions against https://doi.org/10.5194/gmd-2021-88 Preprint. the results of the Multi-angle Imaging SpectroRadiometer (MISR) Plume Height Project (e.g., Kahn et al., 2008). In all of the above mentioned model evaluation studies, the model predictions have agreed well or fairly well with the observations.

150
The plume rise and dispersion treatments of the model contain three physical parameters, the values of which have been determined based on a comparison of model predictions and wind tunnel observations (Kukkonen et al., 2000). The values of these parameters have not been changed or adjusted in any way after their initial determination. The model equations therefore do not contain free parameters that could be adjusted according to the measured values. 155

The source term model
The source term model can be used to evaluate the characteristics of the plume generated by a fire, which will be needed for the subsequent computations on the evolution of a buoyant plume. The source term model includes treatments to evaluate the properties of a fire plume just above the flame tips, based on information on the characteristics of the fire. The model has been designed to be used also for operational 160 purposes; we have therefore tried to keep the amount and nature of the input data as limited and simple as possible.
The current model version can be applied for two significant categories of fires, viz. the forest and liquid pool fires. In case of forest fires, the input data of the source term model has been selected to include (i) 165 the area of the forest on fire, (ii) the number of trunks per unit area of forest, (iii) the average height of trees and (iv) the average bole diameter at breast height. The input data in case of liquid pool fires includes (i) the burning substance, and (ii) the surface area of the liquid poo l on fire. The results of the source term model include (i) the mass fluxes of gaseous and particulate matter produced by a fire, the mass flux of entrained air, and (ii) the characteristic scales of radius, temperature and vertical velocity of 170 the plume. The source term model does not include an evaluation of the propagation of the actual fire in the terrain; separate models have been developed for this purpose in other studies. The influence of the wind also has not been included in the source term dispersion regime. In case of intensive fires under prevailing 175 light or moderate wind speeds, this is a reasonable assumption. In case of very high wind speeds, allowing for the influence of the wind would increase the dilution of the source term, compared with the present computations. However, the influence of the vertical wind structure in the atmosphere has been allowed for in the treatment of the buoyant plume and naturally, in the subsequent atmospheric dispersion (Kukkonen et al., 2014). The fire is assumed to be in the flaming stage; that is the fire burning regime 180 that produces the highest atmospheric concentrations of pollutants.
In the following, we will first address the modelling of (i) the heat fluxes generated by fires and (ii) the average heights of the flames. These results will subsequently be used for deriving equations for (iii) the radius, velocity, temperature and molar flux of a fire plume. 185

Heat fluxes generated by fires
The heat energy of a fire is mostly convected or radiated from the fire region (e.g., Heskestad, 1984 and. A smaller fraction of the heat will also be conducted into the ground (e.g., Ichoku and Kaufman, 2005;Freeborn et al., 2008;Ichoku et al., 2012). The so-called theoretical heat release rate (i.e., heat energy flux) from the fire can be defined to occur, if the burning material is completely combusted. The 190 theoretical heat release rate (Q) can be evaluated as (Heskestad, 2016) where qm,f is the mass burning rate (i.e., rate of mass burned per time); and Hc is the lower heat of complete combustion (heat energy per burned mass). The lower heat of combustion refers to a situation, 195 in which the fire products are in the state, in which they have been formed (e.g., Drysdale, 2016), i.e., the potential subsequent phase transitions and chemical reactions have not been taken into account. For instance, in this situation any liquid water in the fuel that has been vaporized during the burning process is assumed to be in vapour form. The possible condensation of water therefore does not by definition contribute to the heat released by the fire. 200 The combustion efficiency  is defined as the ratio of the total to the theoretical heat release rate. This efficiency has been found to be close to unity for some fire sources (e.g., methanol and heptane pools), but it may deviate significantly from unity for others (Heskestad, 2016). However, as this efficiency is not commonly known in operational applications, we have assumed for simplicity a complete combustion 205 ( = 1).
The heat flux generated by a fire (Q) is propagated in the form of convection (Qc), radiation (Qr) and by other (Qo) means (e.g., conduction), where c, r and o are the fractions of convective, radiative and other heat fluxes of the total heat flux Q, respectively (by definition, c + r + o = 1). Laboratory experiments on biomass burning have demonstrated values of o  0.35 (Freeborn et al., 2008). For large fires, the radiative fraction r tends to decrease with the increasing size of fire (Heskestad, 2016). 215 The convective heat flux Oc can be written simply as (e.g., Achtemeier et al., 2012;Kukkonen et al., 2014;Heskestad, 2016)

220
where cp is the specific heat capacity of the plume,  is the density of the plume, u is the characteristic velocity of the plume, r is the characteristic radius of the plume and T = T -Ta is the excess temperature of the plume. T is the characteristic temperature of the plume and Ta is the temperature of ambient air.

Mean flame height
The flame intermittency, I(z), is defined as the fraction of time, during which part of the flame is above 225 the height z (Heskestad, 2016). Clearly, the flame intermittency decreases with height; it is equal to unity at the fire source and vanishes at sufficiently large heights. The mean flame height () is defined as the altitude, at which the flame intermittency has declined to a half of its initial value. At the height , most of the combustion reactions have taken place, and at higher altitudes, the plume can therefore be considered to be inert with a fairly good accuracy (Heskestad, 2016). 230 Several experimentally derived correlations have been proposed for . As would physically be expected,  has been found to correlate positively with the fire Froude number, i.e., dimensionless heat release rate Q * and the pool diameter d; which can be written as  ~ d Q* (e.g., Luketa and Blanchat, 2015). In general, Froude number is a dimensionless number defined as the ratio of the flow inertia to the external field; in 235 many applications the external field is gravity. Grove and Quintiere (2002), Dupuy et al. (2003), Luketa andBlanchat (2015) and Heskestad (2016) have presented comparisons of several of the correlations for  in terms of experimental data and against each other.
In this study, we have adopted the correlations of  presented by Zukoski et al. (1985). These correlations 240 performed amongst the best ones in a comparison of predictions and experimental data in large-scale liquid natural gas (LNG) burner experiments that were reported by Luketa and Blanchat (2015). The correlations of Zukoski et al. (1985) also provided physically sensible results in conceivable fire scenarios. These correlations can be written as where d is the diameter of a fire source or an equivalent diameter for a noncircular fire with the same area, and Q * is the Froude number of the fire, defined here as * = ( ) 1 2 ⁄ 2 ,

250
where a is the density of ambient air, cpa is the specific heat capacity of air at constant pressure and g is the acceleration due to gravity. According to Eq. (5), Froude number is large ( * ≫ 1) for fires, in which the energy output is relatively large compared to its physical diameter.

Radius, velocity, temperature and molar flux of a fire plume
The source term model presented in this study is based on the buoyant plume model, which is commonly 255 called the Morton-Taylor-Turner model (Morton et al., 1956; hereafter referred to as MTT). The MTT model is applicable for a steady plume of buoyant material that rises vertically in a calm atmosphere.
The MTT model applies to sources that have a relatively small difference in density compared to ambient air density (Boussinesq approximation), or to a region above the source , in which air entrainment has brought the plume density sufficiently close to the ambient value. 260 The coupled first-order differential equations of the MTT model govern the evolution of the characteristic plume radius (r), velocity (u) and density deficit ( = a -) above a point source. For readability, we have presented these equations and their mathematical solution in Appendix A.

270
The entrainment assumption of the MTT model states that the velocity of entrained air across the plume edge (ue) is proportional to the plume velocity = .
(7a) 275 However, there are many instances, in which plumes can not be modelled according to the Boussinesq approximation. Observations indicate that in the non-Boussinesq case, the entrainment velocity depends on the ratio of the plume density and the ambient density (Ricou and Spalding, 1961). Morton (1965) suggested, based on experimental evidence and theoretical considerations, an additional proportionality of the entrainment velocity, 280 which is usually referred to as the Ricou-Spalding entrainment model; hereafter referred to as RS.
Equation (7b) indicates a reduced entrainment into lighter plumes, in comparison with entrainment to plumes of near ambient air density. 285 For extending the model to other than point sources, a concept of virtual source can be introduced. The virtual source is located below the actual area source and is accounted for by replacing the height z by zzvs, where zvs is the height of the virtual source. In addition, to accommodate large density deficiencies, Morton's extension of the MTT model results in that the plume radius in the right-hand side 290 of Eq. (6a) has to be multiplied by a factor ( ⁄ ) 1 2 ⁄ and the ratio ∆ ⁄ in Eq. (6c) has to be replaced by ∆ ⁄ (Morton, 1965).
The comparison of model predictions and measurements of fire plumes above the flames have to a large extent supported the use of the above described theory (Heskestad, 1984 and1998). However, the 295 prediction accuracy can be improved by using experimentally adjustable coefficients. The plume radius ( ∆ ), and the mean values of velocity (u 0 ) and excess temperature (T 0 = T 0 -Ta) at the plume centre line have been found to obey the following relations (Heskestad, 1984 and1998):

300
where ∆ is the plume radius at the point where the excess temperature is 0.5Δ T 0 , T 0 is the mean temperature in the centre line of the plume, and C 1 = 0.12, C 2 = 3.4 and C 3 = 9.1 are experimental coefficients. The above values of the coefficients Ci (i = 1,3) are based on experimental investigations of heated air jets and large-scale pool fires (George Jr. et al., 1977;Kung and Stavrianidis, 1982). For burn experiments in rack storage fires, Kung et al. (1988) and Tamanini (2010) have recommended slightly 305 different values of Ci (i = 1,3).
For fire sources, which do not have substantial in-depth combustion, the height of the virtual source zvs can be estimated based on the experimental relation of Heskestad (1984), A fire source does not have substantial in-depth combustion, if a major fraction of the volatiles released (2/3 or higher) undergoes combustion above the fuel array (Heskestad, 1984). Fire sources with substantial in-depth combustion include, e.g., very openly constructed, or well-ventilated wood cribs.

315
Assuming ideal gas behaviour, the molar flux of gaseous material (qn) of the plume is where pa is the atmospheric pressure (pressure within the plume is assumed to be equal to the ambient value) and Rg is the molar gas constant. The flux qn comprises of the molar fluxes of air (qn,a) and gaseous 320 combustion products (qn,c), Detailed modelling of the two selected application areas has been presented in Appendix C.

Interface of the source term and plume rise regimes 325
The mean velocity (u 0 ) and excess temperature (T 0 ) in Eqs. (8b-c) are values at the centre line of a fire plume. As the radial distance from the centre line increases, u and T approach their ambient values. We have assumed Gaussian distributions for the radial distributions of velocity and temperature in the source term regime.
However, the subsequent modelling phase in the plume rise regime assumes that the plume is described by a uniform (top-hat) distribution of physical quantities. If the centre line values of the Gaussian distributions (u 0 , T 0 ) would be used to represent a top hat profile, the convective heat would not in general be conserved at the boundary of the source term model and the plume rise model.

335
We have therefore presented a solution in Appendix B, in which the convective heat is conserved at the boundary of the two modelling regimes. According to this solution, the top hat velocity (u) and excess The model is not applicable for 1 2 2 3 < 1, as required by the conservation of heat energy, as presented in detail in Appendix B.

345
Detailed modelling of the fluxes of compounds originated from the two selected fire types have been presented in Appendix C. For readability, we present a brief description of the experiments and measured data, which were used in this study. For a more detailed description of the experiments, the reader is referred to a special issue of 375 the International Journal of Wildland Fire, which was aimed to document the RxCADRE study (Peterson and Hardy, 2016). In particular, Ottmar et al. (2016a and has presented overviews of the RxCADRE study.

Summary of the use of the model equations and the computer model
The RxCADRE measurement campaign consisted of six smaller (less than 10 ha) and ten larger 380 (10 -900 ha) prescribed fires (Ottmar et al., 2016a). Measured data that is sufficient for dispersion model evaluation is available for three larger fires: two grass fires (named by the authors as L1G and L2G, in which L presumably refers to larger experiments and G to grass) and one sub-canopy forest fire (named as L2F, in which F refers to forest) (Clements et al., 2016;Dickinson et al., 2016;Strand et al., 2016). In case of these three experiments, data is available regarding the meteorological conditions, fire emissions 385 and airborne concentrations.
As the main focus of the present study was on forest fires, we have selected the L2F fire for model evaluation. This experiment was conducted over a burn block of size 151 ha at Eglin Air Force Base,  (Urbanski, 2014b;Strand et al., 2016). However, there were neither aircraft measurements on particulate matter nor ground-based measurements in the publicly available dataset. We have therefore conducted the model evaluation using the concentration data on gaseous substances from aircraft measurements in the experiment L2F.

400
The accuracy of the meteorological measurements within the L2F burn was evaluated by Clements (2015a and2015b). They concluded that the data are qualitatively reliable. Urbanski (2014b) analyzed the accuracy of the airborne concentrations, regarding the positioning by GPS and the applied spectroscopic methods (Cavity ring-down spectroscopy, CRDS). The estimated analytical uncertainty of the CRDS measurements varied from 1 % to 1.5 % for CO 2 and CH 4 , and from 2 % to 15 % for CO. 405

The selected experiment, the sub-canopy forest fire (L2F)
The radiative heat fluxes of the fire were evaluated using long-wave infrared (LWIR) measurements (Dickinson et al., 2016). The measurements were done onboard a twin -engine Piper Navajo aircraft, which was used to make repeated passes at about three-minute intervals (Dickinson et al., 2016).

410
The aircraft measurements yielded data of fresh emissions, vertical profile of the smoke, plume height and the dispersion of smoke (Urbanski, 2014b). Measurements were conducted at distances of up to 25 kilometres downwind from the source. Measurements were done as so-called parking garage and corkscrew flight profiles. Parking garage vertical profiles involved short (approximately 10 km) horizontal transects, roughly perpendicular to the axis of the smoke plume, taken at multiple elevations. 415 Corkscrew profiles were centered on the plume downwind from the burn unit. Parking garage and corkscrew manoeuvres were designed for measuring the horizontal and vertical concentration distributions, respectively.
The observations of the heat fluxes were used for deriving the fire radiative power (FRP) against time. 420 The evaluated temporal evolutions of the FRP and the burning area are presented in Fig. 2. We applied the meteorological measurements prior to the ignition of the L2F burn. Ten-minute averages 450 of measured wind speed, wind direction, ambient temperature and ground level air pressure were applied for further processing. The Monin-Obukhov length (L) was estimated by fitting the atmospheric vertical profiles used in the BUOYANT model (Kukkonen et al., 2014) to the averaged temperature and wind measurements, using the method presented by Nieuwstadt (1978). The two -layered thermal structure above the atmospheric boundary layer (ABL) was evaluated by applying the measured temperature 455 profile, according to the method by Fochesatto (2015).
Prevailing wind direction was evaluated to be south-easterly at 132°. The modelled wind speed and ambient temperature at the altitude of 10 m were evaluated to be 3.2 m s -1 and 24 °C, respectively. The atmospheric stability was estimated to be moderately stable (L -1 = 0.0011 m -1 ). Based on the observed 460 temperature profile, the height of the ABL was estimated to be 2.2 km. The gradients of potential temperature were estimated to be 0.0193 K m -1 and 0.0094 K m -1 within the inversion and upper layers, respectively. Wind speed above the ABL was assumed to be constant and equal to the modelled value at the top of the ABL, 14 m s -1 .

Evaluation of the properties of the fire source 465
The properties of the fire source term could in principle be evaluated using the source term model that has been presented in this article. However, for this particular measurement campaign, it is better to use directly the values that were reported in the database and the relevant publications regarding the experiment. The source term model presented in section 2.2. was therefore not applied in the following evaluation exercise. Instead, the evaluation has been made for the BUOYANT model, by selecting the 470 option of not using the fire source term module.
We have assumed a steady state of the fire in the modelling, as the model is not capable of treating time-

Results of the model evaluation 485
We consider in the following observed and modelled excess concentrations, i.e., concentrations subtracted by background concentration. These represent the contributions of the fire. We focus on the comparison of the measured and predicted spatial concentration distributions, in the horizontal and vertical directions. The aircraft measurements were specifically designed to measure such distributions.

490
The model does not contain any free parameters and was not adjusted in any way to the measured data.
Measured and modelled vertical excess concentrations of CO 2 are presented in Fig. 3, for three parking garage and the second corkscrew flight paths. The flight durations varied from 6 (PG #1) to 16 minutes (PG #2). 495 Concentration peaks tend to be overestimated by the model (also seen in Fig. 3), while the widths of the 520 plume are slightly underestimated. First, the modelling assumed a steady state of the fire. This implies that the fire intensity was assumed to be temporally constant throughout the whole duration of the experiment. The model therefore tends to underpredict the concentrations in the initial stages of the fire (PG #1 and CS #2) and overpredict these in the later stages (PG #2 and #3). However, the temporal agreement of the measured and modelled 535 highest concentrations was good.
For the modelling of fire plumes, it is crucial to evaluate sufficiently accurately the vertical structure of the atmosphere, especially the potentially existing temperature inversions. Kukkonen et al. (2014) previously compared the predictions of the BUOYANT model against the measurements in two other 540 field measurement campaigns. They commented that, e.g., evaluating the meteorological conditions in the SCAR-C experiments (Kaufman et al., 1996;Hobbs et al., 1996;Gassó and Hegg, 1998 ) using two different meteorological methods, resulted in substantially different meteorological input data values for the model.

545
In case of the RxCADRE measurements, the relevant meteorological parameters have been carefully measured and well reported. However, the application of such datasets in determining vertical atmospheric profiles of the relevant quantities, and the atmospheric stability conditions, will result in some degree of inaccuracy in the dispersion modelling.

The operational version of the model 550
We have compiled an operational version of the model. The main aim of developing this version was to provide a user-friendly tool of assessment for various emergency response personnel. The operational model can be used for emergency contingency planning and for the training of emergency personnel, in case of forest and pool fires.

Overview and functioning of the operational model 555
The operational program has been named as FLARE (Fire pLume model for Atmospheric concentrations, plume Rise and Emissions). An overview of the model structure has been presented in Fig. 5. This model contains the program BUOYANT for conducting the physical and chemical computations, a graphical user interface, and various modules for processing the input and output data of the model. The operational version can be used remotely via an internet connection. 560 The model addresses forest fires and liquid pool fires. The user needs to specify as input values only the 575 location and time of the fire event, the estimated area of the fire, and in case of a pool fire, the released substance. In addition, the model pre-processes and provides for the computations three main types of input data: meteorological parameters, forest information and geographic maps.
The program will subsequently check that all the user-specified input data values and their combinations 580 are physically reasonable. The program will also check that the computations address cases, which are within the applicability range of the model. In case of unrealistic or unreasonable input values, the program will either request the user to confirm the value or to input a more realistic value.
However, the current version of the operational model can be used only for locations that are situated in 585 Finland, or in close vicinity of the country. The operational model could be extended to function also in other countries and regions, by expanding especially the cartographic and forest inventory datasets. In case of missing input datasets, the model could also be modified to skip some of the input processor modules and ask the user to input the corresponding values. For instance, if there would n ot be a suitable forest inventory available for the considered domain, the user would be asked to supply the required 590 information on the characteristics of the forest.
The user can archive descriptions of fire events, which contain input data for a range of potential fire scenarios. These cases can then be retrieved, edited as necessary, and used for further computation.

595
The program presents the numerical results as pollutant iso-concentration curves on maps. The current operational version presents spatial concentration distributions of carbon dioxide (CO 2 ) near the ground level.

The pre-processing of the input datasets
The functioning of the operational model has been made as user-friendly as possible, by an automatic 600 pre-processing of several input materials. The meteorological parameters and forest inventory data will be extracted and preprocessed for the spatial coordinates and the time specified by the model user. The model also presents the results on geographic maps, for the domain selected by the u ser.
The automatic on-line use of weather and forest data makes the use of the model substantially quicker 605 and simpler. This will also reduce potential human errors. For non-expert users, the determination of the required meteorological variables would otherwise be a very challenging task. In case of long-term fires, the user can also use the forecasted meteorological values, for forecasting the spread of the fire plumes for up to two days ahead in time.

610
These input datasets and pre-processors are briefly described in the following. Additional information regarding the meteorological data is presented in Appendix D.

Meteorological data
The program can use either real-time or forecasted meteorological data produced by the numerical weather prediction (NWP) model HARMONIE; this model is run operationally at the FMI. The acronym 615 HARMONIE has been derived from "HIRLAM ALADIN Research on Mesoscale Operational Numerical weather processing In Euro-Mediterranean Partnership" (Nielsen et al., 2014). The modelling domain includes Fennoscandia, the Baltic Countries and the surrounding regions in the eastern Atlantic, Northern Central Europe and Russia.

620
The HARMONIE model was selected for three main reasons. First, this model has been thoroughly evaluated against experimental data and it is known to provide accurate, high-resolution weather forecasts for the whole of Fennoscandia. The treatments of this model have been specifically 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 operational program and the 625 NWP model. Third, the HARMONIE system will be the operational weather forecasting model in Finland in the near future (replacing the older HIRLAM NWP model).
However, the output data of this NWP model does not directly include all the input values required by the BUOYANT model. We have therefore constructed a continuously functioning pre -processor model, 630 which evaluates the required meteorological parameters based on the output of the HARMONIE model. These parameters are the ambient temperature and pressure, the lateral wind components, the inverse Monin-Obukhov length, the height of the atmospheric boundary layer and the vertical profiles of temperature and wind speed above ABL.

635
The meteorological roughness length, which will be needed in the atmospheric dispersion computations, is evaluated based on the CORINE (CoORrdination of INformation on the Environment) land cover information in 2012, using in addition the weighting coefficients modified by Venäläinen et al. (2017).

Forest information
In case of forest fires, the amount of burnt material is evaluated based on a national inventory of forests 640 (Mäkisara et al., 2019). This inventory has been compiled by the Natural Resources Institute Finland, and it is called Multi-Source National Forest Inventory of Finland (MS-NFI). The methods and results of this inventory have been presented by Tomppo and Halme (2004). The inventory is publicly available.

Geographic map information
The model provides as output spatial concentration distributions near the ground level, presented on 655 digital maps. The model uses open-access maps provided by the National Land Survey of Finland. The user can specify the location of the accident by simply clicking that point in the map, by specifying the geographic coordinates or by writing the street address. The accident location will then be searched and the coordinates will be extracted from the Finnish National Geoportal (National Land Survey of Finland, 2021). This location will then be automatically placed on a map. For efficient functioning, this service 660 has been adapted to the computer facilities at the Finnish Meteorological Institute (FMI).

Conclusions
We have presented a refined version of a mathematical model, BUOYANT, which has been designed for analyzing the formation and dispersion of plumes originating from major fires. The model addresses the cross-plume integrated properties of a rising plume in a vertically varying atmosphere; the model also 665 takes into account the impacts on plume rise of possibly occurring inversion layers (Kukkonen et al., 2014). In the present study, the BUOYANT model has been extended to include a more detailed description of the early development of the fire plume. This generalization also made it possible to compile an operational model version, which can be used in a much more straightforward way, compared with the use of the original research model. 670 The developed source term model can be used to evaluate the characteristics of the fire plume, which can be used as input for the subsequent computations on the evolution of a buoyant plume. The source term model uses as input the information on the characteristics of the fire, and it is used to evaluate the properties of a fire plume just above the flame tips. The current version of the source term model can be 675 applied for two significant categories of fires, viz. the forest and liquid pool fires. In future work, the source term model could be generalized to address also other fire types. The main model structure of the source term model is based on the differential equations for releases of buoyant material, which govern the evolution of the plume radius, velocity and density difference. However, the model can be considered to be semi-empirical, as it also relies on various experimental results on fire plumes. 680 We have compared the predictions of the refined BUOYANT model against the experimental field-scale data of the RxCADRE campaign (Prescribed Fire Combustion and Atmospheric Dynamics Research Experiment). These experiments were designed to collect extensive quantitative data regarding the burning of prescribed fires. These datasets have provided accurate measurements of various aspects of 685 the fires, including meteorology, the evolution of fires, their energy, the emissions and airborne concentrations. The predicted concentrations of CO 2 agreed fairly well with the aircraft measurements of the RxCADRE campaign. For instance, the model captured well the observed vertical excess concentration distributions of CO 2 during the parking garage flight manoeuvres, for most of the highest concentrations. However, the model tended to moderately overpredict the highest concentrations, 690 whereas the widths of the plume were slightly underestimated.
There are several reasons for the differences of the measured and predicted concentrations. Previous comparisons of the predictions of plume rise models and experimental field -scale data have illustrated several major challenges in determining accurately the source properties and the meteorological 695 conditions (e.g., Kukkonen et al., 2014). These are also major sources of uncertainty in the present comparisons against the RxCADRE data.
An important limitation of the present modelling is that it has assumed a steady state of the fire. This implies that the fire intensity was assumed to be temporally co nstant throughout the whole duration of 700 the experiment. Clearly, the influenc of this assumption could to some extent be taken into account by conducting several computations with the model, using various values of the fire intensity. The influence of the steady state assumption is that the model tended to underpredict the concentrations in the initial stages of the fire and to overpredict these in the later stages. The model currently assumes that the burned material consists solely of standing tree trunks. Clearly, also other kind s of plant material contribute to 705 the burning in a forest.
Another source of uncertainties in the modelling is the evaluation of the relevant meteorological parameters. The meteorological measurements in the RxCADRE campaign have been carefully measured and reported. However, the application of such data for determining vertical atmospheric 710 profiles of the relevant quantities, and the atmospheric stability conditions, will inherently result in some inaccuracies.
We have also compiled an operational version of the model. The operational model is a user-friendly tool of assessment that can be used by various emergency response and rescue personnel. This model 715 can be used for emergency contingency planning and for the training of emergency personnel, in case of forest and pool fires. The model has been used by Finnish rescue authorities up to date . However, it would be possible to use both the original research model and its operational application also worldwide, the latter after some adjustment of the processing of the model input datasets. This could potentially result in improved preparedness and better, knowledge-based rescue actions in case of major fires. 720

Appendix A. The Morton-Taylor-Turner model for a buoyant plume
Let us consider a plume from a point source, assuming no momentum flux at the source, uniform ambient 725 air density and the Boussinesq approximation. The conservation of mass, momentum and buoyancy can be written as (Morton et al., 1956) ( 2 ) = 2 , (volume/mass), ( 2 2 ) = 2 ( − ), (momentum) and (A1b) where z is the height above ground; r is the radius of the plume; u is the vertical velocity of the plume; 730 ue is rate of entrained air across the plume edge (entrainment velocity ); g is the acceleration due to gravitation; a is the density of ambient air; and  is the density of the plume.
The entrainment velocity is assumed to be proportional to some characteristic velocity at height z (Morton et al., 1956) where  is an experimentally defined proportionality constant (the entrainment constant) relating the entrainment velocity to the vertical velocity within the plume. Equation (A2) is often referred to as the Morton-Taylor-Turner entrainment model. 740 The solution of Eqs. (A1a-A1c) is (Morton et al., 1956) = 6 5 , = 5 6 ( 9 10 ) where the constant buoyancy flux B is 745 Assuming ideal gas behaviour, the buoyancy flux can be written in terms of convective heat flux (Qc) as (Heskestad, 2016) 750 = .
(A5) Appendix B. Centre line properties of a fire plume in the source term regime and the equivalent

top-hat profiles
The mean velocity (u 0 ) and excess temperature (T 0 ) at the centre line of a fire plume in the source term 755 flow regime have been presented in Eqs. (8b) and (8c). These values approach their ambient values, as the radial distance from the plume centre line increases. We present in the following a model for determining the equivalent mean velocity and excess temperature for uniform (i.e., top -hat) profiles of the plume cross-sections, under the condition that convective heat energy is conserved.

760
We assume Gaussian radial profiles of the excess temperature, T(r), and the mean velocity, u(r) (Heskestad, 2016), where r is the radial distance measured from the centre line of the plume; and  T and u are the measures 765 of the plume width corresponding to the radial distributions of excess temperature and velocity, respectively. The density of the plume is assumed to have a constant value within each cross-section of the plume, equal to the centre line value ( 0 ).
The radius of the plume, ∆ , has been defined in terms of ∆ 0 (Eq. (8a)). A velocity radius (ru) can be 770 defined correspondingly: let ru be the plume radius at the point, in which the gas velocity has declined to 0.5u 0 (Heskestad, 2016). The temperature and velocity profiles have in general differing scales, i.e., According to Heskestad (2016), an optimal value is a = 1.1, based on the most reliable measurements 775 (George Jr. et al., 1977).
where the subscript x is either u or T; and b  1.201.
The equivalent top-hat excess temperature (T) and velocity (u) of the plume can derived by integrating Eqs. (B1a) and (B1b), and using the relations (B3), 785 where R is a radial distance from the centre of the plume.
Equations (B4a-b) can be written more simply in terms of the error function (erf), defined as 790 Therefore For C < 1, the conservation of convective heat energy cannot be achieved applying the presented method.

Appendix C. Detailed modelling of the selected application areas
Semi-empirical models of mass burning rate are presented in the following for liquid pool and forest fires. Mass fluxes of the emitted chemical compounds (e.g. CO, CO 2 ) from the fire are, by definition, 845 determined employing the modelled mass burning rate and emission factors. Hottel (1959) suggested how to analyze liquid pool burning, according to heat transfer principles (Babrauskas, 1983). According to Hottel (1959), the mass burning rate is governed by the heat exchange between the flames and the pool surface. The heat exchange mechanisms are (a) radiant flux from the 850 flames into the pool, (b) convective flux from the flames into the pool, (c) re-radiant heat loss (Qrr) due to high temperature of the pool and (d) conduction losses and non-steady terms (Qmisc).

C1. Mass fluxes of pollutants originated from liquid pool fires
The term Qrr is commonly small, and quantitative expressions for Qmisc are usually not available (Babrauskas, 1983). For simplicity, the terms Qrr and Qmisc are therefore customarily ignored 855 (Babrauskas, 1983). Hottel (1959) analyzed the experimental data of Blinov and Khudiakov (1957), concluding that two burning regimes are possible: radiatively dominated burning for larger pools and convectively dominated burning for smaller pools. The distinction between the two regimes can be drawn at the pool diameter of approximately 0.2 m (Hottel 1959;Babrauskas, 1983;Chatris et al., 2001). For the purposes of fire hazard analysis, liquid pool fires will rarely be significantly dangerous, if these are 860 smaller than about 0.2 m in diameter (Babrauskas, 2016). It is therefore commonly necessary to treat only pool burning in the radiative regime. Zabetakis and Burgess (1961) suggested (cf. Babrauskas, 1983;Chatris et al., 2001;Brambilla and Manca, 2009), based on the work of Hottel (1959), the following relationship to represent the mass 865 burning rate (qm,f) of a liquid pool in the radiative dominated regime: where ,∞ is the mass burning rate per unit area of an infinite-diameter pool; A is the surface area of the burning liquid pool; k is the extinction coefficient of the flame; β is a mean-beam-length corrector; 870 and d is the diameter of pool. For small d the flames are said to be optically thin, while for larger d the flames become optically thick (Babrauskas, 1983). For optically thick flames, a further increase in d does not result into a corresponding increase in back radiation into the pool. Such attenuation is accounted for by the coefficient k (Brambilla and Manca, 2009).

875
Values of the empirical coefficients ,∞ and kβ for a variety of fuels have been proposed by, for instance, Babrauskas (1983), Rew et al. (1997) and Chatris et al. (2001). The surface area of the burning liquid pool and the name of the liquid fuel have to be provided as input data for the model.
The total heat generated by a liquid pool fire is here assumed to be propagated only through convective and radiative processes, i.e. o = 0  c = 1 -r. Radiometer measurements from large fire experiments involving different combustible liquids (such as crude oil, heptane and kerosene) suggest that the 885 radiative fraction (r) decreases with increasing fire diameter (d), according to (McGrattan et al., 2000) where max = 0.35; and k = 0.05 m -1 .
where mw ,i is the molecular weight of species i; and qn ,i is the molar flux of species i. Experimental values 895 of yields under well-ventilated fire conditions have been listed by, for instance, Ross et al. (1996) and Hurley (2016).  (10) and (11). 900 Examples of fuel property data are shown in Table C1. Table C1. Examples of the fuel property data (Babrauskas, 1983, Hurley, 2016

C2. Mass fluxes of pollutants originated from forest fires
McAllister and Finney (2016a and 2016b) have evaluated the mass burning rate of wildland fires. Wood 910 cribs, such as the one presented in Fig. (C1), have been used in fire testing. Block (1971) developed a theoretical model of the crib burning rate. Heskestad (1973) combined the experimental results of Gross (1962) and Block (1971) with the theoretical findings of Block (1971). This resulted in a relation of the mass burning rate (qm,f) on the porosity () of the crib (see also McAllister and Finney, 2016a), 915 10 3 , where As is the exposed surface area of the sticks in the crib; and b is the thickness of the sticks (defined in Fig. C1). The functional form of f was found to be approximately (McAllister and Finney, 2016a) 920 For well-ventilated cribs or loosely-packed porous burning cribs,  is large and f() approaches unity.
However, according to Tang (2017), both well-ventilated and under-ventilated fires can exist in forested regions. We have assumed in this study, for simplicity, that the fuel beds are porous ( ( ) = 1).
Let us define the diameter of a tree trunk at human breast height, dbh. Commonly dbh is measured 930 approximately at a height of 1.3 m. Assuming that dbh is a representative value of b, we can approximate the mass burning rate of a porous wildland fire to be Assuming that all of the trees from the ground to treetop are on fire, the exposed surface area of one tree 935 is equal to ℎ ℎ , where ht is the average height of the burning trees. The exposed area of all the trees can therefore be approximated by where nt is the number of burning trunks per unit forest area burning; and A is the area of forest on fire. 940 Heat generated by a forest fire is estimated from Eqs. (1) and (C6). The lower heat of combustion (Hc) of woody fuel ranges typically from 17.8 to 20.4 MJ kg -1 (e.g. Trentmann et al. 2006;Hurley, 2016). We have therefore applied the middle value within this range (Hc = 19.1 MJ kg -1 ).

945
The fraction of total energy released by combustion that is available for convection depends on the ambient and fuel conditions (Trentmann et al., 2006;Freitas et al., 2010;Kukkonen et al., 2014).
Laboratory experiments of biomass burning (Freeborn et al., 2008) have indicated a mean convective fraction of 51.8  9.0 % (determined in terms of higher heat of combustion, i.e., including latent heat released during condensation of water vapour generated by the fire). We have assumed that 55 % of the 950 total heat generated by a forest fire is available for convection (c = 0.55). This is simply in the middle of the commonly accepted range from 0.4 to 0.8 (Trentmann et al., 2002;Freitas et al., 2010).
The emission factor can be defined to be the amount of the chemical species released per mass of dry biomass burned (e.g., Andreae and Merlet, 2001). Therefore, emission factor is equa l to the yield of 955 combustion products (yi). Data on emission factors for various types of biomass burning have been presented by, for instance, Lemieux et al. (2004), Akagi et al. (2011), Kaiser et al. (2012) and Urbanski (2014a. The current model version applies the emission factors, which are applicable for land cover class of extratropical forest, presented by Kaiser et al. (2012). The extratropical forest class includes forest types typically found in the northern hemisphere (Kaiser et al., 2012). 960 For simplicity, particles formed in a forest or a liquid pool fire are assumed to be spherical. Further, they are assumed to be 2.5 m in aerodynamical diameter having the density of water, i.e. the unit density = 1 kg dm -3 .

965
Appendix D. The extraction and pre-processing of meteorological data.
The program can use either real-time or forecasted meteorological data produced by the numerical weather prediction (NWP) model HARMONIE. HARMONIE is a state-of-the-art NWP model, which has been widely used and developed in Europe. The main limitation of the HARMONIE model, as applied in the present study, is its fairly limited geographic domain . 970 HARMONIE is a non-hydrostatic convection-permitting NWP model. The horizontal grid spacing of the model is 0.022° (approximately 2.5 km). The vertical grid consists of 65 vertical hybrid levels. In this study, we applied the HARMONIE version cy40h11, which is in operational use at the FMI.
Meteorological forecasts are continuously produced four times a day, with a temporal resolution of 975 one hour and a forecast length of about two days ahead in time.
Most of the meteorological variables required by the BUOYANT model are directly available from HARMONIE forecasts. The vertical structure of the atmosphere in the BUOYANT model assumed to comprise three distinct layers: the atmospheric boundary layer (ABL), capping inversion layer and upper 980 layer (Kukkonen et al., 2014). Variables that are readily available in HARMONIE include height of the ABL and the vertical profiles of the temperature, pressure and wind speed.
In the ABL the vertical variations of wind speed and temperature are in this study assumed to be described with profiles based on the Monin-Obukhov similarity theory, as presented by Kukkonen et al. (2014). 985 Monin-Obukhov length is estimated based on the values of the turbulent momentum stress near the ground surface, as forecasted by HARMONIE. The two-layered thermal structure above ABL (inversion and upper layers) is evaluated by applying the HARMONIE predictions with a method modified from Fochesatto (2015).

990
In the upper layer (above the inversion layer), the wind speed is assumed to be constant (representing the geostrophic flow), whereas within the inversion layer the wind speed is assumed to change with constant gradient from its value at the top of the ABL to the geostrophic value. The constant geostrophic wind speed was assumed to be equal to the arithmetic mean of HARMONIE forecasts between the top of the inversion layer and the height of 5 km. 995

Code and data availability
The code and the relevant data are available in Zenodo at https://doi.org/10.5281/zenodo.4744300 (Kukkonen et al., 2021). These contain the source code of the BUOYANT model (v4.20), the technical reference of the model, the user manual of the model and the model input data corresponding to the work 1000 described in this paper. The model code, documentation and the input data are published under the license of Creative Commons Attribution 4.0 International.
The experimental data of the RxCADRE campaign used in this paper can be downloaded from the Research Data Archive of the U.S. Department of Agriculture Jimenez andButler, 1005 2016;Clements, 2015a and2015b;Urbanski, 2014b).

Author contribution
The research version of the BUOYANT model, including the source term module, has been developed by Juha Nikmo, Jaakko Kukkonen and Kari Riikonen. All the authors have contributed to the development of the operational model version. Ilmo Westerholm, Pekko Ilvessalo, Tuomo Bergman and 1010 Klaus Haikarainen performed most of the research and coding that was necessary for the functioning of the operational model. Juha Nikmo performed the model computations for the evaluation of the model. Jaakko Kukkonen was the leader and coordinator of the project, in which this work was performed. Jaakko Kukkonen and Juha Nikmo prepared the manuscript, with contributions from all co-authors.

Conflicting interests 1015
The authors declare that they have no conflict of interest.