When is the fire spreading and when it travels? – Numerical simulations of compartments with wood crib fire loads

The requirements for the occurrence of travelling fires and their distinction from previously recognized stages of fire development are not well understood. Here we present a CFD method that couples a fuel-area correction to a single-step wood pyrolysis model. The model was validated using SP medium-scale tunnel, demonstrating 19% under-prediction in peak HRR, and BST/FRS large-scale compartment where the predicted temperatures were within ±6% of the measurements. Detailed simulations of compartment geometries identified three modes of fire propagation: (1) With high ceiling heights, the ignition was followed by a fuel-controlled spread/travelling at speeds <1 cm/s (2) With lower ceiling heights, we observed a rapid (1–9 cm/s) spread towards the opening, driven by the smoke layer radiation and leading to a ventilation -controlled fire. (3) Finally, fire travelled back to the compartment interior at speed 0.1–0.9 cm/s, driven by large flames and controlled by the fuel-burnout. The analytical travelling fire model is designed to describe the first mode, but could also be used for modelling the back-travelling stage. Comparison between the CFD simulation and the analytical model indicates that further development of the analytical model is needed to account for the compartment’s heating history and tunnel -like geometries.


Introduction
The issue of fire spread modelling and its implications on the response of structures is an active research topic that has recently gained immense interest. Especially after multiple fire accidents (the Interstate bank fire in Los Angeles, One Meridian Plaza fire in Philadelphia, World trade centre fires, Windsor tower in Madrid and Faculty of Architecture at TU Delft) where the fires had been observed to move along the floor and also vertically across different floors, the prescriptive approach of the Eurocodes [1] was found to be insufficient to ensure fire safety in large complex structures. These fire scenarios are being classified as travelling fires based on the travelling fire model that was proposed by Stern-Gottfried and Rein in 2012 [2]. Such fires are likely to occur in large compartments (floor area above 100 m 2 ) and create a heterogeneous temperature field which can be more detrimental to the structure. This non-homogeneous nature of fire has also been observed during several experiments in the past ( [3][4][5][6][7]) and has also been demonstrated recently in the Malveira fire tests [8].
Fire spread is a process where a previously unburned fuel surface is heated up by hot combustion products and consequently ignited. As the fire spreads to the new fuel surface, it may increase in terms of burning surface area or power. In case of local fuel burnout, ignition and extinction fronts are formed, and a semi-steady state propagation of a fire line of constant width can be observed. Other factors, such as ventilation, may limit the power even when the area increases.
The first attempt at a detailed analysis of fire spread in compartments was made by Clifton in 1996 [9]. In a report for the Heavy Engineering Research Association (HERA) of New Zealand, Clifton first proposed that in large compartments (firecells) without any partitions, fires could move across the floor. He suggested that this could occur in compartments with floor areas above 100 m 2 and the typical height of a firecell is between 3 and 4 m. Each firecell was divided into four zones -fire, pre-heat, smoke logged and burned-out. Clifton mentions the lack of experimental data as a reason for the numerous assumptions needed in his model. In 2012, the travelling fire methodology (TFM) was proposed to describe the dynamics of fire spread in large open compartments and to support performance based design of large structures [2,10]. TFM assumes a uniformly distributed fire load and the fire spread to be fuel controlled with constant rate of propagation. The assumed fire size can vary from 1% to 99% of the floor area. It considers that a fire in a compartment has two regions: near field, flame dominated high temperature region and a far field, cooler smoke dominated regions. The Alpert's ceiling jet correlation is used to calculate the far-field temperature. Heat release rate (HRR) of the fire is calculated from the heat release rate per unit area and local burning time is obtained from the HRR and fire load density, indicating a semi-steady state fire spread. Some of the uncertainties associated with the original model were later improved by excluding unrealistic fire sizes and behaviour by accounting for the localised fire dynamics and is reported as the improved travelling fire methodology (iTFM) [11]. A limitation was proposed for the fire spread rates, based on data from a few large-scale fire experiments. Flame fluctuations was also taken into account, lowering the near field temperatures to 800 • C-1200 • C range. In 2017, Dai et al. compiled a comprehensive review describing the existing travelling fire experiments and models [12]. Additionally, they also propose another model to analyse travelling fire behaviour based on the Hasemi's localised fire model combined with a FIRM zone model calculation for the far field fuel [12,13]. This model is capable of producing spatially and temporally varying thermal fields. As it is combined with a zone model, it also responds to ventilation controlled burning. This model was also coupled with Opensees for thermo-mechanical analyses of a structure.
These models are significant steps towards improving engineering calculations of fire spread, but from the research perspective they make strong simplifying assumptions. All of them assume a constant rate of fire spread and continuous fire load distribution. The assumption of a constant spread rate is not always applicable to real fire scenarios, and gaps between the fire loads may have significant implications on fire spreading. This observation regarding the constant spread rate was also made in a recent review of the Tisova fire test by Rush et al. [14]. Additionally, other unaccounted factors that can influence the fire spread behaviour and subsequently the structural response are the geometry of the structure, radiative heating of far field fuel, pre-heating times and the influence of ventilation. In the initial formulation of travelling fire methodology [2], it was argued based on a work by Majdalani and Torero that fuel controlled scenarios would be more likely and would have a greater thermal impact on the structures than ventilation-controlled scenarios. But the same authors [15], in 2014, suggest that a mixture of fuel and ventilation-controlled burning is the more likely scenario in large compartments pointing to the experiments conducted during the Edinburgh tall building fire tests (ETFT) at the Building Research Establishment (BRE), UK. Harmathy, in his extensive review of compartment fires, also emphasizes that the ventilation-contolled fires are more severe than fuel controlled fires because they exhibit longer fully developed fire stage [16]. Harmathy reports that the highest effective heat flux was observed at the edge of the ventilation and fuel controlled regimes suggesting that to reduce fire losses, it is essential to provide large ventilation openings such that the fire is almost always fuel controlled. Maluk et al. report that in the ETFT burns for different types of fires, more than 50% energy is lost from the openings [17]. These findings imply that ventilation-controlled fires will likely have significantly longer burn times (low burning rates) and can possibly impose higher thermal loads on the structures than fuel-controlled fires, which have shorter burn durations (higher burning rates) for the same fire load density. The ventilation controlled fires also tend to produce the back and forth (travelling) aspect of the fire leading to varying heating and cooling cycles which are known to be more detrimental to structures. Dai et al. [13] also provide an excellent review and discussion about the influence of ventilation on the travelling/spreading behaviour of fire.
The influence of compartment geometry on the occurrence of travelling fires is not very well known. Drysdale describes three sources for radiative heating along with the vertical flames from the burning object/ s: hot upper surfaces, flames under the ceiling and hot gas layer. These sources can also be expected to influence the occurrence of travelling fires [18]. The relative importance of each is dependent on the compartment geometry. Most experiments have been performed in long, narrow tunnel-like geometries where the flow of hot gases is aided in one direction, enhancing the fire spread. The ceiling heights have usually been about 3 m, with the exception of Tisova test, where the ceiling height was 4.4 m. Rush et al. report that the expected fire spread behaviour was not obtained in the Tisova test [14].
Understanding the influence of the above-mentioned factors on the behaviour of travelling fires is crucial for proper performance based design of buildings. Computational Fluid Dynamics (CFD) -based fire modelling allows us to explore the dependence of various parameters more closely albeit at the expense of high computational cost. A pioneering study on fire spread modelling was performed by Koutlas [19] where he opted to model wood cribs as a thick layer of wood instead of using a stick-by-stick approach. Another initial attempt at modelling a travelling fire was by Sandström et al., in 2011 using Fire Dynamics Simulator (FDS) v5.5 [20] where the fire was modelled using multiple 2-dimensional fire sources that were turned on/off to mimic a travelling fire behaviour. In 2013, Horova et al. simulated the Veseli fire test using the ignition temperature -based simple pyrolysis model in FDS 5 [5]. The wood cribs were modelled as 5 cm thick sticks and an ignition temperature of 260 • C was used. The peak gas temperatures obtained matched the experiments but temperatures were under-predicted during the other stages of the fire. Anderson et al. used FDS to perform a-priori and a-posteriori simulations of a travelling fire test in an elongated structure (18 m × 6 m × 3 m). Six diesel pools with prescribed HRRPUA were used to model the travelling fire scenario. A grid dependence study was also conducted and it was reported that the 10 cm grids yielded satisfactory results [21]. Recently, Dai [22]. Wood boxes instead of sticks were used as fuel and a simple pyrolysis model with an ignition temperature of 200 • C was applied to the wooden surfaces. Compartment dimensions in the Cardington test were 22.8 m × 5.6 m × 2.75 m. The qualitative behaviour and temporal evolution of the predicted fire agreed with the experimental results. The ceiling temperatures matched experimental values when the fire moved towards the opening but were consistently under-predicted when it moved from the opening to the rear of the compartment. In a previous publication [23], the authors of the current study also used a calibrated ignition temperature based pyrolysis model to simulate the SP tunnel fire spread experiments conducted by Hansen and Ingason [24]. The model under-predicted the peak experimental HRRs by 3% on average. Later, the model was used to perform predictive simulations of uniformly distributed fire load in a large compartment and this was used to address the issue of grid dependence in fire modelling, suggesting the use of a surface area correction parameter to improve the predictions in coarse grids, thus allowing computational times that were feasible for engineering applications. All these studies mostly tried to calibrate the HRR of the fuel to reproduce the fire scenario.
Understanding the behaviour of travelling fires, or fire spread in general, would need models that take into account some of the complex fire dynamics involved in spreading fires. Such a model needs to account for the difficulty/ease of fuel ignition and burning by connecting it to the thermal environment which requires the coupling of burning rate to the thermal environment and the ventilation conditions. In a tunnel like geometry, it may be possible to prescribe the fire's ventilation response (front to rear movement) even with an ignition temperature based pyrolysis model, as there is only a single path the fire can propagate, but the finer details, such as the change in the local burning rate as the fire becomes under-ventilated, would not be captured. To capture the different stages of the fire development, a more complex pyrolysis model is required. These aspects become more important when the geometry of the domain is complex and fire has multiple propagation paths.
The objective of this study is to provide a reliable modelling approach for the CFD -based predictive fire spread simulations and also to explore the geometrical and ventilation effects on fire spread. The single-step, reaction kinetics based pyrolysis model used here is a significant improvement over the ignition temperature-based pyrolysis model. It helps us to solve the issue of ventilation dependent burning behaviour by coupling the mass loss rate to the temperature and ventilation conditions. The numerical model is first validated using the SP tunnel fire experiments [24] and the BST/FRS large-scale fire tests [25]. We then use the method to investigate the effects of compartment ceiling height and ventilation arrangement. Results from the simulations are also compared against the analytical travelling fire model [11].

Numerical method
Fire Dynamics Simulator (FDS) is a Large Eddy Simulation (LES) -based Computational Fluid Dynamics code that solves low Mach number combustion equations on a rectilinear grid over time. The simulations in this study were performed using FDS versions 6.7.1 and 6.7.3. Combustion was modelled as a single-step reaction of C 3.4 H 6.2 O 2.5 using the Eddy dissipation concept method and a soot yield of 0.01. Turbulence was modelled using the default Very Large Eddy Simulation (VLES) mode where eddy viscosity was modelled using Deardoff turbulence model. Turbulence near the walls is modelled using the Smagorinsky model with Van Driest damping to model the near-wall eddy viscosity [26]. The default radiation model was used in the simulations which computes the radiation intensity over 104 solid angles using gray absorption coefficient and a radiative fraction of the local heat release rate equal to 0.35.
The external boundaries in the simulations are specific to each case and the details are specified in their respective sections. In the simulations, the domain was extended outside the compartment to minimise the influence of the constant pressure assumption used at the open boundary and also to account for the energy produced by the flames that burn outside the compartment. Mesh resolution was constant within the domain and cell size doubled in the extended domain.
Within the solid, heat transfer is calculated using the onedimensional heat conduction equation where ρ s , k s , c s are the density, conductivity, specific heat of the solid and q ′ ′ ′ s is the heat source term due to the heats of reaction of the pyrolysis and evaporation reactions. The boundary conditions are the convective and radiative heat flux at the solid surfaces [27]. Two pyrolysis models were used in this work: ignition temperature based pyrolysis (ITP) and single-step pyrolysis (SSP) models. The drying and charring processes were modelled using the SSP model with two parallel reactions where y f is the fuel yield of the wood pyrolysis reaction. The reaction rates are modelled using the Arrhenius equations where α is the material component i.e. fuel or moisture, ρ s (0) is the initial density of the material, X O2 is the oxygen concentration, A is the pre-exponential factor, E is the activation energy and n is the reaction order. The chemical source term depends on the reaction rates as where ΔH r,α is the latent heat of vaporization for equation (2) and the heat of reaction for equation (3).
The kinetic parameters and the product yields were taken from the work by Rinta-Paavola and Hostikka [31] where they estimate and optimize the parameters using GPyro. The thermal properties were calibrated manually by reproducing the cone calorimeter experiment at 50 kW/m 2 radiation level. The experimental cone calorimeter data was provided to the authors through a private communication [29]. The kinetic and thermal parameters are summarized in Table 1. Fig. 1a shows a comparison of the measured and predicted mass in TGA at four different heating rates. Fig. 1b shows a comparison of the measured and predicted heat release rates in the cone calorimeter experiment.
Charred wood cribs are assumed to lose their shape and collapse when they lose their mechanical strength. To simulate this, the wood obstructions were removed from the simulation when their mass loss reaches the available user-defined bulk density (BD). In FDS, the feature is activated using 'BURN AWAY' and 'BULK DENSITY' parameters. The bulk density depends on the required/reported fire load density (FLD expressed in terms of combustible energy per floor area) and is calculated as: (6) where A floor is the floor area of the compartment, V solid is the volume of solid fuel, y char is the char yield, y m is the wood moisture content, ΔH c is the heat of combustion (MJ/kg), φ is the appropriate surface area correction factor for the crib, LD is the contribution from Layer Divide factor and n cribs is the number of cribs in the compartment. For calculating the bulk density of a single crib the value of A floor becomes the area covered by the crib. If fire load density is not available then the bulk density can also be calculated based on the mass of individual cribs. The purpose of the LD factor, defined using 'LAYER DIVIDE' input parameter, is to ensure the conservation of mass in a uniform heating scenario by controlling the thickness of the solid from which fuel vapours are assigned to each surface. This method is used because the pyrolysis solver in its current form does not account for the internal pressure distribution and gas transfer. The default value of LD is 0.5, i.e. the division is symmetrical as both sides are equally likely to produce pyrolysis products. In the present simulations, except the SP tunnel fire tests, the fuel objects are heated mainly from above. Therefore, the fuel generated from the entire thickness is expected to be released from the same side, i.e. LD = 1. The mass conservation is ensured by the bulk density parameter and burn-away feature, regardless of the LD value.
The area adjust parameter (φ) is the same as the one described in Ref. [23]. It is a surface area correction factor used to compensate for the decrease in total surface area of a wood crib when it is modelled using a coarse mesh.
where A total, exp and A total, model are the experimental and modelled surface areas, respectively. In addition to the SSP model, the ITP model was used in the simulations for the initial fire to reduce the sensitivity to the ignition source properties. ITP model solves the heat conduction equation (1) without the chemical source term, and starts to release fuel gas at a prescribed rate once the surface temperature reaches the user-defined ignition temperature [23].

Validation: SP tunnel fire tests
The first set of experiments used for the validation study were performed by Hansen and Ingason [24] at SP (now known as RISE) in Sweden to understand the influence of ventilation on HRR of wooden pallets placed at varying distances. The parameters tested were the distance between piles of wood pallets and longitudinal ventilation velocities (u c ) within a 10 m × 0.6 m × 0.4 m domain. The models and experimental setups are identical to those used in Ref. [23] except for the pyrolysis model. In this study, we simulated only tests 3 and 6 where u c = 0.6 m/s and the distances between the four piles were as follows: Test 3: 0.4 m-0.7 m -0.6 m; Test 6: 0.5 m-0.8 m -0.9 m.
The gas phase cell sizes were 2.5 cm × 2.5 cm × 1.8 cm in the sections of the tunnel containing wood piles, and 4 cm × 4 cm × 3.6 cm elsewhere. The chosen mesh sizes were the best combination to model the geometry of the wood pallets as per the experiments and also to reduce the computation time. The surface area correction factor was not used in these simulations and 'LAYER DIVIDE' was 0.5. Inlet boundary condition was modelled as a supply vent with enforced constant velocity of 0.6 m/s. The outlet was modelled as a HVAC duct connected to ambient as the performance parameters of the exhaust system used in the experiments was not available. The walls, floor and ceiling of the tunnel were modelled according to the material properties listed in Ref. [24].
Three different combinations of ITP and SSP models were used to simulate the fire spread, as shown in Table 2. The wood piles modelled with the ITP model had a HRRPUA of 260 kW/m 2 and a HRR ramp to control the burning and SSP model was used to specify the burning behaviour of the wood pallet piles. The bulk density parameter was used to specify the combustible mass of each pile such that the simulated pile mass is equal to the experimental mass. A vent with a specified HRRPUA of 750 kW/m 2 was used as the ignition source. Further details of this model are provided in Ref. [23]. The mass loss rate of the cribs in the simulation was obtained by recording the integral of the burning rate over the entire crib surface area.

Validation: BST/FRS large-scale compartment fire tests
The second set of experiments used for validation was the Test 2 from the British Steel technical and Fire Research Station (BST/FRS) largescale fire tests conducted at the BRE experimental facility in Cardington in 1999 by Kirby et al. to assess the applicability of parametric fire curves proposed in the Eurocodes [25]. This was one of the earliest large-scale fire spread (or travelling fire) experiments conducted.
Eleven rows of wood cribs were arranged discretely in a 22.85 m × Fig. 1. Comparison of the experimental and simulated TGA results [31], and experimental and simulated heat release rate per unit area in a cone calorimeter under a radiation of 50 kW/m 2 for the pine sample [29].
) was reported as 1.4795. Temperature measurements were made 300 mm below the ceiling along crib lines 2, 6 and 10 using K-type thermocouples. The simulation model of the experimental domain and the arrangement of cribs is shown in Fig. 2. The compartment domain was connected to an external, 2.0 m × 5.6 m × 5.6 m domain with open end and top boundaries. The compartment lining properties were set as mentioned in Ref. [25]. The wood cribs were spaced 1.0 m apart along the compartment length and 0.9 m apart along the width. They were modelled stick-by-stick, and arranged in four alternately perpendicular layers with five sticks in each layer. The simulations were carried out with a 10 cm resolution as this was the computationally feasible option. The sticks were thus represented by 10 cm thick obstructions, for which the area correction factors were calculated using the area-accurate, 5.0 cm resolution models as a reference. The details of the correction calculation are shown in Table 3. In the heat conduction solver connected to the stick surfaces, the stick thickness was 5.0 cm.
A combination of ITP and SSP models was used for the wood cribs. The first line of cribs (shown using red colour in Fig. 2) was modelled using the ignition temperature of 300 • C and a specified HRRPUA of 320 kW/m 2 . The HRRPUA was obtained through calibration to match temperature profile along crib line 2 for the first 12 min. The ITP cribs were ignited by a 0.3 m 2 vent on the floor, under these cribs, with HRRPUA increasing to 500 kW/m 2 in 60 s, staying constant to 220 s, and stopping burning at 300 s. The remaining cribs were modelled using the SSP model. Based on the work of G.M.E Cooke, as reported by Dai et al. [22], the heat of combustion was set to 19500 kJ/kg. A bulk density of 277 kg/m 3 was used to enforce the total mass of 78 kg per crib. This value accounts for the char produced, the calculated surface area correction factor and the Layer Divide = 1.

Simplified tunnel simulations
To enable more detailed investigation of the details related to the longitudinal fire spread, we defined an artificial tunnel scenario. The domain is 30 m long, 1.0 m wide, and open to the ambient at one end. The ceiling height was either 3.0 m or 6.0 m, to explore its effect on the fire spread behaviour. The compartment boundaries were modelled as generic concrete walls with ρ = 2200 kg/m 3 , c p = 0.75 kJ/kgK and k s = 1.7 W/mK [18].
Wood sticks were arranged in the tunnel as uniformly distributed cribs with a fire load density of 600 MJ/m 2 . The heat of combustion of wood volatiles was 16000 kJ/kg [18]. The cribs were modelled stick-by-stick and arranged in mutually perpendicular layers. Fire spread was initiated using a 2.0 m 2 burner with maximum HRRPUA = 110 kW/m 2 , growing at ultra fast rate, and reacing the maximum value at 180 s. Beyond the prescribed region, fire spread was predicted using SSP model.
The simulations were performed using FDS 6.7.3 at mesh resolutions of 5.0 cm, 10 cm and 15 cm. The fuel surface area correction factors were calculated using the 5.0 cm resolution as a reference. The calculation details are shown in Table 4.

Square compartment with (a)symmetric ventilation
To investigate the role of the ventilation openings in a realistic, nontunnel -like compartment, a 3D model of a square compartment was built with symmetric and asymmetric natural ventilation. The compartment dimensions are 30 m × 30 m × 3.6 m, similar to the model in Ref. [23]. The compartment boundaries were modelled as 30 cm thick concrete layers with ρ = 2200 kg/m 3 , c p = 0.75 kJ/kg.K and k s = 1.7 W/(m.K) [18].
The fire load was made of a uniform layer of 0.8 m high wood cribs. With the 20 cm resolution of the simulations, four layers of sticks were arranged in mutually perpendicular arrangement. A combination of prescribed and predicted fires was used. The prescribed fire consisted of a 16 m 2 region of wood crib at the centre of the compartment, with ignition from the middle and a specified spread rate of 0.00693 m/s outwards. As a result, the prescribed HRR developed at ultra-fast rate, and the fire front reached the edge of the prescribed crib at 290 s. Around the prescribed crib, the SSP model was used. Similar to the BST/ FRS simulations, the FDS default Layer Divide value (0.5) did not produce sustained fire spread, and hence its value was increased to 1.0. The Burn Away parameter was not used in this case to demonstrate the influence of char oxidation on the temperatures and heat fluxes. This  Table 3 Calculation of the surface area correction factor for a 10 cm -resolution model of a 1 m 2 crib in the BST/FRS compartment scenario.  Table 4 Calculation of surface area correction factors for the 10 cm and 15 cm resolution models of a 1 m 2 crib in the simplified tunnel model. means that the fire load density depends on the specified wood density which is 393 kg/m 3 in this case. The surface area correction factor used was 1.55. The fire is naturally ventilated through two or four vertical, 20 m × 2.2 m openings, as shown in Fig. 3. The amount of the ventilation openings can be evaluated using the inverse ventilation factor where A T is the total surface area excluding the floor area, A w is the vent area, and H is the vent height. According to Thomas and Heselden [18], fires with ψ < 10m 1/2 can be classified fuel controlled. In the current geometry, the asymmetric two vent -geometry has ψ = 9.6 m 1/2 , and the four vents give ψ = 4.5 m 1/2 , indicating that the symmetric, four vent -scenario should exhibit a more fuel-controlled behaviour. Dai et al. discuss the inverse ventilation factor and ventilation openings in detail in Ref. [13].

Validation: SP tunnel fire test
The SP tunnel simulations were performed on the CSC-Taito supercomputer using FDS 6.7.1. The experimental tunnel was modelled as per dimensions and the domain was discretized into 8 meshes, which were assigned to separate MPI processes for the parallel computation. The average computational time for each simulation was around 30 h. Fig. 4 shows the measured and predicted heat release rates (HRR) for the different pyrolysis model combinations in tests 3 and 6. The heat release rates (HRR) are under-predicted by 38.2% for Test 3 and 52% for Test 6 in the first set of simulations where only the first wood pile was modelled using the ITP model. In the second set, the first two piles modelled using the ITP model and they provide sufficient energy to support burning of the other piles. In the third set, the first three wood piles were modelled using an ITP model and we move closer to the modelling method used in Ref. [23]. Using three ITP piles instead of two slightly increases the growth rate but does not have a noticeable effect on the peak HRR. Overall, the under-prediction of HRR falls within the 20% range of uncertainty in the set 2 and set 3. HRR in the decay phase is under-predicted in the second and third sets whereas simulation predictions of the first set match the experimental HRR. In set 1, the simulated Mass Loss Rate (MLR) in Fig. 5 curves indicate that the piles in the decay phase are still burning. MLR curves of the second and third set show that the combustibles are almost completely consumed by 600 s. Even though the simulations specify the same mass as reported in the experiments, the HRR behaviour implies that oxidation of char plays an important role during the extinction phase.
The simulated Mass Loss Rates (MLR) of each pile in Test 6 are shown in Fig. 5. In the ITP1-SSP3 simulation, piles 2-4 ignite about 15 s later than in the other simulations. The low burning rate of pile 2 leads to reduced heat transfer to piles 3 and 4 resulting in a lower HRR. These piles are still burning during the decay stage of the overall HRR and this contributes to the good agreement in HRR observed during the decay stage in the ITP1-SSP3 model. MLR curves of the ITP2-SSP2 and ITP3-SSP1 simulations show that the piles are almost completely consumed by 600 s and subsequently show a clear under-prediction of overall HRR. Even though the simulations reproduce the experimental mass losses, the HRR of the decay stage seems to include a significant contribution from the char oxidation -a phenomenon that we did not include in our model.
Overall, the different inter-pile distances in tests 3 and 6 affect the HRR in only the first model combination (ITP1-SSP3). The small changes in pile heat transfer propagate to the kinetics-based pyrolysis rates of the SSP model, but are not affecting the ignition times and subsequent prescribed HRR of the ITP models.
These validation simulations showed that starting the fire spread process using the SSP model is difficult and, from the modelling perspective, it is better to combine the SSP model with a more robust and easily ignitable pyrolysis model. This could either be a prescribed fire for a certain period of time or an ITP model along with an ignition source.

Qualitative analysis and validation
The simulations of the BST/FRS large-scale compartment fire test were performed using FDS 6.7.1 on a personal workstation. The simulations were performed with mesh resolutions of 10 cm and 20 cm to investigate the sensitivity of the model against mesh size. The domain was divided into 12 meshes for parallel computing, assigning each to a separate MPI process. The simulation uncertainty is evaluated based on the measured averaged gas temperatures above the crib-lines 2, 6 and 10.
The results from the 10 cm mesh resolution are presented here for comparison against the experimental results. To build a qualitative understanding of the fire development, we first investigate the predicted temperature, velocity and gas concentration fields. Figs. 6 and 7 show the vector fields of temperature and horizontal u-velocity component (uvelocity) in the middle of the compartment at different times. Figs. 8 and 9 illustrate the oxygen and fuel concentrations.
12 min The flames form a burning front along the first row of cribs, and a classical ceiling jet flow is established, with hot gases (400-700 • C) flowing towards the opening, and ambient air entraining into the plume. Only the first row of cribs is producing fuel, as shown in Fig. 9. The maximum u-velocity of the hot gases is over 3 m/s and u-velocity of the entrained air is about 1.3 m/s. The spacings between the cribs produce small disturbances in the entrained flow as seen in Fig. 7. Fig. 8 shows that oxygen is mainly consumed along the first row, but a small reduction in the oxygen concentrations can be also observed along rows 2-6 due to moisture evaporation from the cribs.
13.5 min Flaming starts on the second row of cribs. The tunnel -like geometry (low ceiling height) promotes drying of the cribs which then allows for the rapid fire propagation. The velocity data also shows how the flow turbulence is generated due to the interaction between the plumes of from the two rows.
15 min Burning is more prominent along the second row and the hot gas layer thickness has also increased. Temperature along the third row begins to increase and the oxygen levels along the first two rows start to decrease.
15.5 min At 15.5 min, the burning front starts to propagate rapidly towards the opening, occupying the third and fourth rows. Small flames are also seen along the fifth and sixth rows.   but the cribs still continue to release fuel vapours, as seen in Fig. 9.

min
The flame front is in-between rows 7 and 8, but flaming is observed in the remaining rows as well. The flow field is symmetric about the flame front, with hot gases flowing away from the plume, and colder air at the bottom flowing into the plume.
17.5 min By this time, the burning front has reached the opening, and burning mainly occurs along rows 10 and 11. The propagation speed of the fire during this rapid spread stage was 9.2 cm/s. Most of the hot gases flow out of the compartment, and the air is entrained into the plume at over −2.5 m/s. The gas temperature in the closed end varies from 350 • C near the cribs to 800 • C near the ceiling. 20.5 min Fire has started its propagation back towards the closed end. The cribs in row 11 have lost about 60% of their specified combustible mass. Fig. 8 shows that the oxygen levels at this time is almost zero in the entire domain, but a small pocket of oxygen is seen at row 10. Fig. 7 shows that the flow field in the closed end exhibits three velocity layers instead of the two seen before: At the ceiling, the gases flow from the burning front towards the closed end, and the velocity ranges from −1.2 m/s to −2.5 m/s. Above the cribs, the gases circulate back towards the burning front with velocities between 1.5 m/s and 2.0 m/s. In-between the cribs, the gases flow away from the burning front with velocities of −0.04 m/s to −0.15 m/s. It seems that some of the entrained air is flowing towards the closed end.
28.5 min The flame front has reached the eighth row. A very broad temperature range 20-1200 • C is observed towards the opening, and the flow behaviour is the same as few min before. The temperature around the cribs is above the level required for pyrolysis, and it can be seen from Fig. 9 that the cribs produce fuel vapours. Oxygen levels in the closed end cannot support combustion at this point. The three layered flow behaviour is seen in the closed end until 54 min, after which the usual two-layer flow occurs as the burning stops. The propagation speed during this steady backward burning stage was 0.78 cm/s. Fig. 10 shows a comparison of the experimental and simulated temperatures along crib lines 2, 6 and 10. Three distinct stages of   temperature development can be observed: First, a growth stage between 0 and 18 min, where the temperatures increase at all measurement points. During this stage, according to Figs. 6-9, the fire spreads from the ignition source to the opening. A semi-steady stage between 18 and roughly 60 min follows, where temperatures of the crib lines 2 and 6 first drop, as the fire front passes them on the way towards the opening, and then gradually increase towards the second peak. At crib line 2, the second peak is observed in the end of the semi-steady stage, but crib line 6 has a peak at about 40 min. At crib line 10, the temperature does not drop at 18 min, but reaches the peak right away and starts to gradually decrease. This second stage corresponds to the steady burning of the last cribs and the subsequent propagation back to the closed end, as the combustible mass of the cribs is consumed. This stage best corresponds to the travelling fire behaviour. The third stage of the time-temperature curve is the decay from 1000 • C range to about 300 • C at the end of the experiment.
Generally, the simulations reproduce the observed stages both in temperature and in duration with both the mesh sizes. In the 10 cm mesh simulations, the accuracy of predictions at crib line 2 during the first minutes is, however, a consequence of the ITP model calibration. After the first 18 min, the first peak of temperature is seen and reproduced precisely, but during the travelling stage, these closed-end temperatures are under-predicted more than 20%. Also, the second temperature peak is delayed by 6.0 min and under-predicted by 9%.
At crib line 6, the first temperature peak is over-predicted by 5%. The second temperature peak is predicted precisely both temporally and in magnitude. The drop of the predicted temperatures across crib-lines 2 and 6 around 20 min is more prominent than in the experiment. At this instance, the temperatures near the crib top surfaces are around 350 • C (Fig. 6) which is above the level required for pyrolyzing the wood. The fuel vapours released from the cribs do not burn due to lack of oxygen but are mixed with the hot gases, reducing the gas temperatures. The strongest fuel production happens at crib line 1, as this was modelled using ITP model, and the crib is consumed before the fire travels back to the end of the compartment. At crib line 10, the correspondence is excellent during the first 60 min, i.e. from ignition to the time at which the backwards-travelling fire front passes the point of measurement. The temperatures during the decay stage are under-predicted, mainly due to the absence of glowing char in the simulations. Another factor would be the premature burn-away of the first row of cribs. As these cribs were modelled using ITP model, their MLR is decoupled from the temperature after ignition, and they continue to produce fuel vapours once ignited, causing them to be consumed too soon. In the experiments, cribs in the first row were the last remaining fire load in the compartment; whereas in the simulation, cribs in the second row are the last ones to burn. This is an important shortcoming of the modelling choice of combining the ITP and the SSP model.
In the 20 cm mesh simulation, similar peak temperature values are reproduced along the three locations as in the 10 cm mesh. The mesh size mainly affects the temporal development. At crib line 2, the first peak is delayed by about 3.5 min and the second peak by 9.5 min with respect to experiments. The temperatures during the travelling stage are under-predicted by 27%-38% in comparison to the experiments. At crib line 6, the first temperature peak is delayed by 5 min and the second peak by 15 min. Temperature during the travelling stage is underpredicted more than 30%, and the second temperature peak is overpredicted by 2%. At crib line 10, the predicted temperature value is within the experimental uncertainty, but the entire temperature development is delayed by about 5 min. The temperatures during the decay stage are under-predicted due to the absence of glowing char as in the 10 cm mesh simulations. Overall, a temporal delay and a significant under-prediction of closed-end temperatures during the travelling stage is observed in the 20 cm mesh simulations. One reason for the temporal delay is the reduced porosity and heat transfer inside the crib. Nonetheless, the 20 cm mesh simulation still manages to capture the main qualitative and quantitative behaviour of the fire.

Comparison to the analytical iTFM model
To investigate how well the analytical travelling fire model iTFM [11] can describe the observed behaviour, we next compare the simulated temperatures against temperature estimates from the iTFM model. The iTFM model allows the user to vary the compartment dimensions, spatial and temporal resolution of the model, fire length, heat release rate per unit area, fire load density, room temperature, peak fire temperature and the flapping angle of the fire. The spatial resolution and time step of the iTFM calculation was 0.1 m and 0.1 s, respectively. Fig. 11 presents a comparison between the temperatures from the FDS and iTFM models at two different locations along the BST/FRS compartment. The comparison is limited to the period when the fire travels back towards the closed end of the compartment, and the locations are indicated as a distance from the opening. Multiple fire sizes (3.25%, 4.3% and 6.3% of floor area) were used in the iTFM model. The fire size of 4.3% corresponds to the crib size in the experiments, 3.25% corresponds to a fraction of the first crib size burning at the moment when the second crib ignites, and 6.3% corresponds to 1.5 times the crib size. The compartment dimensions were 23.0 m × 5.6 m × 2.8 m and the reported experimental fire load density of 380 MJ/m 2 was used. The HRRPUA, 1800 kW/m 2 , was calibrated to produce the best agreement with the simulation data.
The results show that, with an appropriate value of the HRRPUA, the iTFM model produced temperature curves with similar peak values as the simulations, but the temperatures ahead the burning front were under-predicted. This is obviously caused by the fact that the simulated fire had already developed for 20 min and pre-heated the region. The temperature behaviours during the decay stage were similar to the simulated ones.

Effect of increased ceiling height
In the qualitative analysis of the BST/FRS large-scale compartment simulations, the fire development in the early spread stage was found to be driven by the hot gas flow below the ceiling, and the flow was characterized as 'tunnel-like'. To investigate the effect of the ceiling height on the fire development, two additional simulations were performed where the ceiling height was increased from 2.8 m to 3.6 m and 5.6 m, corresponding to 29% and 100% increase, respectively. All the other boundary conditions were kept unchanged.
The influence of the increased ceiling height on the average gas temperature above crib-line 2 and the incident heat flux on the center crib of crib-line 2 is shown in Fig. 12. With the ceiling height of 2.8 m, ignition occurred when the heat flux exceeded 20 kW/m 2 (Waterman criterion) [18]. When the ceiling height is increased to 3.6 m or 5.6 m, the heat flux decreases by 94% and 96%. As a consequence, the fire does not spread, and the temperature histories reflect the ignition, burning and burnout of the first line of cribs, modelled using ITP. The temperature values are reduced by 56% and 72%, respectively.
Figs. 13 and 14 show the comparison of temperature and oxygen concentrations along the middle of the compartment for a ceiling height of 2.8 m and 3.6 m at 860 s. Increasing the domain height seems to decrease temperatures along the downwind side which results in a lack of fire spread. Hinkley and Wraight [30] suggested that when the ceiling height is low, flaming occurs away from the ceiling where the oxygen concentration is higher. The opposite is observed when ceiling height is increased: The flaming occurs closer the ceiling, as more oxygen is available. Fig. 14 shows how the flaming occurs away from the ceiling where oxygen rich air is available.
Figs. 13 and 14 show the instantaneous temperature and oxygen concentration fields from the 2.8 m, 3.6 m and 5.6 m ceiling height simulations at 860 s (14.3 min) after ignition. In the 2.8 m compartment, the flame reaches the ceiling, at least occasionally, and the plume bends over the neighbouring crib. In the taller rooms, the fire plume is more prominently pushed towards the back wall, and the combustion mainly occurs close to the wall, as more ambient air is entrained into the plume. As a result, the effective radiating region of high temperatures and soot concentration appears at a smaller solid angle to the second crib than in the 2.8 m room. From the viewpoint of the fire spreading, the gap between the cribs plays a crucial role. Fire might spread also with a higher ceiling if the fire load was arranged in continuous pattern instead of the discrete crib arrangement. Fig. 11. Comparison of predicted gas temperatures from iTFM model [11] and FDS along the midline of the BST/FRS large-scale test no. 2 [25].

Effect of ceiling height and grid resolution
The simplified tunnel simulations were performed using FDS 6.7.3 on two computing clusters: Puhti cluster, maintained by CSC -IT centre for science Ltd., Finland, and the Aalto Science-IT cluster, Finland. The computational domain was decomposed into 61, 31 and 11 meshes in 5 cm, 10 cm and 15 cm grid simulations, respectively. The 5 cm grid simulations were run for 1500 s, taking 3 weeks of computational time. The average computational times for the 3600 s simulations with 10 cm and 15 cm grids were about 52 h and 32 h, respectively. The 5 cm -simulation is used as a reference for evaluating the accuracy of the 10 cm and 15 cm -simulations.
The overall fire development can be analysed using the instantaneous oxygen concentrations, shown in Fig. 15 for the 3.0 m high tunnel simulation with 15 cm mesh. The burning starts in the predicted region and forms a smoke layer with low oxygen concentrations. Rapid fire spread towards the opening follows at speed 5 cm/s. During the spread stage, the fire load is not entirely consumed. The fire reaches the opening at about 540 s, and for the following 1000 s, burning occurs predominantly at the opening. After 1500 s, fire starts to move back towards the closed end at speed 0.14 cm/s, consuming practially all the fuel on the way. From the viewpoint of fuel burnout -driven, semi-steady propagation, this stage can be called a travelling fire stage, although it appears after a fully developed stage which in smaller compartment follows flashover.
Interestingly, the oxygen concentrations on the close end -side of the fire front seem to have very different values in the different stages of the fire. During the spread stage and in the middle of the travelling stage, they are about 7 vol % or less, but at times 1600 s and 5756 s, the average concentration is clearly above 10 vol%. This is likely due to the continuous fuel bed acting as a channel entraining air into the   compartment. Spatially, the oxygen field is almost uniform, except for the shallow layer inside the cribs where the concentrations are below the gas phase -average. Fig. 16 shows the cellulose concentrations inside the cribs at 1600 s. A small amount of fuel is entrained towards the closedend through the wood cribs. As the temperatures within the fuel bed are too low to support combustion, the likely explanation for the reduced oxiygen concentrations is the mixing with the fuel vapours. Fig. 17 gives more insight to the gas dynamics by plotting the gas temperatures (15 cm below the ceiling) and oxygen volume fractions (vertically middle of the domain) at every 5.0 m along the tunnel. The temperature peak of 600 • C at 100-500 s represents the rapid fire spread towards the opening. The fire then burns at the opening until 1500 s, as seen in Fig. 15, and the highest temperature is observed near the opening at 30 m. As the fire travels back towards the closed end, the location of the peak temperatures also moves. During the spread stage, the temperature increase at each location is followed by a drop in oxygen concentration from the initial value to about 7%. From about 1500 s onwards, we can observe periodic fluctuations with cycle length of 1900 s ± 100 s. In oxygen concentrations, the fluctuations are seen in the region between the closed end and the burning front, and in ceiling jet temperature they are seen in the region between the burning front and the opening. This indicates that the oxygen level on the closed side of the tunnel is slightly replenished by the entrained air. Comparing Figs. 15 and 17 reveals the high closed side oxygen levels in Fig. 15 were due to the fluctuations. Fig. 18 shows the temperatures within the domain. At 800 s, the fire burns at the opening and creates a one-sided ceiling jet along the tunnel. Considerable amount of the heat flows into the tunnel, and not out from opening, as one would expect based on the classical picture of the ventilation-controlled compartment fire. Temperatures near the closed end of the tunnel are over 100 • C. At 6000 s, the heterogeneous nature of the temperature field is obvious, best described as a travelling fire plume.
More detailed understanding of the flow field can be obtained from the horizontal velocity components, shown in Fig. 19. During the growth stage, a hot gas layer, approximately 60 cm thick, is formed under the ceiling maximum speed of 4.0 m/s. The large-scale flow field is close to symmetric about the burning front (fire plume) at all times: Hot gases under the ceiling flow away from the plume, and colder gas entrains into the plume from below the smoke layer. But only at the vent side does the entraining air provide oxygen. The flow produces a lot of eddies leading to significant mixing and dissipation of the temperature gradients.
With a ceiling height of 6.0 m, burning starts in the predicted region, but after the prescribed period of burning, the combustion (low oxygen) and high temperatures are not observed under the ceiling, as with the 3.0 m ceiling height (Fig. 20). The fire spreads slowly towards the opening as a counterflow flame spread process. Spread rate is approximately 0.13 cm/s. Theoretically, the length of the burning region depends on the fire load, and a constant spread rate can be achieved. This stage of fire spread can be classified as a travelling fire. Fig. 21 shows a comparison of the heat release rates obtained with different grid resolutions and ceiling heights. With the 3.0 m ceiling height, the HRR curves show three distinct stages: rapid growth stage, controlled by the prescribed fire; reduction stage when the fire spreads from the prescribed region to the opening; and a semi-steady stage when the fire travels back from the opening. The growth stage lasts about 150 s and is consistent across all the resolutions. During the reduction stage, spurious oscillations are observed in the 5.0 cm and 15 cm -resolution simulations. They are likely due to small numerical inconsistencies in gas concentrations across mesh boundaries. The reduction stage is more unstable than the other two stages, but the duration of the simulation with 5.0 cm resolution is too short to conclude about the stability of fire dynamics. The three stages are very similar to the BST/FRS large-scale compartment.
The simulations with the 6.0 m ceiling height (Fig. 21 right) show the same three stages as the 3.0 m ceiling height, but they correspond to different physical behaviours: The growth stage is again dominated by the prescribed fire, but the peak HRR achieved at the end of the growth stage is higher. This indicates that the availability of oxygen through the tunnel entrainment has an effect on the peak HRR. The reduction stage represents a transition from the ignition source -dominated, thermally assisted burning to the natural spreading, and the semi-steady stage represents slow counterflow spreading. Unlike the BST/FRS large-scale compartment with increased ceiling, the fire continues spreading. This is mainly due to the continuous fuel arrangement. Table 5 shows the relative differences between the lower and higherresolution predictions of the HRR maxima and semi-steady stage averages for the ceiling heights 3.0 m and 6.0 m. In engineering applications, it would be difficult to use a 5 cm grid and therefore the reported differences can be used as the expected uncertainties for the 10 cm and 15 cm grid simulations.

Comparison to the analytical iTFM model in 3.0 high tunnel
Similar to the BST/FRS large-scale compartment, we compare the simulated temperatures (20 cm below the ceiling) from the 10 cm and 15 cm resolution simulations of the 3.0 high tunnel against the temperature estimates from the iTFM model. The parameters of the iTFM model are chosen to match the simulated burning conditions closely. As the fire size is not constant in the simulation, three different fire sizes (1.6%, 2.5% and 3.4% of floor area) are used to estimate the temperatures in the iTFM model. The HRRPUA of 800 kW/m 2 and peak fire temperature of 500 • C are estimated based on the simulations. The flapping angle used is 10 • . The spatial resolution and time step of the iTFM calculation are 0.1 m and 0.1 s, respectively.
The first important difference between the FDS and iTFM results is related to the spread rates. The rapid fire spread in the beginning of the fire, observed in the CFD simulation, is not predicted by the analytical model which assumes a fuel-controlled steady spread process. While FDS predicts spread rates of 5.0 cm/s and 0.14 cm/s in the forward and backward directions, respectively, the analytical model with 3.4% floor area produces a steady fire spread at rate 0.11 cm/s, i.e. 21 % less than the semi-steady state back-propagation rate of FDS. Comparison between the FDS and iTFM temperatures, in Fig. 22, therefore starts at the moment when the FDS fire starts moving from the opening towards the back of the tunnel, and the Location-axis indicates the distance from the opening. The first set of plots shows predicted temperatures as a function of location at different times between 10 min and 120 min. As noted, the high-temperature region propagates at similar speed, reaching the distance of 10-15 m in 180 min. The second set of plots on the right shows temperatures at specific locations as a function of time. The results from the 10 cm resolution are not shown in some plots as these instances exceed the simulated time of 3600 s.
As we saw in Fig. 18, the simulated temperature field does correspond to the phenomenological assumption behind the travelling fire models, but quantitative differences are important. Namely, the predicted far-field temperatures are approximately 200 • C higher in FDS. Increasing the iTFM flapping angle increases the analytical far-field temperatures slightly but also reduces the peak temperature. It is clear that the calculating the far-field temperatures with the Alpert's ceiling jet model, which assumes a symmetrical spreading of flow and energy in all directions, leads to greatly under-estimated temperatures in tunnellike scenarios, where the energy is accumulated in a limited domain.

Square compartment with (a)symmetric ventilation
The simulations were performed on a personal workstation with 18 processor cores and 64 GB memory. The domain was decomposed into 10 meshes and each process was assigned to a separate MPI process with two OpenMP threads each. Computational time for the 7200 s simulations was around 42-60 h.
We will first look at the fire development in the compartment with two vents. The predicted flame positions at different times are visualized in Fig. 23 that also shows oxygen volume fractions at the coloured slices. At 2.5 min after ignition, the flame is mainly in the middle of the compartment but tilts towards the closed corner due to the ventilation flow, creating a thicker smoke layer in that region with about 16% oxygen concentration. At 4 min, the fire has increased in size, the smoke layer temperature is above 700 • C (Fig. 24) and the flame moves towards the closed end due to the strong negative velocity at height 0.8 m (Fig. 25). At 6 min, the low oxygen level in the closed corner drives the growing flame towards the openings. The flame movement towards the openings is slow, and it reaches the opening at 16 min. At this point, the highest smoke layer temperatures close to the vents are 800 • C, negative velocity components (u and v) exist only at the openings, and positive velocities of the same magnitude are observed inside the compartment. The positive velocity components result from the circulating flow that entrains back to the plume. Burning at the opening continues until 90 min, after which the fire starts to slowly move back into the compartment. At 120 min, the gas temperatures continue increasing, being 100-200 • C higher than the values at 16 min. The fire continues spreading towards the closed corner, where it ends at 250 min.
In the compartment with symmetric ventilation through four vents, the fire spread process is not as complicated as in the asymmetric situation. Fig. 26 illustrates the flame locations and oxygen concentrations. Due to the symmetric flow conditions, the flame remains above the ignition region longer, showing >700 • C smoke layer temperature (Fig. 27). At 6.0 min, the fire has reached the edge of the prescribed region, and large flames bend over the unburned wood cribs. At 9.0 min, the fire starts to spread outwards against the entrainment gas velocity of about 1.5 m/s (Fig. 28). The temperature contours in Fig. 27 show a preferential fire spread towards the x = 30 m boundary, which is likely due to the channelling effect caused by the arrangement of the wood cribs. By 16 min, the maximum temperatures and flow velocities are observed mainly at the openings, with room centre temperatures of about 500 • C and oxygen concentration of 5%. The fire burns at the openings until the fire loads next to them are consumed, after about 65 min, and the fire is seen to be travelling back towards the compartment centre at 82.5 min. The compartment centre is clearly under-ventilated (Fig. 26), and the temperature is about 800 • C. At 93.5 min, the temperature within the compartment increases as the fire travels back, and by 99.5 min, the fire is at the compartment centre. Temperatures along the corners are reduced and oxygen levels increased close to ambient.
The simulations illustrate how the different physical phenomena contribute to the flame positioning inside a compartment: In the beginning, the flame is affected by entrainment velocity with a dynamic   time scale of few seconds. As soon as the room's internal entrainment cannot feed the required oxygen, the flame starts to lean towards the vents within few tens of seconds. These two phenomena are observed in the asymmetric ventilation scenario, but not so clearly present in the symmetric case. The propagation towards the vents depends on the heat transfer from the flame and smoke layer to the unburned fuel heat transfer, and takes place within about 10 min. Finally, the travelling back is governed by fuel burnout, for which the time scale mainly depends on the fire load density, being in the order of 1 h in the current scenario. The last two phenomena are quite similar in the two ventilation scenarios. Fig. 29 shows the simulated HRR, average temperature and average heat flux to fuel surface within the compartment. The HRR, fire duration and thermal environment are heavily dependent on the ventilation conditions. The HRR curves show three distinct stages: a growth stage (≈15 min), semi-steady burning, and a decay stage. In the four-vent compartment (ψ = 4.50 m 1/2 ), the peak HRR is obtained at the end of the semi-steady stage, after which the decay starts at about 80 min, and the fuel is completely consumed by 102 min. In the simulation of the two-vent compartment (φ = 9.6 m 1/2 ), the peak HRR is 50% lower and it is obtained in the middle of the semi-steady stage. The semi-steady burning period occurs from 20 min to 166 min and the fire lasts 2.5 times longer than in the four-vent compartment. The horizontal lines in the HRR-plot of Fig. 29a show the classical values of the HRR at the onset of flashover (HRR FO = 750A w H 1/2 ) [18]. The simulated HRR barely exceed these values, although the flame shapes and the burning of the entire fuel bed indicate a flashover-like combustion situation.
Average temperatures within the compartment increase rapidly during the growth stage and continue to rise slowly during the semisteady burning stage. The peak values are reached at the end of the semi-steady stage when the flames are in the middle of the room, and the compartment boundaries experience the longest heating. Beyond the temporal difference, the average temperature in the two ventilation scenarios are quite similar. Averaged heat fluxes, in turn, show more drastic differences, with peak values of 96 kW/m 2 and 60 kW/m 2 for the four and two-vent scenarios, respectively.

FDS model
The method of surface area correction was introduced in Ref. [23] to enable large-scale CFD analyses of fire spread scenarios where the fuel consists of objects with length scale that is similar to or finer than the numerical gas phase resolution. The method compensates for the changing fuel surface area, thus allowing the use of coarser grids, but does not ensure grid independent solutions. Besides the fuel area conservation, reaching grid independent fire spread predictions requires that the thermal feedback from the flames and smoke to the surfaces does not vary with resolution.
In FDS, the primary method for producing grid independent flame heat flux is the use of radiative fraction, which ensures that a specified fraction of the volumetric heat release rate is assigned to the emission source term of the radiative transfer equation. In Ref. [23], we investigated the grid dependency in the SP tunnel and square-shaped compartment scenarios. We showed -for the ignition temperature pyrolysis model (ITP) -that the size of the prescribed fire must be scaled according to the domain, so that the prescribed fire can be resolved by the CFD grid but it does not dominate the entire fire behaviour. For instance, in the 30 m × 30 m compartment, it was necessary to prescribe an initial fire that covered 16 m 2 of the floor area. The same prescribed area was found suitable in the current work too, where the pyrolysis was based on the single-step model (SSP) with reaction kinetics. In the simplified tunnel simulations of the current article, the initial HRR was 220 kW, i.e. about one fourth of the predicted steady-state HRR. In the current validation scenarios, the fire load was modelled as a combination of ITP and SSP, and it was possible to use prescribed fires that are comparable to the actual ignition burners used in experiments. However, the attempts to model the entire fire load using SSP (not reported here) required substantially larger prescribed fires, which then started to dominate the thermal environment. The dependency of the predicted HRR on the mesh resolution was demonstrated and quantified in the simplified tunnel scenario.

Influence of fuel geometry
Besides the fuel material and surface area, fire spreading on solid fuels also depends on the geometrical configuration of the fuel. In the square compartment, presented here and in Ref. [23], we found that the convective heat transfer within the cribs affected the burning behaviour. In Ref. [23], we modelled the fire spread using the ITP model with 10 cm and 20 cm resolutions. The results were not very sensitive to the fuel geometry, as a consistent fire spread was observed by maintaining the crib height (0.4 m) across the two mesh sizes. With the current model, using the same fuel height did not produce a consistent fire spread, and the height of the wood cribs (subsequently the number of layers) had to be artificially increased to capture the internal, convective heat transfer processes. This affects the ignition and fire spread on the crib's top surface and reduces the accuracy of modelling. More researchexperimental and numericalis needed for understanding the wood cribs' internal fire dynamics and to propose methods for modelling them with limited mesh resolutions.

Fire spread behaviour
The fire spread scenarios simulated in this work attempt to investigate the different modes of fire propagation on wood cribs and highlight some influencing parameters. The propagation speeds were obtained from the simulations by monitoring the heat release rate per unit volume just above the fuel bed, and calculated as the propagated distance divided by the propagation time. See Table 6 for a summary.  When the ceiling heights were higher than 3.6 m, the fire spread was controlled by the local radiative heating, analogous to the counterflow flame spread. In the simplified tunnel, a travelling fire was formed with a spread rate of 0.13 cm/s, but the BST/FRS compartment with increased ceiling height did not show spread at all. With lower ceiling heights, the smoke layer radiation pre-heated the fuel surface, leading to rapid spread over the entire fuel in flashover manner. In tunnel -like scenarios the spread rates were 5-9 cm/s and in square compartment 2 cm/s. As usual, the post-flashover flames concentrated at the vent, and the fuel deeper in the room stopped pyrolyzing. The last stage of fire in the low ceiling height scenarios was a backward travelling stage. During this stage, the flame geometry and temperature field were controlled by ventilation, while the steady propagation was controlled by the fuel burnout. The speed of backwards travelling varied between 0.1 and 0.8 cm/s. The lack of fire spread in the high-ceiling version of the BST/FRS compartment was most likely caused by the non-uniform fuel loading. The gaps between the piles reduced the level of near-field radiation reaching the unignited fuel. On the other hand, the isolated wood piles made the flow field more turbulent, and identifying the exact position of the fire front became more difficult.
In the square compartment, a two-fold variation of ventilation openings and the simultaneous change of the ventilation arrangement from asymmetric to symmetric had a strong effect on the fire duration (factor 2.5). The forward spread rates of the two ventilation scenarios were the same (see Table 6), and the speeds of backward travelling were within a factor of two (0.2-0.4 cm/s). The averaged temperatures and   heat fluxes were slightly higher in the shorter fire with stronger ventilation. We realise that a much more detailed analysis is necessary to obtain a well-rounded understanding of the factors influencing the spread behaviour. From the viewpoint of adopting the travelling fire model as a design fire describing the temperature field in large open-space compartments, it is important to understand that the rapid fire spread mode, discussed above, was not governed by the burnout of fire load. Nevertheless, the burnout -controlled propagation was observed in these simulations too, but only after the fire had first reached the ventilation opening and stayed there for some time. The simulated temperatures in the tunnellike scenarios were clearly higher than the far-field temperatures predicted by the analytical travelling fire model. The Alpert's ceiling jet equations, used in the analytical model, are not suitable for such geometries, as it was developed for axi-symmetric jets under flat ceilings and cold ambient. In addition to the geometrical discrepancy, a postflashover back-travelling occurs in a significantly pre-heated environment.
The studied scenarios are much smaller than real compartments proposed for spaces such as open-plan offices, which are typically in the range of 1000-5000 m 2 . This is an important limitation in our capability to generalize the conclusions about the observed spread modes. Furthermore, the predictions of the counterflow -type fire spread belong to the most challenging class of fire CFD. Further research is needed to confirm the results on the qualitative spread behaviours and to validate the corresponding predictions of propagation speed. The validations within the current work mainly cover the spread modes under hot smoke layer.
The sensitivity of spread results to modelling choices, such as the combustion model, radiation model parameters and alternate pyrolysis models, was not tested. As the strongest stages of the simulated fires involved clearly under-ventilated conditions, significantly higher yields of CO and soot could be expected, and the conditions of radiation exchange could be affected. The most obvious, uncaptured pyrolysis phenomenon was the lack of char oxidation, which would change the temperature and oxygen concentration of the entrainment flow during the back-travelling stages.

Conclusions
This paper attempts to provide a reliable and robust modelling approach for the CFD-based, predictive fire spread simulations and to explore the effects of ceiling height and fuel and ventilation arrangements on fire spread. At the same time, we want to use the simulations to improve our understanding of fire spread process in general, and to examine the applicability of the travelling fire -concept in describing the simulated spread modes in particular.
Combinations of the ignition temperature and single-step kinetics -based pyrolysis models were used in the validation simulations. The previously developed method for correcting the fuel surface area on coarse meshes was found to work well with the kinetics -based pyrolysis model, in addition to the previously used, ignition temperature model. Consistent HRR results were obtained with three different resolutions, with roughly 20% difference in the semi-steady state HRR between the simulations with 15 cm and 10 cm resolutions. Two different scenarios were used for the validation: In the SP tunnel fire, we found that using an ignition temperature model for the first two cribs produced the best results, with 19% under-prediction of the peak heat release rates. In the BST/FRS large-scale simulations, we used the ignition temperature model for the first row of cribs and the single step pyrolysis model for the rest. The overall burning behaviour and point-wise gas temperatures were in excellent agreement with the experiments, with relative errors in range −6 to +6%. The results show that the modelling approach is robust and can be reliably applied to model fire spread scenarios.
Fine details of the gas velocity, temperature and oxygen volume fractions were reported for the BST/FRS compartment as well as the two artificial scenarios, where the modelling was based on the combination of sufficiently strong prescribed fire and kinetics -based pyrolysis model. Based on the observations, we can identify three different modes of fire spread in compartments with wood crib fire loads: 1. When the ceiling height was high, the heat transfer was mainly local, and the fire started to spread or travel at speed < 1 cm/s, or stopped completely. The spreading process was analogous to the counterflow flame spread. In this mode, vent flow and entrainment provided enough oxygen to support flaming above the pyrolysis zone. 2. The second observed mode was a relatively fast (~ 10 cm/s) spread of flames from the ignition source towards the ventilation opening. It was observed with low ceiling heights when the heat flux to the unburned fuel was high and the oxygen entrainment limited. 3. The third mode of spreading was observed when the fire had reached the ventilation opening and consumed the fire load over there. The rate of spread towards the compartment interior was controlled by the high heat fluxes from the large flame and the rate of fire load burnout. The spread rates were ~ 0.8 cm/s for the 380 MJ/m 2 fire load density of the BST/FRS large-scale compartment, and 0.1 cm/s for the 600 MJ/m 2 fire load of the simplified tunnel.
A three-layer flow field was observed on the closed-end side of the compartment, exhibiting (1) a classical ceiling jet flow under the ceiling, (2) a re-circulation flow towards the fire plume, and (3) a weak but noticeable flow of oxygen-containing gas passing below the fire front, inside the wood cribs. The role of the internal gas flow in the combustion dynamics of burning wood cribs would deserve a more detailed examination, along with the observed difference in the turbulence levels between the uniform and discrete fire load arrangements.
Regarding the application of the travelling fire concept in modelling the observed fire behaviours, we can conclude that it fundamentally corresponds to the first mode of fire spread, but could also be used to describe the third mode. The temperatures given by the analytical iTFM model and the CFD were compared in the back-propagation stages of the tunnel -like scenarios. In the 3.0 m high tunnel, the analytical model gave a 21% slower propagation speed than the CFD model with the same fire load density, an arbitrary 3.4% burning region and a calibrated HRRPUA. Further calibration could lead to even better agreement between the models. The near-field temperatures, as given by the analytical model, were within few percent from those predicted by the CFD, but far from the burning zone, the analytical model showed only few degrees of temperature increase, while the CFD predictions were in the range 100-200 • C. An application of the travelling fire concept to the third, back-propagation mode would thus require consideration of the pre-heating prior to the back -propagation stage. Additionally, the accumulation of heat in narrow, tunnel -like compartments results in much higher far-field temperature than what is given by the axisymmetric ceiling jet correlation used in the analytical model. Nevertheless, the results show that CFD simulations can provide valuable insights into the fire spread phenomenon, and serve to further the development of analytical models intended for every-day design purposes.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.