A Mixed Convection Model for Estimating the Critical Velocity to Prevent Smoke Backlayering in Tunnels

: A novel mathematical model for the critical ventilation velocity to prevent smoke backlayering in tunnels is presented, addressing limitations of prior approaches. The basis of the model is a rigorous characterisation of the physical processes by the characteristic quantities. Empirical parameters within the new model are determined, to align with results from both full-size and small-scale tunnel experiments. Data from numerical simulations (CFD, Computational Fluid Dynamics), validated by known test data, are then used to estimate the effects of tunnel slope and other parameters on the critical velocity. The model is seen to approximate the critical velocity well, following all trends identified by test data and CFD parameter studies. The empirically calibrated equation permits prediction of the critical velocity beyond the narrow range of tunnel geometries where known results already give an answer. The resulting equation has practical application for tunnel design.


Maximum nondimensional critical velocity value for different non-dimensional trend lines (-)
Specific isobaric heat capacity based on average temperature between approaching air and mean temperature just downstream of the fire, according to equation (4).
(J/kg K) Empirical constant used in fire length function.See Table 2 for pool fires.

𝑓 ℎ
Empirical factor relating   to the length of that part of the fire that influences backlayering propensity.See Table 2 for pool fires. (-) Empirical constant used in fire length function.See Table 2 [6] and in Annex D of NFPA 502 version 2023 [48].Height of duct or tunnel at fire site in Annex D of NFPA 502 version 2014 [44] (m) Tunnel height of a scaled tunnel (m)  𝑓𝑖𝑟𝑒 Fire intensity (heat release rate per unit area) (e.g.pool fire with oil fuel 2.25 MW/m 2 [27])

𝐾 𝑔
Grade factor, see "Critical velocity" is that speed of oncoming tunnel air which just prevents upstream backlayering of smoke in a tunnel fire.Achieving critical velocity as a minimum airspeed in fire response is often a design criterion for tunnel ventilation, particularly in the US, Australia and some countries in Asia.There is generally no practical upper limit in those jurisdictions.Other jurisdictions do not see this critical velocity approach as the best or recommended approach [50].European approaches [38][39][40][41][42][43] for example encourage more moderate air speeds, acknowledging that risks from over-ventilating are also significant.The US standard NFPA 502 historically relied on the critical velocity philosophy, but recently moved away from preventing backlayering to controlling backlayering [48], allowing some extent of smoke backlayering.However, even if absolute prevention of backlayering from a tunnel fire should not always be the first choice during the self-rescue phase, it can still be a useful design or benchmark criterion.Therefore, there is motivation to have a reliable tool for estimating critical velocity.
In the published literature, there have been three approaches to get to the answer: (i) conducting full-scale fire tests (expensive and hence rare), (ii) using a Computational Fluid Dynamics (CFD) model validated against full-scale tests, or (iii) calculating the critical velocity using accepted semi-empirical equations.
Over the last decades, many full-scale fire tests were carried out, but few of them were designed to accurately determine critical velocity in tunnels for large fires.The best-known full-scale data for critical velocity with a wide and relevant range of heat release rates (HRRs) are the Memorial Tunnel fire tests [27][28][29].These tests are well documented and thus well suited for validation.Several authors have compared different CFD models and approaches with the Memorial Tunnel full-scale tests previously (see Karki et.al. [35], Rhodes [55] and Galdo Vega [56]).For the purpose of the present study, a new CFD model validated by data from the Memorial Tunnel fire tests was developed and recently presented by Beyer and Stacey [51].Conducting CFD simulations can be time consuming and requires knowledge about the software and the choice of physical models and boundary conditions but, when validated appropriately, is able to provide reasonably accurate answers for a specific tunnel project with tunnel geometry and fire scenario different from the known full-scale data.It also allows a "confinement velocity", i.e., the air speed allowing for a designated length of hot smoke backlayering under the tunnel ceiling, to be determined.
The simplest and fastest way for estimating a critical velocity is the use of a wellvalidated semi-empirical equation.There are many semi-empirical approaches based on dimensional or non-dimensional models .The results and answers of those approaches can be very different.However, only a few of them have been adopted in versions of the US standard NFPA 502, which confers contractual relevance for real road tunnel projects.Therefore, more focus is laid on these calculation methods for critical velocity.Some of the references noted above  are summarised and discussed in Section 2. For a more detailed explanation of the different methods, reference is made to other work [22][23][24][25][26].
For many years, Annex D of NFPA 502 [44] proposed an equation introduced by Kennedy [6].Even though the equation proposed by Kennedy was mainly based on the Memorial Tunnel fire test results, the equations in the Annex D of NFPA 502 were first amended [45] and then replaced [46] to better represent small-scale test results published by Li et al. [18][19][20][21].
As discussed above, there has been conflicting and confusing guidance from NFPA 502 for designers on estimating critical velocity of smoke backlayering in tunnel fires.The 2023 edition of NFPA 502 [48] did not achieve the required clarity.One recommendation to designers was to look at the original test data for values appropriate to a particular tunnel (Beyer et. al. [50]).However, even with real test data available, it is still helpful to have a physically based mathematical relation allowing the critical velocity to be estimated for tunnel geometries or cases not well represented by known tests.
The present work reports such a mathematical formulation for critical velocity and its validation.Section 2 reviews the current methodology.Section 3 presents the proposed physical model, and Section 4 provides the CFD validation and a systematic parameter variation study to test the derived physical model.In Section 5, the model is validated against real test data.As the equation in Section 3 is perhaps too complex for design purposes, Section 6 provides a simpler, practical estimation of the critical velocity for design purposes.Conclusions of the work are summarized in Section 7.

Current Methodology
The present section reviews frequently discussed methods for calculating the critical velocity.Several articles that provide a comprehensive overview are already available [22][23][24][25][26] and they are not repeated here.The methods are based on characteristic numbers of fluid mechanics, which are briefly reviewed.The models essentially yield the required critical velocity as a dimensional result.

Froude Number
The Froude number  is a non-dimensional number representing the ratio of the inertial force to the gravity force; where  is the gravitational acceleration,  a characteristic velocity and  a characteristic length scale of the flow field.In fluid mechanics, the Froude number is significant for flows with free surface effects, like water waves in open channel flows, and for any isothermal fluid flow, which can be significantly influenced by the body force due to gravity.In the tunnel community and in fire engineering, a different form and definition of the Froude number, sometimes called the 'critical Froude number', can be found.The distinction is discussed in the following Section 2.2.

Dimensional Model
One of the first calculation models for the critical velocity was introduced by Thomas [1][2].He concluded that the upstream spread of hot smoke in tunnel fires can be described by the ratio between the inertial force (velocity head) and the buoyancy force (buoyancy head).It was assumed that, at the critical condition where backlayering of hot smoke is just prevented, the stagnation pressure (velocity) head will be of the same order of magnitude as the buoyancy head.The inertial force,   and the buoyancy force,   , as defined by Thomas [1][2] are: (2) where  is a characteristic length scale in the problem,  is the isobaric thermal expansion coefficient of the gas, and ∆ is the temperature difference relevant for the buoyancy force.The approach taken by Thomas determines ∆ in the buoyancy as the temperature increase of a heat capacitance flow rate      caused by the heat release rate (HRR)  ̇, as it follows from the First Law of thermodynamics for an open system.
Here,  is the tunnel cross-sectional area and   the specific isobaric heat capacity of the gas.By expressing the temperature rise as a function of air speed, the equation above connects the temperature increase, and hence the buoyancy force, to inertia effects.Substituting this equation into equation (3) and forming the ratio of inertia force to buoyancy force (  /  ) gives the Froude number as defined by Thomas [1].
The number above is also named the 'critical Froude number', but in fluid dynamics, the number is a Richardson number  (i.e., a ratio of Grashof number  to Reynolds number squared  2 ) rather than a Froude number as given by equation (1).In this approach there is only one characteristic length scale, other than that embodied in the tunnel cross-sectional area A. As a result, the method does not distinguish between tunnel height or hydraulic diameter, or the fire plume height.The length scale here is taken to be the tunnel height,  = .Thomas then concluded that the critical Froude number is expected to be in the order of unity at the flow situation where backlayering of hot smoke is just prevented, arriving at a critical velocity formula of where  has been replaced by the inverse of the temperature of the hot gas, 1/  , treating the gas as ideal.
A similar equation was proposed by Danziger and Kennedy [4,6].In contrast to equation (6), the total HRR  ̇ was replaced by the convective HRR,  ̇(1 − ), with the fractional convective heat reduction  due to radiation.A correction for the tunnel slope was also introduced.Furthermore, the temperature   was replaced by the bulk temperature downstream of the fire,   .This led to the critical velocity with the grade correction factor for downgrade slopes  given as The proposed grade correction factor was derived from studies on the control of methane layers in coal mines.It was suggested that the grade correction should only be used for ventilating downgrade.For uphill grades,  , is 1.0.
The temperature   in equation ( 7) is the bulk temperature downstream of the fire and is expressed as where   is the approaching ambient air temperature.For the critical Froude number, Kennedy [6] used a conservative value of 4.5 instead of unity as Thomas proposed.This conclusion resulted from the observations by [5].
Based on that value, Kennedy [6] validated the derived formula against the longitudinal ventilation test results of the Memorial Tunnel Fire Ventilation Test Program (MTFVTP) [27][28] and concluded an under-prediction of the critical velocity for low fire heat release rates (<20 MW) and an over-prediction for high heat release rates.The proposed formula was also adopted into Annex D of NFPA 502 [44].Attention needs to be paid on the definition of   in the different quotations of this equation.The proposed formula by Kennedy [6] uses for   the distance from the base of the fire source to the highest point in the ceiling.The 2014 NFPA 502 Annex D uses the height of duct or tunnel at fire site ( ) for   .From a design perspective, it makes sense to use the full tunnel height for   as the height of a design fire is usually not specified, and it is a more onerous case to assume that the base of the fire is at the roadway.The 2023 version of Annex D [48] also references the 2014 Annex D equation, but with   inexplicably redefined as being 'height from the base of the fire to the tunnel ceiling at the fire site (not tunnel height  )' without providing guidance on a typical fire height for a design fire.

Non-dimensional Froude Model
Non-dimensional Froude Models refer to the Froude number as defined by equation (1) and not to the 'critical Froude number' or the Richardson number.The former neglects density variation and the latter two include density variation.
Non-dimensional Froude models are based on a non-dimensional critical velocity   * and HRR  ̇ * , which are defined as ̇ * =  ċ    2 √   (11) Both numbers use the reference velocity √ for the non-dimensionalisation.In analogy to the definition of the Froude number in equation (1), equation (10) defines the non-dimensional velocity, and the velocity  in equation ( 4) is also replaced by √.The tunnel cross-sectional area  is represented by  2 , and the temperature difference ∆ by the ambient temperature   .Relations between the nondimensional critical velocity and fire heat release rate are exclusively investigated in model experiments at a scale as small as 1:20, as demonstrated by Oka and Atkinson [13], Wu and Bakar [15] and Li et al. [18].In all these investigations, Froude similarity for scaling the velocity and fire HRR was assumed as follows where the index  refers to the values in the scaled tunnel model.In the work of Oka and Atkinson [13], as well as in Li et al. [18],  is the duct (tunnel) height, but in Wu and Bakar [15], tunnel height  is replaced by the hydraulic diameter of the tunnel  ̅ .The scaling as given by equation ( 12) leads to the following relationship for velocity and heat release rate: In the definition of the Froude number used by those authors, buoyancy is not included.Consequently, this approach inherently presumes that the temperature difference, and therefore the density deficit ∆ between the scaled and the real tunnel, is the same.Moreover, it also assumes that the geometric proportions between area, hydraulic diameter and tunnel height are the same between the scaled and the real tunnel.For example, if the width to height ratio of the real tunnel is not the same as in the scaled duct, the mass flow rate, and therefore the gas temperature after the fire, are not scaled in the right manner.The same applies for the ratio between hydraulic diameter and tunnel height, which are relevant for the ratio of inertia to buoyancy force.
Beyond the velocity and the HRR, the fire intensity (heat release rate per unit area) also needs to be scaled in the right manner to avoid altering the plume density deficit between the real and the scaled tunnel.Failure to attend to the fire intensity would violate the assumption that allows Froude scaling to even be considered.This also includes the need to scale the flame length and plume volume in the right proportions.
The test setups of the different small-scale tests [13][14][15]18] are very similar.Wu and Bakar [15], as well as Li et al. [18], used a fixed burner size for all tests in the same test tunnel.The fuel flow rate was increased to increase the fire HRR, which automatically increased the heat release rate by unit area (fire intensity).Oka and Atkinson [13] also investigated different burner sizes and arrangement and derived a set of equations to account for the different results.The trendline for the critical velocity dependence on the HRR is different for the individual tests performed, but always has a similar characteristic as given below.
The parameters ,  and  for the different results are listed in Table 1.As discussed above, and detailed by Stacey and Beyer [36], it is unlikely that all essential similarities are preserved in the scaled model, and hence unlikely that the trend lines provide relevant numbers when scaled to real tunnels.
Oka and Atkinson [13] and Wu and Bakar [15] used a water spray device to cool the tunnel wall near the fire source, which may have altered the heat loss from the smoke layer.
Any approach that does not include the density difference cannot be regarded as a physically reasonable approach for approximating density-driven flow behaviours with the confidence needed for design.

Practical physical problems with prior approaches
There are several observations to be made on the limitations of models previously applied.
The first is the dependence of the critical velocity on the heat release rate.That was a problem for Kennedy et al. [6] in matching the Memorial Tunnel data [27][28][29][30], and has been a problem since.Looking at the data relevant to road tunnel design, such a dependence could not be identified.Above a certain fire heat release rate, there is no further increase in critical velocity (see Figure 16).
Related to the problem of the heat release dependence is the second problem, which is the influence of the fire intensity.Fire intensity depends on the fuel type (heat of combustion) and the ability of the fuel to mix and burn before being advected downstream (burning rate).If the burning takes place while advecting downstream, the fire area has been increased and the intensity reduced.So, for a given fire type, an increase in the fire heat release rate involves an increase in the physical extent of the fire.For example, hydrocarbon pool fires that burn at around 2.25 MW/m 2 [27][28][29][30] do not obligingly burn at a higher rate (without changing the fuel type).To double the heat release rate, the pool area (area of the fuel source) also needs to be doubled.That has the effect of keeping the temperature in the hot plume relevant for the backlayering nearly constant with increasing HRR.A more intense fire (greater heat release rate per unit area) would result in a hotter plume relevant for the backlayering and greater propensity to backlayer.Experiments with a fixed burner area and varying gas flow are actually varying the intensity of the fire, introducing a variation that is not there for large real fires, yet the resulting formulae do not control for fire intensity.
The third problem is the extent of the fire, and of course, it is related to intensity.For a 20 m long fire (e.g.20 m long truck on fire or 20 m long pool fire), the most upstream part of the fire will have the most chance of contributing to backlayering, with subsequent parts being progressively less influential.It seems, from CFD investigations, that heat released more than one tunnel height or so beyond the fire front has less additional influence on the backlayering propensity (see Figure 15).So, a 20 m long fire (e.g.truck or pool fire) uniformly burning at 100 MW would have a similar critical velocity value as a 10 m long fire burning at 50 MW, and probably not too much different to that for a 5 m long fire burning at 25 MW, when assuming that the width of the fire stays the same.Especially for wide fires, it seems that the downstream part of the fire after the first few metres does not contribute much to the required critical velocity value (as is seen in the Memorial Tunnel test results [27][28][29][30]).The heat release rate in those tests was increased from 10 to 100 MW by keeping the width of the fire the same but adding fire pans in the longitudinal direction.Consequently, the critical velocity value did not change very much from 25 MW to 100 MW (see Figure 16).
The fourth problem is the effect of tunnel aspect ratio.There are a couple of proposed corrections to formulae to account for the ratio of the tunnel width to the height [12], [19][20].The reason that such corrections are required is clear.The recent (since about 1990) prior formulae all inherently use a temperature rise calculated by diluting the heat released by the fire into the entire flow along the tunnel (e.g.see equations ( 6), ( 7) and ( 9)).It is obvious that for most tunnel fires, only a part of the tunnel flow enters the fire plume in the first few metres of the fire (where backlayering will be determined).As a tunnel gets wider, that fraction of the flow entering the plume goes down (see Section 4.4).Starting the formulation with an inappropriate temperature rise and then seeking to correct the answer empirically based on geometry does not seem the best approach.As part of this work, a more realistic basis for buoyancy, roughly modelling the plume area, is sought.Natural break points then enter the formulation when the tunnel is so narrow that the plume fills most of the cross section.
The next section seeks to develop a new, physically reasonable model, which addresses the above shortcomings.

The Present Mixed-Convection Model
All the transport processes in tunnel fire events are governed by the equations of motion of continuum fluid mechanics.These are the balance equations for mass, momentum and thermal energy.The fluid motion contributes convective transport of the conserved quantities.Tunnel fires inherently include both forced and natural convection, where the latter is due to Archimedian buoyancy caused by the temperature (or rather density) differences in the flow field.The occurrence of both types of convection at similar strengths produces what is called mixed convection.This is particularly the case when analysing the equilibrium where the buoyancy-driven convection upstream of the fire (natural convection) interacts with the inertial force of forced convection.The governing non-dimensional equations of motion are formulated in a Cartesian coordinate system, with x1 (x) aligned with the tunnel axis in the flow direction, x2 (y) horizontal, and x3 (z) pointing in the direction of the tunnel height.For the present flow regime (steady state, Newtonian fluid and Mach number M < 0.3), the continuity equation for the incompressible fluid, the momentum equation and the thermal energy equation are given as respectively.The non-dimensional variables and unknowns are defined with the length scale , reference velocity  , the coordinates   and velocities   .The non-dimensional pressure  * = /( 2 ), and the non-dimensional temperature is defined by Θ * = ( −   ) ∆ ⁄ .The thermal buoyancy caused by density differences owing to temperature gradients   ∆ is expressed via the body force in the momentum equation.This replaces the Froude number by the ratio    2  ⁄ .That is, the Froude number does not exist in the governing equations for mixed convection and hence should not have been central to scaling of such flow regimes as discussed in Section 2.3.The non-dimensional numbers ,   ,  and  are the Reynolds, Grashof, Eckert and Prandtl numbers, respectively, defined for this problem as The Reynolds number is the ratio of inertial to viscous forces, where  is the characteristic velocity,  the characteristic length of the flow field and   the kinematic viscosity of ambient air.In the current flow problem, the hydraulic tunnel diameter  ℎ will be adopted as the characteristic length  in the Reynolds number.The Grashof number is the ratio of buoyancy to viscous forces.As the buoyancy is driven by the hot plume, the kinematic viscosity of ambient air   becomes the kinematic viscosity of the plume   .For an ideal gas,  = 1/.The characteristic length  is defined here as being the tunnel height from the base of the fire to highest point of the ceiling   .The term   represents the gravitational acceleration component in the Cartesian coordinate direction .For tunnels with a moderate (typical) longitudinal slope  and a transverse slope   ,   can be approximated as: That is,   in the x and y directions becomes very small so that the body force in these coordinate directions can be neglected for the current problem (strength of forced convection dominates the strength of natural convection).In a further simplification,   in the z direction will be declared just as .
The Eckert number is the ratio of the kinetic energy to the enthalpy difference due to the temperature increase, and the Prandtl number represents the ratio of the diffusivities for momentum and thermal energy of the fluid.In flow with strong enthalpy differences as compared to kinetic energy, such as the present one, the Eckert number is not important for the physical analysis of the flow.Yet, the Prandtl number appears in a product with the Reynolds number in front of the heat conduction term in equation (18).In flow dominated by convection, heat conduction is of minor importance, so that the Prandtl number does not appear as a quantity relevant for the flow modelling either.
When assuming that the Reynolds number is high (say Re > 10 5 ), the extra stress term with 1/Re in front must be modelled for turbulent flow.The ratio of Gr and Re 2 expresses the strength of natural convection relative to the strength of forced convection.It may be of the order of unity for mixed convection processes.The ratio   2  ⁄ is also known as the Richardson number Ri as per As noted above (Equation ( 5)) this has been referred to by prior authors as a "critical Froude number" (not to be confused with the Froude number that has been used in attempting to scale mixed convection problems), but is here referred to in standard fluid mechanics terminology.
The gas temperature  in the relation  = 1/ can be expressed as the sum of the approaching air temperature and the temperature increase due to the fire, i.e.

𝑇 = 𝑇 𝑎 + ∆𝑇 (25)
Based on the relations above, the hydraulic diameter of the tunnel cross section is selected as the length scale L in the Reynolds number.The physical flow processes, however, are influenced by other length scales also, which enter via boundary conditions.For the modelling of the mixed convective flow effects of the present flow, therefore, the Grashof number for the z direction is modified by incorporating a length scale ratio and a grade factor for sloped tunnels.The Reynolds and modified Grashof numbers are: The grade factor for sloped tunnels   is defined as being the ratio of the critical velocity in a sloped tunnel to the critical velocity in a flat tunnel, with the fire scenario and other geometrical tunnel parameters kept the same (see Table 5).
Formulating the Richardson number with the Reynolds number and the modified Grashof number as defined by equations ( 26) and ( 27) and solving for the velocity  gives: The velocity  is related to smoke propagation, but is not critical velocity until  and any empirical calibration is resolved for the equilibrium of the smoke layer at the front of the fire site.Equation ( 29) suggests a dependence of the smoke control related velocity on the kinematic viscosity of the gas.Since the flow regime analysed is turbulent, the kinematic viscosity relevant for the flow is the turbulent or eddy viscosity, which is orders of magnitude higher than the molecular kinematic viscosity of the gas.Even if there is a big difference in the dynamic viscosities of the plume and the approaching air, the effective viscosity (i.e., the sum of the molecular and turbulent viscosities) of the two regimes is very similar.The viscosity ratio in equation ( 29) is therefore dominated by the eddy viscosities, with similar values of   and   .This being the case, the viscosity ratio in (29) will be close to unity, thus simplifying the equation.To check this simplifying approximation, Figure 1 shows results from simulations (specified in Section 4) carried out with either a constant dynamic viscosity or a temperature-dependent dynamic viscosity.The results regarding smoke propagation are, as expected, very similar and confirm the conclusion.
(a) (b) As already discussed, it is presumed that, at the critical velocity, the strength of forced convection is of the same order as the strength of natural convection.According to (24), this requires that the Richardson number is close to unity.
The frontal area of the rising plume represents an obstacle for the flow past the fire.If the plume frontal area is small compared to the tunnel area, a big fraction of the flow goes around the plume, so that the momentum acting onto the plume is lower than it would be for a plume frontal area that is nearly as big as the tunnel area (wide fires).The plume in the latter case becomes more deflected by the upstream flow, and more mixing into the plume occurs, so that the upstream velocity required to prevent backlayering goes down.The influence of the plume width on critical velocity can be considered by the factor   defined here as: This correlation is a function of the plume frontal area     in relation to the tunnel cross-sectional area .The term 'min(1,     / )' makes sure that the plume frontal area does not exceed the tunnel cross-sectional area.The parameters   and   are empirical constants that will be calibrated to measurements and a validated CFD model.
It is assumed that the effective temperature difference ∆ is a function of the fire intensity (heat release rate by unit area), and the cooling due to the approaching air.The higher the approaching air mass flow and the bigger the plume frontal area, the more the hot plume gases are mixed and cooled down.The higher the fire intensity, the higher is the density deficit and thus the buoyancy force opposing the forced convection at the initial plume regime (fire front).It is also acknowledged that an increase in the heat release rate caused by extending the fire source in the longitudinal direction contributes less than the first part of the fire to the strength of natural convection at the front of the fire.This reasoning leads to the following representation of the effective temperature difference: where  ̇ is the total fire heat release rate,   is the fire intensity (total fire heat release rate per unit area),   the width of the fire,   the tunnel height from the base of the fire to highest point of the ceiling,  is the heat reduction due to imperfect combustion and/or radiation,  is the tunnel cross-sectional area,  is the air velocity onto the fire, and  ℎ   is describing the fraction of the length of the fire (as a proportion of   ) that contributes to the natural convection at the front of the fire.The chosen function for the effective fire length   is based on the ratio of the longitudinal extent of the fire   =  ̇/(     ) (see equation (34) below) to the tunnel height from the base of the fire   and the empirical constants   ,   and   as follows: It was chosen simply as it has the desired form, and is able to be calibrated to observations.
Taking the minimum of  ̇ and   times the effective fire area      ℎ   ensures that the modelled effective heat release rate does not exceed the actual total fire heat release rate if    ℎ   is larger than   .The plume area interacting with the approaching air is taken as   (  +   ).The basis for this is an assumption that plume shapes will be geometrically similar, growing in width with grade  as the plume rises.As the plume width increases with height, the plume frontal area interacting with the ambient air will then vary with height from the base of the fire (conservatively tunnel height for design purposes).If a tunnel is wider than the mean plume width, the additional air passing the fire to the side of the fire does not dilute the effective temperature difference and therefore does not decrease the buoyancy force.However, the plume frontal area cannot be bigger than the tunnel area itself, which requires that the minimum of  and   (  +   ) be taken as representing the plume frontal area interacting with the ambient air.Deluge or water mist systems reduce the convective heat release rate and therefore the effective fire intensity, which can be considered via the parameter .
With the definition of the fire intensity the longitudinal extent of the fire can also be expressed as where  ̇ is the total fire heat release rate.
The definitions and assumptions above, finally lead, together with equations ( 29) and (32), to the following system of implicit equations for the critical velocity For flat tunnels, the grade factor   is assumed to be unity.The characteristic length scale   is defined as the height from the bottom of the heat source to the tunnel ceiling (conservatively the tunnel height for design purposes),  ℎ is the tunnel hydraulic diameter, and   is the fire width (most likely between 2.0 and 3.5 m for road tunnels).The fire intensity   depends on the type of the fire and the fire scenario.A typical value for a pool fire is 2.25 MW/m 2 [27] and is considered as a conservatively intense fire scenario.A higher fire intensity might be possible with different fuels or, in the experimental case, with burners introducing fuel at high rates in small areas.The empirical factors and constants in the semiempirical correlation for calculating critical velocity ( ℎ ,   ,   ,   ,   ,   and  ) might be variable, depending on the fire characteristics.For pool fires, it has been found that the values provided in Table 2 accurately represent the physical phenomena (see Section 5).These values were first selected as fitting the Memorial Tunnel test data and results from a validated CFD model (see Section 4), and were subsequently found to require no alteration to fit with the available small scale test data.The validity of the parameters listed in Table 2 is limited to the cases analysed or validated against (see Section 4 and 5).However, with no test data conflicting with the empirical results, they likely have wider applicability.How wide, is subject to further investigations.If water mist or deluge systems are not considered,  = 0. Values for  usually lie between 0.2 and 0.4 [54].

CFD Validation, Grade Factor and Parameter Variation
To properly calibrate equations ( 35) and (36) together with ( 29) and ( 32) against the Memorial Tunnel data [27][28][29] it is required to understand the grade factor, as the Memorial Tunnel had a 3.2% downgrade.There are no critical velocity data from full size tunnels with varying grade, and scale models are not considered appropriate for the purpose.Further, the available data in the literature [6,14] provide very different answers regarding the slope effect.It cannot be reliably said if one or none of them are correct.To reference the results back to a flat tunnel requires an analysis of the grade factor   between a flat tunnel and the 3.2% downgrade.Clearly the same information on grade effects will also be useful in design of tunnel ventilation.For this, a CFD study using the simulation software ANSYS Fluent version 2021 R1 [31][32][33] has been carried out.The first task was to carefully validate the CFD technique against Memorial Tunnel data, to be confident that simulations of slope effect would be meaningful.That validation has already been done and was documented in a separate paper by Beyer and Stacey [51].A short summary of the different fire cases and results will be provided here (Section 4.1 and 4.2).Table 3 provides an overview of the simulation methodology and boundary conditions used.More comprehensive description of the Memorial Tunnel tests, the CFD model including all boundary conditions, and the methodology used, can be found in the work by Beyer and Stacey [51].The CFD model validated for the Memorial Tunnel was then used to analyse the same model but with different tunnel slopes (Section 4.3).In a further step, a few parameters like width of fire pan, shape of tunnel profile and area, longitudinal extent of fire, and fire intensity were adapted and the results analysed to test the derived equations for a broader range of parameters (Section 4.4).Changes to the original validated CFD model, and the results of the additional CFD study, are documented here.

Validation Cases
The Memorial Tunnel has a total length of about 850 m and a typical tunnel profile with a curved ceiling (see Figure 4).Jet fans were installed in banks of three, up-and downstream of the fire, to adjust the tunnel air velocity during the fire tests.Instrument trees ("loops") were located adjacent to the fire and further upand downstream to measure velocity and temperature profiles and further flow parameters along the tunnel.The measurement equipment and data acquisition units created a fair blockage and were therefore implemented in the CFD model, based on sketches and the reported blockage area.Location of the measurement loops, as well as of the individual sensors in the tunnel profile, are shown in Figure 3.
The Memorial Tunnel fire tests included pool fires with nominal heat release rates (HRR) of 10, 20, 50 and 100 MW.Each test used the appropriate fire pan, or a combination of the fire pans, as illustrated in Figure 2 and listed in Table 4.To validate the CFD model for low as well as high heat release rates, one test with a nominal HRR between 10 MW and 20 MW, one at 50 MW, and a third one at 100 MW were selected.The CFD methodology and model boundary conditions are listed in Table 3. Essential case parameters are listed in Table 2.The flow parameters and measured HRR values listed are averages over a period when the flow and smoke layer were fairly constant.This allows the comparison of the measured values with a steady-state simulation in a reasonable way.More details about the geometry, test setup and simulation methodology are provided by Beyer and Stacey [51].

Thermal effects
Incompressible ideal gas law.Piecewise-polynomial function for the temperature dependent specific heat capacity for the individual species.Thermal dependency of the thermal conductivity and the dynamic viscosity for the individual species were estimated based on Sutherland's law with three coefficients, in combination with a massweighted mixing law.

Representation of fire
Eddy-dissipation combustion model.

Thermal radiation
Not simulated.Instead, a radiative fraction approach as explained in the literature referenced [35,51] has been adopted, capturing the loss of thermal energy in a gross sense.

Fuel
Ideal combustion of fuel oil with a chemical formula of C19H30.Additional mass source term for ideal combustion and the radiative fraction approach (see [51]) are incorporated.

Tunnel inlet
Mass flow inlet with air temperature that was measured upstream of fire (~Loop 207).Mass flow is based on Loop 214 (mass fraction of O2 = 0.23, N2 = 0.77 and 0 for all other species).See Table 4

Fire pan
Mass flow inlet with ambient temperature (~Loop 207) on surface located on top of fire pan.Mass fraction of O2 and N2 was 0, mass fraction of fuel and other species according to adopted fire HRR, radiative fraction and combustion efficiency.Combustion efficiency for extra mass source is assumed to be 95% [27,28].See Table 4 and Beyer & Stacey [51] for further details.Radiative fraction 2 (-) 0.21 0.24 0.17 1 Measured heat release rate corrected by the CO to CO2 ratio to account for non-ideal combustion [27], [28] 2 Ratio between measured bulk temperature and the theoretical bulk temperature rise resulting from the total HRR is assumed to be the radiative fraction. 3Value as stated in Section 8.8.4 in the comprehensive test report [27], and in the paper from Kile & Gonzalez [29].The fraction of the tunnel area blocked is about 17%, as also reported in the same references. 4Tunnel height and free tunnel cross sectional area are calculated assuming that the removed ceiling had been flat.It may in fact have had a slight camber (see Figure 4) such that the true height and area inferred from the dimensions were slightly higher than tabulated (possibly 1%).

Validation Results and Discussion
For the three simulated validation cases listed in Table 4, Figure 5 to Figure 7 compare the predicted bulk temperature along the tunnel axis and the temperature profiles over the tunnel height against the corresponding test results [27][28].Multiple thermocouple readings at a given elevation are averaged in the temperature profiles at the different measurement loops.The representative elevation area is a horizontal slice of the tunnel cross-section around the elevation instruments as specified in the test reports [27][28].Identical horizontal slices have been used for evaluating the CFD simulations.For a better understanding of the temperature profiles, a contour plot of the temperature through the middle of the tunnel around the vicinity of the fire is also given for each case.
In all cases, the temperature profiles downstream of the fire were predicted with reasonable accuracy with the applied simulation methodology.The higher temperature readings of the test data at the upper part of the temperature profile close to the fire (Loops 304 and 303) might be related to the influence of thermal radiation on the sensors (see Figure 5 and Figure 6) as noted in the test reports [27][28].
Besides the deviation of the bulk temperature close to the fire (influence due to thermal radiation on the temperature sensors during the tests), the bulk temperature at Loop 202 also shows a slight deviation, especially at higher downstream temperatures (higher HRR).Loop 202 is located in the smaller cross section near the south portal with the reduced ceiling height [51].The step change in the crosssection area is just 1.5 m upstream of Loop 202.This step will have created a separation zone and recirculation close to the ceiling.The difference in bulk temperature seen at Loop 202 may be related to the difficulty in turning point measurements of velocity and temperature into bulk averages, for such a non-uniform flow.
The velocity in the tunnel during the Memorial Tunnel tests was adjusted using the jet fans installed in the tunnel upstream of the fire.These jet fans were not speed controlled.To change the air speed in the tunnel, a jet fan could be switched on or off, making the speed control fairly coarse.The buoyancy force and flow resistance of the 11 MW fire are quite low, certainly much lower than for a 50 MW or 100 MW fire.With the coarse control on jet fan forcing, the low resistance made it difficult to achieve a steady velocity where backlayering was just prevented.For the tests with low heat release rates, two active jet fans were insufficient to stop smoke propagation upstream of the fire, and three jet fans were too strong, so that the approaching air velocity was higher than required for preventing backlayering.The point at which backlayering was just prevented was analysed during the velocity change after turning the third jet fan on or off.This made it difficult to accurately determine the critical velocity for low heat release rates.The scatter in the reported data at low HRR reflects this (see Figure 16).
As discussed in detail in Beyer and Stacey [51], the results below are evaluated at times when smoke upstream of the fire was recently moved downstream and reached almost steady state conditions for a short period.Some of the temperature readings upstream of the fire show slightly increased values which may possibly be explained by the thermal inertia of the temperature sensors only recently out of the smoke layer.Especially for the lower HRR cases, Loop 304 just a few metres downstream of the fire pans gives higher temperature readings in the upper part of the tunnel than predicted.This may be related to sensor hardware thermal inertia and a transient behaviour of the flames and hot plume moving up and down and heating up the sensors in the upper part of the tunnel.A more comprehensive discussion about the smoke propagation during the tests and the deviations of the temperature profiles for validation cases is given in Beyer and Stacey [51].Based on the observations provided by Beyer and Stacey [51] and above, and acknowledging a few uncertainties in the temperature readings, the smoke propagation and temperature distribution can be predicted within a reasonable accuracy by the simulation techniques and procedures applied.Most importantly for the present purpose, the smoke extent matches the tests very well for the same conditions.However, it is important to check other parameters, to see that the critical velocity is not correct by an 'accident' of compensating errors that may not compensate each other in other scenarios.The temperature profiles away from radiation effects are also predicted quite well, giving confidence that the analysis is a reasonable representation of the plume and smoke layer dynamics.That is; the technique is a reasonable tool for exploring trends in critical velocity, around the validated cases.
The trend of interest here is the effect of tunnel slope on the critical velocity, and as noted at the start of the section, the CFD work reported above was all with the primary purpose to obtain a validated simulation methodology for exploring the tunnel slope influence, with a secondary objective to look at the influence of other geometric and fire parameters.The outcomes are described in the following sections.

Grade Factor
The validated CFD-model described by Beyer & Stacey [51] and summarized earlier was used to identify the grade factor   as defined in Section 3.For this purpose, the jet fans and all the instrument trees and other obstacles have been removed from the model, first to reduce the complexity and computational expense, and second, to simplify assessment of when critical velocity is achieved, which is made more difficult with the obstructions.As the models with different slopes are compared in a relative sense, these geometric modifications are not thought to be significant, and assessing critical velocity consistently across all tunnel slopes was seen as far more important.Grade factors have been analysed for a 50 MW fire (validation case 2), which is believed to be the most relevant fire scenario for road tunnels.
To accurately reproduce the point where upstream smoke propagation is just prevented, a simple feedback controller for adjusting the upstream velocity was implemented in the simulation model.Upstream velocity is automatically increased or decreased so that the tip of an upstream smoke layer is confined at the front edge of the fire pan.This was done by tracking the maximum occurring temperature at the tunnel cross sectional area in front of the fire front ('tracking' area) and adjusting the approaching air velocity in increments proportional to the temperature ratio (between maximum temperature at the 'tracking' area and the approaching air temperature).Once the temperature at the tracking area is equal to the ambient air temperature, the approaching air velocity was decreased by small, fixed increments.The inlet velocity stays constant once the temperature difference between the 'tracking' area and ambient air is >10 K and <40 K.By allowing sufficient computational time, the same criteria could be applied precisely to each slope case, eliminating any human judgement based on CFD images, and hence refining the relative grade effect.The outcome of that CFD study is summarized in Table 5 below.The conclusion from Table 5 is that, for most design purposes, the grade factor can be ignored.For the practical range of road tunnel longitudinal slopes, the effect of slope appears to be much smaller than tolerances in calculating the required ventilation to achieve critical velocity, or in measuring the result on commissioning.That is; for practical purposes, the factor   may be deleted from equation (35).Eliminating grade factor is a significant departure from prior practice, but the basis of the above assessment that grade effects are not significant for typical road tunnels is considered more robust than the prior trends based on very small models or tests with methane layering [6,14].It seems likely that no grade factor is required for transverse gradient (crossfall), but that has not been confirmed.

Parameter Variation
Based on the validated CFD model and the presented simulation methodology [51], geometric parameters like width of the fire, tunnel height, tunnel profile, longitudinal extent of the fire, as well as the fire intensity, were varied to better understand the effect of those parameters on the critical velocity.For most scenarios, the same velocity/backlayering controller described in Section 4.3 was used for the parameter variation study.The results in turn were used to verify the presented physical model (equations ( 29), (32) (35) and (36)) and the empirical parameters as listed in Table 2.

Fire Width
The first parameter to be varied was the width of the fire pan.Fire width effects were studied based on the simplified Memorial Tunnel model that was used for the grade factor study.The tunnel slope of -3.2% was kept the same for the different scenarios.The original 50 MW fire pan according to Figure 2 was varied in width as illustrated in Figure 8, maintaining the fire pan length of 6.1 m and the fire intensity of 2.43 MW/m 2 according to validation case 2 (see Table 4).The upstream air temperature and therefore the wall temperature in the simulation was 6.13°C.Further parameters and results of the simulations are listed in Table 6.
As already discussed, the wider the fire, the greater is the interaction with the approaching air.This has two effects.First, the increased plume frontal area causes more mixing and secondly, it provides less space for the upstream air to go around the plume (higher momentum onto the plume frontal area).Both have the effect of decreasing the required critical velocity.As an outcome of that study, it also seems that the critical velocity does not change if the fire width (fire blockage) goes below a specific value (here 2.5 m) when the fire intensity is kept the same.
The proposed physical model shows a reasonably good agreement with the simulation results but tends to underestimate critical velocity for very narrow fires by about 3%.Table 6 and Figure 9 provide a comparison between the simulation results and the semi-empirical model.29), (32) ( 35) and (36).

Tunnel Height and Aspect Ratio
The influences of the tunnel height and the aspect ratio were explored for a flat rectangular 3-lane tunnel with a width of 12 m.This was done by using the simplified Memorial Tunnel model parameters, replacing the Memorial Tunnel profile by a rectangular profile as depicted in Figure 10.The tunnel slope was also removed.The rectangular pool fire surface with a length of 20 m and a width of 2 m was on the floor and located 245 m downstream of the entrance portal.Fire intensity was 1.25 MW/m 2 which resulted together with the fire surface in a total HRR of 50 MW1 .The upstream air temperature and wall temperature in the simulation was set to 21.4°C.Only the tunnel height varied from 3 m up to 12 m but all other parameters were kept the same.Further parameters and results of the simulations are listed in Table 7 below.
The simulation showed that the critical velocity increases from a tunnel height of 3 m (aspect ratio of 1:4) to 12 m (aspect ratio of 1:1) by 55%.
The proposed physical model accurately predicts the trend of the critical velocity with tunnel height and aspect ratio, but underestimates the values by up to 14%.This may be related to the rectangular tunnel profile and the reduced velocity at the corners.The momentum onto the fire and potential backlayering at those areas is reduced so that the smoke starts to propagate upstream along the corners.That physical effect related to a rectangular tunnel profile is not captured in the proposed equations.For practical purposes, it is suggested either to increase the resulting values by 7% to 15% for rectangular tunnel profiles, or simply to accept some backlayering in the corners.Table 7 and Figure 11 provide a comparison between the simulation results and the semi-empirical model.29), (32) ( 35) and ( 36).

Variation of Fire Dimensions
For this study, the Memorial Tunnel profile of the simplified validation model used for the grade correction has been replaced with a similar profile but with 3 lanes as shown in Figure 12.The total tunnel length was shortened to 550 m and is downgrade at 5%.For the pool fire, a rectangular surface was defined and placed on the floor.The fire front is located 300 m downstream of the entrance portal.To understand the influence of the fire dimensions, the width and length of the fire source area was varied by keeping the fire intensity of 2.5 MW/m 2 constant (fuel mass flow per m 2 was constant).Depending on the fire source area, the constant fire intensity resulted in different total fire heat release rates.In a further step, some of the cases were analyzed with a downgrade slope of 1.7%.The upstream air temperature and wall temperature in all simulation cases was set to 13.85°C.Further parameters and results of the simulations are listed in Table 8 below.
The results for the bigger tunnel profile suggest that there is no difference in critical velocity between 1.7% and 5% downgrade.These results support the conclusion made in Section 4.3, that for practical purposes, the factor   may be taken to be unity.These simulations confirm that for the same fire intensity, a wider fire reduces the required critical velocity but a fire more extended in longitudinal direction increases the critical velocity as indicated by equation (32).
The effect of the variation of the fire dimensions is also captured accurately by the derived physical model, especially if the grade correction is neglected (  = 1) as advocated.Only the case where the fire is wider than it is long shows an overprediction (of 8%).That overprediction can be accepted as such a design fire scenario may be an edge case anyway.Table 8 and Figure 13 provide a comparison between the simulation results and the semi-empirical model.29), (32) ( 35) and (36).

Variation of Fire Intensity and Fire Length
In a final step, the validated simulation methodology [51] was applied on a TBM (Tunnel Boring Machine) tunnel profile with 3 lanes and a false ceiling (to understand the effects of varying fire intensity and fire length on critical velocity.In this case, the tunnel is flat and has a total length of 530 m.The rectangular fire source with a constant width of 2.0 m was placed on the floor, with the front of the fire source located 300 m downstream of the entrance portal.Fire intensity was varied by keeping the total fire heat release rate constant with varying fire length.In a second step, fire intensity was varied while keeping the fire length constant.Note, this variation acknowledges that the critical velocity mainly depends on the physical dimensions of the fire and the fire intensity and not on the fire heat release rate itself.The upstream air temperature and wall temperature in all simulation cases was 20.0°C.Further parameters and results of the simulations are listed in Table 9 below. The results for this tunnel profile confirm that critical velocity increases strongly with fire intensity for the same total heat release rate.As before, the critical velocity increases slightly with the fire length (as modelled in equation ( 34)).
That behavior was also replicated with the proposed physical model, as listed in Table 9 and illustrated in Figure 15.The maximum deviation between the CFD and the proposed model is 3%.29), (32) (35) and (36).

Validation of Mixed-Convection Model Against Real Test Data
The behaviour of the proposed physical model with varying parameters has been shown above to agree quite closely with results obtained from a validated CFD model.This section documents the agreement of the proposed model with critical velocity values obtained from longitudinal ventilation full-scale test results of the Memorial Tunnel Fire Ventilation Test Program [27], [28][29][30], and also the agreement with values obtained from small-scale tests as recorded by several researchers [13], [15], [18].This validation against the small-scale tests was done without attempting any scaling, but simply by using the real dimensions of the test ducts, as well as the real measured velocities and heat release rates.As discussed earlier, there is no approach to scaling the formula that does not involve major assumptions about at least some of the important parameters.Using the real parameters ensures that the physical phenomena are retained, and not 'assumed away' and therefore allows a verification that the mathematical model captures the important physics appropriately.

Memorial Tunnel Fire Ventilation Test Program
Figure 16 compares the critical velocity values from the Memorial Tunnel tests [27], [28][29][30] with the values given by the derived physical model (equations ( 29), (32) (35) and (36)).The scatter in the plotted data points produced from the proposed equations is due to the fire intensity being calculated from the reported HRR and the actual pan area.That is; they include some element of the experimental scatter within the original Memorial Tunnel data.For the proposed model, a fixed radiative fraction of 0.21 was adopted, as being the average observed from the validation cases (see Table 4).The grade factor was set according to Table 5.
Apparently, there are two versions of the ASHRAE paper from Kile & Gonzalez published in 1997 [29], [30].As it is not clear which version presents the most reliable critical velocity values, both versions are plotted in Figure 16 ( [29]: triangles, [30]: circles).The data points from the first version [29] are the same as published in the comprehensive test report [27].However, the data points from the second version [30], seem to be more aligned with the results presented by Kennedy [6].The difference in the critical velocity values for high HRRs (>50 MW) is minor and probably less of a concern, but the difference is quite big for low HRRs (8 to 16 MW).
As noted earlier, the air velocity in the Memorial Tunnel tests was adjusted by jet fans without variable speed drives.The change in air velocity caused by starting an additional jet fan or stopping a running jet fan was too large to allow an accurate determination of the velocity where backlayering of hot smoke was just prevented at low heat release rates.This is likely to be the reason for the high variation of the reported critical velocity values at low heat release rates.The outlier at low heat release rates (critical velocity > 3.2 m/s) was re-evaluated by looking at the raw data and 3 m/s may be a more reasonable value (see squares in Figure 16).
With those notes on the Memorial Tunnel data, we return to the comparison against the validated CFD model.The Memorial Tunnel results show that for the 50 MW fire (Case 2 in Table 4) an air velocity of 3 m/s onto the fire was sufficient to prevent backlayering (see Figure 16).That is also consistent with the grade factor study as given by Table 5.The results of the 100 MW case (Case 3 in Table 4) indicate that a velocity of 2.85 m/s was high enough to control upstream smoke propagation (Figure 16) but was perhaps too low to completely avoid backlayering (see Figure 7).A value of 3.0 to 3.1 m/s may be sufficient to prevent backlayering for that 100 MW case.Also, the validation case 1 (Table 4 and Figure 5) shows that the velocity of 2.9 m/s onto the 10 MW fire was higher than critical velocity.Consequently, a value around that or slightly less for a 10 MW fire may be sufficient.That is; across the range of HRR of the Memorial Tunnel tests, the proposed mixed convection model is consistent with the original test results and the CFD validated against those results.
When acknowledging some uncertainties in the published data for low HRR, as discussed above, the agreement of the proposed physical model with the Memorial Tunnel data in Figure 16 can be seen as acceptable.
Figure 16 includes the predicted critical velocity according to Kennedy [6] and those versions of the 2014 NFPA 502 equation published in 2014 [44] and in 2023 [48] for comparison.The difference between the various versions is discussed in Section 2.2.

Small-Scale Test Data
All but one of the test rigs giving results used for comparison below, used circular burners.The current formula was based on rectangular heat sources, with the Memorial Tunnel tests in mind.For the comparison purposes, the circular smallscale burners have been taken as equivalent rectangular burners as follows.From the proposed model, the width is the most important dimension of the fire.Provided that the maximum width (the diameter) falls within the effective length, which it does in all cases, the fire width can be taken as the burner diameter.In order to preserve both the fire intensity and the HRR, the burner area needs to be unchanged, so the fire length is then taken as the burner area divided by the width.For circular burners of diameter   , this makes the fire width   , and the fire length   =   /4. Figure 17 and Figure 18 compare the critical velocity values from different small-scale tests with the values predicted by our mixed convection model.As previously discussed, and also elaborated by Stacey and Beyer [36], the fire intensity is a measure for the density deficit, with the buoyancy force increasing with the fire intensity.The burner size was kept the same during the referenced model tests [13], [15], [18] (for the individual test setup), which means that the fire intensity (MW/m²) in these models increases in proportion to the HRR.This is the reason why the critical velocity in the scale model tunnels climbs steeply with increasing HRR.
As discussed by Wu & Bakar [15], with further increase in HRR, the plume gets elongated in the downstream direction and the critical velocity become independent of the HRR once the intermittent flame zone reaches the tunnel ceiling.That is seen by the kink in the measured data, which can be observed earlier in the smaller and more confined test tunnels.If the tunnel gets higher and wider, the position of the kink moves towards higher heat release rates (see Figure 17 (a) vs. Figure 17 (b)).In all small-scale tests considered here, the kink occurs when the average temperature after the fire exceeds 330°C to 370°C.This behaviour could be reproduced with the proposed model by restricting the effective temperature difference ∆ in equation ( 35) once the average temperature after the fire exceeded 330°C to 370°C (see green dashed lines in Figure 17 and Figure 18).Whatever the precise location of the kink is, the mechanism is believed to be that the fire intensity becomes nearly constant, with the effective area of combustion extended downstream in order to get all the injected fuel burnt.If that is so, as indicated by Wu and Bakar [15], the proposed model predict a constant critical velocity.However, it appears that this is an artefact of the burner types used in the very confined space of small-scale tunnels, and therefore not so relevant for real tunnel fires.
Oka and Atkinson [13] as well as Wu and Bakar [15] used a water spray device to cool the tunnel walls near the fire source.As noted by Oka and Atkinson [13], the wall temperature has significant impact on the critical velocity.This can also be observed when comparing the test results of similar test tunnels (Tunnel A in Figure 17 and Tunnel B in Figure 18).The cooling effect of the water spray system could be reproduced with the proposed equations by considering a reduction of the convective heat release rate ( = 30%) and hence further restricting the temperature difference ∆ in equation (35) (see blue dotted lines in Figure 18).The value of  = 30% was derived by comparing Tunnel A in Figure 17 with Tunnel B in Figure 18 and then kept the same for the wider Tunnels C (Figure 18 (c)) and D (Figure 18, (d)).For the small-scale tests conducted by Oka and Atkinson [13], a value of  = 10% gave agreement when applied in the proposed equations.
Acknowledging the uncertainties around thermal restriction due to flame/burner behaviour and the use of water cooling, our proposed model shows a reasonably good agreement with small scale test results, as illustrated in Figure 17 to Figure 18, and therefore also represents the physical phenomena in small-scale tunnels adequately.

Practical estimation of critical velocity
The proposed physical model is seen above to accurately model the Memorial Tunnel outcomes, and also accurately follows small-scale test data until the flame/burner behaviour in the very confined space affects the results.However, equations ( 35) and ( 36) require knowledge of the real fire dimensions (at least the fire width) as well as the heat release rate or fire intensity.A tunnel design cannot be based on a particular fire geometry, but must consider a more general case, so there is a need to consider nominal inputs that would lead to reasonably 'safe' values of critical velocity, without unnecessary over-prediction which also increases risk and cost if used as a design basis.

Design assumptions
It is unlikely that a truck will become completely sideways in the tunnel.The width of a truck is 2.5 m.Spilt liquid fuel will run to the gutter following the cross-fall in combination with any longitudinal gradient component.It would start as a splash zone and perhaps grow somewhat as it flows, depending on the road surface slope.It is felt that 2.5 m would cover that also, noting that wider fires reduce the critical velocity.
The fire intensity can be taken as that of a hydrocarbon pool fire, adopting 2.25 MW/m 2 as the figure.As shown by the Memorial Tunnel tests (see Table 4), 0.2 seems to be a reasonable value for ε (radiative heat fraction), and η (heat fraction lost to water suppression) is taken as zero (no active suppression).Section 4.3 demonstrates that grade factor   can also be taken as 1.0 for tunnels with typical gradients.The fire can conservatively be taken as starting at roadway level, and so the buoyant convection scale   can be taken as equal to the tunnel height  .The empirical numbers can be adopted as listed in Table 2, as it was shown in Section 5 that those values provide good approximation to real tests.
Making a sea level assumption, air density is 1.2 kg/m 3 with an ambient temperature of 294 K.

Outcome of design assumptions for typical tunnel profiles
To explore the variation of the critical velocity for typical two to three lane tunnels, the equations ( 35) and (36) were evaluated for different tunnel profiles with the assumptions as noted in Section 6.1.
Figure 19 and Figure 20 show the range of critical velocity values for a wide variation of TBM and rectangular tunnel profiles.Table 10 lists the geometric parameters for the TBM tunnel profiles analysed.
When considering that the design HRR for road tunnels typically lies in the range of 30 MW to 60 MW, the 'design' critical velocity in the TBM profiles varies between 2.82 m/s and 3.63 m/s.For the same HRR range, the critical velocity in the rectangular profiles varies between 2.47 m/s and 3.50 m/s.As discussed in Section 4.4, the critical velocity values in rectangular tunnels are underestimated by about 10%.Of course, in locking in the adoption of a 2.5 m wide fire and assuming that the base of the fire is at the floor, the model results will no longer reproduce the Memorial Tunnel tests, for which the fire pans were elevated and 3.66 m wide.

Simplified equation for design purposes
The design HRR for road tunnels is typically >10 MW, which allows the strong HRR dependency of critical velocity at low HRR to be ignored, and therefore the "min()" function in equation (36) to be deleted.The specific heat capacity can be taken as 1.007 kJ/(kg•K) over the range of interest.With those assumptions and the substitutions according to Section 6.1 and Table 2, the critical velocity as given by equations ( 35) and (36) only depends on the total HRR and some geometric parameters (tunnel height, tunnel area and hydraulic diameter) as shown by the simplified equations ( 37) and (38) In equations ( 37) and ( 38), the variables have the units noted in the nomenclature section.Several parameters also have units, as tabulated below by their numerical value.
Table 11.Units of parameters used in equations ( 37) and (38).The parameters resulted from a reduction of ( 35) and ( 36) after substituting the noted design assumptions and the empirical parameters listed in Table 2. Therefore, the physical meaning of the parameters listed is related to the substitutions and combinations of the different numbers.

Parameter
Unit Location 1.476 -First factor in eq. ( 37) 2.64 m First factor in eq. ( 37) The applicability of the simplified equations ( 37) and ( 38) is limited to the real test data and simulations on which the validity of the empirical parameters listed in Table 2 relies.However, with the broad agreement between the full equations and a very wide range of experiments, and the fact that the simplification of the equations involved conservative parameter choices, the simplified equations are seen as useful for road tunnels.

Conclusions
A fresh approach to accommodating the physics of the zero-backlayering tunnel fire flow field in a mathematical model has yielded equations with reasonable predictive capability for critical velocity.The model is based on non-dimensional numbers characterizing the gas flow field in the tunnel with mixed convection.The model was validated against CFD and experimental measurement data.The work has produced the following outcomes, the first three of which are novel, with the last being reconfirmed.
• Explained by geometry and fire intensity considerations, heat release rate beyond a certain value (about 10 to 15 MW for realistic fires in real road tunnels) causes no further significant increase in critical velocity, matching the Memorial Tunnel results.• Tunnel aspect ratio (width to height ratio) no longer requires an additional adjustment, being handled natively by the current formula.• Corrections for typical tunnel grades are small enough to be ignored, such that no grade correction factor is required in design.• It is clear why scaling up from small models by using Froude number (ignoring density, fire intensity etc.) leads to unreliable numbers.
Recommendations on design parameters are provided to reduce the complexity of the proposed equation.Based on those substitutions, it seems that a value between 2.8 m/s and 3.8 m/s is sufficient to prevent backlayering in typical road tunnels for typical design fires with heat release rates between 30 MW and 60 MW.Within the limitations noted, the proposed equations are ready for application in tunnel designs, with the cautions, as before, that critical velocity may not be so critical to design, and that application to specific geometries and fire scenarios requires thought.
Further comparison of the proposed model with any full-scale data and relevant test data would be welcome.The determination of reasonable η values to account for the effect of fixed firefighting systems on the critical velocity value may be the subject of future work.

Figure 1 .
Figure 1.Simulation results of validation case 2 according to Table 1 in [51], with constant viscosity (top) and temperature-dependent viscosity (bottom).(a) Contour plot of dynamic viscosity in a vertical plane through the tunnel centreline.(b) Contour plot of gas temperature (in Kelvin) in a vertical plane through the tunnel centreline.Temperature and viscosity are clipped to 400°C and 4×10 -5 Pa•s, respectively, to better visualise the differences in the upstream temperature layer.The clipped area is shown in magenta.

Figure 2 .
Figure 2. Computational mesh around the fire pans and adjacent instrument loops.

Figure 3 .
Figure 3. Instrument loop locations, with instrument types at each loop.Sketch from the Memorial Tunnel report [27-28].

Figure 4 .
Figure 4. Memorial Tunnel cross section with dimensions, looking north [27-28].The dimensioning did not include what may have been a slight camber in the ceiling, which if included, would increase tunnel height and area by perhaps 1% compared to the values used here.

Figure 5 .
Figure 5. 11 MW validation (Case 1 from Table 4) with combustion model.Comparison of simulation results with test results for test 606A [51]: (a) Bulk temperature along the tunnel axis, (b) Temperature profiles over the tunnel height for different loops in the vicinity of the fire, (c) Contour plot of the temperature in °C through the middle of the tunnel and fire pans.Temperature is clipped to 150°C for a better visualisation of the temperature layering.

Figure 6 .
Figure 6.54.3 MW validation (Case 2 from Table 4).Comparison of simulation results with test results for test 610 [51]: (a) Bulk temperature along the tunnel axis, (b) Temperature profiles over the tunnel height for different loops in the vicinity of the fire, (c) Contour plot of the temperature in °C through the middle of the tunnel and fire pans.Temperature is clipped to 400°C for a better visualisation of the temperature layering.

Figure 7 .
Figure 7. 105 MW validation (Case 3 from Table 4).Comparison of simulation results with test results for test 615B [51]: (a) Bulk temperature along the tunnel axis, (b) Temperature profiles over the tunnel height for different loops in the vicinity of the fire, (c) Contour plot of the temperature in °C through the middle of the tunnel and fire pans.Temperature is clipped to 800°C for a better visualisation of the temperature layer.

Figure 8 .
Figure 8. Schematic of Memorial Tunnel profile used for fire pan width variation.

Figure 9 .
Figure 9. Critical velocity values in the Memorial Tunnel profile for different fire pan widths, but constant fire intensity and fire pan length obtained from the validated CFD model and the semi-empirical equations (29), (32) (35) and(36).

Figure 10 .Table 7 .Figure 11 .
Figure 10.Schematic of the tunnel profile used for the tunnel height and aspect ratio variation.

Figure 12 .
Figure 12.Schematic of the tunnel profile used for the fire width and length variation.

Figure 13 .
Figure 13.Critical velocity values in a 3-lane TBM tunnel profile for different fire widths and lengths as well as tunnel slopes, but constant fire intensity.Values were obtained from the validated CFD model and from the semi-empirical equations (29), (32) (35) and(36).

Figure 14 .
Figure 14.Schematic of the tunnel profile used for the fire intensity and fire length variation.

Figure 15 .
Figure 15.Critical velocity values in a 3-lane TBM tunnel profile with a false ceiling, for different fire intensities and fire lengths but constant fire width, obtained from the validated CFD model and the semi-empirical equations (29), (32) (35) and (36).

Figure 16 .
Figure 16.Critical velocity determined from the Memorial Tunnel tests (triangles, circles and squares) compared to the mixed convection model (orange dots) and the equation proposed by Kennedy and two 'versions' of NFPA 502 2014 [44, 48] (dotted and dashed lines).

Figure 17 .Figure 18 .
Figure 17.Critical velocity values measured in small scale tunnels as published by Li et al. [18] compared to the proposed physical model and the trend line proposed by Li et al. [18]: (a) Test tunnel A; (b) Test tunnel B.

Figure 19 .
Figure 19.Critical velocity values for different TBM tunnel profiles and total fire heat release rates based on proposed design assumptions.

Figure 20 .
Figure 20.Critical velocity values for different rectangular tunnel profiles and total fire heat release rates based on proposed design assumptions.
Exponent in exponential function of eq.(38) for pool fires.

Table 1 .
Parameters for non-dimensional trend lines from different scale model tests.

Table 2 .
Empirical factors and constants used in the derived equations for calculating critical velocity.Values are derived from validations against pool fires and results from a validated CFD model (see Section 4 and 5).

Table 3 .
CFD methodology and boundary conditions for validation cases.

Table 5 .
Grade factor for different tunnel slopes based on a 50 MW fire.

Table 6 .
Case parameter and simulation results for different fire pan widths.

Table 8 .
Case parameter and simulation results for different fire widths and lengths.

Table 9 .
Case parameters and simulation results for different fire intensities and fire lengths.

Table 10 .
Geometric parameters for typical TBM or arched tunnel profiles used for Figure19. below.