Mostafa Taha

Song Zhao

Aymeric Lamorlette

Jean-Louis Consalvi

Pierre Boivin
email: pierre.boivin@univ-amu.fr

Large eddy simulation of fire-induced flows using Lattice-Boltzmann methods

Large-eddy simulations (LES) of the near-field region of large-scale fire plumes are performed for the first time with a pressure-based Lattice Boltzmann method (LBM) with low-Mach number approximation. Two scenarios are considered: the large-scale non-reactive helium plume and the 1 m methane pool fire, both investigated experimentally at Sandia. In the second scenario, a simplified modeling of the combustion and radiation processes is introduced involving a one-step irreversible reaction eddydissipation concept-based combustion model and a radiant fraction model, respectively. In both scenarios, a quantitative agreement is observed with the experimental data and model predictions are consistent with previouslypublished numerical studies. Our simulations demonstrate the computational efficiency of the proposed LBM solver to tackle fire-induced flows, suggesting that LBMs are a good alternative candidate for the modeling of fire-related problems.

Introduction

Computational Fluid Dynamics models play today an increasing role in fire safety engineering with the development of LES-based fire simulators such as FDS by the NIST and FireFoam by FM Global [START_REF] Mcgrattan | Dynamics simulator technical reference guide[END_REF][START_REF] Firefoam | [END_REF]. The near field of well-controlled fire plumes offers a canonical scenario for the validation of such simulators since it involves most of the coupled physical processes encountered in fire problems, namely, buoyancy-controlled flows, buoyancy-induced turbulence, turbulent combustion, thermal radiation, soot generation, and, in the case of pool fires, burning rate [START_REF] Brown | Proceedings of the first workshop organized by the iafss working group on measurement and computation of fire phenomena (macfp)[END_REF].

LES of non-sooting and sooting fire plumes were reported in the literature with different levels of sophistication in the modeling of subrid-scale turbulence [START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF][START_REF] Desjardin | Modeling of conditional dissipation rate for flamelet models with application to large eddy simulation of fire plumes[END_REF][START_REF] Ferraris | Large-eddy simulation of a large-scale methane pool fire[END_REF][START_REF] Maragkos | Application of fds and firefoam in large eddy simulations of a turbulent buoyant helium plume[END_REF][START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF][START_REF] Maragkos | Large eddy simulations of ch4 fire plumes[END_REF][START_REF] Maragkos | Towards predictive simulations of gaseous pool fires[END_REF][START_REF] Maragkos | On the use of dynamic turbulence modelling in fire applications[END_REF][START_REF] Ahmed | Large eddy simulation of the unstable flame structure and gas-to-liquid thermal feedback in a medium-scale methanol pool fire[END_REF][START_REF] Ma | Large eddy simulation of medium-scale methanol pool fires -effects of pool boundary conditions[END_REF][START_REF] Ma | Exploring subgrid-scale variance models in les of lab-scale methane fire plumes[END_REF], turbulent combustion [START_REF] Ma | Exploring subgrid-scale variance models in les of lab-scale methane fire plumes[END_REF][START_REF] Wang | Large eddy simulation of fire plumes[END_REF][START_REF] Chen | Large eddy simulation of a medium-scale methanol pool fire using the extended eddy dissipation concept[END_REF], radiative heat transfer [START_REF] Chatterjee | A model for soot radiation in buoyant diffusion flames[END_REF][START_REF] Sikic | Non-grey radiative heat transfer modelling in les-cfd simulated methanol pool fires[END_REF][START_REF] Nmira | Influence of gas radiative property models on large eddy simulation of 1 m methanol pool fires[END_REF], and soot modeling [START_REF] Chen | Large eddy simulation of a medium-scale methanol pool fire using the extended eddy dissipation concept[END_REF][START_REF] Chatterjee | Application of a subgrid sootradiation model in the numerical simulation of a heptane pool fire[END_REF][START_REF] Nmira | Large-eddy simulation of lab-scale ethylene buoyant diffusion flames: Effects of subgrid turbulence/soot production interaction and radiation models[END_REF]. One of the difficulty in the modeling of the near field of fire plumes is that the flow, and the resulting air entrainment that controls the combustion process, is governed by the formation and growth of the flame base non-dissipative laminar instability near the edge of the pool that develops periodically to form energy containing large-scale toroidal vortices [START_REF] Cetegen | Experiments on the periodic instability of buoyant plumes and pool fires[END_REF][START_REF] Tieszen | Transport phenomena that affect heat transfer in fully turbulent fires[END_REF]. A consequence is that these instabilities have to be spatially-resolved which requires a grid-resolution of the order of 2 mm [START_REF] Ahmed | Large eddy simulation of the unstable flame structure and gas-to-liquid thermal feedback in a medium-scale methanol pool fire[END_REF][START_REF] Ma | Large eddy simulation of medium-scale methanol pool fires -effects of pool boundary conditions[END_REF][START_REF] Ma | Exploring subgrid-scale variance models in les of lab-scale methane fire plumes[END_REF].

The LES of fire plumes reported in the literature were performed by using classical Navier-Stokes solvers. On the other hand, Lattice Boltzmann Methods (LBM) appear as an attractive alternative to the classical Navier-Stokes solvers due to their lower computational cost, their capacities for massively parallel computing and the ease to deal with complex geometries using multi-level Cartesian grids [START_REF] Krüger | The lattice boltzmann method[END_REF]. These methods have been successfully applied to different classes of problems, including compressible flows [START_REF] Feng | A three dimensional lattice model for thermal compressible flow on standard lattices[END_REF][START_REF] Feng | A compressible lattice boltzmann finite volume model for high subsonic and transonic flows on regular lattices[END_REF][START_REF] Guo | An efficient lattice boltzmann method for compressible aerodynamics on d3q19 lattice[END_REF][START_REF] Farag | A pressurebased regularized lattice-boltzmann method for the simulation of compressible flows[END_REF][START_REF] Farag | A unified hybrid lattice-boltzmann method for compressible flows: bridging between pressure-based and density-based methods[END_REF][START_REF] Coratger | Large-eddy lattice-boltzmann modelling of transonic flows[END_REF], atmospheric flows [START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows[END_REF][START_REF] Cheylan | Compressible pressure-based latticeboltzmann applied to humid air with phase change[END_REF][START_REF] Jacob | Wind comfort assessment by means of large eddy simulation with lattice boltzmann method in full scale city area[END_REF][START_REF] Merlier | Lattice-boltzmann large-eddy simulation of pollutant dispersion in street canyons including tree planting effects[END_REF][START_REF] Merlier | Lattice-boltzmann large-eddy simulation of pollutant dispersion in complex urban environment with dense gas effect: Model evaluation and flow analysis[END_REF], reactive flows [START_REF] Tayyab | Experimental and numerical lattice-boltzmann investigation of the darrieus-landau instability[END_REF][START_REF] Tayyab | Hybrid regularized lattice-boltzmann modelling of premixed and non-premixed combustion processes[END_REF][START_REF] Tayyab | Lattice-boltzmann modelling of a turbulent bluffbody stabilized flame[END_REF][START_REF] Bhairapurada | A lattice-boltzmann study of premixed flames thermo-acoustic instabilities[END_REF] and, recently, to buoyancydriven flows encountered in the far field of plumes [START_REF] Taha | Lattice-boltzmann modeling of buoyancy-driven turbulent flows[END_REF].

The objective of this article is to explore for the first time the capacity of a pressure-based LBM solver with low-Mach number approximation [START_REF] Wang | A new hybrid lattice-boltzmann method for thermal flow simulations in low-mach number approximation[END_REF] to handle fire-related flows. An hybrid strategy is considered where continuity and momentum are solved with the LBM equation and transport equations for other scalars (enthalpy, chemical species,...) are solved by classical finite difference methods. The numerical model is applied to provide LES of two scenarios representative of unwanted fires. The first involves the large-scale helium plume investigated at Sandia [START_REF] Tieszen | Spatial and temporal resolution of fluid flows: Ldrd final report[END_REF][START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. The main advantage of this configuration is that it mimics the dynamics and structure of large-scale fire plumes while avoiding the complexities associated with combustion and radiation. The second configuration is relative to the 1 m methane fire plumes also investigated at Sandia [START_REF] Tieszen | Experimental study of the flow field in and around a one meter diameter methane fire[END_REF][START_REF] Tieszen | Experimental study of the effect of fuel mass flux on a 1-m-diameter methane fire and comparison with a hydrogen fire[END_REF]. This scenario enhances the complexity of the first as it involves turbulent combustion and thermal radiation.

The article is organized as follows. A detailed description of the experimental setup is first given. Section 2 presents the numerical method, including the LBM strategy, the different turbulence, combustion and radiation models. The results for the helium plume scenario and the methane pool fire are presented in Sections 3 and 4, respectively. Conclusions are finally drawn to highlight the findings of this study and suggestions for future work.

Experimental setup

Comparisons are made in this study with experiments performed in the Fire Laboratory for Accreditation of Models and Experiments (FLAME) at Sandia National Laboratories in Albuquerque, New Mexico [START_REF] Tieszen | Spatial and temporal resolution of fluid flows: Ldrd final report[END_REF]. The considered experiments involve either a non-reactive helium plume scenario or a reactive scenario with methane as fuel. The central FLAME chamber consists of a 6.1 m cubical enclosure with a 2.4 m diameter chimney.

The plume source is 1 m in diameter surrounded by a 0.51 m wide floor (i.e. the ground plane). Further details about the facility can be found in Tieszen et al. [START_REF] Tieszen | Spatial and temporal resolution of fluid flows: Ldrd final report[END_REF].

For the non-reactive case, a gas mixture, composed of 96.4% helium, 1.7% acetone, and 1.9% oxygen by volume, with a molecular weight of 5.45 g/mol, was injected at an average velocity of 0.325 m/s. Measurements were obtained using particle image velocimetry (PIV) for velocity field measurements and planar laser-induced fluorescence (PLIF) for helium mass fraction [START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. The experimental uncertainty on the measured velocities, concentrations, turbulent statistics, and concentration fluctuations is estimated to be around ±20, ±18, ±30, and ±21%, respectively.

For the methane pool fires, Tieszen et al. [START_REF] Tieszen | Experimental study of the flow field in and around a one meter diameter methane fire[END_REF][START_REF] Tieszen | Experimental study of the effect of fuel mass flux on a 1-m-diameter methane fire and comparison with a hydrogen fire[END_REF] measured the velocity field for methane fire plumes with different heat release rates ranging from 1.56 to 2.61 MW. Nonetheless, no measurements for temperature were performed. Numerous tests were conducted during the experimental study [START_REF] Tieszen | Experimental study of the flow field in and around a one meter diameter methane fire[END_REF][START_REF] Tieszen | Experimental study of the effect of fuel mass flux on a 1-m-diameter methane fire and comparison with a hydrogen fire[END_REF] where the fuel inlet conditions were varied to cover a wide range of fire regimes. TEST-24 is chosen as a reference for our numerical study and whose inlet and boundary conditions are summarized in Tab. 1.

Lattice-Boltzmann model for Low-Mach flows

Macroscopic governing equations

The leading order flow mass, momentum, and energy conservation equations are introduced as [START_REF] Wang | A new hybrid lattice-boltzmann method for thermal flow simulations in low-mach number approximation[END_REF]

∂ρ ∂t + ∂ρu i ∂x i = 0 (1)
∂ρu i ∂t + ∂(ρu i u j + δ ij p h ) ∂x j = ∂Π ij ∂x j + (ρ -ρ ∞ )g i , (i = 1, 2, 3) (2) ∇P = 0 (3)
ρ ∂h ∂t + ρu i ∂h ∂x i = dP dt - ∂q i ∂x i + Π ij ∂u i ∂x j , (4)
ρ ∂Y k ∂t + ρu i ∂Y k ∂x i = ∂ ∂x i (ρY k V k,i ) + ωk , (5)
where ρ is the density, u i is the velocity vector, p h is the hydrodynamic pressure, P is the thermodynamic pressure, g i is the gravitational acceleration, δ ij as the Kronecker delta symbol, ρ ∞ is the ambient density. The viscous stress tensor Π ij in Eqs. [START_REF] Firefoam | [END_REF][START_REF] Brown | Proceedings of the first workshop organized by the iafss working group on measurement and computation of fire phenomena (macfp)[END_REF][START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF] reads:

Π ij = µ ∂u i ∂x j + ∂u j ∂x i -δ ij 2 3 ∂u k ∂x k , (6)
with µ the dynamic viscosity. In enthalpy equation (Eqn. 4), h represents the mass enthalpy. Throughout this work, we set the temporal variation of the thermodynamic pressure to zero, i.e., dP/dt = 0, to meet the open boundary conditions [START_REF] Wang | A new hybrid lattice-boltzmann method for thermal flow simulations in low-mach number approximation[END_REF][START_REF] Poinsot | Theoretical and numerical combustion[END_REF]. The heat flux q i in the enthalpy equation reads:

q i = -λ ∂T ∂x i + Nsp k=1 ρh k Y k V k,i , (7)
where T is the temperature, λ the heat conductivity, obtained assuming constant Prandtl number (Pr):

Pr = c p µ λ . ( 8
)
h k is the partial enthalpy of species k, Y k is its mass fraction, and V k,i is its diffusion velocity vector. The term Y k V k,i constitutes the species mass flux. The species mass flux is modeled by Hirschfelder and Curtiss approximation equipped with a correction term ensuring global mass conservation:

Y k V k,i = -D k W k W ∂X k ∂x i + V c i Y k , (9)
where D k is the molecular diffusion coefficient of species k, W k is its molecular weight and X k is its molar fraction. W is the average molecular weight and V c α the correction velocity introduced in order to ensure the conservation of total mass (i.e., ensuring Nsp k=1 Y k V k,α = 0) which can be evaluated by:

V c α = Nsp k=1 D k W k W ∂X k ∂x α . ( 10
)
D k is linked to the viscosity through Schmidt number (Sc) as:

Sc = µ ρD k . (11)
The system of Eqs (1-5) is fully closed by the choice of an equation of state:

P = ρ.r.T , h = Nsp k=1 h k Y k = Nsp k=1 Y k ( T T 0 c p,k (θ)dθ + ∆h 0 k ) (12)
with P the thermodynamic pressure, c p,k and ∆h 0 k the constant pressure heat capacity and the formation enthalpy of species k respectively, and r = R/W the gas constant where R = 8.314J K -1 is the universal gas constant and its molecular weight is W.

Turbulence modeling

The filtered expressions for Eqs. 1-5 are widely reported in the literature (see, e.g., Ref. [START_REF] Poinsot | Theoretical and numerical combustion[END_REF]) and not recalled here. Subgrid-scale models (i.e. turbulence models) are essential to estimate the subgrid terms appearing in the filtered equations. Applying the subgrid-scale model numerically comes down to modifying the viscosity µ through the addition of a turbulent viscosity µ t . Turbulent models aimed to estimate the correct value of µ t that compensate the non-resolved structures filtered by the mesh (i.e. smaller than the mesh size).

Smagorinsky model

The first model to be introduced is the Smagorinsky model [START_REF] Smagorinsky | General circulation experiments with the primitive equations: I. the basic experiment[END_REF]. It is constructed of a length scale term, a time scale length, and a dimensionless constant C s called Smagorinsky constant. The turbulent viscosity µ t is calculated through:

µ t = ρ(C s ∆ m ) S , (13)
where ∆ m is the local mesh size, and S is the filtered rate of stress tensor written as:

S = 2 S ij S ij , (14)
with S ij reads:

S ij = 1 2 ( ∂ u i ∂x j + ∂ u j ∂x i ), (15)
Care should be taken when using the Smagorinsky model because it has many drawbacks listed below:

• C s is a predefined input to the simulation, this constant can not represent correctly various turbulent flows.

• The eddy viscosity does not vanish for laminar flows.

• The back-scatter of energy from small scale to large scale is prevented completely since:

(C s ∆ 2 m ) 2S ij S ij ≥ 0 (16)
• Smagorinsky model introduces too much diffusion to the flow.

Vreman model

Inspired from the Smagorinsky model, Vreman subgrid model was introduced [START_REF] Vreman | An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications[END_REF], where the µ t is described as:

µ t = ρC B β α ij α ij , (17)
with,

α ij = ∂u j ∂x i , (18)
β ij = ∆ 2 m α mi α mj , (19)
B β = β 11 β 22 -β 2 12 + β 11 β 33 -β 2 13 + β 22 β 33 -β 2 23 , (20) 6
The constant C is related to the Smagorinsky constant C s as C = 2.5C 2 s . Vreman model was developed in a try to fix the downsides of Smagorinsky model. Both Smagorinsky and Vreman models are simple to implement as they only require the local filter width (i.e. mesh size) and the first order derivatives of the velocity field.

Combustion modeling

The combustion modeling aims to compute the species production term ωk (i.e. combustion source term) in the species conservation equation (Eqn. 5). The eddy dissipation concept (EDC), which is commonly incorporated in fire simulators such as FDS [START_REF] Mcgrattan | Fire dynamics simulator user's guide[END_REF] and FireFoam [START_REF] Firefoam | [END_REF], is adopted in the present study [START_REF] Magnussen | On mathematical modeling of turbulent combustion with special emphasis on soot formation and combustion[END_REF]:

ωF = C EDC ρ 1 τ t min Y F , Y O s ( 21
)
where ωF is the fuel mass consumption rate per unit volume (i.e. Kg m -3 s -1 ), C EDC is the model constant, ρ is the averaged density, τ t is the turbulent/mixing time scale, Y F and Y F are the filtered mass fraction of fuel and oxidizer respectively and finally s is the stoichiometric oxygen-to-fuel mass ratio. The main parameter that controls and defines the dynamics of the EDC model is the turbulent/mixing time scale τ t . Following the work of Yaga et al. [START_REF] Yaga | Modeling of eddy characteristic time in les for calculating turbulent diffusion flame[END_REF], the C EDC is set to 4.0 whereas the eddy characteristic time τ t is estimated by considering the Kolomogorov scale as:

τ t = ν ϵ 1 2 , ( 22
)
where ϵ is the dissipation rate of turbulent kinetic energy and ν is the kinematic molecular viscosity. Assuming local equilibrium between production and dissipation of turbulent kinetic energy (same hypothesis adopted while deriving Smagorinsky SGS model), eddy dissipation rate ϵ is written as follows:

ϵ = 2ν t S ij : S ij , (23)
where ν t is kinematic turbulent viscosity and S ij is the strain rate, both described in Sec. 2.2. The eddy characteristic time scale will then be:

τ t = ν 2ν t S ij : S ij 1 2 . ( 24
)
Note the term S ij : S ij is an inner product and should calculated as follows:

S ij : S ij = S 11 2 + S 22 2 + S 33 2 + 2 S 12 2 + S 13 2 + S 23 2 . (25)
A simple radiant fraction-based model is introduced to account for radiative loss. In this approach, commonly used in fire simulators [START_REF] Mcgrattan | Dynamics simulator technical reference guide[END_REF], a specified part of the local heat release rate is assumed to be lost by radiation [START_REF] Mcgrattan | Dynamics simulator technical reference guide[END_REF][START_REF] Firefoam | [END_REF].

LBM description for Low-Mach number approximation

For a complete description of the numerical method, the reader is referred to Farag et al. [START_REF] Farag | A pressurebased regularized lattice-boltzmann method for the simulation of compressible flows[END_REF][START_REF] Farag | A unified hybrid lattice-boltzmann method for compressible flows: bridging between pressure-based and density-based methods[END_REF]. Lattice-Boltzmann methods are derived from a space, time, and velocity discretization of the Boltzmann equation [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF]. In the present model, the probability density function f i is solved at each point x via the Boltzmann equation discretized as:

f i (x + c i δt, t + δt) = f eq i (x, t) + 1 - δt τ f neq i (x, t) + δt 2 F E i (x, t) (26)
where δt is the time-step, and c i is the i th discrete velocity of the D3Q19 lattice [START_REF] Krüger | The lattice boltzmann method[END_REF], and F E i is a volume force including gravity and correcting terms, as defined in Appendix A. The equilibrium and off-equilibrium populations (f eq i , f neq i ) are to be defined in Eqs [START_REF] Farag | A pressurebased regularized lattice-boltzmann method for the simulation of compressible flows[END_REF][START_REF] Farag | A unified hybrid lattice-boltzmann method for compressible flows: bridging between pressure-based and density-based methods[END_REF]. In Eq. ( 26), the relaxation time τ is related to the dynamic viscosity as:

τ = µ ρc 2 s + δt 2 , (27)
where c s = δx/( √ 3δt) is the characteristic velocity of the D3Q19 lattice [START_REF] Krüger | The lattice boltzmann method[END_REF].

The equilibrium function is obtained as

f eq i = ω i ρθ + H (1) i,α c 2 s ρu α + H (2) i,αβ 2c 4 s ρu α u β + a (3) + a (4) , (28)
where ω i is the D3Q19 weight of discrete velocity c i , H i are the discrete Hermite polynomials, defined in Appendix A, and a (3) and a (4) are third and forth order terms also provided in Appendix A.

The off-equilibrium counterpart f neq is obtained as:

f neq i = ω i H (2) i,αβ 2c 4 s a (2),neq αβ + H (3r) i,γ 6c 6 s a (3r),neq γ , (29) with a
(2),neq αβ

≡ a * (2),neq αβ - δ αβ 3 a * (2),neq γγ , (30)
a * (2),neq αβ = i f i -f eq i + δt 2 F E i H (2) i,αβ . (31)
and the third-order contribution defined in Appendix A.

Following the work introduced recently by Wang et al. [START_REF] Wang | A new hybrid lattice-boltzmann method for thermal flow simulations in low-mach number approximation[END_REF], the compressible description of LBM can be adjusted to follow the low-Mach number approximation (LMNA). This approach divides the pressure field into the thermodynamic part, p th , and the hydrodynamic part, P . The thermodynamic pressure is linked to the equation of state and is spatial independent. Thus, only the hydrodynamic pressure appears in the momentum conservation equation. In practice, the density field is updated using the thermodynamic pressure and the temperature field, i.e.

ρ(x, t + δt) = P (t + δt) rT (x, t + δt) . (32)
For the simulations in this paper, the thermodynamic pressure P is assumed to be constant, considering the open boundary conditions. The hydrodynamic pressure cooperates with the LBM system through the reduced temperature

θ ≡ p h ρc 2 s , (33)
and p h will be updated following

p h (x, t+δt) = p h (x, t)+c 2 s i f col i (x -c i δt, t) -f col i (x, t) -(ρ(x, t + δt) -ρ(x, t)) .
(34) The velocity fields are calculated through the first-order moment of the LBM population

ρu i (x, t + δt) = i c i f i (x, t + δt) + δt 2 F g i , (35)
where F g i is the gravity force term defined in Appendix A. The enthalpy equation ( 4) is solved at the same time using a finite difference discretization under the non-conservative form, exactly as presented by Tayyab et al. [START_REF] Tayyab | Experimental and numerical lattice-boltzmann investigation of the darrieus-landau instability[END_REF][START_REF] Tayyab | Hybrid regularized lattice-boltzmann modelling of premixed and non-premixed combustion processes[END_REF]. Second-order consistency to the macroscopic equations (1-4) can be shown via Chapman-Enskog [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF], or Taylor [START_REF] Farag | Consistency study of lattice-boltzmann schemes macroscopic limit[END_REF][START_REF] Wissocq | Hydrodynamic limits and numerical errors of isothermal lattice boltzmann schemes[END_REF] expansions.

Sandia's 1-m diameter helium plume

Numerical setup

The simulation is performed on a Cartesian mesh for a cubic domain sizing 8 × 8 × 4 m 3 depicted in Fig. 1. The choice of the domain size as well as the grid resolution is consistent with previous numerical studies [START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF][START_REF] Maragkos | Application of fds and firefoam in large eddy simulations of a turbulent buoyant helium plume[END_REF][START_REF] Chung | Buoyancy-corrected k-ε models and large eddy simulation applied to a large axisymmetric helium plume[END_REF][START_REF] Ma | Verification and validation of a variable-density solver for fire safety applications[END_REF].

Three meshes were included in our study, and refinements zones were used to reduce the calculation cost. The zones of refinements (Fig. 1) are denoted by ZoneI, ZoneII, ZoneIII and ZoneIV with grid size reduced in half between two adjacent zones. The coarse mesh contains ZoneI, ZoneII and ZoneIII corresponding to grid size of 8 cm, 4 cm and 2 cm, respectively, resulting in a 2 million-cell mesh. This coarse grid will be referred to as ∆x = 2 cm hereafter. In the medium and fine meshes, an additional refinement zone is added close to the plume source, indicated by ZoneVI in Fig. 1, to end up with four levels of refinements. For the medium grid, the largest and smallest grid size are of 8 cm and 1 cm, respectively, resulting in a 4 million-cell mesh. This medium grid will be referred to as ∆x = 1 cm hereafter. For the fine, the mesh sizes are deduced from the medium grid by applying a scaling factor of 0.8 in each zone. A consequence is that the the largest and smallest grid size are of 6.4cm and 0.8 cm, respectively, resulting in a 7.8 million-cell mesh. This fine grid will be referred to as ∆x = 0.8 cm hereafter.

Boundary conditions are given in Fig. 1. Helium is injected at the center of the bottom plane through a 1 m diameter inlet with a uniform vertical velocity U He inlet = 0.325 m/s. The helium inlet is surrounded by an annulus solid plate of 0.5 m width beyond which a co-flow of air is injected at a much lower speed U Air cof low = 0.01 m/s. At the upper outlet plane a Dirichlet boundary condition is used for pressure and Neumann for other variables, with a clip on the streamwise velocity to prevent any backflow from the outlet. The sides are set to a typical In/Out flow.

At every point in the domain, the composition of the fluid corresponds to a mixture of helium and air. The air was treated as a whole specie without separating Oxygen and Nitrogen with an average molecular weight of the air is taken W air = 28.9g mol -1 . For experimental visualization reasons, the injected helium was not pure, but rather a mixture of 96.4% helium (He), 1.9% oxygen (O 2 ) and 1.7% acetone (CH 3 COCH 3 ) resulting in an average molecular weight W He = 5.45g mol -1 . The ambient (also inlet) temperature and pressure are T ∞ = 285 K and p ∞ = 80900Pa, respectively, in order to match the experimental setup.

The Prandtl number Pr was set to 0.7 while the Schmidt number of helium Sc He is set to 0.2. The dynamic viscosity µ has a constant value of 1.8877 × 10 -5 kg m -1 s -1 for both air and helium. Resulting in a Reynolds number Re ≃ 3220 and a Richardson number Ri ≃ 75.4, matching the experiment.

In the case of homogeneous isotropic turbulence, the model constant is estimated as C s = 0.2 [START_REF] Poinsot | Theoretical and numerical combustion[END_REF]. However, C s depends on the flow configuration with values of C s ≈ 0.1 -0.2 often used [START_REF] Garnier | Large eddy simulation for compressible flows[END_REF]. For this simulation, the Smagorinsky model is used (Sec. 2.2) with a constant C s = 0.1 as suggested by Maragkos et al. [START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF].

The un-resolved SGS heat and species diffusion fluxes are estimated using the turbulent Prandtl Pr t and the turbulent Schmidt Sc t , respectively. Both are set to 0.5 as Maragkos et al. [START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF] and Chung et al. [START_REF] Chung | Buoyancy-corrected k-ε models and large eddy simulation applied to a large axisymmetric helium plume[END_REF].

The large-eddy simulation is run for 26 s of physical time. The first 13 s are to ensure the evacuation of the initial field and to reach statistically stationary flow conditions, and the other 13 s are used to construct the mean and RMS quantities of the plume. The simulations are performed at a constant CFL number equaling to 0.6 such that the time step would depend on the local mesh size. The CFL number is defined as:

CFL u = (u + c black )δt δx (36)
where c black denotes the fictitious reduced speed of sound following the low-Mach number approximation. The time-step is local (it is halved/doubled at each mesh transition), and is δt = 38 µs on the smallest mesh region (where ∆x = 0.8cm). In other words, the CFL number ( 36) is kept almost constant in the whole domain. The simulations are parallelized on 120 cores on a Dell PowerEdge C6420 server with 4× 32-core Intel Xeon Gold 6142, 2.6 GHz, and 96 GB RAM. The computational cost of the simulation with the finest mesh is compared to other studies in Table 2. The corresponding reduced computational time (cost per iteration per grid point) is close to 2.4µs, an excellent score when comparing with available computational benchmarks [START_REF] Boivin | Benchmarking a lattice-boltzmann solver for reactive flows: Is the method worth the effort for combustion?[END_REF][START_REF] Abdelsamie | The taylor-green vortex as a benchmark for high-fidelity combustion simulations using low-mach solvers[END_REF][START_REF] Suss | Comprehensive comparison between the lattice boltzmann and navier-stokes methods for aerodynamic and aeroacoustic applications[END_REF]. Another advantage of the present approach is that complex geometry can be easily tackled at no extra cost [START_REF] Zhao | Lattice-boltzmann modelling of the quiet and unstable preccinsta burner modes[END_REF]. Finally, 12

the considered Lattice-Boltzmann regularized kernel is well suited for crossplatform cpu-gpu usage [START_REF] Latt | Cross-platform programming model for many-core lattice boltzmann simulations[END_REF], allowing to envision significant speedups in the future.

Reference

Results and discussions 3.2.1. Turbulence

Figure 2 shows the energy spectrum based on the axial velocity at distance z = 0.5 m. The -5/3 power law, predicted by the the theory of Kolmogorov [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large reynolds numbers[END_REF], is well reproduced.

Instability modes and vortex dynamics

Figure 3 shows a typical puffing cycle. A puffing cycle can be divided into four distinct phases: In the first phase, instabilities form near the edge of the plume as observed in Fig. 3(a). The trigger of these instabilities is the misalignment of pressure and density gradients which generates a localized torque (baroclinic torque). Those instabilities grow and entrain large quantity of the surrounding fluid, forming into a toroidal vortex as shown in Fig. 3(b). While this vortex moves upwards more fluid is pumped at the center of the plume causing a large streamwise velocity at the centerline as shown in Fig. 3(c). The increase in streamwise velocity at the centerline causes a sharp increase in cross-stream velocity near the base of the plume, as a result of mass conservation, providing a perturbation of the mixing layer for the next cycle (Fig. 3(d)).

Strong instabilities are observed at the helium-air interface near the plume source, and it is important to capture those small dynamics as discussed in the introduction [START_REF] Ahmed | Large eddy simulation of the unstable flame structure and gas-to-liquid thermal feedback in a medium-scale methanol pool fire[END_REF][START_REF] Ma | Large eddy simulation of medium-scale methanol pool fires -effects of pool boundary conditions[END_REF][START_REF] Ma | Exploring subgrid-scale variance models in les of lab-scale methane fire plumes[END_REF]. These instabilities are generated by buoyancy-driven (gravitational and baroclinic) vorticity generation. Thus, Rayleigh-Taylor and Kelvin-Helmholtz instabilities are triggered, near the base of the plume, and at last they will form the toroidal structures [START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF][START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF] as shown in Fig. 3(d).

A better understanding can be made from the analysis of the vorticity equation:

Dω Dt = (ω • ∇)u vortex stretching -ω(∇ • u) dilatation term + 1 ρ 2 (∇ρ × ∇p) baroclinic torque + ρ ∞ ρ 2 (∇ρ × g) gravitational torque + ∇ × ( 1 ρ ∇ • τ ) viscous diffusion (37)
On the right hand side of Eq. ( 37), five different physical mechanisms affect the vorticity transport. Those five terms are identified as vortex stretching, dilatation term, baroclinic torque, gravitational torque and viscous diffusion [START_REF] Jiang | Spatial direct numerical simulation of the large vortical structures in forced plumes[END_REF]. In incompressible flows, only the first term prevails. Vortex stretching represents the enhancement of vorticity by stretching, it is an essential mechanism by which turbulent energy is transferred to smaller scales. The dilatation term represents the effects of vorticity field expansion which causes a decrease in vorticity magnitude. Gravitational torque generates vorticity due to the misalignment of gravity and density gradients, while baroclinic torque generates vorticity as a result of non-aligned pressure and density gradients. In buoyancy-driven flows, both the baroclinic and gravitational torque terms are the principal mechanisms generating flow vorticity [START_REF] Jiang | Spatial direct numerical simulation of the large vortical structures in forced plumes[END_REF][START_REF] Jiang | Spatial dns of flow transition of a rectangular buoyant reacting free-jet[END_REF]. The DNS [START_REF] Jiang | Dynamics and structure of transitional buoyant jet diffusion flames with side-wall effects[END_REF] demonstrated that the gravitational torque is fundamental mechanism promoting cross-stream vorticity. LES studies [START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF][START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF] observed that during typical puffing cycle, the maximum grav-itational torque is located at the base of the plume. The same behavior can be observed in the current study, as seen in Fig. 3 where the maximum value of gravitational torque is observed at the base of the plume. The contribution of the baroclinic torque has similar importance. Maximum values can be located in zones where large pressure and density gradients exist near the base of the plume. Both torques trigger the instabilities at the base and then nurture its growth. For large plumes and pool fires, this vortex destabilizes rapidly forming secondary azimuthal, or "finger-like" instabilities that can be detected in Fig. 3, which were also observed experimentally [START_REF] Cetegen | Experiments on the periodic instability of buoyant plumes and pool fires[END_REF][START_REF] Weckman | The oscillatory behaviour of medium-scale pool fires[END_REF]. Those secondary instabilities create streamwise vorticity that promote the breakdown of large scale toroidal structures and improves the local mixing eventually. Capturing these instabilities and mechanisms is therefore a prerequisite for pool fire simulation because the combustion process for this type of flows is controlled mainly by the mixing process of fuel and oxidizer.

Puffing frequency

Figure 4 shows both the temporal signal of the centerline streamwise velocity at z = 0.5 m and the corresponding power density spectrum obtained from Fast Fourier transform (FFT). The signal extends over 5 s of physical time, from 13 to 18 s, in order to be compared with the experimental one. It is worth mentioning that results revealed a small phase shift compared to the experiment, probably because we do not have the exact same initialization as in the experiment. Thus, the simulation time signals in Fig. 4 have been shifted along the time axis to match the first peak of the experimental signal for the sake of clearer comparison.

Around 9 cycles are identified for the coarse mesh and 7 cycles for both medium and fine meshes, versus 7 cycles found in the experiment [START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. In addition, the predicted puffing frequency is 1.69 Hz for the coarse grid and 1.39 Hz for both the medium and fine grids whereas the observed experimental frequency was 1.37 ± 0.1 Hz. These results suggest that a converged solution is obtained for the medium grid. This will be further demonstrated when discussing the profiles of velocities and helium mass fraction. The predicted puffing frequency on both the medium and fine meshes is also consistent with the experimental correlation proposed by Cetegen and Kasper [START_REF] Cetegen | Experiments on the oscillatory behavior of buoyant plumes of helium and helium-air mixtures[END_REF] f = 0.8Ri 0.38 U p /D p = 1.34 Hz from measurements of heliumair plumes for Ri < 100. In this expression, Ri is a modified Richardson number defined as Ri 1 = (ρ ∞ -ρ)gD p /(ρ ∞ U 2 p ) ≃ 76. Furthermore, the results are also inline with the experimental correlation given by Cetegen and Ahmed [START_REF] Cetegen | Experiments on the periodic instability of buoyant plumes and pool fires[END_REF] for buoyant diffusion flames (i.e. pool fire) of various fuels f = 1.5/ D p = 1.5 Hz. The disparity from this latter correlation may be attributed to the fact that it was suggested for diffusion flames. These results highlight that the finer mesh enables the correct prediction of the instabilities generated near the inlet which is the main trigger of the puffing motion. The same conclusion was drawn by Maragkos et al. [START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF], DesJardin et al. [START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF] and Ma et al. [START_REF] Ma | Verification and validation of a variable-density solver for fire safety applications[END_REF]. The sampling time ts to construct the power spectrum is taken to be each 10 time steps that is ts coarse = 9.6 × 10 -4 sec, ts medium = 4.8 × 10 -4 sec and ts f ine = 3.8 × 10 -4 sec for the coarse, medium and fine mesh, respectively.

Statistics of velocity and mass fraction

The mean values and high order statistics (i.e. rms) are discussed in this section. In all the figures the data with bars represents the experimental results of O'Hern et al. [START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF] including the experimental uncertainties as explained in Sec. 1. Furthermore, we added the results of the LES presented by Marakgos et al. [START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF], DesJardin et al. [START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF] and Ma et al. [START_REF] Ma | Verification and validation of a variable-density solver for fire safety applications[END_REF] for further comparisons.

Figure 5 shows the centerline mean and rms profiles of the streamwise velocity up to z = 0.8 m above the helium source. Both profiles predicted by the medium (∆x = 1 cm) and fine (∆x = 0.8 cm) grids are in close agreement and within the experimental uncertainties at all locations along the plume axis. The rms are more difficult to predict which is consistent with the findings of Maragkos et al. [START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF].

The mean and rms profiles of helium mass fraction along the centerline up to a height of z = 0.8 m are displayed in Fig. 6. The mean Y He The experimental data includes the respective uncertainties reported in the experimental study [START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. The numerical results of Maragkos et al. [START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF] are also added.

profile follows that predicted by Maragkos et al. [START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF] and decays in a slower manner as compared to the experiment beyond z = 0.3 m. Consistently with the results of Marakgos et al. [START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF], the rms Y ′ He profile is significantly over-predicted over the entire axis. It can be observed that the predictions obtained with both the medium (∆x = 1 cm) and fine (∆x = 0.8 cm) grids are very similar, confirming that converged solutions are obtained with the medium grid. Hereafter, only the results obtained with medium grid will be shown. The experimental data includes the respective uncertainties reported in the experimental study [START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. The numerical results of Maragkos et al. [START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF] are also added.

The radial profiles of mean streamwise velocity at several downstream locations (z = 0.2 m, 0.4 m and 0.6 m above the inlet) are presented in Fig. 7.

Our simulation shows a satisfactory agreement with the experimental data and falls within the experimental uncertainties. The results are also inline with the numerical results [START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF][START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF][START_REF] Ma | Verification and validation of a variable-density solver for fire safety applications[END_REF]. We should mention that all the radial profiles are azimuthally averaged all around the plume. The experimental data includes the respective uncertainties reported in the experimental study [START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. Numerical data are also presented [START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF][START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF][START_REF] Ma | Verification and validation of a variable-density solver for fire safety applications[END_REF].

Figure 8 shows RMS values of the streamwise velocity at different heights. The results reproduce in overall well the experimental data. Nevertheless, at z = 0.2 m, the RMS profile exhibits an M-shaped profile, whereas the experimental profile is flatter. A similar behavior is observed by Maragkos et al. [START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF]. On the other hand, Desjardin et al. [START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF] overpredicted the centerline RMS values at all heights. The experimental data includes the respective uncertainties reported in the experimental study [START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. Numerical data are also presented [START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF][START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF][START_REF] Ma | Verification and validation of a variable-density solver for fire safety applications[END_REF].

The cross-stream velocity is presented in Fig. 9. Our results show a good overall agreement with the experimental data. Maragkos et al. [START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF], DesJardin et al. [START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF] and Ma et al. [START_REF] Ma | Verification and validation of a variable-density solver for fire safety applications[END_REF] reported an overestimation on the left hand side of the plume as seen in Fig. 9. While keeping in mind that the experimental data are not perfectly symmetric, we have the same behavior at z = 0.2 m but the disparities in our results are lower. Afterwards farther downstream, the agreement becomes better compared to experimental data. The application of SGS model did not increase the cross-stream velocity as reported by Maragkos et al. [START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF].

A rise of the cross-stream velocity indicates an increased entrainment from the surrounding fluid which will result in a surge in the streamwise velocity following the conservation of mass. So accurate predictions of the cross-stream velocities is crucial, as in these type of flows, entertainment controls the mixing (as elaborated in a previous work [START_REF] Taha | Lattice-boltzmann modeling of buoyancy-driven turbulent flows[END_REF]), which is a fundamental parameter in pool fires where the combustion process is mainly controlled by mixing (i.e., diffusion flame) as will be seen in Sec. 4.

-0.5 0 The experimental data includes the respective uncertainties reported in the experimental study [START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. Numerical data are also presented [START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF][START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF][START_REF] Ma | Verification and validation of a variable-density solver for fire safety applications[END_REF].

Figure 10 shows the RMS values of the cross-stream velocity at different heights. The results agree with the experiment falling within the uncertainties except at some radial locations where they exhibit a slight overestimation. Our simulation shows an M-shaped for the RMS of the cross-stream velocities, coherent with the observations of Chung et al. [START_REF] Chung | Buoyancy-corrected k-ε models and large eddy simulation applied to a large axisymmetric helium plume[END_REF], DesJardin et al. [START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF] and Maragkos et al. [START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF]. This M-shaped transits progressively to a bell-shaped profile as the height increases.

The radial profiles of mean helium mass fraction at different heights are presented in Fig. 11. Slight overestimation can be observed, especially at z = 0.2 m at the left side of the plume. The discrepancies of the centerline value increase with height consistently with the overestimation in the axial profile seen in Fig. 6. The profiles align better with the experiment farther from the centerline. The LES of Maragkos et al. [START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF] is in better agreement in comparison with the experimental data. The experimental data includes the respective uncertainties reported in the experimental study [START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. Numerical data are also presented [START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF][START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF][START_REF] Ma | Verification and validation of a variable-density solver for fire safety applications[END_REF]]. The experimental data includes the respective uncertainties reported in the experimental study [START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. Numerical data are also presented [START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF][START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF][START_REF] Ma | Verification and validation of a variable-density solver for fire safety applications[END_REF].

Figure 12 shows the radial profiles of the RMS values of helium mass fraction. Globally, a significant over-prediction from the experimental data is observed close to the plume inlet for all the simulation. However, the discrepancies reduce farther downstream from the source. The results are consistent with those reported in previous numerical studies [START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF][START_REF] Maragkos | Application of fds and firefoam in large eddy simulations of a turbulent buoyant helium plume[END_REF][START_REF] Ma | Verification and validation of a variable-density solver for fire safety applications[END_REF].

Sandia's 1-m diameter methane fire

Numerical Setup

The computational domain is a 4×4×7 m 3 box, on which a non-uniform Cartesian grid was used with refinement zones for cost reduction. The largest and smallest mesh sizes in the computational domain are δx max = 4 cm and δx min = 2 cm, resulting in a 6-million cell mesh. The configuration of the 1 m diameter methane fire was simplified such that the fuel (methane) enters the domain through a 1 m diameter inlet surrounded by a 0.51 m wide steel plate (representing the ground plane). 21 The experimental data includes the respective uncertainties reported in the experimental study [START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. Numerical data are also presented [START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF][START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF][START_REF] Ma | Verification and validation of a variable-density solver for fire safety applications[END_REF].

Methane is injected at the ambient temperature and pressure, T inlet = 285 K and P inlet = 81.0 kPa, respectively. Uniform bulk inlet velocity for methane is U CH 4 inlet = 0.097 m/s corresponding to TEST-24 conditions (see Table 1). In addition, a co-flow air stream is injected outside the annular plate with a uniform velocity of U Air cof low = 0.14 m/s. The value of the coflow velocity was taken following the numerical studies of Han et al. [START_REF] Han | Les/pdf modelling of a one-meter diameter methane fire plume[END_REF] and Koo et al. [START_REF] Koo | Model sensitivities in les predictions of buoyant methane fire plumes[END_REF]. Boundary conditions of the sides and the outlet are identical to the ones used in the helium study in Sec. 3. Velocities with negative values were forced to zero at the outlet to prevent flow rentering the domain (reverse flow), which might trigger numerical instabilities and pollute the results. The steel plate was modeled as an adiabatic non-slip wall.

A single-step irreversible chemical reaction of methane is considered:

CH 4 + 2(O 2 + 3.76N 2 ) → CO 2 + 2H 2 O + 7.52N 2 (38)
The unresolved turbulent stress is modeled using the Vreman model (Sec. 2.2) with C s = 0.1 [START_REF] Maragkos | Large eddy simulations of ch4 fire plumes[END_REF]. In order to resolve the filtered heat and species fluxes, a value of 0.7 is assigned to both turbulent Prandtl (Pr t ) and Schmidt (Sc t ) numbers [START_REF] Maragkos | Large eddy simulations of ch4 fire plumes[END_REF]. The molecular Prandtl number (Pr) is set to 0.7, while the molecular Schmidt numbers of the species are summarized in Table 3. The molecular viscosity follows a power law because the simulation will exhibit high levels of temperature due to combustion, which impacts the viscosity significantly. The radiant fraction is 25% of the total heat release.

The simulation is run over 27 physical seconds. The first 7 seconds are thought sufficient to evacuate the initial solution then the remaining 20 seconds are used to construct the statistics for quantitative comparisons. The time step is local and depends on the local mesh size while keeping a 22 Figure 13 shows the energy spectrum based on the axial velocity at z = 0.5 m. The spectrum exhibits the theoretical energy cascade characteristic reported by the theory of Kolmogorov [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large reynolds numbers[END_REF] and that we have the correct power law of -5/3 in the inertial range, also the dissipation range was detected at higher frequencies.

Pool fire dynamics: Instability modes

Large-scale pool fires exhibit a particular behavior where a periodic shedding of large toroidal vortices pulsates at a specific frequency and is proportional to the square root of the diameter [START_REF] Cetegen | Experiments on the periodic instability of buoyant plumes and pool fires[END_REF]. The "puffing" mechanism is fundamental for fires as it controls the mixing between fuel and oxidizer, which conducts the combustion process. Figure 14 shows the different stages of a puffing cycle (four distinct stages); at first, instabilities are generated near the edges of the plume at the fire base due to baroclinic and gravitational torques; afterwards, those instabilities grow towards the center of the source; then, a large toroidal vortex is formed due to Rayleigh-Taylor instability which self-propagates and entrains a large amount of surrounding air; finally, the destruction of this toroidal vortex because of the creation of secondary instabilities and Kelvin-Helmholtz instabilities that grow causing a non-linear breakdown of the toroidal vortex. These stages are then repeated at every cycle, and the frequency of this motion was quantified experimentally, and numerical correlations were proposed.

Extensive experiments covering different scenarios and regimes of fires concluded that the puffing frequency depends only on the fire source diameter regardless of the fuel type. Correlations can be found in the literature that predicts the puffing frequency, such as the one given by Cetegen and Ahmed [START_REF] Cetegen | Experiments on the periodic instability of buoyant plumes and pool fires[END_REF]:

f = 1.5 √ D Hz, (39)
and the one given by Zukoski [START_REF] Zukoski | Properties of fire plumes[END_REF]:

f = 0.5 g D 1 2 Hz ( 40
)
Figure 15 shows the Fast Fourier transform (FFT) performed on the time signal of axial velocity at z = 0.5 m in order to identify the principal frequency of puffing. Our simulation predicts a puffing frequency of ∼ 1.32 Hz which underestimates the experimental values of 1.57 Hz, while the correlation of Cetegen gives 1.5 Hz and the correlation of Zukoski estimates 1.57 Hz.

Mean and RMS profiles

The axial profiles of axial velocity and temperature are displayed in Fig. 16. The velocity U z in the near-field region is in reasonable agreement with the experiment. On the other, comparisons with the McCaffrey correlations suggest that it is underestimated in the far-field region. The temperature decrease rate in the far field is consistent with the correlation of McCaffrey. However, the quantitative comparison suggests that the temperature is overestimated. Fig. 17 shows the radial profiles of axial velocity at different elevations downstream the inlet z = 0.3, 0.5, 0.9 m. A decent agreement can be observed with the experimental data, the predicted profiles remaining within the experimental uncertainties.

Fig. 18 depicts the radial profiles of horizontal velocity, which may be 25

TKE = 1 2 u ′2 x + u ′2 y + u ′2 z (41)
Fig. 19 shows the turbulent kinetic energy at different heights. The agreement is reasonable, although noticeable discrepancies are observed further downstream the fire base (i.e., at z = 0.9 m).

Flame height

A quantitative criterion to define the average flame height is to identify the centerline location where the difference between the mean gas temper-27 Comparison is done with experimental and numerical data [START_REF] Maragkos | Large eddy simulations of ch4 fire plumes[END_REF][START_REF] Tieszen | Experimental study of the effect of fuel mass flux on a 1-m-diameter methane fire and comparison with a hydrogen fire[END_REF] ature and the ambient reaches 500 -600 K [START_REF] Drysdale | An introduction to fire dynamics[END_REF]. In our study, a threshold of 550 K is selected following Maragkos and Merci [START_REF] Maragkos | Large eddy simulations of ch4 fire plumes[END_REF]. Our LES estimates a height of 5.2 meters which is slightly higher than the reported value of 4.8 m [START_REF] Marjanovic | On the evolution of the plume function and entrainment in the near-source region of lazy plumes[END_REF].

Concluding remarks

LES of fire plume scenarios representative of unwanted fire were performed for the first time by using a pressure-based LBM solver with low-Mach number approximation. Two scenarios were selected involving largescale helium and methane pool fires. The capability of the LBM-based LES to reproduce quantitatively the near-field dynamics of purely buoyant flows was clearly demonstrated, and numerical results were found similar to predictions reported with classical Navier-Stokes solvers at a lower computational cost. These results establish that LBM solvers are good candidates to serve as a basis for developing fire simulators. The next natural steps of this work will be to incorporate state-of-the-art turbulence and radiation sub-models to model more complex configurations. stance, the equilibrium distribution in equation ( 28) reads f eq,19r i = ω i a (0),eq + H It is worth noting that the the forth-order coefficients (a (4),eq ) are added to improve the isotropicity of the lattice, which could be quite important considering the round jet simulation in the current study.

The third-order off-equilibrium terms are reconstructed recursively from the second-order non-equilibrium tensor as Depending on the order of the Gauss-Hermite quadrature [START_REF] Krüger | The lattice boltzmann method[END_REF] used in the LB model, an adequate forcing term should be added to achieve a correct viscous stress tensor: where all the differential operations are performed using first order upwind FD except for the divergence operator for which a second order centered FD scheme was employed. The final expression of the forcing term is then

F E i = ω i 2c 4 s H (2)
i,αβ a F E αβ + F g i , (A.30)

Figure 1 :

 1 Figure 1: Schematic of the computational domain including refinement zones and boundary conditions.

Figure 2 :

 2 Figure 2: Temporal energy spectrum at z = 0.5 m for axial velocity. Dashed line indicate the expected characteristic slope of -5/3 [63].

Figure 3 :

 3 Figure 3: Instantaneous iso-volume of density colored by gravitational torque magnitude showing a complete puffing cycle. From left to right 4 different instances were captured which extend over around 0.7 sec which is the periodic time of a puffing cycle.

Figure 4 :

 4 Figure 4: (left) Time series of centerline streamwise velocity at z/D = 0.5 in a window of 5 s. (right) Power spectrum of the centerline streamwise velocity U z . The principle frequency (i.e. puffing frequency) is indicated on the figure.The sampling time ts to construct the power spectrum is taken to be each 10 time steps that is ts coarse = 9.6 × 10 -4 sec, ts medium = 4.8 × 10 -4 sec and ts f ine = 3.8 × 10 -4 sec for the coarse, medium and fine mesh, respectively.

Figure 5 :

 5 Figure 5: Comparison of mean (left) and rms (right) of centerline streamwise velocity.The experimental data includes the respective uncertainties reported in the experimental study[START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. The numerical results of Maragkos et al.[START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF] are also added.

Figure 6 :

 6 Figure 6: Comparison of mean (left) and rms (right) of centerline helium mass fraction.The experimental data includes the respective uncertainties reported in the experimental study[START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. The numerical results of Maragkos et al.[START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF] are also added.

Figure 7 :

 7 Figure 7: Comparison of mean radial profiles of streamwise velocity for different resolutions and different turbulence models at heights z = 0.2 m (left), z = 0.4 m (middle) and z = 0.6 m (right).The experimental data includes the respective uncertainties reported in the experimental study[START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. Numerical data are also presented[START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF][START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF][START_REF] Ma | Verification and validation of a variable-density solver for fire safety applications[END_REF].

Figure 8 :

 8 Figure 8: Comparison of rms radial profiles of streamwise velocity for different resolutions and different turbulence models at heights z = 0.2 m (left), z = 0.4 m (middle) and z = 0.6 m (right).The experimental data includes the respective uncertainties reported in the experimental study[START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. Numerical data are also presented[START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF][START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF][START_REF] Ma | Verification and validation of a variable-density solver for fire safety applications[END_REF].

Figure 9 :

 9 Figure 9: Comparison of mean radial profiles of cross-stream velocity for different resolutions and different turbulence models at heights z = 0.2 m (left), z = 0.4 m (middle) and z = 0.6 m (right).The experimental data includes the respective uncertainties reported in the experimental study[START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. Numerical data are also presented[START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF][START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF][START_REF] Ma | Verification and validation of a variable-density solver for fire safety applications[END_REF].

20 -Figure 10 :

 2010 Figure 10: Comparison of rms radial profiles of cross-stream velocity for different resolutions and different turbulence models at heights z = 0.2 m (left), z = 0.4 m (middle) and z = 0.6 m (right).The experimental data includes the respective uncertainties reported in the experimental study[START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. Numerical data are also presented[START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF][START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF][START_REF] Ma | Verification and validation of a variable-density solver for fire safety applications[END_REF].

Figure 11 :

 11 Figure 11: Comparison of mean radial profiles of helium's mass fraction for different resolutions and different turbulence models at heights z = 0.2 m (left), z = 0.4 m (middle) and z = 0.6 m (right).The experimental data includes the respective uncertainties reported in the experimental study[START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. Numerical data are also presented[START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF][START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF][START_REF] Ma | Verification and validation of a variable-density solver for fire safety applications[END_REF].

Figure 12 :

 12 Figure 12: Comparison of rms radial profiles of helium's mass fraction for different resolutions and different turbulence models at heights z = 0.2 m (left), z = 0.4 m (middle) and z = 0.6 m (right).The experimental data includes the respective uncertainties reported in the experimental study[START_REF] O'hern | Experimental study of a turbulent buoyant helium plume[END_REF]. Numerical data are also presented[START_REF] Desjardin | Large eddy simulation and experimental measurements of the near-field of a large turbulent helium plume[END_REF][START_REF] Maragkos | Large eddy simulations of the flow in the near-field region of a turbulent buoyant helium plume[END_REF][START_REF] Ma | Verification and validation of a variable-density solver for fire safety applications[END_REF].

Figure 13 :

 13 Figure 13: Temporal energy spectrum at z = 0.5 m for axial velocity. The dashed line indicates the expected characteristic slope of -5/3 [63].

Figure 14 :

 14 Figure 14: Snapshots of heat release rate iso-volume throughout 1 second depicting a complete puffing cycle.

Figure 15 :Figure 16 :Figure 17 :Figure 18 :

 15161718 Figure 15: Fast Fourier transform (FFT) based on the axial velocity at z = 0.5 m showing the puffing frequency.

Figure 19 :

 19 Figure19: Radial profiles of turbulent kinetic energy (TKE) at z = 0.3 m (left), z = 0.5 m (middle) and z = 0.9 m (right). Comparison is done with experimental and numerical data[START_REF] Maragkos | Large eddy simulations of ch4 fire plumes[END_REF][START_REF] Tieszen | Experimental study of the effect of fuel mass flux on a 1-m-diameter methane fire and comparison with a hydrogen fire[END_REF]

1 = 3 ( 2 = 3 ( 3 = 3 ( 4 = 5 = 6 =

 132333456 ρu x u x u y + ρu y u z u z ) , ρu x u z u z + ρu x u y u y ) , (A.17) a (3r),eq ρu y u y u z + ρu x u x u z ) , (A.18) a (3r),eq ρu x u x u y -ρu y u z u z , (A.[START_REF] Nmira | Influence of gas radiative property models on large eddy simulation of 1 m methanol pool fires[END_REF]) a (3r),eq ρu x u z u z -ρu x u y u y , ρu y u y u z -ρu x u x u z ,

  αβ is a correction tensor due to the deflection of second order moments of the population introduced by the modification of the mass equation, which can be evaluated asa cor αβ ≡ c 2 s δ αβ ∂(ρ(1 -θ)) ∂t , (A.28)which can be discretized using a backward Euler operator and a F D αβ the correction tensor due to the defect of the lattice at third ordera F D αβ = -  (ρu 3 x ) ,x (ρu x u y u z ) ,z (ρu x u y u z ) ,y (ρu x u y u z ) ,z(ρu 3 y ) ,y (ρu x u y u z ) ,x (ρu x u y u z ) ,y (ρu x u y u z ) ,x (ρu 3 z )

Table 1 :

 1 Initial and boundary conditions of TEST-24 in the experiments[START_REF] Tieszen | Experimental study of the flow field in and around a one meter diameter methane fire[END_REF][START_REF] Tieszen | Experimental study of the effect of fuel mass flux on a 1-m-diameter methane fire and comparison with a hydrogen fire[END_REF]

	]:

Table 2 :

 2 Cost comparison between different numerical studies for the finest mesh.

Table 3 :

 3 Molecular Schmidt number Sc for different species constant CFL number (CFL ≈ 0.7). This gives a minimum time step in the finest zone δt min ≃ 2.0 × 10 -4 sec. The simulation is parallelized on 128 cores on a Dell PowerEdge C6420 server with 4× 32-core Intel Xeon Gold 6142, 2.6 GHz, and 96 GB RAM producing a calculation time of ∼ 1050 CPU hours.

	4.2. Results and discussion
	4.2.1. Turbulence: Energy Cascade

Acknowledgments

Centre de Calcul Intensif d'Aix-Marseille and GENCI-TGCC/CINES (Grant 2022-A0132B11951) are acknowledged for granting access to their high-performance computing resources.

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Appendix A. Expressions for the LBM solver

In regularized LBM, distribution functions will be constructed using an orthogonal polynomial basis. The basis of the D3Q19r lattice used in the current study consists of 19 polynomials, read [START_REF] Farag | A unified hybrid lattice-boltzmann method for compressible flows: bridging between pressure-based and density-based methods[END_REF]

i,yzz , (A.2)

i,xyy , (A.3)

i,xxz , (A.4)

i,yzz , (A.5)

i,xyy , (A.6)

i,yyz -H

i,xxz , (A.7)

i,zyx , (A.9)

where

i,γγ , (A.12)

Any distribution function in the D3Q19r lattice can be written as a (weighted) sum of the contributions from each base polynomial. For in-where F g i is the gravity force term defined as

where g α is the gravity acceleration in the direction α.

Appendix B. Vorticity equation

Helmholtz equation is the transport equation for vorticity, which can be deduced as follows. Starting from writing the momentum equation of incompressible flow, that is ∇ • u = 0, but with a variable density ρ.

The convective term u • ∇u can be written in the following form:

The following vector identity will be used:

where ∇ • ω ≡ 0 (vorticity is solenoidal by definition), the vorticity equation can be rearranged as:

4) The pressure term can be written as:

Further more for a homentropic flow, the pressure is a function of only of density resulting in the alignment of the two vector ∇p and ∇ρ and consequently the so called baroclinic torque -(∇ρ × ∇p) /ρ 2 will vanish. By the same manner, the buoyancy term will be treated : 38


Ibrahim Hatoum

Nassim Bouchoul

Mélissandre Richard

Christophe Dujardin

Revealing origin of hydrogen-carbonate species in CO oxidation over Pt/Al2O3: a SSITKA-IR study

Keywords: CO oxidation, SSITKA, Operando IR spectroscopy, Hydrogen-carbonate, Mechanism

The catalytic oxidation of CO by O2 to form CO2 over Pt based catalysis is one of the most studied catalytic reaction leading to controversial mechanism descriptions. Steady State Isotopic Transient Kinetic Analysis (SSITKA) coupled with InfraRed (IR) spectroscopy is a powerful technique to study heterogeneous reaction mechanisms combining both the observation of adsorbed species on the catalyst surface and kinetic measurements. In this paper, SSITKA-IR technique was applied to study the CO oxidation reaction in presence and absence of CO2 in order to distinguish between active and spectator intermediates formed during the reaction. Linear carbonyl and bridged carbonyl type species adsorbed on metallic Pt0 were clearly identified as active intermediates at low temperature (131 °C). In contrast, the hydrogen-carbonate species formed during CO oxidation reaction were proven to be inactive species in CO2 formation but rather due to the re-adsorption of CO2 product itself on alumina surface.

Introduction

Air pollution, along with climate change, is one of the most serious environmental threats to human health responsible for 7 million premature deaths per year and loss of millions healthy life years [START_REF]New WHO global air quality guidelines aim to save millions of lives from air pollution[END_REF].

Transportation industry is the largest contributor to air pollution with emissions of carbon monoxide (CO), nitrogen oxides (NOx), unburned hydrocarbons (HC) and particulate matter [START_REF] Casapu | Origin of the normal and inverse hysteresis behavior during CO oxidation over Pt/Al2O3[END_REF]. In Europe, legislation imposing restrictions on the emissions of pollutants from vehicle exhaust has been introduced since 90's and the next coming Euro 7 regulations are expected to reduce drastically the acceptable emissions thresholds. Thus, one of the greatest challenges for current catalytic converter technologies will be to meet the major issue of cold-start (warm up) emissions, when catalyst has not yet reached its light-off temperature [START_REF] Gao | Review of thermal management of catalytic converters to decrease engine emissions during[END_REF]. In that context, the catalytic CO oxidation at low temperature has gained increasing attention [START_REF] Carlsson | Low-temperature oxidation of carbon monoxide and methane over alumina and ceria supported platinum catalysts[END_REF][START_REF] Heck | Automobile exhaust catalysts[END_REF][START_REF] Shelef | Twenty-five years after introduction of automotive catalysts: what next?[END_REF]. This reaction is one of the most extensively studied in the history of heterogeneous catalysis. Despite the fact that CO oxidation appears to be a relatively straightforward reaction, the fundamental understanding of the process over supported noble metal particles is still a matter of debate. Over 20 different mechanistic steps have been proposed [START_REF] Soubaihi | Critical review of low-temperature CO oxidation and hysteresis phenomenon on heterogeneous catalysts[END_REF]. The Langmuir-Hinshelwood (LH) dual-site mechanism, in which the reaction occurs between CO and O2 after both molecules have been adsorbed on the surface of a catalyst, is one of the most widely accepted pathway for low-temperature CO oxidation [8][START_REF] Salomons | CO and H2 oxidation on a platinum monolith diesel oxidation catalyst[END_REF][START_REF] Carlsson | A transient in situ FTIR and XANES study of CO oxidation over Pt/ Al2O3 catalysts[END_REF]. Recently, Newton et al. proposed a carbonatemediated mechanism for CO oxidation at low temperature over a highly dispersed γ-alumina-supported Pt catalyst [START_REF] Newton | Room-temperature carbon monoxide oxidation by oxygen over Pt/Al2O3 mediated by reactive platinum carbonates[END_REF][START_REF] Newton | Time resolved operando X-ray techniques in catalysis, a case study: CO oxidation by O2 over Pt surfaces and alumina supported Pt catalysts[END_REF][START_REF] Newton | Kinetic studies of the Pt carbonate-mediated, room-temperature oxidation of carbon monoxide by oxygen over Pt/Al2O3 using combined, timeresolved XAFS, DRIFTS, and mass spectrometry[END_REF]. However, the participation of carbonate and hydrogen-carbonate species is not explicitly demonstrated in the literature. The Steady State Isotopic Transient Kinetic Analysis (SSITKA) coupled with operando InfraRed spectroscopy (IR) is one of the most powerful system to investigate catalytic reaction mechanisms. The combination of steady state and isotopic transient conditions highlights the advantages of both techniques as obtaining kinetic information and identifying the structural properties of the surface-adsorbed active or inactive species at molecular level under the realistic steadystate conditions of the reaction [START_REF] Burch | A steady-state isotopic transient kinetic analysis of the NO/O2/H2 reaction over Pt/SiO2 catalysts[END_REF][START_REF] Ledesma | Recent approaches in mechanistic and kinetic studies of catalytic reactions using SSITKA Technique[END_REF]. In this work, combined SSITKA-transmission IR experiments will be applied for the first time to study the oxidation of CO on Pt/Al2O3 catalyst unraveling the true role of hydrogen-carbonate species during the reaction at low temperature.

Experimental

Catalyst Preparation and Characterization

The catalyst was prepared by incipient wetness impregnation of ɣ-Al2O3 support (supplied by the Centre de Recherches de SOLAIZE TOTALEnergies) with platinum (II) acetylacetonate (Fluka) precursor salt to obtain 1 wt% metal loading. The sample was dried overnight at 85 °C, calcined 2 h under air flow at 500 °C (ramp 10 °C/min) and reduced with dihydrogen 2.5 h at 400 °C (ramp 5 °C/min). Final Pt loading was determined by ICP-OES 5110 (Agilent Technologies). Nitrogen physisorption measurements were performed on a Micromeritics TriStar II instrument. The experiments were conducted using 0.2 g of powder sample, previously outgazed at 150 °C under vacuum for 2 h. Brunauer-Emmett-Teller (BET) theory was applied to calculate the specific surface area. Pore size distributions were obtained from the complete isotherms based on the Barret-Joyner-Halenda (BJH) theory. H2 chemisorption measurement was carried out at 30 °C in a Micromeritics Autochem II instrument using a pulse technique. Prior to the adsorption, the catalyst (50 mg) was preheated in pure H2 at 400 °C and outgassed in flowing Argon at 400 °C. Then the temperature of the catalytic bed was reduced to 30 °C in order to perform the measurement in which H2 pulses (5 vol% H2 in Ar) were injected until saturation. Metal dispersion was calculated assuming H/Pt = 1.

SSITKA-IR Experiments

The Pt/Al2O3 catalyst was finely crushed and pressed into a 5.9 mg square pellet of 1.21 cm2 under a pressure of 100 MPa for operando SSITKA-IR measurements. The pellet was loaded in the homemade IR cell with an internal dead volume inferior to 0.5 cm3 and equipped with KBr windows. The cell was connected upstream to an automatic regulator for gas flow and switch valves controls. Prior to the SSITKA-IR study, the catalytic material was pre-treated at 300 °C during 5 h under He (35 mL/min).

Then, 2060 ppm CO/10% O2/ He gas flow mixture was sent to the catalyst at a total flow rate of 12 mL/min (122 L h-1 g-1). The temperature was increased from 50 °C up to 131 °C (1 °C/min) corresponding to a CO conversion of 11%. After 30 min of isotherm equilibration, the SSITKA-IR experiment was performed by switching the reactant from 12CO to its labelled counterpart 13CO at the same concentration, pressure and temperature. This 12CO → 13CO gas switch step was repeated three times to record the concentration evolution of reversibly adsorbed CO and active Cx intermediates but also to check reproducibility (true steady-state achieved). 5000 ppm of Kr and 7900 ppm of CH4 were used as inert gas for kinetic analysis (i.e. SSITKA gas switching procedure: 12CO/ O2/He → 13CO/O2/He/Kr/CH4). The comparison of Kr and CH4 mass spectrometry (MS) signals is presented in supplementary information and confirms the absence of time delay between the two gas responses (Fig. S1). The use of methane allows direct comparison of MS (gas phase) and IR (adsorbed species) signals.

Finally, the catalyst was heated up to 300 °C (1 °C/min) before cooling it down to 50 °C. The SSITKA-IR experiment was subsequently repeated at the same CO conversion (11% -136 °C). The same experimental protocol was carried out for the test in the presence of CO2 by adding 0.2 vol% CO2 in both gas feeds: 12CO/12CO2/O2/He → 13CO/12CO2/O2/ He/Kr/CH4. The gas flow composition from IR cell inlet/outlet was analyzed by micro-gas chromatography (Varian 490) and a quadrupole mass The evolution of the adsorbed species on the catalytic material was analyzed by a Thermo Scientific Nicolet 6700 FTIR spectrometer (rapid-scan mode, 1 scan per spectra and resolution 2 cm-1) equipped with a MCT detector. CO conversion (χCO) was calculated using the following Eq. (1):

𝜒 (%) = [ ] [ ] × 100 (1)
The concentration of reversibly adsorbed CO (NCO), active C-containing intermediates (NC) leading to the formation of CO2, and their respective surface coverage (θ) were calculated using the following Eqs.

(2)-( 4) based on material balances:

𝑁 (µ𝑚𝑜𝑙 𝑔 ) = ( ) ∫ [𝑍 (𝑡) . -𝑍 (𝑡)] 𝑑𝑡 (2)
𝑁 (µ𝑚𝑜𝑙 𝑔 ) = ∫ [𝑍 (𝑡) . -𝑍 (𝑡)] 𝑑𝑡 (3)
𝜃 = 𝑁 𝑁 , (𝑖 = 𝐶, 𝐶𝑂) (4)
Where FT is the total molar flow rate (μmols-1) of the gas feed, Yf CO is the mole fraction of CO in the feed, χCO is the CO conversion, W is the mass of catalyst (g), Y13CO2 is the mole fraction of 13CO2 at the outlet of the reactor in the new steady-state obtained, Z is the dimensionless concentration of a given gas-phase species, ts. s is the time at which the new steady-state is obtained under the 13CO gas mixture, and NPt, surf is the molar amount of surface platinum per gram of catalyst (μmol Pt gcat -1) obtained based on the dispersion and loading values of Pt metal. The surface residence times (τ) of surface adsorbed species were calculated using Eqs. ( 5) and ( 6):

𝜏 (𝑠) = ∫ [𝑍 (𝑡) . -𝑍 (𝑡)] 𝑑𝑡 (5)
𝜏 (𝑠) = ∫ [𝑍 (𝑡) . -𝑍 (𝑡)] 𝑑𝑡 (6)
The turnover frequencies, namely TOFCO, ITK and TOFCO2, ITK (s-1), given by Eqs. ( 7) and ( 8) are based on the concentrations of active reaction intermediates found in the reaction path of CO2 formation, NCO and NC, respectively.

𝑇𝑂𝐹 , = (7)
𝑇𝑂𝐹 , = (8)
where rCO and rCO2 are the steady-state rates of CO conversion and CO2 formation in μmol g-1 s-1, respectively.

3 Results and Discussions

Initial Catalyst Characterization

The catalyst exhibits a specific surface area of 232 m2g-1 with a pore volume of 0.79 cm3g-1 comparable to the area and pore volume obtained for the ɣ-Al2O3 support itself, 225 m2g-1 and 0.76 cm3g-1, respectively. A loading of 0.82 wt% of Pt on alumina was determined using inductive coupled plasma optical emission spectroscopy (ICP-OES) analysis. Finally, H2 chemisorption resulted in a dispersion of 68% corresponding to a mean Pt particle size of about 1.65 nm (Table S1 in supplementary material). or more adsorbed reaction intermediates after that of 13CO before it forms 13CO2. Based on the SSITKA theory (Eqs. 5 and 6), the integrations of observed delays between the response curves of CH4 and 13CO, then 13CO and 13CO2 provide the mean surface residence times τ of reversibly adsorbed CO (τCO = 5.0 s) and active intermediates (τC = 29.2 s), respectively. The latter values allow calculating the surface concentrations and coverages of each species using the Eq. 2 to 4. A concentration of 12.9 μmol of CO surface-adsorbed species (NCO) per gram of catalyst is reached corresponding to a Pt surface coverage of about 45% (θCO = 0.45 in Table 1). In addition, a concentration of 9.6 μmolg-1 of Cx-adsorbed intermediate is determined corresponding to a surface coverage of 33% (θC = 0.33 in Table 1). The black line in Fig. 2a, corresponding to the spectrum recorded in presence of 12CO and O2 in the carbonyl wavenumber region (2200 -1700 cm-1), shows a strong IR band at 2058 cm-1. It corresponds to the stretching C ≡ O vibration of the 12CO species linearly adsorbed on metallic Pt atom [START_REF] Bourane | Heats of adsorption of the linear CO species on a Pt/Al2O3 catalyst using FTIR spectroscopy: comparison between TPD and adsorption equilibrium procedures[END_REF][START_REF] Bourane | Oxidation of CO on a Pt/Al2O3 catalyst: from the surface elementary steps to light-off tests: IV. Kinetic study of the reduction by CO of strongly adsorbed oxygen species[END_REF]. A red shift of this band upon 13CO exchange is observed, giving rise to a band at 2010 cm-1 corresponding to the 13 C-labelled linear carbonyls adsorbed on Pt0. Carbonyls adsorbed on oxidized Ptδ+ are also observed on both spectra in Fig. 2a as a smaller IR bands around 2102 cm-1 [START_REF] Bourane | Impact of Pt dispersion on the elementary steps of CO oxidation by O2 over Pt/ Al2O3 catalysts[END_REF]. No isotopic shift is shown in that case. The presence of bridged 12CO carbonyl species adsorbed on two metallic platinum is also observed at 1837 cm-1 along with an isotopic 13CO shift after the transient exchange at 1797 cm-1. The carbonate region represented in Fig. 2b (1800 -1300 cm-1) can be decomposed with the bands at 1434 and 1657 cm-1 corresponding to the hydrogencarbonate species ν(OCO) asymmetric and ν(OCO) symmetric stretching vibrations, respectively [START_REF] Ojala | Isotopic oxygen exchange over Pd/Al2O3 catalyst: study on C18O2 and 18O2 exchange[END_REF]. A band near 1215 cm-1 corresponding to the bending δ(OH) vibrations as well as the stretching ν(OH) band at higher wavenumber would be expected for hydrogen-carbonates species. However, due to the low energy in these ranges, we were not The TOF,ITK were also calculated based on the concentration of reversibly adsorbed CO ( NCO) and the sum of the concentration NCO+ NC for TOFCO and TOFCO2, respectively (Eqs. 7 and 8, Table 1). It is worth noting that this calculation takes into account the active sites defined by SSITKA measurements (concentration of all the active species that truly participate in the formation of a given reaction product)

SSITKA-IR Studies for CO Oxidation

and not the amount of Pt exposed surface sites determined by H2-chemisorption as conventionally adopted for TOF determination [START_REF] Vasiliades | The effect of preparation route of commercial Co/γ-Al2O3 catalyst on important Fischer-Tropsch kinetic parameters studied by SSITKA and CO-DRIFTS transient hydrogenation techniques[END_REF]. As a direct result of the decrease of Cx-adsorbed intermediate 13 CO 2 13 CO CH 4 12 CO concentration (NC), the TOFCO2, ITK increases by about 45% in the presence of CO2 (from 0.015 s-1 to 0.027 s-1) while the values of TOFCO, ITK remain similar in both cases. IR spectra analysis has also been performed to observe the evolution of adsorbed species on the catalytic surface after reactant exchange with 13CO in the presence of CO2. Figure 4 shows the IR spectra recorded during the 12CO+O2+12CO2 and 13CO+O2+12CO2 flow stages at the steady state (1 min after isotopic reactant switching in case of 13CO) at 136 °C.

0.00E+00 The Fig. 4a shows the shift of the band at 2058 cm-1 in the carbonyl wavenumber region corresponding to 12CO linear species adsorbed on metallic Pt towards the band of its labelled counterpart 13CO

adsorbed linear species at 2010 cm-1 as seen previously in Fig. 2. Once again, CO-Ptδ+ carbonyls are observed at 2102 cm-1 in both feeds without any isotopic shift effect while bridged carbonyls ( Pt2CO) are shifted from 1837 cm-1 to 1797 cm-1 regarding to 12CO→13CO reactant switch. Interestingly, in the presence of CO2, the bands at 1434 cm-1 and 1657 cm-1 observed in the carbonate region (Fig. 4b) and attributed to O = C = O asymmetric and symmetric stretching vibrations of hydrogen-carbonate species are not shifted after reactant exchange to 13CO. Other low-intensity non-shifted bands appear at 1594 cm-1, 1395 cm-1 and 1372 cm-1 which can be assigned to OCO asymmetric stretching, C-H bending and OCO symmetric stretching vibrations of formate species formed on alumina surface, respectively [START_REF] Lorito | Determination of formate decomposition rates and relation to product formation during CO hydrogenation over supported cobalt[END_REF][START_REF] Mcinroy | The application of diffuse reflectance infrared spectroscopy and temperatureprogrammed desorption to investigate the interaction of methanol on η-Alumina[END_REF].

General Discussion

The interest for SSITKA-IR methodology is illustrated in this study for the discrimination of active and spectator species in a catalytic reaction. This combines both the observation of adsorbed species and the kinetic measurements during CO oxidation reaction at low conversion level (i.e. in chemical regime).

The CO adsorption on Pt/Al2O3 mainly forms linear carbonyl type species adsorbed on metallic Pt0 as observed at 2058 cm-1 (12CO) and 2010 cm-1 (13CO) with operando IR spectroscopy and in lower extent bridged carbonyl species on metallic Pt0 (Fig. 2a). Moreover, after 1 min of the isotopic exchange, a complete vanish of the bands at 2058 cm-1 and 1837 cm-1 is observed which can be related to a high exchange rate and so high reactivity of CO adsorbed on Pt0 during the experiment. This kinetic result is furthermore coherent with the gas phase evolution of 13CO observed by mass spectrometry (Fig. 1) that reaches its maximum at around 60 s. In the meantime, no isotopic exchange during the 12CO→13CO SSITKA experiment is observed for carbonyls adsorbed on Ptδ+ suggesting a low reactivity and so their spectator nature in CO oxidation reaction at low temperature. An important lag was observed between the normalized 13CO and 13CO2 concentration curves indicating in theory the formation and participation of C-adsorbed reactive intermediates in the CO2 production. According to the IR results (Fig. 2b) hydrogen-carbonate species seems to be the most serious contender to play the role of active intermediates in the CO oxidation reaction due to a high 12 C→13 C exchange rate at around 130 °C comparable to the rate of CO2 product formation. The SSITKA-IR experiment in the presence of CO2

in the feed gives additional information concerning the role of hydrogen-carbonate during CO oxidation.

As a matter of fact, the simultaneous saturation of the catalyst surface with adsorbed CO2 reduces by a factor of seven the amount of C-containing intermediate species ( NC) for the same CO conversion and reaction temperature when CO2 is present in the reaction mixture. This result proves that the lag firstly observed between 13CO and 13CO2 during the SSITKA experiment without CO2 is due in reality to the re-adsorption of the CO2, produced by the CO oxidation, on the surface of the catalyst rather than any intermediate formation or reactivity to form CO2. IR spectra recorded in the presence of CO2 also confirm this trend in which IR bands of hydrogen-carbonate species don't show any 12C→13C isotopic shift anymore. This observation clearly indicates that their formation is induced in that case by the adsorption of 12CO2 already present in the gas feed rather than the conjugation of the reactants 13CO+O2 with surface hydroxyls from alumina. The intermediate formation of hydrogen-carbonate during CO oxidation can therefore be definitely ruled out.

Conclusion

The formation of carbonyl and hydrogen-carbonate species was investigated on 0.82wt% Pt/Al2O3 Supplementary Information : https:// doi. org/ 10. 1007/ s11244-022-01722-2.

  spectrometer (Pfeiffer QMS 200, 70 eV electronic impact) by recording He+ (m/z = 4), CH4 + (m/z = 15), 12CO+ (m/z = 28), 13CO+ (m/z = 29), 12CO2 + (m/z = 44), 13CO2 + (m/z = 45) and Kr+ (m/z = 84).

3. 2 . 1

 21 CO Oxidation Without CO2 in the Reaction MixtureFirstly, the oxidation of CO in an excess of O2 has been conventionally carried out without the presence of CO2 in the reaction feed. The normalized MS signals of 13CO, 13CO2 and CH4 are presented in Fig.1as a function of time after applying a reactant isotopic switching from 12CO to 13CO at low conversion (χCO = 11% at 131 °C). It is worth noting that the 13CO transient response curve lags behind that of the tracer (CH4) due to the formation of reversibly adsorbed CO during oxidation reaction. The 13CO2 transient response curve lags behind that of 13CO indicating that, in theory, 13CO2 passes through one

Fig. 1 ̶Fig. 2 ̶

 12 Fig. 1 ̶ Normalized concentrations of CH4, 13 CO and 13 CO2 obtained during the SSITKA transient ( 12 CO/O2/He 13 CO/O2/Kr/CH4/He) on Pt/Al2O3 at 131 °C (χCO = 11%)

  during CO oxidation in lean conditions. The present work has provided fundamental insight into the role of hydrogen-carbonate species in CO oxidation reaction. SSITKA-IR technique has been applied for the first time to study this reaction. A fast shift of IR signals related to linear carbonyl and bridged carbonyl type species adsorbed on metallic Pt0 was observed during the 12CO+O2 and 13CO+O2 transient exchange in absence and presence of CO2 in the reaction mixture, confirming their high reactivity toward CO2 formation contrary to carbonyls adsorbed on oxidized Ptδ+. The hydrogencarbonate species formed during CO oxidation reaction were finally proven to be inactive adsorbed species. The lag between the 13CO and 13CO2 observed on SSITKA normalized concentration curves was in fact explained by the readsorption of CO2 product, formed from the CO oxidation, on the surface of the catalyst. This result was supported by the non-shifted IR bands of the hydrogen-carbonate species under 12CO→13CO transient exchange in presence of CO2 in the gas feed.

  able to distinguish them. As for linear carbonyl species, a red isotopic shift is clearly spotted over the hydrogen-carbonates O = C = O bands. Upon exchange with 13CO, the bands at 1434 and 1657 cm-1 shifts towards the bands at 1395 and 1613 cm-1, respectively. 11% at 136 °C). In the presence of CO2, the mean surface residence times of reversibly adsorbed CO (τCO) and active intermediates (τC) appear to be equal to 5.2 s and 3.6 s, respectively. While no significant effect of CO2 is observed on τCO, the area between the 13CO and 13CO2 curves vanished quasi-completely leading to a drastic decrease of the τC value (Table1). Consequently, new calculations of the surface concentration and coverage of CO on catalytic surface from MS signals appear to be quite similar to the value determined without the presence of CO2, i.e. NCO = 12.8 μmolg-1 and θCO = 0.45, unlike those obtained for Cx-adsorbed intermediate with NC and θC dropping from 9.2 μmolg-1 and 0.33 to 1.4 μmolg-1 and 0.05 in presence of CO2, respectively.

			a)		2058	2010		CO + O 2
								13 CO + O 2
	A=0.02 3.2.2 CO Oxidation in Presence of CO2
	The same experiment of CO oxidation with O2 on Pt/Al2O3 catalyst has been repeated by adding, in that
	case,	a	significant	2102 amount	of	12CO2	in	both	SSITKA	feeds	(12CO/
	1837 1797 O2/12CO2/He➙13CO/O2/12CO2/Kr/CH4/He). The normalized IR signals of 13CO, 13CO2 and CH4 are
	represented as a function of time in Fig. 3 after applying the reactant isotopic switch at low conversion
	1300 1700 12 CO + O 2 1400 1800 13 CO + O 2 1395 1434 Wavenumber (cm -1 ) 1500 1600 1900 2000 Wavenumber (cm -1 ) 1657 1700 2100 1613 A=0.02 1800 b) (χCO = Fig. 3 ̶ Normalized concentrations of CH4, 13 CO and 13 CO2 obtained after the SSITKA transient in the
	presence of CO2 in the feed ( 12 CO/O2/ 12 CO2/He 13 CO/O2// 12 CO2/Kr/CH4/He) on Pt/Al2O3 at 136 °C (χCO =
	11%)

Table 1 -

 1 Mean residence time τ, concentration N, surface coverage θ, and TOF of adsorbed CO reactant and active intermediates in CO2 formation during SSITKA CO oxidation experiment with and without the presence of CO2 in the feed on Pt/Al2O3

	COads	τCO (s)	NCO (µmol/g)	θCO	TOFCO,ITK (s -1 )
	Without CO2	5.0	12.9	0.45	0.20
	With CO2	5.2	12.8	0.45	0.19
	Cx intermediates	τC (s)	NC (µmol/g)	θC	TOFCO2,ITK (s -1 )
	Without CO2	29.2	9.6	0.33	0.015
	With CO2	3.6	1.4	0.05	0.027

  Fig.4̶ IR spectra recorded during CO oxidation at 136 °C in the presence of12 CO/O2/ 12 CO2/He (black) and13 CO/O2/ 12 CO2/Kr/CH4/He (red): (a) carbonyl region and (b) carbonate region

	a)	2058			1.75E-01 12 CO + 12 CO 2 + O 2
			2010			1.50E-01 13 CO + 12 CO 2 + O 2
							1.25E-01
	A=0.025					1.00E-01
							7.50E-02
		2102		1837 1797	2.50E-02
							0.00E+00
	2200	2100	2000	1900		1800	1700
		Wavenumber (cm -1 )
	b)						1.75E-01
							1.50E-01
						1372
				1434	1395	1.25E-01
	A=0.025					1.00E-01
		1657	1594				7.50E-02
						5.00E-02 12 CO + 12 CO 2 + O 2
						2.50E-02 13 CO + 12 CO 2 + O 2
	1800	1700	1600	1500	1400	1300
		Wavenumber (cm -1 )

CO

Acknowledgements The authors thank the Region "Hauts-de-France" (Grant STaRS 2019), Centrale Lille Institute, the Ministere de l'Enseignement Superieur et de la Recherche (CPER IRENE and CPER ECRIN) and the European Fund for Regional Economic Development for their financial support. ICP-OES analyses were performed in the « Spectrometrie par torche a plasma » platform of the Research Federation Michel-Eugene Chevreul hosted by the LASIRE laboratory.

Author Contributions All authors contributed to the study conception and design. Material preparation and data collection were performed by IH and NB. All authors participated to the methodology development and data analysis. The first draft of the manuscript was written by IH and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. Data Availabilty The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Declarations Conflict of interest The authors have no relevant financial or non-financial interests to disclose.