Guillard

E Drinks

Celia Lepeytre

C Lorentz

F Dappozze

Chantal Guillard
email: chantal.guillard@ircelyon.univ-lyon1.fr

UV-A photocatalytic degradation of the radionuclide complexants tributylphosphate and dibutylphosphate

Keywords: Tributylphosphate, dibutylphosphate, photocatalysis, 31P NMR, carboxylic acid

UV-A-a photocatalytic degradation of the radionuclide complexants tributylphosphate and dibutylphosphat

Introduction

The approach most commonly used to reprocess spent nuclear fuel involves three steps: a mechanical step, a dissolution step, and then liquid/liquid separation using a specific extractant [START_REF] Wei | Progress and Discussion on Chemical Separation Technologies for Nuclear Fuel Reprocessing Developed Abroad[END_REF]. In the conventional PUREX (Plutonium and Uranium Refining by EXtraction) chemical treatment process, uranium, plutonium, other actinides and most of the fission products are dissolved in a nitric acid solution. Uranium and plutonium are then separated and purified using tributylphosphate (TBP) as an extractant. This extractant is diluted using an organic compound, hydrogenated tetrapropylene in France [START_REF] Dodi | Improved determination of tributyl phosphate degradation products (mono-and dibutyl phosphates) by ion chromatography[END_REF].

The cycles of extraction and back extraction steps yield purified uranium and plutonium nitrates and after conversion, the resulting plutonium oxide is used to prepare uranium and plutonium MOX (mixed oxide) fuel. The nitric acid, extractant and solvent used in these processes have to be recycled and the liquid and gaseous effluents and solid waste have to be managed. Furthermore, the irradiation of TBP during the extraction process by the actinides and fission products, leads to the formation of radiolysis products, namely dibutylphosphate (DBP), monobutylphosphate, organic acids and butanol [START_REF] Dicholkar | Direct determination of tri-n-butyl phosphate by HPLC and GC methods[END_REF][START_REF] Gao | Studies on Pu-238 induced alpha radiolysis of the solvent TBP[END_REF]. Therefore, the secondary liquid waste could contain residual TBP and DBP, which can cause emulsification and act as actinide complexants. Tributylphosphate is also used in other industries, as an extractant for rare earth elements, in aircraft hydraulic fluid [START_REF] Ouyang | Research on Air Content Estimation of Tributyl Phosphate Hydraulic Fluids: A Novel Approach Based on the Vacuum Method[END_REF][START_REF] Marklund | Traffic as a source of organophosphorus flame retardants and plasticizers in snow[END_REF], as a plasticizer and flame retardant for vinyl resins [START_REF] Katsoulis | The combined effect of epoxy nanocomposites and phosphorus flame retardant additives on thermal and fire reaction properties of fiberreinforced composites[END_REF], as a wetting agent in the textile and paper industries, and as an anti-foaming agent.

The decontamination of liquid waste containing TBP and its degradation species is therefore a common industrial problem.

Tributylphosphate is a refractory molecule particularly difficult to degrade. The most promising results have been obtained using ultraviolet (UV) C radiation and TiO 2 nanoparticles in the presence of H 2 O 2 [START_REF] Watts | Photooxidation and subsequent biodegradability of recalcitrant tri-alkyl phosphates TCEP and TBP in water[END_REF][START_REF] Seshadri | Photocatalytic performance of combustion synthesized beta-Ga2O3 for the degradation of tri-n-butyl phosphate in aqueous solution[END_REF][START_REF] Seshadri | Synthesis, characterization and evaluation of nanocrystalline anatase titania for the degradation of dibutyl phosphate[END_REF]. This photocatalytic process follows first order kinetics [START_REF] Watts | Photooxidation and subsequent biodegradability of recalcitrant tri-alkyl phosphates TCEP and TBP in water[END_REF][START_REF] Seshadri | Photocatalytic performance of combustion synthesized beta-Ga2O3 for the degradation of tri-n-butyl phosphate in aqueous solution[END_REF][START_REF] Seshadri | Synthesis, characterization and evaluation of nanocrystalline anatase titania for the degradation of dibutyl phosphate[END_REF] and allows a 400 mg/L solution of TBP to be totally degraded in 70 min, and a 240 mg/L solution of DBP in roughly 45 min. The present study investigated the degradation of TBP and DBP by TiO 2 photocatalysis using UV-A irradiation. The degradation process was followed by recording the solution concentrations of residual carbon and phosphate ions, measured respectively by total organic carbon (TOC) analysis and ion chromatography. The appearance and disappearance of carboxylic acids, intermediate compounds in the reaction, was also monitored. To our knowledge, this article is the first to demonstrate TBP and DBP degradation using UV-A radiation and the first to provide a mechanistic explanation of the mode of degradation.

Experimental details

Materials

Photocatalysts

Photocatalytic experiments were performed using commercial (Degussa) P-25 TiO 2 nanoparticles with an 80/20 anatase/rutile ratio (diameter, 20-30 nm; density, 3.8 g/cm 3 ; specific surface area, about 50 m 2 /g).

Chemical products

Tri-and dibutylphosphate (99 and 97% purity, respectively) were purchased from Sigma-Aldrich and used as received. All the aqueous solutions used in this study were prepared using ultra-pure water with a resistivity of 18.2 MΩ•cm.

Photoreactor and light sources

All experiments were carried out in a 1 L cylindrical glass reactor containing a central 125 W Philips HPK lamp inside a Pyrex tube. A jacket with circulating water was placed around the Pyrex tube to absorb IR radiation and thus prevent any heating of the solutions.

The use of Pyrex glass ensured that < 290 nm radiation was blocked. The solutions were mixed constantly using a magnetic stirrer placed at the bottom of the reactor.

Procedure

The TiO 2 nanoparticles (375 mg) were added to 750 mL of aqueous solutions of TBP or DBP at different concentrations. The TiO 2 concentration (0.5 g/L) was chosen to ensure all the incident photons were absorbed. All reactions were conducted at room temperature, using ambient air as oxidant. The suspensions were kept under magnetic stirring to ensure adequate mixing and contact between the nanoparticles and the organic molecules. Before irradiation, the solutions were kept in the dark for 2 h under stirring to reach adsorption equilibrium between the pollutant and the photocatalyst. Samples were regularly withdrawn: 8 mL for gas chromatography-mass spectrometry (GC/MS), TOC, and carboxylic acid analyses, 5 mL for ion chromatography. The 8 mL samples were directly filtered through 0.45 µm Millipore filters to remove the TiO 2 particles. The 5 mL samples were basified at pH 13 using NaOH before filtration to desorb the phosphate ions from the surface of the photocatalyst.

Analytical Measurements

The overall efficiency of the process was determined by measuring the decrease in the concentration of total organic carbon using a Shimadzu TOC-VCPN analyzer equipped with an auto-sampler. The concentration of TBP was also measured by GC/MS using a DB-5MS capillary GC column (length, 30 m; inner diameter, 0.25 mm × 0.25 µm). The temperature conditions were 250°C for injection, 280°C for the detector, and the oven was initially kept at 50°C for 1 min, raised to 160°C at 10°C/min, then raised up to 270°C and kept at this temperature for 5 min. The helium flow was 1.3 mL/min and the mass spectrometer was operated in electron impact ionization mode at 70 eV.

The formation of phosphate ions was followed by ion chromatography using a Dionex DX-120 device and an AS14A (250 mm x 4 mm) column. The flow rate was 1.2 mL/min. The mobile phase consisted of 84 mg of NaHCO 3 and 848 mg of Na 2 CO 3 per liter of pure water.

The injection volume was 50 µl.

The formation of carboxylic acids was monitored by high-performance liquid chromatography with ultraviolet detection using a Waters 600 pump, a Waters 486 UV detector (detection at 210 nm), and a Transgenomic COREGel 87H3 7.8 mm (300 mm × 7.8 mm) column. The flow rate was 0.7 mL/min. The injection volume was 50 µL and the mobile phase was H 2 SO 4 (5 × 10 -3 mol/L).

Liquid-state 31 P-NMR analyses of the photocatalytic degradation of aqueous solutions of 100 ppm of TBP or DBM were performed by mixing 400 ml of solution irradiated for different times in 200 ml of D 2 O. The experiments were performed on a Bruker Avance HD 400 spectrometer equipped with a 5 mm Bruker BBFO probe, operating at a 31 P Larmor frequency of 161.97 MHz. The initial delay was 2 s and the 90° pulse length was 7.85 μs.

Results and discussion

Photocatalytic degradation of TBP

Disappearance

We first confirmed that the concentration of TBP in the solution did not decrease under direct irradiation at wavelengths above 300 nm without TiO 2 (photolysis, data not shown). Fig. 1 shows that TBP disappears completely from the solutions when both TiO 2 and of UV irradiation are present, indicating a purely photocatalytic process. The initial disappearance rate of TBP under these conditions was found to be 6 µmol/L/min.

`

The photonic efficiency, ρ = r/ ϕ, is the ratio of the reaction rate, r, the number of reactant molecules converted per second, and the incident flux of efficient photons, ϕ (i.e. the ones absorbable by the TiO 2 nanoparticles). Considering the number of photons emitted per second (290 nm (0.8 mW/cm2) + 365 nm (3.2 mW/cm2) and assuming all the photons are absorbed by the photocatalyst, this gives a photonic efficiency of 0.6%. This value is of the same order as those reported in the literature for the photocatalytic degradation of other aqueous pollutants [START_REF] Gregori | Hybrid sol-gel porous nanocomposites as efficient photocatalyticcoatings: Insights in the structure/reactivity relationships[END_REF].

Variation of the pH during the degradation of TBP

Fate of the phosphorus atoms

Fig. 3 shows that phosphate ions do not immediately appear following the disappearance of TBP, indicating that the formation of phosphate ions is not the first step in the degradation process. After disappearance of TBP only small amounts of phosphorus are transformed into phosphate ions indicating that the majority of organic compounds still contain phosphorus. After about 15 h, more than 95% of the phosphorus atoms are mineralized into phosphate ions (Fig. 3). The DBP peak appears after 1 h of irradiation along with trace peaks from H 2 PO 4 -and the unidentified organic phosphate. Subsequently, the latter two peaks increase in intensity while the DBP peak weakens. At 3 h of irradiation, new peaks appear next to the one from the unidentified compound and remain until 5 h of irradiation. Beyond 12 h of irradiation, only the H 2 PO 4 -peak is observed. mol/L for instance, more than 95% of the organic carbon is still present after 5 h of irradiation. This indicates that the first steps of the degradation involve either the oxidation of the carbonated chain and/or the breaking up of TBP. The appearance of DBP (Fig. 4) suggests that TBP does indeed break up. At the lowest concentration of TBP, the TOC concentration decreases more rapidly, suggesting that the surface of the TiO 2 nanoparticles was not saturated. Trace amounts of propionic acid are only detected after 20 h of irradiation. This late arrival suggests that it is a product of the decarboxylation of 2-ketobutyric acid:

Mineralization

Formation of carboxylic acid

CH 3 -CH 2 -CO-COOH + ½ O2 CH 3 CH 2 COOH + CO 2 (1)
Since the concentrations of lactic acid, pyruvic acid and acetic acid peak at different times, our interpretation is that lactic acid is progressively degraded by oxidation to form pyruvic acid, which is then decarboxylated to form acetic acid:

CH 3 -CHOH-COOH + ½ O2  CH 3 -CO-COOH  CH 3 COOH + CO 2 (2)
However, the amount of acetic acid formed during the reaction (Fig. 6a) indicates that this is not the only mechanism through which it is formed. The fact that acetic acid is the carboxylic acid that reaches the highest concentration during the reaction may be explained by its reactivity toward hydroxyl radicals (1.6 × 10 7 L/mol/s), which is about 100 times less than those of other carboxylic acids [START_REF] Buxton | Critical Review of rate constants for reactions of hydrated electrons, hydrogen atoms and hydroxyl radical in aqueous solution[END_REF]. Fig. 7 shows the concentrations of organic carbon in the form of carboxylic acid and the total concentration in the solution at different irradiation times. At 5 h, carboxylic acids account for about a third of the TOC present in solution, and after 20 h of irradiation of the TBP solution, all the organic carbon is in the form of carboxylic acids.

Comparison of the photocatalytic degradation of TBP and DBP

0,E+00 The photocatalytic degradation of alkylphosphate compounds in general was investigated by comparing those of TBP and DBP. This was achieved by comparing the disappearance of organic carbon, the fate of phosphate ions by 31 P NMR, and the formation of carboxylic acids.

Comparison of the mineralization of TBP and DBP

Fig. 8a shows the time required to photocatalytically degrade more than 99% of the carbon in solutions of TBP and DBP with different initial concentrations. Fig. 8b shows the corresponding average rates of carbon elimination. The time required to degrade more than 99% of the TOC increases linearly with the initial concentration of the pollutant and does not depend on its nature. At the same TOC concentration, the time required to remove more than 99% of the TOC is similar for TBP and DBP. Fig. 8b shows that the average mineralization rate increases linearly with the initial carbon concentration before reaching a plateau at concentrations greater than 4 × 10 -3 mol/L. This curve is typical of a Langmuir-Hinshelwood process, in which the reaction rate is given by: where k is the rate constant of the reaction and K is the equilibrium (or adsorption) constant.

The rate constants and adsorption constants calculated using these data for TBP and DBP are identical, within experimental errors (k ≈ 2.2 × 10 -4 mol/L/h and K ≈ 0.45 × 10 3 L/mol). In summary, this 31 P NMR analysis shows that during the degradation of DBP, an unidentified organic compound containing a phosphorus atom (resonating at 0.55 ppm) forms that is not detected during the degradation of TBP. The other unidentified peaks at 3.74 ppm and the weak shoulder at 3.74 ppm are also observed during the degradation of TBP. Further work is required to identify these two compounds. Acetic acid is the carboxylic acid that reaches the highest concentration during the degradation of both alkylphosphates, because of its low reactivity with hydroxyl radicals [START_REF] Buxton | Critical Review of rate constants for reactions of hydrated electrons, hydrogen atoms and hydroxyl radical in aqueous solution[END_REF].

Comparison of the fate of phosphate ions during TBP and DBP degradation.

Comparison of the carboxylic acids formed during the degradation of TBP and DBP.

The maximum concentrations measured for all the carboxylic acids are lower for the DBP solution and they disappear sooner than from the TBP solution, which is consistent with the lower carbon content of DBP.

Conclusions

The present study demonstrates for the first time that TBP and DBP can be mineralized using non-damaging radiation (UV-A), with an irradiance similar to that of solar light (4 mW/cm 2 ) and without adding an oxidant.

Our results show that less than 10% of the TOC of the solution has disappeared at the time all the TBP has degraded, indicating that the first steps in this process are the oxidation of the alkyl chain and/or the breaking of the bond between the phosphate and the alkyl group. The The same carboxylic acids are formed during the degradation of TBP and DBP, namely 2ketobutyric acid, lactic acid, pyruvic acid, formic acid, and acetic acid, the latter being the one found in the highest concentration because of its low reactivity with hydroxyl radicals. These carboxylic acids only appear once the alkylphosphates have completely disappeared. Pyruvic acid arises from the oxidation of lactic acid and is then transformed into acetic acid. In keeping with the smaller number of carbon atoms in DBP than in TBP, smaller amounts of carboxylic acids are formed during the degradation of DBP. Finally, our data also show that before being evacuated as CO 2 , all the organic carbon in the solution is in the form of carboxylic acids.

Fig. 1 .

 1 Fig. 1. Photocatalytic degradation of a solution of tributylphosphate (TBP) at an initial concentration of 3.9 × 10 -4 mol/L under UV-A irradiation and in the presence of 0.5 g/L TiO 2

Fig. 2

 2 Fig.2shows that the pH of the solution decreases rapidly in the initial stages of the

Fig. 2 .

 2 Fig. 2. Variation of pH as a function of UV-A irradiation time in an aqueous solution of

Fig. 3 .Fig. 4 .

 34 Fig. 3. Concentration of tributylphosphate and phosphate ions as a function of time in a solution with an initial concentration of 3.8 × 10 -4 mol/L tributylphosphate (100 ppm) under

Fig. 5 Fig. 5 .

 55 Fig.5shows the decrease in the TOC concentration as a function of time for three solutions with initial concentrations of between 4 × 10 -4 mol/L and 5 × 10 -5 mol/L TBP, highlighting

Fig. 6 Fig. 6 .

 66 Fig.6shows the concentrations of different carboxylic acids measured by HPLC

Fig. 7 .

 7 Fig. 7. Total concentrations of organic carbon in the solution and concentration in the form

Fig. 8 .

 8 Fig. 8. (a) Time taken to remove 99% of the total organic carbon from solutions of tri-and

Fig. 9

 9 Fig. 9 shows 31 P NMR spectra acquired at different times during the photocatalytic degradation of a DBP solution with a carbon content (4.8 × 10 -3 mol/L) similar to that of the

Fig. 9 .

 9 Fig. 9. 31 P NMR spectra acquired at different stages of the degradation of a solution of 5.63 × 10 -4 mol/L (100 ppm) dibutylphosphate, corresponding to 4.51 × 10 -3 mol/L carbon, under

Fig. 10 shows

 10 Fig. 10 shows how the concentrations of different carboxylic acids vary as a function of irradiation time during the degradation of a 4.1 × 10 -4 mol/L solution of DBP, and during

Fig. 10 .

 10 Fig. 10. Comparison of the time evolutions of the concentrations of carboxylic acids during the photocatalytic degradation of 4.1 × 10 -4 mol/L solutions of dibutylphosphate (filled

  by 31 P NMR spectra showing that DBP is formed during the degradation of TBP. The NMR data also show that similar phosphorus compounds form during the degradation of DBP and TBP, one of which is probably monobutylphosphate.Analysis of the TOCs of solutions of TBP and DBP at different concentrations suggest that both mineralization process have identical Langmuir-Hinshelwood kinetics. The average carbon disappearance rate and the time required to remove 99% of the organic carbon are directly proportional to the initial carbon content of the solution.

Acknowledgments

The authors thank ORANO inc. for financial support and Stephane Mangematin for help with the GC/MS analysis.


M Tayyab

S Zhao

Y Feng

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

Hybrid regularized Lattice

Keywords: lattice Boltzmann, reacting flows, combustion, detailed kinetics

published or not. The documents may come

Introduction

Motivated by the rapid development of Lattice-Boltzmann (LB) methods in the field of low-Mach external aerodynamics and aeroacoustics [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF], and particularly its potential as an engineering tool [START_REF] D'hooge | Application of real-world wind conditions for assessing aerodynamic drag for on-road range prediction[END_REF][START_REF] Gleason | Comparison of computational simulation of automotive spinning wheel flow field with full width moving belt wind tunnel results[END_REF][START_REF] Khorrami | Computational evaluation of airframe noise reduction concepts at full scale[END_REF][START_REF] Khorrami | Simulation-based airframe noise prediction of a full-scale, full aircraft[END_REF][START_REF] Casalino | Turbofan broadband noise prediction using the lattice boltzmann method[END_REF][START_REF] Romani | Rotorcraft blade-vortex interaction noise prediction using the lattice-boltzmann method[END_REF] we recently proposed a hybrid LB framework able to tackle combustion in low-Mach flows [START_REF] Feng | A lattice-boltzmann model for low-mach reactive flows[END_REF] in line with the quest of extending the LB capabilities to reactive flows [START_REF] Filippova | A Novel Lattice BGK Approach for Low Mach Number Combustion[END_REF][START_REF] Yamamoto | LB simulation on soot combustion in porous media[END_REF][START_REF] Chiavazzo | Coupling of the model reduction technique with the lattice Boltzmann method for combustion simulations[END_REF][START_REF] Xu | Multiple-relaxation-time lattice boltzmann kinetic model for combustion[END_REF][START_REF] Sun | A diffuse interface method for simulating the dynamics of premixed flames[END_REF][START_REF] Ashna | Extended lattice boltzmann scheme for droplet combustion[END_REF][START_REF] Hosseini | Mass-conserving advectiondiffusion lattice boltzmann model for multi-species reacting flows[END_REF].

Based on nearest neighbor lattices, the model previously proposed in [START_REF] Feng | A lattice-boltzmann model for low-mach reactive flows[END_REF] is easy to implement, but was derived under constant heat capacity and singlestep chemistry assumptions, a stark limitation compared to state-of-the-art combustion softwares. A second limitation resides in the simple collision model, merely derived from the BGK model, known to have serious limitations in complex flows [START_REF] Latt | Lattice boltzmann method with regularized precollision distribution functions[END_REF][START_REF] Malaspinas | General regularized boundary condition for multi-speed lattice boltzmann models[END_REF][START_REF] Feng | Regularized thermal lattice boltzmann method for natural convection with large temperature differences[END_REF][START_REF] Freitas | Analysis of latticeboltzmann methods for internal flows[END_REF][START_REF] Gehrke | Scrutinizing lattice boltzmann methods for direct numerical simulations of turbulent channel flows[END_REF][START_REF] Karlin | Gibbs' principle for the latticekinetic theory of fluid dynamics[END_REF][START_REF] Geier | The cumulant lattice Boltzmann equation in three dimensions: Theory and validation[END_REF].

This contribution aims at lifting these two limitations as to make the model fully functional for the simulation of reactive flows, including detailed chemistry description.

Lifting the first limitation is attained through integration of the classical thermodynamic closure based on NASA polynomial coefficients [START_REF] Esch | Thermodynamic properties in polynomial form for carbon, hydrogen, nitrogen, and oxygen systems from[END_REF], as used in Chemkin [START_REF] Kee | The chemkin thermodynamic data base[END_REF], or Cantera [START_REF] Goodwin | Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes[END_REF]. Detailed chemistry is also accounted for, through a 12-step mechanism for H 2 -air combustion [START_REF] Boivin | An explicit reduced mechanism for H 2 -air combustion[END_REF] derived from the San Diego mechanism [START_REF] Williams | Chemical-kinetic mechanisms for combustion applications[END_REF]. Each of the nine species is assigned a Schmidt number relating its diffusion properties with a temperature-dependent viscosity coefficient, and thermal diffusion is defined via a constant Prandtl number as in [START_REF] Schönfeld | Steady and unsteady flows simulations using the hybrid flow solver AVBP[END_REF].

The second limitation is equally lifted, after derivation of a new collision model compatible with thermal, multicomponent flows. This was achieved following Jacob et al [START_REF] Jacob | A new hybrid recursive regularised bhatnagar-gross-krook collision model for lattice boltzmann methodbased large eddy simulation[END_REF], who recently proposed a robust collision model in the context of single-component athermal LBM showing promising results for a wide range of applications [START_REF] Feng | Hybrid recursive regularized lattice boltzmann simulation of humid air with application to meteorological flows[END_REF][START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows[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][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] Wilhelm | An explicit power-law-based wall model for lattice boltzmann method-reynolds-averaged numerical simulations of the flow around airfoils[END_REF]. As further discussed in this study, the collision kernel adopts a regularization strategy [START_REF] Malaspinas | Increasing stability and accuracy of the lattice boltzmann scheme: recursivity and regularization[END_REF] in which the nonhydrodynamic moments of the distribution function (or ghost modes) relax infinitely fast to equilibrium. This strategy is similar to the popular Multiple Relaxation Time (MRT) models [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF][START_REF] Dubois | Towards higher order lattice Boltzmann schemes[END_REF][START_REF] Premnath | Three-dimensional multi-relaxation time (MRT) lattice-Boltzmann models for multiphase flow[END_REF][START_REF] Xu | Sensitivity analysis and determination of free relaxation parameters for the weakly-compressible MRT-LBM schemes[END_REF], in which these modes relax at a finite rate.

Validation is carried out by comparisons with Cantera computations [START_REF] Goodwin | Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes[END_REF] for the planar premixed flame propagation. To further demonstrate the capability of the model in premixed configurations, the Darrieus-Landau instability is simulated, and the associated growth-rate is compared with asymptotic descriptions [START_REF] Denet | Numerical study of thermal-diffusive instability of premixed flames[END_REF][START_REF] Denet | A numerical study of premixed flames darrieus-landau instability[END_REF][START_REF] Matalon | Intrinsic flame instabilities in premixed and nonpremixed combustion[END_REF], showing good agreement for the linear onset of the instability. To conclude the validation regarding premixed combustion, vortex-flame interactions are simulated, allowing to recover three regimes classically observed experimentally and numerically [START_REF] Rutland | Simulations of Flame-Vortex Interactions[END_REF][START_REF] Kadowaki | Numerical simulation of dynamics of premixed flames: flame instability and vortex-flame interaction[END_REF][START_REF] Poinsot | Quenching processes and premixed turbulent combustion diagrams[END_REF][START_REF] Roberts | Images of the quenching of a flame by a vortex-to quantify regimes of turbulent combustion[END_REF][START_REF] Mueller | Generation and attenuation of vorticity by flames: measured vorticity field time evolution during a premixed flame-vortex interaction[END_REF].

Next, the capabilities of the model are presented through comparison with Direct Numerical Simulations (DNS) on a double periodic shear-layer including auto-ignition and establishment of a diffusion flame in between the streams. That test case -besides being highly challenging numericallyhighlights the accuracy of the proposed model at a reasonable CPU cost.

The paper is organized as follows. First, the LB combustion model is presented, with extensive implementation details. The model is then validated, in premixed and non-premixed combustion test cases. Following the validation, cost and stability issues are discussed prior to drawing conclusions.

Governing equations: continuous formulation

Mass and momentum conservation: Lattice Boltzmann approach

The Boltzmann equation models the gas kinetics with the particle velocity distribution function (VDF) f (x, ξ, t), which presents the probability density of finding (virtual) gas particles at position x with velocity ξ at time t. Using the single-relaxation-time collision model of Bhatnagar-Gross-Krook (BGK) [START_REF] Bhatnagar | A model for collision processes in gases. i. small amplitude processes in charged and neutral onecomponent systems[END_REF], the evolution of VDF through phase space and time can be written as

∂f ∂t + ξ α ∂f ∂x α = - 1 τ (f -f eq ) , ( 1
)
where τ is the relaxation time associated to the dynamic viscosity of the fluid and f eq is the Maxwell-Boltzmann equilibrium distribution function associated to the local thermodynamic state

f eq = ρ (2πrT ) D/2 exp - |ξ -u| 2 2rT , (2)
with ρ and u being respectively the mixture mass volume and velocity, D the spatial dimension, T the temperature and r the gas constant of the mixture

r = R/W . (3)
R = 8.3145 J/mol/K is the universal gas constant and W is the molecular weight. In principle, ( 2) is valid for mono-constituent gas, for which the molecular weight is constant. Here, however, we'll consider an identical equilibrium distribution function for multi-component gas mixtures introducing the mean molecular weight

W = 1 k Y k /W k , (4)
where Y k is the k-th component mass fraction and W k its molecular weight.

It will be shown further that the multi-component macroscopic conservation equations are nonetheless recovered.

Under this formalism, the macroscopic variables ρ and u are obtained from the VDF through

ρ = f dξ , ρu = f ξdξ . (5)
The Chapman-Enskog (CE) expansion [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF] in the low-Mach number limit leads to the following macroscopic mass and momentum equations

∂ρ ∂t + ∂ ∂x α (ρu α ) = 0, ∂ρu α ∂t + ∂ ∂x β (ρu α u β ) = - ∂p ∂x α + ∂T αβ ∂x β , (6)
where

p = ρrT (7)
is the classical ideal gas thermodynamic closure, and T αβ is the viscous stress tensor

T αβ = 2µS αβ = µ ∂u β ∂x α + ∂u α ∂x β , (8)
with µ the dynamic viscosity. Note that the bulk viscosity is neglected, as often in the combustion community [START_REF] Poinsot | Theoretical and numerical combustion[END_REF][START_REF] Kuo | Principles of combustion[END_REF].

The CE expansion, provided in AppendixA, differs from the classical one [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF] in that the perfect gas equation of state is recovered instead of the athermal EOS. It also differs from our recent study [START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows[END_REF] as (i) the gas is multi-component, e.g. r in (3) is not a constant and (ii) additional equations for gas components and energy are required to close the system, to account for non-unity Prandtl and Lewis numbers, as well as source terms between the gas components. AppendixA should be read along with Section 3, as it provides most technical details as well as a detailed discussion.

Energy and species conservation

Mass conservation of species k is considered following

ρ ∂Y k ∂t + ρu α ∂Y k ∂x α = ∂ ∂x α (-ρV k,α ) + ωk , (9)
where ωk is the net chemical production rate of species k, and V k,α is its diffusion velocity [START_REF] Poinsot | Theoretical and numerical combustion[END_REF].

The energy conservation is considered by the balance equation of the enthalpy. Following the multi-component ideal gas thermodynamic closure, the enthalpy of a gas mixture consisting of N species is defined as

h = N k=1 h k Y k , h k = T T 0 C p,k (T )dT + ∆h 0 f,k , (10)
where C p,k (T ) is the constant pressure heat capacity of species k at temperature T , and ∆h 0 f,k its formation enthalpy [START_REF] Poinsot | Theoretical and numerical combustion[END_REF]. The conservation law of enthalpy reads

ρ ∂h ∂t + ρu α ∂h ∂x α = Dp Dt - ∂q α ∂x α + T αβ ∂u α ∂x β . ( 11
)
where Dp Dt = ∂p ∂t + u α ∂p ∂xα is the substantial derivative of pressure, usually neglected in low-Mach combustion applications. The heat flux q α reads

q α = -λ ∂T ∂x α + ρ N k=1 h k Y k V k,α , (12)
with λ the thermal conductivity.

Model implementation

Discretization on standard lattices

Following classical LB models, the particle velocity ξ in Eq. ( 1) is discretized as a finite dimension velocity space. Our choice lies with nearest neighbors lattices (D2Q9, D3Q19) for their low associated cost. Such discretization strategy being found throughout the LBM literature (e.g. [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF]), let us simply indicate that the choice of lattice DnQm set defines, besides the discrete velocities c i (i = 1, • • • , m), the lattice sound speed c s and w i , the Gaussian weights associated to each discrete velocity c i .

With these definitions,

f i (x, t) = w i ω(c i ) f (x, c i , t), (13)
and Eq. ( 1) can be written at the accuracy of the quadrature as

∂f i ∂t + c iα ∂f i ∂x α = - 1 τ (f i -f eq i ) . ( 14
)
This equation is further discretized in space and time with second-order accuracy [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF] as

f i (x α + c iα δt, t + δt) = f i (x α , t) - δt τ [f i (x α , t) -f eq i (x α , t)], (15)
where the time step δt is classically linked to the spatial discretization through definition of the set of lattices DnQm [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF]. Note that introduction of f i and τ in [START_REF] Hosseini | Mass-conserving advectiondiffusion lattice boltzmann model for multi-species reacting flows[END_REF] as

f i = f i - δt 2 Ω i = f i - δt 2 - 1 τ (f i -f eq i ) , (16)
and

τ = τ + δt 2 , ( 17
)
is a necessary change of variable to ensure the second-order accuracy in time [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF]. In practice, all physical quantities are normalized in the LBM code using the grid size δx for length, the time step δt for time and a reference density ρ 0 (combined with δx) for mass.

Hybrid regularized collision model

Single relaxation time models often lead to numerical instabilities in shear-flows [START_REF] Hosseini | Theoretical and numerical analysis of the lattice kinetic scheme for complex-flow simulations[END_REF], and are unable to tackle the test case of Fig. 5. To increase numerical stability in such flows, a possibility is to include multiple relaxation times models (see, e.g. [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF][START_REF] Dubois | Towards higher order lattice Boltzmann schemes[END_REF][START_REF] Premnath | Three-dimensional multi-relaxation time (MRT) lattice-Boltzmann models for multiphase flow[END_REF][START_REF] Xu | Sensitivity analysis and determination of free relaxation parameters for the weakly-compressible MRT-LBM schemes[END_REF]), in which relaxation times for the higher-order moments (neither corresponding to the mass conservation nor the momentum) are carefully tuned to optimize the stability properties. On the same idea, it was later proposed [START_REF] Malaspinas | General regularized boundary condition for multi-speed lattice boltzmann models[END_REF] to simply suppress the non-equilibrium part of these modes, in so-called regularized collision kernels. Here a regularization step is adapted from [START_REF] Latt | Lattice boltzmann method with regularized precollision distribution functions[END_REF][START_REF] Feng | Regularized thermal lattice boltzmann method for natural convection with large temperature differences[END_REF][START_REF] Jacob | A new hybrid recursive regularised bhatnagar-gross-krook collision model for lattice boltzmann methodbased large eddy simulation[END_REF] to the multi-component flow model. The regularization strategy reconstructs f i before collision by considering

f i = f eq i + f neq i , (18)
with both the equilibrium (f eq i ) and off-equilibrium part (f neq i ) of the VDF confined in the Hermite basis. Following [START_REF] Feng | A lattice-boltzmann model for low-mach reactive flows[END_REF][START_REF] Feng | Regularized thermal lattice boltzmann method for natural convection with large temperature differences[END_REF][START_REF] Jacob | A new hybrid recursive regularised bhatnagar-gross-krook collision model for lattice boltzmann methodbased large eddy simulation[END_REF], the equilibrium VDF are truncated up to third order to reduce the non-Galilean defect

f eq i = w i a (0),eq + a (1),eq α c 2 s H (1) α (c i ) + a (2),eq αβ 2c 4 s H (2)
αβ (c i ) + a (3),eq αβγ 6c 6 s H (3) αβγ (c i ) , (19)
with the Hermite polynomial tensors associated to each discrete velocity c i defined as

H (1) α (c i ) = c iα , H (2)
αβ (c i ) = c iα c iβ -δ αβ c 2 s , H (3)
αβγ (c i ) = c iα c iβ c iγ -c 2 s (c i δ) αβγ , (20)
and the associated equilibrium coefficients a (i),eq reading a (0),eq = ρ, (21a)

a (1),eq α = ρu α , (21b) a
(2),eq αβ

= ρ u α u β + c s 2 (θ -1) δ αβ , (21c) a
(3),eq

αβγ = ρ u α u β u γ + c s 2 (θ -1) (uδ) αβγ , ( 21d
)
where θ is the normalized temperature, defined as in our initial model [START_REF] Feng | A lattice-boltzmann model for low-mach reactive flows[END_REF]

θ = rT c s 2 , (22)
and the tensor uδ reads

(uδ) αβγ = δ αβ u γ + δ γβ u α + δ αγ u β . (23)
Note that, consistently with the Hermite space, the lattice sound speed is linked to the space-time discretization as c s ≡ δx/δt/ √ 3 [START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF].

The non-equilibrium part of VDF are truncated up to second order to recover a correct viscous tensor [START_REF] Jacob | A new hybrid recursive regularised bhatnagar-gross-krook collision model for lattice boltzmann methodbased large eddy simulation[END_REF]

f neq i = w i a (2),neq αβ 2c 4 s H (2) αβ (c i ). (24)
In the present regularized collision model, the off-equilibrium coefficients a

(2),neq αβ are evaluated through two different approximations, hence its given name: hybrid regularized collision model [START_REF] Jacob | A new hybrid recursive regularised bhatnagar-gross-krook collision model for lattice boltzmann methodbased large eddy simulation[END_REF].

The first approximation arises from direct projection of

f neq i as a (2),neq,PR αβ = m i=1 f i -f eq i H (2) αβ (c i ) (25)
while the second approximation is computed as to best approximate the viscous tensor S αβ in Eq. ( 8). That term is assessed via finite difference method (hence the FD subscript)

a (2),neq,FD αβ = -2pτ S FD αβ . (26)
As shown by Jacob et al. [START_REF] Jacob | A new hybrid recursive regularised bhatnagar-gross-krook collision model for lattice boltzmann methodbased large eddy simulation[END_REF] and further demonstrated in this article, excellent numerical stability of the collision model is obtained through combination of ( 25) and [START_REF] Boivin | An explicit reduced mechanism for H 2 -air combustion[END_REF] as

a (2),neq αβ = σa (2),neq,PR αβ + (1 -σ)a (2),neq,FD αβ (27)
where σ ∈ [0, 1] is a free parameter.

Correction term and final model

Last, let us introduce an additional forcing term to be included in the LBM equation [START_REF] Hosseini | Mass-conserving advectiondiffusion lattice boltzmann model for multi-species reacting flows[END_REF]. That correction corresponds to the deviation due to defect of symmetry of lattices (D2Q9, D3Q19, D3Q27) on the third order moment [START_REF] Feng | A lattice-boltzmann model for low-mach reactive flows[END_REF][START_REF] Feng | A three dimensional lattice model for thermal compressible flow on standard lattices[END_REF], and reads

s i = w i 2c 4 s H (2) αα (c i )∂ α ρu α (3c 2 s -3c 2 s θ -u 2 α ) , (28)
affecting the LB equation ( 15) as

f i (x α + c iα δt, t + δt) = f eq i (x α , t) + 1 - δt τ f neq i + δt 2 s i . (29)
The correction term's first two moments being mass and momentum preserving, the reconstruction of the macroscopic variables from the distribution functions remain unaffected:

ρ = m i=1 f i , ρu α = m i=1 c iα f i . (30)
Note that, as shown in [START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows[END_REF], additional correction terms are necessary for high-Mach flows with the D3Q19 stencil. They are neglected here as the Mach numbers encountered remain moderate.

FD solver specifics

In parallel of the LB solver, we use a FD solver to solve the species (9) and energy [START_REF] Chiavazzo | Coupling of the model reduction technique with the lattice Boltzmann method for combustion simulations[END_REF] conservation equations. The flowchart of the code is identical to that presented in our initial model [START_REF] Feng | A lattice-boltzmann model for low-mach reactive flows[END_REF]. Through systematic testing of the different finite differences methods, we retained the method presented in this Section.

The advection-diffusion part of ( 9) and ( 11) (e.g. without source term) is solved using the following strategy, arguably the simplest and most robust one could think of:

1. The LB and FD solver use the same time-step, and first order forward Euler integration is used. It is important to note that strong-stabilitypreserving Runge-Kutta (SSPRK) methods [START_REF] Gottlieb | Strong stability-preserving highorder time discretization methods[END_REF] up to fourth order have been implemented, with no apparent effect on the validations carried out.

2. Second order central difference operators are used for convection and diffusion. Note that introducing 10% of first-order up-winding on the convective operator was found to enhance numerical stability in high shear flows (as in Section 5)

The chemical source term is applied through a Strang-method splitting [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF][START_REF] Strang | On the construction and comparison of difference schemes[END_REF], allowing the use of sub time-steps for the chemical integration when necessary (as for an implicit chemical solver). Similarly, the diffusion term may be easily split from the advection step, should it become necessary.

Boundary conditions

This model being new, future work is required to present implementation of complex boundary conditions. Following our initial strategy [START_REF] Feng | A lattice-boltzmann model for low-mach reactive flows[END_REF], we picked validation test cases requiring rather simple boundary conditions. Besides the trivial periodic boundary condition of the double shear layer test case, only basic velocity inlet and zero-gradient outlet boundary conditions were necessary for this work.

A non-equilibrium bounce back scheme with a regularization procedure [START_REF] Malaspinas | General regularized boundary condition for multi-speed lattice boltzmann models[END_REF] is employed at the inlets for the LB solver, allowing to set all velocity components and the pressure (or density). (h, Y k ) are hard-coded at inlets for the FD solver accordingly with the target temperature and composition.

Outlets are treated as zero-gradient boundaries on microscopic distribution functions f i as well as h and Y k variables.

Premixed combustion validation test cases

Let us now present validations and discussions of the hybrid model presented above. To that aim, the model was implemented in the ProLB solver [START_REF] Feng | Regularized thermal lattice boltzmann method for natural convection with large temperature differences[END_REF][START_REF] Jacob | A new hybrid recursive regularised bhatnagar-gross-krook collision model for lattice boltzmann methodbased large eddy simulation[END_REF][START_REF] Feng | A three dimensional lattice model for thermal compressible flow on standard lattices[END_REF][START_REF] Marié | Comparison between lattice boltzmann method and navier-stokes high order schemes for computational aeroacoustics[END_REF] considering the following aspects:

Thermodynamic closure is ensured via the classical NASA polynomial formulation for each gaseous constituent [START_REF] Esch | Thermodynamic properties in polynomial form for carbon, hydrogen, nitrogen, and oxygen systems from[END_REF], as in Chemkin [START_REF] Kee | The chemkin thermodynamic data base[END_REF].

Viscosity is set via a temperature power-law (B.2).

Transport of heat and species are defined via constant Prandtl (B.1) and individual Schmidt numbers (B.4), as done for instance in the AVBP solver [START_REF] Schönfeld | Steady and unsteady flows simulations using the hybrid flow solver AVBP[END_REF].

Chemistry for hydrogen-air is accounted for via a 12-step skeletal mechanism [START_REF] Boivin | An explicit reduced mechanism for H 2 -air combustion[END_REF] involving nine species (see Tab. B.4).

Details about these four points are gathered in AppendixB. Since they were not necessary in the development of Sections 2 and 3, they may be easily changed in future works.

The solver being based on the D3Q19 lattice, 1D results presented below use periodic conditions along the y and z axis (with a single cell in either direction), and 2D results use periodic conditions along the z axis. 3D results will be presented in future publications.

1D premixed flame

Let us now consider the standard freely propagating flame configuration, following the setup of [START_REF] Feng | A lattice-boltzmann model for low-mach reactive flows[END_REF].

The computational domain is pseudo one-dimensional of length L = 10 mm, with grid size of 10 -5 m. At the left of the domain, the velocity is set to an arbitrary value U f , and the right boundary is left open (zero gradient).

Periodic boundary conditions are applied to the other boundaries (corresponding to the y and z directions). The time step is δ t = 7.21 × 10 -9 s.

Initial conditions consist of two half-domains initially at u = U f , and a sharp transition at L/2, separating the fresh and burnt gases. The left side of the domain corresponds to the fresh gases, whereas the right side corresponds to the burnt gases initialized at the corresponding Canteracomputed thermo-chemical equilibrium -provided in Table B.5 of AppendixB.

As in our previous contribution [START_REF] Feng | A lattice-boltzmann model for low-mach reactive flows[END_REF], we measure the flame speed as

S L = lim t→∞ U f - ρ f U f -ρ b U b ρ f -ρ b , (31)
where the subscript f represents fresh and b burnt state of the gases, respectively the first and last cell of the computational domain. History of this expression ( 31) is monitored until it converges to the flame velocity, independently of the inlet velocity U f [START_REF] Feng | A lattice-boltzmann model for low-mach reactive flows[END_REF].

The results obtained for U f = 0 are compared with Cantera reference computations in Fig. 1.a. The temperature and mass fractions profiles show an excellent agreement, including for the minor species H and HO 2 , whose profiles in the flame appear indistinguishable from the reference.

The dependence of S L with equivalence ratio is shown in Fig. 1.b, also showing a good agreement -below 10 -3 relative error at ϕ = 1 -besides the important effective Lewis number variation for this range of equivalence ratio [START_REF] Fernández-Galisteo | The differential diffusion effect of the intermediate species on the stability of premixed flames propagating in microchannels[END_REF].

2D Darrieus-Landau instability

To domain with grid size of 10 -5 m. At t = 0, the 2D domain is initialized with the profiles obtained from the 1D computation reported in Fig. 1.a. The front position, defined as the maximum temperature gradient position is then perturbed with a wavenumber (k) of 600m -1 . The evolution of the front position is reported in Fig. 2.a.

Qualitatively, it is interesting to note that the evolution produces the expected behavior, and that the coupling between the LB and FD solvers is robust to more complex flame front shapes. More quantitatively, Figure 2.b concludes the validation of this test case by comparing the growth rate in the linear regime with the analytical solution [START_REF] Denet | Numerical study of thermal-diffusive instability of premixed flames[END_REF][START_REF] Denet | A numerical study of premixed flames darrieus-landau instability[END_REF][START_REF] Matalon | Intrinsic flame instabilities in premixed and nonpremixed combustion[END_REF]. As expected, the amplitude A of the perturbation grows exponentially from its initial value A 0 according to

A = A 0 e ωt , (32)
with ω the perturbation growth rate. Theory indicates that, in the linear a) regime [START_REF] Denet | Numerical study of thermal-diffusive instability of premixed flames[END_REF][START_REF] Denet | A numerical study of premixed flames darrieus-landau instability[END_REF][START_REF] Matalon | Intrinsic flame instabilities in premixed and nonpremixed combustion[END_REF],

ω = S L k -σ + √ σ 3 + σ 2 -σ σ + 1 , (33)
with k the wave number and σ = ρ f /ρ b = 6.822 the ratio between the fresh and burnt gases density.

Flame-vortex interaction

Let us conclude the validation of the model for premixed combustion by testing the classical flame-vortex interaction.

For this study, we will consider the propane-air single-step chemistry

C 3 H 8 + 5O 2 → 3CO 2 + 4H 2
O detailed in our previous study [START_REF] Feng | A lattice-boltzmann model for low-mach reactive flows[END_REF]. Parameters for the thermodynamic closure are the classical NASA polynomials [START_REF] Linstrom | Nist chemistry webbook[END_REF], and all species are assigned unity Lewis numbers

Sc k = P r = 0.682, ∀k. (34)
A two-dimensional 4mm×2mm domain is discretized (δx = 10 -5 m) and initialized with the results of a planar premixed flame computation at equivalence ratio 0.8, corresponding to an expansion ratio of

       u x = U 0 -ǫ y-y 0 R 0 exp -(x-x 0 ) 2 +(y-y 0 ) 2 2R 2 0 , u y = ǫ x-x 0 R 0 exp -(x-x 0 ) 2 +(y-y 0 ) 2 2R 2 0 , ( 35
)
where ǫ is the vortex strength, R 0 the characteristic radius and (x 0 , y 0 ) are the coordinates of the vortex center. Three values of the vortex strength ǫ = (5, 30, 65) were adopted, corresponding to weak, moderate and strong interactions -referred hereafter as cases A, B and C. The initial distance of the vortex center to flame is 0.5mm (= 2δ L ) at t = 0. A summary of the initialization setup is provided in the sketch of Fig. 3. To characterize the three regimes (A, B and C), we measure the following parameters at the instant the vortices start interacting with the flame:

• Vm s L , with V m the effective vortex strength,

• σ δ L , with σ the effective vortex radius.

Following [START_REF] Kadowaki | Numerical simulation of dynamics of premixed flames: flame instability and vortex-flame interaction[END_REF], the Karlowitz number is defined as the ratio of these quantities

Ka = V m δ L s L σ . ( 36
)
The effective vortex strength, radius and Karlowitz number for the three cases are reported in Table 1. Case A Flame wrinkling appears and decays as the vortex intensity decreases.

V m /s L σ/δ L Ka Case A 3.
Case B Moderate interaction leads to the formation of an unburnt gas pocket, which detaches and burns at a later time.

Case C For strong vortices, a significant elongation of the unburnt region appears into the burnt gases, before eventually separating and burning.

These observations are in line with the results reported in the literature [START_REF] Rutland | Simulations of Flame-Vortex Interactions[END_REF][START_REF] Kadowaki | Numerical simulation of dynamics of premixed flames: flame instability and vortex-flame interaction[END_REF][START_REF] Poinsot | Quenching processes and premixed turbulent combustion diagrams[END_REF][START_REF] Roberts | Images of the quenching of a flame by a vortex-to quantify regimes of turbulent combustion[END_REF][START_REF] Mueller | Generation and attenuation of vorticity by flames: measured vorticity field time evolution during a premixed flame-vortex interaction[END_REF].

Non-premixed combustion in a double mixing layer

Design and numerical set-up of the test case

To illustrate the potential of the method, an ad-hoc test case illustrated in Fig. 5 was adapted from the athermal test case presented in [START_REF] Hosseini | Theoretical and numerical analysis of the lattice kinetic scheme for complex-flow simulations[END_REF] to show the stability of the collision model derived in Sec. layer appeared to be a reasonable choice as it is known to be challenging to LB collision kernels [START_REF] Hosseini | Theoretical and numerical analysis of the lattice kinetic scheme for complex-flow simulations[END_REF], the classical BGK model almost always leading to failure.

For the results to resemble the incompressible test case presented in [START_REF] Hosseini | Theoretical and numerical analysis of the lattice kinetic scheme for complex-flow simulations[END_REF],

we decided to keep the same initial velocity field, Reynolds number (30000) and momentum ratio between the streams. Next, we set the initial temperature and composition of the two streams so that a diffusion flame would auto-ignite and stabilize in between (e.g. one oxidizer and one fuel stream, as shown in Fig. 5). This test case was imagined to be as close as can be to a jet diffusion flame, whilst keeping simple (periodic) boundary conditions.

Under these considerations, we considered a [0, L] × [0, L] 2D domain, with L = 5cm, in which initial conditions are fully defined through

                                         u x (x, y) = U 0 tanh λ 1 4 -| y L -1 4 | u y (x, y) = ǫU 0 sin 2π x L + 1 4 , p(x, y) = p 0 ρ(x, y) = ρ 0 , Y H 2 (x, y) = Y F H 2 1 + tanh λ 1 4 -| y L -1 4 | /2, Y O 2 (x, y) = Y O O 2 1 -tanh λ 1 4 -| y L -1 4 | /2, Y N 2 (x, y) = 1 -Y H 2 (x, y) -Y O 2 (x, y), (37)
corresponding to the smooth profiles illustrated in the first contour of Fig. 5, definition proposed in [START_REF] Hosseini | Theoretical and numerical analysis of the lattice kinetic scheme for complex-flow simulations[END_REF], the Reynolds number can be obtained as

using
Re = ρ 0 U 0 L µ = 30 000, (38)
using the fuel stream parameters from Table 2, and the viscosity obtained from (B.2) at T F . Note that the alternative value of 37 300 may be obtained based on the less viscous oxidizer stream. It is also useful for the discussion to define the convective time t c as

t c = L U 0 , (39)
The flow evolution was simulated up to 2t c on a 1024 2 cartesian grid, using both our hybrid LB solver and the DNS solver Ntmix [START_REF] Baum | Ntmix-ckemkin release 2, user's guide version 1.0[END_REF], leading to the temperature contours displayed in the second and third plots of Fig. 5.

Ntmix is an in-house parallel DNS flow solver developed at CERFACS, Toulouse. Space is discretized through 6 th or 8 th order centered finite differences schemes, and time integration is performed with a 3 rd order Runge-Kutta. Ntmix is coupled with a Chemkin version for the computation of transport, thermodynamics and kinetics properties. The simulations conducted for this work use the 8 th order discretization, mixture-averaged transport properties [START_REF] Kee | CHEMKIN-III: A FOR-TRAN chemical kinetics package for the analysis of gas-phase chemical and plasma kinetics[END_REF], the 12-step H 2 -air chemistry [START_REF] Boivin | An explicit reduced mechanism for H 2 -air combustion[END_REF], and neglect the Soret effect.

For the results to be comparable with both codes, the pressure term in [START_REF] Chiavazzo | Coupling of the model reduction technique with the lattice Boltzmann method for combustion simulations[END_REF] Dp Dt was approximated as ∂p ∂t , as to obtain a good agreement on the pressure and temperature fields in this confined configuration.

Validation of the velocity field

In order to validate the flow-field, we successively simulated the coldflow (e.g. without kinetic source term) and the reacting flow. Figure 6 presents vorticity contours after one convective time t = t c , showing excellent qualitative agreement. Let us now confirm this validation quantitatively by plotting the profiles along the diagonal lines of the domain identified in Fig. 5. This is done in Fig. 7: for each simulation, we plotted both velocity components (u x , u y ) -easily distinguishable given the flow configurationboth for the cold and hot test case, and at two successive convective times (t = t c , t = 2t c ). Superimposed in light gray is the domain diagonal along which the profile is plotted. Besides small discrepancies appearing at 2t c , the agreement is overall excellent, hereby confirming the results of Fig. 6.

Reacting flow validation

Let us now focus on the reactive case. Figures 8 and9 present profiles of temperature, as well as H 2 , H and H 2 O mass fractions along the two domain diagonals. Again, agreement along the diagonal passing through the two vortices (Fig. 8) is seen to be excellent after both one and two convective times. The profiles on the second diagonal (Fig. 9) are equally good after one 0. times.

3. Transport models are different in both solver: constant Lewis numbers are assumed in the LB solver, whereas Ntmix solver assumes mixtureaveraged properties.

Factor 1, in our opinion, is responsible for most of the departure appearing in Fig. 9 after two convective times. The agreement remains, however, very good for the Lattice-Boltzmann simulation of such a complex reactive flow configuration.

Robustness, accuracy and computational cost

Let us conclude this study by a discussion on robustness and accuracy.

δ L = T b -T f max ∂T ∂x , (40)
with T f and T b the temperature of the fresh and the burnt gases, respectively, as given in Tab. B.5. The reference for the error is the velocity as obtained on an over-refined 2000 points grid -which we found to be less than 0.1% off the Cantera result presented in Fig. 1.a. Note that, following this definition, the profiles of Fig. 1.a corresponds to δ L δ x = 16 points per flame thickness.

As identified by the line, the convergence is second-order with the spatial discretization. This result is surprising, since the convection term is a simple second-order central difference scheme including 10% of upwinding (which should degrade the order of accuracy), and that the solvers coupling is firstorder. Note, however, that most test cases do not require any upwinding, and that the diffusion scheme is second-order. We have also tested higher-order schemes for the convection term and the FD temporal integration, without appreciable difference in the convergence. This seems to indicate that the obtained solver mostly keeps the low-dissipation feature characteristic of LB methods [START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows[END_REF][START_REF] Marié | Comparison between lattice boltzmann method and navier-stokes high order schemes for computational aeroacoustics[END_REF].

Another interesting result to be noted from Fig. 10 is the robustness of the method: the algorithm does converge even with only 5 points in the flame, an excellent result considering the stiffness of the H 2 -O 2 chemistry -confirmed by the sharpness of the intermediate species in Fig. 1.

To further assess the robustness of the method, we ran again the double shear layer while reducing the number of grid points. assessment of our collision kernel associated hyper-viscosity, responsible for this robustness, will be carefully conducted in future works.

To close this Section, let us add a note on the computational cost involved, which may serve as perspective for practical Large Eddy Simulations applications [START_REF] Löhner | Towards overcoming the LES crisis[END_REF]. This result is encouraging: it corresponds to 5.1 µs per grid point, which, assuming perfect scaling translates to less than 150 000 / cpuh to run 10 6 time-steps on 10 8 grid points.

Use of the hybrid regularized collision model costs a mere 15% more than the conventional single relaxation time (SRT) model, which only provides the required stability for the planar premixed flame propagation presented in Fig. 1. We also measured the relative cost of the cold and hot flow simulation with the athermal simulation (e.g. θ ≡ 1, and no finite difference solver). Figures are presented in Table 3. These results will be improved in the future as no Table 3: Relative computational costs of the method on the double shear layer. Athermal flow : as in [START_REF] Jacob | A new hybrid recursive regularised bhatnagar-gross-krook collision model for lattice boltzmann methodbased large eddy simulation[END_REF]; Cold flow: PR collision with the transport equations for energy and 9 species; Hot flow: cold flow configuration plus chemical sources evaluated using the 12-step chemistry.

optimization of any kind was performed on the FD solver at this stage, but are certainly promising: including a relatively simple chemistry adds about 30% computational cost, indicating a wide margin for progress. Also, the ratio of the reacting simulation cost to the athermal flow is quite classical, indicating that gains in cpu costs obtained in aerodynamics and aeroacoustics [START_REF] D'hooge | Application of real-world wind conditions for assessing aerodynamic drag for on-road range prediction[END_REF][START_REF] Gleason | Comparison of computational simulation of automotive spinning wheel flow field with full width moving belt wind tunnel results[END_REF][START_REF] Khorrami | Computational evaluation of airframe noise reduction concepts at full scale[END_REF][START_REF] Khorrami | Simulation-based airframe noise prediction of a full-scale, full aircraft[END_REF][START_REF] Casalino | Turbofan broadband noise prediction using the lattice boltzmann method[END_REF] may be obtained in the field of combustion in the near future.

Concluding remarks

A new hybrid Lattice-Boltzmann model for low-Mach reactive flows was presented following [START_REF] Feng | A lattice-boltzmann model for low-mach reactive flows[END_REF]. Mass and momentum conservation are addressed within a Lattice-Boltzmann solver, whereas the energy and species conservation are addressed via a classical finite difference solver.

To cope with the limitations of the standard BGK collision model [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF], an ad-hoc model model was presented as an extension to thermal multicomponent flows of Jacob et al's proposal [START_REF] Jacob | A new hybrid recursive regularised bhatnagar-gross-krook collision model for lattice boltzmann methodbased large eddy simulation[END_REF].

We then presented major progress compared to our initial proposal [START_REF] Feng | A lattice-boltzmann model for low-mach reactive flows[END_REF]:

not

AppendixA. Chapman-Enskog expansion

The CE expansion of equation ( 1) is performed in the Hermite polynomial space [START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF]. The Hermite tensor applied in the current study at order n reads

H (n) α 1 •••αn (ξ) ≡ (-r 0 T 0 ) n ω(ξ) (∂ α 1 • • • ∂ αn ω(ξ)) , (A.1)
and the weight function ω(ξ) is defined as

ω(ξ) ≡ 1 (2πc s 2 ) D/2 exp - ξ 2 2c s 2 , (A.2)
with ξ 2 ≡ ξ • ξ and c s the reference speed of sound

c s ≡ r 0 T 0 . (A.3)
where r 0 , T 0 are the reference gas constant and reference temperature respectively. It is worth noting that compared to weight functions used in thermal LBMs [START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF][START_REF] Feng | A three dimensional lattice model for thermal compressible flow on standard lattices[END_REF], a reference gas constant is introduced to overcome the difficulties brought out by the non-constant molecular weights of the gas mixture.

The VDF f can be projected to its Hermite coefficients tensor through

a (n) α 1 •••αn (x, t) ≡ f (x, ξ, t)H (n) α 1 •••αn (ξ)dξ . (A.4)
According to the orthogonality of Hermite polynomials, the VDF can be projected back from the Hermite space via

f (x, ξ, t) = ω(ξ) ∞ n=0 1 (n)!c s 2n a (n) (x, t) : H (n) (ξ), (A.5)
where the operator (:) stands for full contraction of two tensors. Applying the projection (A.4) to the equilibrium distribution (2), the equilibrium coefficient tensors provided in Eq. ( 21) are obtained.

Projecting equation ( 1) to the Hermite space, one has

∂ ∂t a (n) α 1 •••αn + ∂ ∂x j ξ j a (n) α 1 •••αn = - 1 τ a (n) α 1 •••αn -a (n),eq α 1 •••αn (A.6)
Note that the index j repeats from 1 to D (the spatial dimension). Using the Rodrigues recursive relation of the Hermite tensor

ξ j H (n) α 1 •••αn (ξ) = H (n+1) jα 1 •••αn (ξ) + c s 2 n i=1 δ jα i H (n-1) α 1 •••α i-1 α i+1 •••αn (ξ) , (A.7) equation (A.6) becomes ∂ t a (n) α 1 •••αn +∂ j a (n+1) jα 1 •••αn +c s 2 n i=1 ∂ α i a (n-1) α 1 •••α i-1 α i+1 •••αn = - 1 τ a (n) α 1 •••αn -a (n),eq α 1 •••αn . (A.8)
The CE analysis is performed on this equation. The scale separation respect to Knudsen number ǫ is achieved by setting the operator

∂ t = ǫ∂ t 1 + ǫ 2 ∂ t 2 ; ∂ j = ǫ∂ j 1 , (A.9)
and the coefficients tensor

a (n) = a (n),eq + a (n),1 , where a (n),1 /a (n),eq ∼ O(ǫ) ≪ 1 . (A.10)
It is worth noting that a (0),1 = a

(1),1 α = 0 according to equation [START_REF] Khorrami | Simulation-based airframe noise prediction of a full-scale, full aircraft[END_REF].

Inserting (A.9) and (A.10) into (A.6), the system can be written at different orders of ǫ as

ǫ∂ t 1 a (n),eq α 1 •••αn + ǫ∂ j 1 a (n+1),eq jα 1 •••αn + c s 2 n i=1 ǫ∂ α i 1 a (n-1),eq α 1 •••α i-1 α i+1 •••αn = - 1 τ a (n),1 α 1 •••αn ∼ O(ǫ 1 ) (A.11a) ǫ 2 ∂ t 2 a (n),eq α 1 •••αn + ǫ∂ t 1 a (n),1 α 1 •••αn + ǫ∂ j 1 a (n+1),1 jα 1 •••αn + c s 2 n i=1 ǫ∂ α i 1 a (n-1),1 α 1 •••α i-1 α i+1 •••αn = 0 ∼ O(ǫ 2 ) (A.11b) ǫ 2 ∂ t 2 a (n),1 α 1 •••αn = 0 ∼ O(ǫ 3 ) (A.11c)
Summing all the orders together and using the expression of the equilibrium coefficients in equation ( 21), the following relationships are achieved

∂ t ρ + ∂ α (ρu α ) = 0, (A.12a) ∂ t (ρu α ) + ∂ β (ρu α u β + pδ αβ ) = -∂ β a (2),1 αβ , (A.12b) ∂ t (ρu α u β + pδ αβ -ρr 0 T 0 δ αβ ) + ∂ γ [ρu γ u α u β + (p -ρr 0 T 0 )(uδ) γαβ ] + c s 2 [∂ α (ρu β ) + ∂ β (ρu α )] + ∂ t a (2),1 αβ + ∂ γ a (3),1 γαβ = - 1 τ a (2),1 αβ . (A.12c)
Equation (A.12a) gives the mass conservation in [START_REF] Casalino | Turbofan broadband noise prediction using the lattice boltzmann method[END_REF]. The next step is to find an explicit expression of a

(2),1

αβ to fulfill the viscous tensor in the momentum equation (A.12b). Using equations (A.12a) and (A.12b), the transport equation for ρu α u β can be written as then leads to recovery of the viscous tensor in equation [START_REF] Feng | A lattice-boltzmann model for low-mach reactive flows[END_REF]. Note that, unlike with mono-species thermal LBM, the second term in equation (A.15) can not be linked to the bulk viscosity using the energy conservation for mono-species ∂ t θ + u γ ∂ γ θ + 2 D θ∂ γ θ = 0, because this energy equation is generally invalid for multi-species flow.

∂ t (ρu α u β ) + ∂ γ (ρu α u β u γ ) + (∂ α p)u β + (∂ β p)u α = -u β ∂ γ a (2),1 αγ -u α ∂ γ a (2),1 βγ . (A.13) Inserting back to (A.12c), δ αβ [∂ t p + ∂ γ (u γ p)]+2pS αβ -u β ∂ γ a (2),1 αγ -u α ∂ γ a (2),1 βγ +∂ t a (2),1 αβ +∂ γ a (3)

AppendixB. Thermochemistry & transport H 2 -air combustion model

Throughout the paper, we consider the 12-step skeletal mechanism for H 2 -air combustion [START_REF] Boivin | An explicit reduced mechanism for H 2 -air combustion[END_REF] derived from the detailed San Diego mechanism [START_REF] Williams | Chemical-kinetic mechanisms for combustion applications[END_REF].

The mechanism, summarized in Tab.B.4, involves eight reacting species, as well as inert N 2 for combustion with air. The associated thermodynamic data was obtained from the San Diego mechanism website [START_REF] Williams | Chemical-kinetic mechanisms for combustion applications[END_REF].

Thermodynamic closure

The thermodynamic properties required for the thermodynamic closure are specified in the form of the classical NASA polynomials [START_REF] Esch | Thermodynamic properties in polynomial form for carbon, hydrogen, nitrogen, and oxygen systems from[END_REF] species k. Computation of the temperature is done through a Newton iterative procedure with the previous time-step temperature as initial condition.

Considering this thermodynamic closure, the thermodynamic equilibrium, required to initialize the premixed flame presented in Fig. 1, is given in Table B.5.

Simplified transport model

The thermal diffusion coefficient, required in ( 12) is defined as

λ = µ Pr N k=1 Y k C p,k , (B.1)
with µ the dynamic viscosity, assumed here to follow temperature-dependent power-law

µ = µ 0 T T 0 β . (B.2)
The diffusion velocities are evaluated through a Fickian approximation A correction velocity [START_REF] Poinsot | Theoretical and numerical combustion[END_REF] is applied to ensure mass conservation, reading

V k,α = -D k ∂X k ∂x α W k W + V c α Y k , (B.
V c α = N k=1 D k ∂X k ∂x α W k W . (B.

1 )Figure 1 :

 11 Figure 1: Freely propagating flame: a) Temperature profile (thick line), H 2 (△), O 2 (▽), H 2 O ( * ), H ( ), HO 2 (•) mass fractions. Fresh gases are in stoichiometric proportion, at 300K and atmospheric pressure. b) Variation with equivalence ratio of the corresponding flame speed. Cantera reference (plain line), and present model (red dashed line).

Figure 2 : 2 -

 22 Figure 2: 2-D Darrieus Landau instability. a) Representation of the expansion with interval δt = 10 -4 s and k = 600m -1 . The front is initially at x ≈4.5mm. The vertical line marks the end of the linear regime depicted in plot b. b) Comparison of analytical (solid line) and simulated growth rate (dashed line) in the linear region delimited ending at the vertical line in plot a.

ρ f ρ b = 7 . 16 .

 716 The pre-exponential factor is adjusted to obtain a flame speed of s L = 0.5m/s and a flame thickness of δ L = 0.25mm. Next, two counter rotational vortices are added to the flow-field upstream of the planar flame as

Figure 3 :

 3 Figure 3: Numerical set-up for the the premixed flame / vortex interaction simulations.

Figure 4

 4 Figure 4 presents the temperature contours as obtained with our model at two successive instants (with the time normalized by δ L /s L ). Accordingly, each case leads to very different topologies of interaction:

Figure 4 :

 4 Figure 4: Contours of the temperature in a premixed flame -vortex interaction. From top to bottom: case A at t= 2.12 and t=2.83; case B at t=1.42 and t=1.80; case C at t=0.90 t=1.09.

  Fig. 8

Figure 5 :

 5 Fig. 9

2 266

 2 the parameters summarized in Tab. 2. Note that use of these parameters leads to a temperature of T F = 853K in the fuel stream and T O = 1200K in the oxidizer stream. Using the fuel stream conditions, and following the × 10 -3 U 0 96 m s -1 p 0 101 325 Pa ρ 0 0.2926 kg m -3 ǫ 0.05 λ 80

Figure 6 :

 6 Figure 6: Doubly periodic shear layer at Re = 3.10 4 . Vorticity contours (magnitude of the z-component) at t = t c . Ntmix contours (left) compared with the LBM contours (right), for the compressible "cold" flow (top) and the "hot" flow (bottom), e.g. including the chemical source term, on a 1024 × 1024 grid.

Figure 7 :

 7 Figure7: Velocity profiles (m s -1 ) at t = t c (top) and t = 2t c (bottom) as obtained with Ntmix (solid line) and LB solver (dashed line) along the diagonal lines superimposed as illustrated in Fig.5, for the cold (left) and hot (right) case.

Figure 8 : 1 . 2 .

 812 Figure 8: Temperature (K), Y H2 , Y H and Y H2O (from left to right) profiles as obtained with Ntmix (solid line) and LB solver (dashed line) along the diagonal identified in Fig. 5, after one (top) and two (bottom) convective times.

Figure 9 :

 9 Figure 9: Temperature (K), Y H2 , Y H and Y H2O (from left to right) profiles as obtained with Ntmix (solid line) and LB solver (dashed line) along the diagonal identified in Fig. 5, after one (top) and two (bottom) convective times.

Figure 10 showsFigure 10 :

 1010 Figure10shows the relative error obtained for the flame velocity of Fig.1as a function of the number of point per flame thickness. Here, we have defined

1 -Figure 11 :

 111 Figure 11: Vorticity contours in the cold (top) and hot (bottom) double shear layer, for different grids. From left to right: 128×128, 256×256, and 512×512.

  In human time, computing one convective time -or 15000 time steps -of the reactive flow on the 1024 × 1024 grids takes about 1.1 hours on a 20-core Ubuntu desktop (2× Intel Xeon(R) CPU E5-2630 v4 2.20GHz).

  only does the present work lift the constant C p and single step chemistry assumptions, but stability and robustness are tested on challenging two-dimensional flows. The final discussion on accuracy, robustness and cpu cost clearly opens the door to Large Eddy Simulation in that context. Future work will include a detailed account on the model numerical viscosity as to allow the simulation of 3D flames in turbulent flows by inclusion of a turbulence model. A*MIDEX project (ANR-11-IDEX-0001-02), funded by the "Investissements d'Avenir".

2 ) 1 ,

 21 [START_REF] Ashna | Extended lattice boltzmann scheme for droplet combustion[END_REF] whereS αβ ≡ 1 2 (∂ β u α + ∂ α u β )is the strain-rate tensor. Omitting the O(ǫ) terms in the above equation, the viscous tensor appears in the momentum equation (A.12b) asa (αβ = -2pτ S -τ δ αβ [∂ t p + ∂ γ (u γ p)] ,(A.[START_REF] Hosseini | Mass-conserving advectiondiffusion lattice boltzmann model for multi-species reacting flows[END_REF] with the second term being negligible in the low Mach approximation.

  [START_REF] Gleason | Comparison of computational simulation of automotive spinning wheel flow field with full width moving belt wind tunnel results[END_REF] with X k the k-th species mole fraction and D k its diffusion coefficient, which is determined via the component specific Schmidt number Sc k from the viscosity

Table 1 :

 1

		80 1.32 2.88
	Case B	20.00 1.16 17.24
	Case C	52.00 1.12 46.43

Effective vortex strength, radius and Karlowitz number for the three configurations (A, B and C) of vortex -flame interaction.

Table 2 :

 2 Physical parameters for the double shear layer test case.

  for each 1 H + O 2 ⇋ OH + O 7 HO 2 + OH → H 2 O + O 2 2 H 2 + O ⇋ OH + H 8 H + OH + M ⇋ H 2 O + M 3 H 2 + OH ⇋ H 2 O + H 9 2 H + M ⇋ H 2 + M 4 H + O 2 + M → HO 2 + M 10 2 HO 2 → H 2 O 2 + O 2 5 HO 2 + H → 2 OH 11 HO 2 + H 2 → H 2 O 2 + H 6 HO 2 + H ⇋ H 2 + O 2 12 H 2 O 2 + M → 2 OH + M Table B.4:The 12-step skeletal mechanism for the combustion of H 2 -air[START_REF] Boivin | An explicit reduced mechanism for H 2 -air combustion[END_REF]. Up-to-date rates are available[START_REF] Williams | Chemical-kinetic mechanisms for combustion applications[END_REF].

Table B .

 B 5: Initial conditions: 1-D domain is initialized with fresh gases corresponding to (0 : L/2) and burnt gases (L/2 : L) Y H 2 O 2 0 1.444 × 10 -10 Y N 2 7.451 × 10 -1 7.451 × 10 -1

	V ariables fresh gases	burnt gases
	T	300 K	2385 K
	p	1 atm	1 atm
	Y H 2 Y H Y O 2 Y OH Y O Y H 2 O Y HO 2	2.852 × 10 -2 1.145 × 10 -3 0 6.983 × 10 -5 2.264 × 10 -1 7.474 × 10 -3 0 5.458 × 10 -3 0 3.838 × 10 -4 0 2.403 × 10 -1 0 1.074 × 10 -6

  5) The various parameters required for the computation of the viscosity through (B.2) as well as species (B.4) and heat (B.1) diffusion properties are those recommended in Cerfacs' database [64] and validated for this mechanism. They are reported in Tab. B.6. O 0.6000 Sc HO 2 0.8000 Sc H 2 O 2 0.8200 Sc N 2 1.0000 Table B.6: Power-law's viscosity coefficients (SI units), Prandtl number and Schmidt numbers for each species of the 12-step mechanism.

	µ 0 Pr	1.8405 × 10 -5 0.7500	β Sc H 2 0.2100 0.6759
	Sc H	0.1400	Sc O 2 0.8000
	Sc OH 0.5300	Sc O 0.5300
	Sc H 2

Acknowledgements

B. Cuenot is gratefully acknowledged for fruitful discussions and providing the NTMIX code (and A. Pestre for helping setting up the Ntmix computations). J. Jacob provided invaluable help with the implementation. P. Haldenwang is gratefully acknowledged for fruitful discussions on thermodiffusive instabilities.

The French Space agency (CNES) is acknowledged for supporting Song Zhao at M2P2. This work was granted access to the HPC resources of Aix-Marseille Université funded by the project Equip@Meso (ANR-10-EQPX-29-01), and from GENCI-TGCC/CINES (Grant 2018-A0032A07679). We also acknowledge support from Labex MEC (ANR-10-LABX-0092) and the