On the Influence of Kinetic Uncertainties on the Accuracy of Numerical Modeling of an Industrial Flameless Furnace Fired With NH3/H2 Blends: A Numerical and Experimental Study

Ammonia/hydrogen-fueled combustion represents a very promising solution for the future energy scenario. This study aims to shed light and understand the behavior of ammonia/hydrogen blends under flameless conditions. A first-of-its-kind experimental campaign was conducted to test fuel flexibility for different ammonia/hydrogen blends in a flameless burner, varying the air injector and the equivalence ratio. NO emissions increased drastically after injecting a small amount of NH3 in pure hydrogen (10% by volume). An optimum trade-off between NOx emission and ammonia slip was found when working sufficiently close to stoichiometric conditions (ϕ = 0.95). In general, a larger air injector (ID25) reduces the emissions, especially at ϕ = 0.8. A well-stirred reactor network with exhaust recirculation was developed exchanging information with computational fluid dynamics (CFD) simulations, to model chemistry in diluted conditions. Such a simplified system was then used in two ways: 1) to explain the experimental trends of NOx emissions varying the ammonia molar fraction within the fuel blend and 2) to perform an uncertainty quantification study. A sensitivity study coupled with latin hypercube sampling (LHS) was used to evaluate the impact of kinetic uncertainties on NOx prediction in a well-stirred reactor network model. The influence of the identified uncertainties was then tested in more complex numerical models, such as Reynolds-averaged Navier–Stokes (RANS) simulations of the furnace. The major over-predictions of existing kinetic scheme was then alleviated significantly, confirming the crucial role of detailed kinetic mechanisms for accurate predictive simulations of NH3/H2 mixtures in flameless regime.


INTRODUCTION
Global warming and air pollution are pushing to identify new solutions that can reduce CO 2 and greenhouse gas (GHG) emissions (UNFCCC, 2015), through the replacement of fossil fuels with renewable energies (RE). Notwithstanding the intermittent nature of RE sources (wind, solar, etc.), continuous, and on demand, power production can be achieved by exploiting RE surplus in power to fuel techniques (Evans et al., 2012). In this contest, hydrogen is one of the most promising energy carrier, in spite of the challenges related to its storage and transportation. In order to bypass these complications, "green" H 2 may be converted to other molecules with different properties. In this regard, ammonia shows a very high potential (Kobayashi et al., 2018), as it has very high H 2 density and it can be liquefied at pressures higher than 9.9 bar at ambient temperature. From an economic perspective, "green ammonia" was also found to be competitive with natural gas-based ammonia plants ("gray ammonia"), as a result of the recent cost reductions in solar and wind technologies (Philibert, 2017). However, comparing to other liquid fuels (i.e., gasoline), ammonia presents some safety issues, being considered a high health hazard, since it is corrosive to the skin, eyes, and lungs. Once it turns to gas, ammonia is colorless with a sharp, penetrating, intensely irritating odor. Takizawa et al. (2008) and Hayakawa et al. (2015) performed laminar flame speed measurements for ammonia in different conditions, showing its limited reactivity, which may lead to combustion instabilities (Valera-Medina et al., 2017a). To overcome these practical issues, NH 3 combustion doped with hydrogen was tested in internal combustion engines (Frigo and Gentili, 2013), rapid compression machines (Pochet et al., 2017), and swirl burners (Valera-Medina et al., 2017b;Valera-Medina et al., 2019) for utilization within a gas turbine environment, where it was found to have significant NOx emissions. Originally, ammonia kinetics was given prominence for its role in both fuel NOx formation and abatement through selective non-catalytic reduction (SNCR) (see Cremer et al., 2000). The pioneering work from Miller and Bowman (1989) in ammonia kinetic modeling was validated via several species measurements and laminar flame speeds (Maclean and Wagner, 1967;Green and Miller, 1981). Lately, Glarborg et al. (2018) reviewed the nitrogen chemistry in combustion, including the NH 3 sub-mechanism. Nowadays, even other widely validated kinetic models for ammonia oxidation are available in the literature, for instance, the ones from Shrestha et al. (2018), Otomo et al. (2018) and Stagni et al. (2020).
Over the last years, particular attention was paid to flameless combustion, introduced by Wünning and Wünning (1997), which is characterized by preheated and diluted reactants, non-visible flame, and uniform distributed temperatures. These features result in a stable, complete, and efficient combustion, together with a strong reduction of pollutants, such as CO, NOx, and soot. Later on, high temperature air combustion (HiTAC) (Katsuki and Hasegawa, 1998), moderate or intense low oxygen dilution (MILD) (Cavaliere and de Joannon, 2004), and colorless distributed combustion (CDC) (Arghode and Gupta, 2010;Arghode et al., 2012) were also investigated.
In spite of the reasonable number of investigations involving natural gas or methane (see Plessing et al., 1998;Özdemir and Peters, 2001;Dally et al., 2004;Szego et al., 2008;Veríssimo et al., 2011;Veríssimo et al., 2013;Huang, 2018), the amount of detailed studies available for furnaces operating under flameless or MILD conditions using non-conventional fuels is scarce and limited to few operating conditions (see Ayoub et al., 2012;Mosca, 2017;Sabia et al., 2019;Ferrarotti et al., 2018). For the first time in the literature, Sorrentino et al. (2019) investigated the feasibility of pure ammonia combustion using a MILD lab-scale cyclonic burner, varying the mixture equivalence ratio, nominal thermal power, and inlet preheating level. They pointed out the need of operating above 1300 K to ensure both combustion stability and low NOx emissions. A minimum in NOx emissions (less than 100 ppm) was observed for equivalence ratios (ER) between 1 and 1.1, with tolerable ammonia slip. Cavaliere and de Joannon (2004) suggested that, in terms of fluid dynamics, the reactive zone in MILD combustion might evolve as in a well-stirred reactor (WSR). The feasibility of such an approximation was recently demonstrated through systematic comparison with experimental data from attached flames (Chen et al., 2017). Medwell et al. (2016) successfully extended such hypothesis to the Jet in Hot Coflow (JHC) flames, focusing on the chemistry-dominated effects, to distinguish among two different resembling regimes in non-premixed flames, that is, MILD combustion and autoignition lifted flames. Conclusively, Zieba et al. (2019) proposed to adopt a series of plug flow reactors (PFRs) with recirculation, while Rocha et al. (2020) adopted a freely propagating flame model with exhaust gas recirculation (EGR).
The main purpose of this work is to present first-of-their-kind experimental data from an industrial flameless burner fired with NH 3 /H 2 blends over a wide range of operating conditions. Since the main challenges associated with ammonia combustion are NOx emissions and ammonia slip, such campaign aimed at finding a configuration with optimal trade-off between the two, and at testing fuel flexibility in the ULB flameless furnace. The second goal consists in using the data gathered in this study to extend the validation of existing numerical models for turbulence/chemistry interaction (partially stirred reactor-PaSR model), which was extensively validated for other fuels, in traditional (Chomiak and Karlsson, 1996;Hua et al., 2005;Sabelnikov and Fureby, 2013;Duwig and Iudiciani, 2014) and flameless/MILD conditions (Ferrarotti et al., 2018;Ferrarotti et al., 2019;Amaduzzi et al., 2020;Li et al., 2021). The last objective is to verify whether the origin of the deviation between experimental measurements and model predictions for NOx emission is attributable to uncertainty in ammonia kinetics. In these regards, a canonical WSR model with exhaust gases recirculation was adopted to isolate the effect of kinetics, and to identify influential reactive steps for NOx formation, through sensitivity analysis. Subsequently, uncertainty quantification (UQ) techniques were adopted to identify the uncertainty bounds for NOx emissions, due to the identified impactful reactions. A conceptual scheme of the followed methodology is shown in Figure 1.
The remainder of the work is organized as follows. Section 2 presents the proposed methodology. Section 3 describes the aforementioned experimental campaign, exploring the feasibility of flameless/MILD combustion for NH 3 /H 2 mixtures, and the related validation of numerical models, together with the range of NOx predictions induced by the uncertainty factor of influential rate constant in the kinetic sub-model. Finally, conclusions are drawn in Section 4.

METHODOLOGY Experimental Facility and Measurement Techniques
In the following subsections, the experimental facility is presented with particular focus on the description of the operating conditions and the main measurement techniques. The furnace was introduced in Ferrarotti et al. (2018), and it is shown in Figure 2. It is composed of a cubic combustion chamber (1,100 × 1,100 × 1,100 mm) insulated with a 200-mm-thick hightemperature ceramic foam layer, resulting in inner dimensions of 700 × 700 × 700 mm. A commercial WS REKUMAT M150 recuperative Flame-FLOX burner (nominal power of 20 kW) is mounted at the bottom of the combustion chamber. The burner has an integrated metallic finned heat exchanger to extract energy from the flue gases and to preheat the combustion air. The test bench shows a configuration similar to industrial furnaces, allowing to vary geometry, injection system, air excess, and load. The fuel is injected via a centrally located nozzle (inner diameter ID 8 mm) and surrounded by a coaxial air jet, whose dimensions can be varied to adjust the air jet entrainment (ID 16-20-25 mm) ( Figure 2 (bottom right)). The unit is equipped with an air cooling system ( Figure 2) consisting of four cooling tubes [outer diameter (OD) of 80 mm], with a length of 630 mm inside the furnace. Varying the air flow allows the combustion chamber to operate at different stable conditions, thus simulating the effect of a variable load. On each vertical wall of the combustion chamber, an opening is available for measurements. One side is equipped with a 110 × 450-mm quartz window allowing OH* chemiluminescence measurements. The other three openings are blocked using the same insulation material, and they host ports for temperature measurements. Temperature inside the chamber is sampled at different locations. Three sides of the furnace are closed with insulated plates. One of them ( Figure 2) is equipped with twelve equally spaced thermocouple ports, at a related distance of 50 mm. In particular, an air-cooled suction pyrometer equipped with a 0.5-mm diameter B-type thermocouple (platinum rhodium-30%/platinum rhodium-6%) is used to measure the in-flame temperature profiles. It works with a Venturi tube connected to a compressed air circuit at a maximum pressure of 6 bar(g). The thermocouple is protected from chemical attack and from radiation heat exchange with the surrounding walls by two concentric sintered alumina shields (Ferrarotti et al., 2018). According to the specifications of the manufacturer, the associated uncertainty is 0.5% of the reading. Finally, the exhaust gas temperature (before the heat exchanger) is given by a shielded N-type thermocouple (Nicrosil/Nisil) positioned on the central plane and shifted 200 mm with respect to the axis, on the bottom wall ( Figure 2). Other K-type thermocouples measure temperature of the main operating parameters, such as fuel, cooling air, combustion air, and exhaust gases after the heat exchanger temperature. An intensified charge-coupled device (ICCD) camera allows to record chemiluminescence images to analyze the topology of the reaction zones. In the current configuration, images were recorded with an ICCD camera (LaVision 1,392 × 1,040 pixels-16 bits) equipped with a UV 78 mm f/3.8 lens and an interferential filter to collect OH* (310 × 10 nm). Considering the size of the window (110 × 450 mm), the accessible area on the symmetry plane of the furnace goes from 40 to 500 mm from the burner exit. A series of N f 150 frames was taken at a global acquisition frequency of 17 Hz. The images post-processing included 1) noise subtraction recorded with the flame extinguished and 2) averaged of the N f frames. The outlet of the combustion chamber (after the heat exchanger) is equipped with a heated sampling probe to allow flue gas temperature and composition measurements, avoiding condensation. A Fourier transform infrared analyzer (HORIBA MEXA-ONE-FT) is used to measure pollutants (NH 3 , NO 2 , and NO). Since oxygen is transparent to IR sources, a paramagnetic analyser is used, after condensing water from the exhaust gases. With Fourier transform infrared spectroscopy, uncertainty depends on the choice of the concentration ranges, and it has different sources, such as zero noise (1% of full scale (FS)), linearity (1% FS), and water interference (1% FS). The current investigation considered an input thermal power of 15 kW and an equivalence ratio varying between 0.8 and 1. The cooling power subtracted via the air cooling system was set to 5.1 kW to ensure exhaust gas temperature of around 950°C before the heat exchanger. Two different air injector sizes (ID 16 and ID 25 mm) were used, varying therefore the injection velocity and the residence time.
Experiments were performed at steady-state conditions, after a warming period of about 3 h, during which the same burner was used in normal flame conditions, acting on the fluid dynamic of the injection. In the following, "NxHy" term represents the fuel mixture of x %vol. of NH 3 and y %vol. of H 2 .

Computational Fluid Dynamics Model
The Reynolds-averaged Navier-Stokes (RANS) simulations were solved with the commercial code Fluent 19.3 by Ansys, Inc. Turbulence chemistry interactions were handled with the PaSR combustion model (Ferrarotti et al., 2018;Ferrarotti et al., 2019). In the latter model, the chemical time scale is retrieved from the species formation rate, while a static formulation for the mixing time scale is computed as follows: The constant C mix was set equal to 0.5 as suggested by Ferrarotti et al., 2018;Ferrarotti et al., 2019). Standard k-ϵ was used to model turbulence, while two kinetic schemes were used for the H 2 /NH 3 chemistry, namely Stagni et al. (2020), composed by 31 species and 203 reactions, and Glarborg et al. (2018), also containing an ammonia sub-mechanism consisting of 31 species and 211 reactions. Radiation was modeled using a discrete ordinate (DO) approach, combined with the weighted-sum-ofgray-gases (WSGG) model. The computational domain considers an angular sector of 45°due to the symmetry of the problem. Indeed, the presence of the window was not considered in the present study, and the related heat loss was added to the cooling loss. A verification study to ensure grid independence was carried out, following directions from previous studies of the same system (Ferrarotti et al., 2018). The resulting grid consists of 216 k cells. Air and fuel flow rates were set according to the operating conditions of the experiments. The fuel blends were assumed to be fed into the furnace at 343 K, while the inlet air temperature was estimated by solving an energy balance across the heat exchanger. In order to model the cooling surfaces a constant negative heat flux condition was imposed, which also incorporates the energy loss by radiation through the window.

Well Stirred Reactor With Exhaust Gas Recirculation
As recommended by Medwell et al. (2016), an adiabatic, nonisothermal well-stirred reactor, with recirculation, was adopted to model the chemistry of a highly diluted and preheated reactive mixture (i.e., MILD-like conditions). Indeed, as flameless combustion is characterized by a slower chemistry, ignition is likely to take place in flame kernels with premixed fuel and oxidizer, and relatively low strain rate, i.e., distant from the inlet. In this study, the use of WSR with recirculation is intended to map pollutants emissions (i.e., NOx) in different operating conditions, and to qualitatively reproduce experimental trends. Also, this simplified model was used for UQ (see Uncertainty Quantification). Figure 1(top) shows a schematic representation of the adopted network. The mixing unit receives three streams in input, namely fuel, air and exhaust gases. The temperature of each stream is equal to that of experiments. It has to be pointed out that the temperature of exhaust gases was assumed to be equal to that measured at the outlet (before the heat exchanger). This was done to avoid modeling of heat losses along the recirculation region, where the mixture may be assumed to be nonreactive (Zieba et al., 2019). As proposed by Medwell et al. (2016), intermediate species were included in the EGR, since they were found to take part in pre-ignition chemistry in MILD regime (Sidey et al. (2014)). The complete list of recirculated species in the exhaust gas is NH 3 , H 2 , O 2 , N 2 , H, O, H 2 O, OH, HO 2 , and NO. The network consists in a MATLAB script involving the perfectly stirred reactor solver in OPENSMOKE++ by Cuoci et al. (2015).
The residence time of the reactor τ res was estimated from CFD simulations (see Figure 1). For each operating conditions, the reactive region was identified by looking at the contour of OH, which is a well-known marker for reactivity in combustion. Then, the residence time of the reactor was calculated as the average mixing time scale given by the PaSR model (τ res τ mix ), in this region. The first reactor in the network is solved by introducing air as oxidizer. Afterward, the outlet composition is retrieved and assumed to be equal to the EGR for the next simulation. The recirculation degree k v was also estimated using data from RANS simulations, considering the following equation: where _ m mix represents the flow rate of the reactants mixed with exhaust gases (Figure 3). The latter was calculated considering the abovementioned reactive zone, identified using OH as marker. Within this region, three planes were defined at different axial locations. The flow rate passing through a clip of positive velocity (i.e., toward the top wall) was then calculated along these planes. _ m mix was defined as the average between the resulting three values. In the mixer, k v moles of EGR are mixed with air and fuel. The resulting mixture is then fed to the WSR solver. This procedure is iteratively repeated until convergence for both the outlet temperature and NO moles fraction is achieved; that is, residuals are lower than a certain threshold. OPENSMOKE++ by Cuoci et al. (2015) also enables the user to perform sensitivity, and rate of production, analyses. These capabilities were used to identify influential reactions to be further investigated with UQ. The analysis was carried out using the recent mechanism from Stagni et al. (2020).

Uncertainty Quantification
The uncertainty related to a reaction rate constant is often expressed in terms of a factor (see Baulch et al., 2005), which is usually defined as follows: f log 10 k max k nom log 10 k nom k min , where k is the kinetic rate constant for a generic reaction, and the subscripts nom, max, and min represent its nominal, maximum, and minimum value, respectively. According to Eq. 3, the kinetic rate constant, for each reactions in a kinetic mechanism, can be addressed as a random variable with estimated uncertainty (i.e., with a factor of 10 f ). Indeed, deviations from the nominal value of k for the abovementioned sensitive reactions (see Well Stirred Reactor With Exhaust Gas Recirculation) might have a strong impact on the model output, NOx in this case. In order to evaluate the combined effect of those reactions, which were found to be influential for NOx formation, a latin hypercube sampling (LHS) (Iman and Conover, 1980;Florian, 1992) was used. By sampling from the uncertainty range of the aforementioned reactions, several alternative kinetic mechanisms could be obtained, together with the related model responses. Thus, this approach enables the quantification of the variability in a model output, due to inherent uncertainties in kinetics, using forward propagation through the model itself. So, for each combination of operating conditions (i.e., equivalence ratio and fuel composition), NOx emissions can be represented as a region, rather than a curve. Performing a LHS study using RANS simulations would have been prohibitive from a computational point of view. For this reason, the WSR network was used instead. By analyzing the system responses to input variations, Arrhenius parameter combinations corresponding to both maximum and minimum of the NO formation distribution could be identified. Figure 4 shows the 500 samples, which were withdrawn from the three-dimensional space associated with the uncertainty ranges of influential reactions in this work. The deviation of the WSR model responses from the nominal mechanism is reported in terms of color: white points are characterized by null deviation, while black indicates maximum deviation. The two parameter combinations corresponding to the maximum negative (i.e., min NOx) and positive (i.e., max NOx) deviations are also reported. These two mechanisms were then used in numerical simulations (i.e., RANS) to quantify the uncertainty in pollutant predictions, due to kinetics, in more complex models.

Experimental Results
This section focuses on the quantification of NOx emissions and ammonia slip observed while burning different ammonia/ hydrogen mixtures using an industrial flameless burner. To the authors' knowledge, this is the first time such analysis is carried out. Different combinations of air injectors (diameters of 25 and 16 mm hereafter referred to as ID25 and ID16), equivalence ratio (0.8-1), and fuel composition were tested. Rich conditions were excluded from the study to ensure a full conversion of the fuel and to minimize the content of ammoniaslip released in atmosphere. Figure 5 shows the intense yellow color typical of ammonia combustion for different NH 3 -H 2 blends, for the case ID25 and ϕ 1. Figures 6, 7 compare averaged experimental temperature profiles extracted at different axial locations and OH* imaging for the N50H50 mixture, varying the air ID and the ER. Using ID16 and ϕ 0.8, the reaction region is located in a region between 110-160 mm from the nozzle, with a maximum temperature of around 1750 K, at z 150 mm. The OH* contour also appears more spread and less intense than the other cases. Indeed, ID16 ensures a very high injection velocity (∼185 m/s), leading to a high strain rate value close to the burner exit. At a certain axial distance, when the latter weakens, ignition occurs, leading to a noticeable lift-off. Keeping the same injector, but reducing the air excess (ϕ 1), thinner reaction layer shifted toward the burner exit is observed. However, for this case, as well as for ID25, the actual maximum temperature is likely to be located below the first available measurement port (z 100 mm). When the ID25 is employed, the OH* region is shifted even more upstream (between 50 and 80 mm) and the temperatures profiles are smoother for the investigated region. Figure 8 (left) shows the normalized NO (at 3% of O 2 ) and ammonia-slip emissions, varying the ammonia molar fraction and the ER for ID25. Differently from nitrogen-free fuels (i.e., methane and hydrogen), when a fuel blend containing ammonia is used, different pathways are involved. With a small amount of NH 3 (10% in volume), NO emissions grow considerably (from 159 to 827 ppm for ϕ 0.8), reaching a peak between 50% and 60% NH 3 of 3,500 ppm. Moreover, results suggest that the stoichiometry has a major impact on NO formation, confirming literature outcomes from Somarathne et al. (2017) and Sorrentino et al. (2019). As expected, the minimum NO emission levels were obtained close to stoichiometric conditions. Under these conditions, NO is less sensitive to the reaction O + NH 2 H + HNO (R31 in Table 1) due to a lower availability of the radical O. HNO is then converted to NO via the reaction HNO + H → NO + H 2 (more details in next Section). Furthermore, going toward ϕ 1, the peak is shifted progressively toward lower ammonia molar fraction up to 10% NH 3 for ϕ 1 (137 ppm). Very low NO emissions (single digit) can be achieved for this last condition (ϕ 1), for a percentage of ammonia above 50%. The stabilization of pure ammonia combustion was not achieved, since extinction occurred above 80% NH 3 , for all the investigated conditions. In the literature, there are example of pure ammonia burning in MILD regime, for instance, Sorrentino et al. (2019) managed to use pure ammonia in a cyclonic burner, under specific conditions. Further investigations will focus on extending the extinction limit preheating ammonia and/or reducing the thermal power to enhance the reactivity and increase the residence time, respectively.
Finally, reaching conditions close to stoichiometry, unburned ammonia might be found in the exhaust gases (NH 3 -slip). At ϕ 1 (Figure 8 (right)), NH 3slip rapidly increases. reaching values about ∼3,000 ppm, while it is almost zero for lean conditions. An optimal window can be found between ϕ 0.95 and ϕ 1.00 with a strong reduction in NO emission (maximum value 400 ppm) as well as low NH 3 -slip. However, it must be pointed out that it is easier to clean the exhaust gases removing ammonia (i.e., by condensation and adsorption) than adopting techniques to abate NO (i.e., DeNOx).
The effect of the air injector ID is shown in Figures 9 for both NO (left) and NH3-slip (right) emissions. A higher air inlet velocity tends to increase NOx emissions as well. This might be explained considering the following: a higher recirculation ratio k v decreases NO since it increases the level of dilution; however, a reduced residence time (ID16) might not guarantee a sufficient time to convert NO into N 2 . An in-deep explanation is offered in the next section. Analogous trends can be observed varying the equivalence ratio.

Uncertainty Analysis
The WSR analysis is not intended to quantitatively predict the experimental data shown in Figures 8, 9, but to provide qualitative information about NO formation in  hydrogen-ammonia mixtures. Figure 10 shows the NO estimations computed with the WSR network, varying the ammonia molar fraction in the fuel, for the two different ER (i.e., ϕ 0.8, 1). These estimations are referred to as mean values. Their variability due to separate influence of the uncertainties in k v , τ res , and kinetics, is also displayed. Remarkably, the simplified reactor network is capable of qualitatively capturing the NOx dependence on the equivalence ratio, as observed experimentally (see Experimental results). In fact, the following conclusions can be withdrawn by looking at the mean model predictions: • at ϕ 1 (see Figure 10) (left), for ID16, a peak is observed in correspondence of the N10H90 mixture, then emissions diminish as ammonia concentration in the fuel raises; • at ϕ 0.8 (see Figure 10) (right), NO emissions increase with respect to stoichiometric conditions; • at ϕ 0.8, a lower NO production can be achieved using a larger air injector (ID25). This is in line with what was found experimentally (see Figure 8). The lower inlet velocity, due to ID25 configuration, reduces the entrainment of exhaust gas (k v ) and increases the residence time in the reactive zone (τ res ). Globally, this results in reduced emissions.
As the adopted k v and τ res values, for each fuel composition, are average quantities extracted from the reactive zone in RANS simulations (see Well stirred reactor with EGR), a sensitivity analysis was performed by multiplying and dividing them by a factor of 2. The model responses are reported in Figure 10, where a higher sensitivity to the recirculation degree than to residence time is observed. Yet, uncertainties in the kinetic model yield a much greater impact on the NO variance, which is also reported in Figure 10. The displayed uncertainty bands were estimated by performing a LHS (see Uncertainty Quantification) of 500 samples from the cube shaped by three coordinates, namely, the pre-exponential factors of reactions R39, R80, and R85, which are uncertain random variables of the kinetic model. These reactions were selected using sensitivity analysis in OPENSMOKE++. Figure 11 reports the most influential reactions for the mixtures N25H75  and N50H50, at different equivalence ratio, namely, 0.8 and 1.0. In particular, reactions involved in the hydrogen core mechanism were discarded a priori. The sub-mechanism for thermal NOx, involving R89, R90, and R91, was found to be particularly sensitive, especially for stoichiometric mixtures. However, it was not considered for the UQ, as the authors assumed it to be well-characterized (see Baulch et al., 1994). Interestingly, R85 is the most important reaction with negative sensitivity coefficient. It is also important to point out that NO is very sensitive to R31, which forms HNO, for lean conditions, where more oxygen is available. Indeed, HNO is then converted to NO via HNO + H NO + H 2 . This may explain why lean conditions produce higher NO emissions. However, R31 and R76 were found to have high sensitivity coefficient for temperature and consequently ruled out from the UQ study. Regarding other reactions, R80 converts NH into HNO and impact positively the sensitivity, and R39 is only sensitive for higher ammonia content in the fuel than 25%. Indeed, this reaction affects more and more the formation of NO as NH 2 , and NH production increase, due to higher availability of NH 3 as well as lower H and OH radicals' concentrations. As a consequence, only R80, R85, and R39 were selected for the UQ study. Table 1 reports the adopted uncertainty factors for reactions in Figure 11, based on literature information. In addition, a flux analysis was performed for both N25H75 and N50H50 mixtures (see Figure 12), to explain the NO emissions trends.
Ammonia reactivity proceeds along NH 2 → NH → N, and NO is part of its oxidation. In fact, NH 2 mainly forms NH and HNO in R31, which gives NO. The NH intermediate has a crucial role, as not only leads to NO through HNO in R80, but also reacts with it in R85 to form N 2 O, which is almost completely converted to N 2 in the termination step H + N 2 O N 2 + OH. In addition, NH is converted to N, which exhibits an analogous behavior, that is, it produces NO in R89 and R90, but also reacts with it in R91, again as a termination step. Up to N10H90, hydrogen concentration is so high that the radical pool is extremely rich in H and OH, prompting HNO production (R80) and its next conversion to NO, determining an emissions peak. The latter peak is even more pronounced at ϕ 0.8 because of the higher availability of local O radical, prompting the HNO production via R31. As the NH 3 percentage in the fuel increases (i.e., at N50H50), these pathways weaken, and R39 starts competing. The latter reaction offers an alternative path to NH, namely, N 2 H 2 → NNH → N 2 , which tends to reduce the NO formation by subtracting NH and NH 2 from the pool of reactants. So, R39 is part of the reason why richer fuel mixtures in NH 3 show decreased NO emissions. In fact, this reaction shows a positive impact of NO sensitivity in Figures 11 as it competes with NH → NO →N 2 O → N 2 , which represents the preferential way for the system to reduce NO emissions. Finally, one possible explanation for the existence of a shifted peak at ϕ 0.8 (see Figure 8) is the higher oxygen content, which pushes the NO production through HNO in R31 for richer mixtures with respect to stoichiometric conditions, delaying the effect of R39. As previously reported in Figure 4, the LHS study suggests that the maximum NO production in Figure 10 is located nearby the maximum of R85, minimum of R39, and R80. The opposite is true for the minimum NO formation in Figure 10, in agreement with the sensitivity and flux analysis.

Uncertainty Propagation in the CFD Simulation
This section is focused on the influence of existing uncertainties in detailed kinetic mechanisms on NO predictions with more complex numerical models, that is, RANS simulations with PaSR sub-model for turbulence/ chemistry interaction. To do so, the two mechanisms that were found to give maximum and minimum NO output values from the previously described LHS study were tested. A first validation for temperature profiles is presented for case N50H50 ID16 ϕ 1. Here, a value of 0.5 was employed for C mix in the PaSR model, together with standard k-ϵ for turbulence. Figure 13 shows the comparison between measured and predicted temperature profiles along the axis (a) and at different axial positions (b, c, d). The results from two different kinetic mechanisms, developed for ammonia combustion, that is, Stagni et al. (2020) and Glarborg et al. (2018), are also reported. Looking at Figure 7, the reaction region (maximum of OH* counts) is located between 110 and 160 mm from the burner exit. However, the two models predict a late ignition compared to experimental data, as they under-predict the temperature peak at 100 mm away from the inlet. The abovementioned under-prediction corresponds to a 2% and 1% relative error for Stagni et al. (2020) and Glarborg et al. (2018) models, respectively. Even though both mechanisms performed well on temperature profiles, strong differences were detected for pollutant emissions estimates. Regarding NOx, a pronounced overestimation was observed using both models (see Figure 14). In particular, for ϕ 1 (Figure 14 left), Glarborg et al. (2018) predict much higher values, that is, 13580 and 5300 ppm at N25H75 and N50H50, respectively, vs. 9,473 and 2,851 ppm for the mechanism from Stagni et al. (2020). Anyway, even using the latter model, predictions are still far from the experimental data (460 and 160 ppm for N25H75 and N50H50, respectively). For ϕ 0.8 (Figure 14 right), a better agreement with experimental data is observed, even though the model under-predicted the NO emissions for N50H50.
In order to verify the major role of an accurate kinetic submodel in this chemistry-controlled regime (flameless), results from the previously performed LHS study on the WSR network were tested in CFD simulations. Again, this means performing additional simulations using the set of kinetic parameters corresponding to the maximum and minimum NO emissions on the LHS chart ( Figure 4). Figure 14 also shows the uncertainty propagation associated with R39-80-85 on the NO emissions. At ϕ 1, the lower band moves toward the experimental values, allowing a massive reduction compared to the original model from Stagni et al. (2020). This is true especially for N50H50, where NO emissions decrease from 2,851 to 539 ppm. On the contrary predicted, temperature profiles in the furnace (not shown here) remain almost constant, meaning that the effect of kinetic of the three reactions is relevant only for NO. Much better results were achieved at ϕ 0.8, where the uncertainty bounds almost intersect the experimental data region for N10H90 and contains it for N25H75 and N50H50.

CONCLUSION
A novel, interdisciplinary approach to study ammonia/hydrogen blends in flameless combustion was presented in this work. To the authors' knowledge, these mixtures were experimentally investigated for the first time in such regime. The main findings of this work are summarized in the following: • An experimental campaign was performed to confirm the enhanced fuel flexibility of the ULB furnace fired with ammonia/hydrogen blends. Operating configurations in terms of trade-off between NOx emissions and ammonia slip could be identified. In fact, ammonia slip emissions are negligible in lean conditions, while they become relevant close to ϕ 1. The optimal working point was identified for all fuel mixtures at ϕ 0.95, which allows to reduce NO emissions with respect to leaner conditions, while keeping NH 3 -slip below 10 ppm. Additionally, a peak in NO emissions was observed for both tested injectors (ID16 and ID25) in stoichiometric conditions, for the mixture N10H90. Then, emissions decreased up to extinction, which occurred for a volume content of ammonia above 80%. In general, for lean cases (ϕ 0.8), the NO production increased, and a shifted peak toward higher ammonia content in the fuel (N60H40) was observed. Using a bigger air injector (ID25) also help controlling pollutants emissions, as the associated residence time increasing enhances NO conversion to N 2 . • Subsequently, these high-fidelity data were used to extend the validation of existing turbulence/chemistry interaction sub-model (PaSR), for non-conventional fuels and conditions. The agreement between temperature measurements and estimations was found to be satisfactory (see Figure 13), and nearly insensitive to the adopted kinetics, that is, Stagni et al. (2020) and Glarborg et al. (2018). However, in terms of NOx emissions, substantial differences between predictions with the two mechanisms were observed (see Figure 14 (left)). In light of this pronounced difference, the error in pollutant emissions can be partially attributed to the uncertainties intrinsically embedded in the kinetic submodel. Finally, in spite of the large over-prediction at ϕ 1, the mechanism from Stagni et al. (2020) was found to be more accurate. • A WSR network was designed, and fed with boundary conditions from experiments and CFD simulations for the residence time (τ res ), and the averaged recirculation degree (k v ). The WSR allowed isolating the effect of chemistry on NOx production. Sensitivity and rate of production analyses could be performed to identify influential reactions for NO production and consumption, i.e. R39-80-85 in Table 1. An LHS method was adopted to propagate the uncertainty of the abovementioned reactions to NO production/ emissions. Two kinetic mechanisms were developed from Stagni et al. (2020), representing the maximum and the minimum of the NO distribution, reported for both the WSR and CFD simulations, in Figures 10, 14, respectively.
In conclusion, a significant portion of the discrepancy between the pollutant emissions from simulations and experimental observations can be associated to inherent uncertainties in recent kinetic mechanisms for ammonia/hydrogen combustion. In fact, the latter were found to be reliable for temperature prediction, while over-predicting NO emission significantly (see Figure 14). Thus, this work indicates the need to improve existing NH 3 /H 2 models, especially in diluted conditions. In particular, specific reactions ( Table 1) need better characterization, to improve models of practical systems. This is particularly true for stoichiometric conditions, where the discrepancies between numerical model predictions and experiments were found to be most significant.
Future experimental campaigns will focus on reducing the emissions of pollutants, namely, NH 3 and NOx. In this perspective, the effect of fuel preheating and water dilution will be investigated (Ariemma et al., 2020), still need to be characterized. Using a larger injector (ID > 25 mm) could further increase the residence time, thus reducing NO emissions. Finally, a possible option could be to fire the furnace in slightly fuel rich conditions in combination with ammonia slip capture at the furnace outlet. The latter approach would reduce the combustion efficiency though.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article, and further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
All the authors contributed to the conception and design of the work. In addition, MF performed the experimental and numerical analyses, collected and interpreted the data, and drafted the paper. AB performed the kinetic and UQ analysis and drafted the paper. RA contributed to the experimental and numerical analyses. AP supervised the work, provided feedback about the methods and results, and contributed to the writing (review and editing).

FUNDING
This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement no. 714605. The first, second, and third authors wish to thank the Fonds de la Recherche Scientifique (FNRS) Belgium for financing their research.