# Bayesian calibration of a finite-rate nitridation model from molecular beam and plasma wind tunnel experiments

Michele Capriati, Anabel del Val, Thomas E Schwatzentruber, Timothy K Minton, Pietro Marco Congedo, Thierry E Magin

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Michele Capriati, Anabel del Val, Thomas E Schwatzentruber, Timothy K Minton, Pietro Marco Congedo, et al.. Bayesian calibration of a finite-rate nitridation model from molecular beam and plasma wind tunnel experiments. Aerospace Europe Conference 2023 – Joint 10th EUCASS – 9th CEAS Cp,fere,ce, Jul 2023, Lausanne, Switzerland. 10.13009/EUCASS2023-897. hal-04381531

HAL Id: hal-04381531

https://inria.hal.science/hal-04381531v1

Submitted on 9 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Bayesian calibration of a finite-rate nitridation model from molecular beam and plasma wind tunnel experiments

Michele Capriati<fnref n="1" />,<fnref n="2" />,†, Anabel del Val<fnref n="2" />, Thomas E. Schwartzentruber<fnref n="3" />,
Timothy K. Minton<fnref n="4" />, Pietro M. Congedo<fnref n="1" /> and Thierry E. Magin<fnref n="2" />,<fnref n="5" />

## Abstract

The accurate modeling of gas-surface interaction phenomena is crucial for predicting the heat flux and the mass loss experienced by atmospheric entry bodies. Gas-surface interactions refer to the phenomena occurring between the reacting gas and the material. An important part of the modeling is the description of the surface chemical reactions. We propose to stochastically calibrate the rates of the elementary reactions between a nitrogen gas and a carbon surface by using molecular beam data. The parameters' joint posterior is then propagated through a CFD model to reproduce plasma wind tunnel experiments. The predictive quantities exhibit good agreement with the experimental counterparts.

## 1. Introduction

Gas-Surface Interaction (GSI) describes the set of phenomena occurring when a chemically reacting gas interacts with a solid surface. For example, vehicles traveling at hypersonic speeds, such as atmospheric entry spacecraft, experience important energy transfer from the gas to their surfaces. Furthermore, the highly excited post-shock gas promotes the development of a chemically reacting Boundary Layer (BL) that interacts with the surface. Vehicles are thus protected by using Thermal Protection Materials (TPMs) for a safe landing. Two different classes exist: catalytic and ablative. Catalytic TPMs dissipate energy mostly by radiation. However, the exothermic nature of the catalytic reactions increases the heat flux experienced by the vehicle. On the other hand, ablative TPMs dissipate energy mostly by degradation and they experience mass loss. In both cases, the GSI phenomena need to be accurately modeled to predict the in-flight heat flux and mass loss.

One crucial aspect of the GSI modeling is the description of the surface chemical reactions. From an experimental standpoint, Inductively-Coupled Plasma (ICP) facilities<fnref n="4" />,¹¹ and arc jets<fnref n="1" />,⁹,²⁹ are widely used for material characterization, as they can duplicate the reacting BL in front of a hypersonic vehicle. Nevertheless, along with the GSI phenomena that one is interested in characterizing, phenomena such as gas chemistry and gas diffusion are also present. The assumptions contained in the modeling of these processes directly affect the characterization of the surface reactions.⁶ On the other hand, molecular beam-surface scattering experiments²¹⁻²³ isolate the GSI phenomena by bombarding a sample placed in a high vacuum chamber with a beam of molecules or atoms at energies thought to be relevant to hypersonic flight. While these experiments manage to isolate GSI phenomena, we cannot fully trust that they are representative of hypersonic flight conditions.

From a modeling point of view, surface chemical reactions can be computed using phenomenological approaches<fnref n="2" />,<fnref n="3" />,¹²,²⁶ or Finite-Rate Chemistry (FRC) models.¹⁸,²⁷,²⁸ Phenomenological approaches are widely used in hypersonics due to their straightforward implementation. Such approaches prescribe an efficiency to each macroscopic reaction, while FRC models describe the elementary surface reactions.

Recently, Prata et al. derived an Air Carbon Ablation (ACA) model,²⁸ which accurately captures the trends of the reaction efficiencies observed both in molecular beam and in plasma wind tunnel experiments. The model includes 20 elementary reactions involving both oxygen and nitrogen. Further insights into the surface chemical processes are

<footnote n="1">¹Inria, Centre de Mathematiques Appliquees, Ecole Polytechnique, IPP, Route de Saclay, 91120 Palaiseau, France</footnote>

<footnote n="2">²von Karman Institute for Fluid Dynamics, Chaussee de Waterloo 72, 1640 Rhode-Saint-Genese, Belgium</footnote>

<footnote n="3">³Aerospace Engineering and Mechanics Department, University of Minnesota,
110 Union Street SE, 55455 Minneapolis, Minnesota (USA)</footnote>

<footnote n="4">⁴Aerospace Engineering Sciences, University of Colorado, 3775 Discovery Dr, 80303, Boulder, Colorado (USA)</footnote>

<footnote n="5">⁵Aero-Thermo-Mechanics Laboratory, Ecole Polytechnique de Bruxelles, Universite Libre de Bruxelles,
91050 Brussels, Belgium
michele.capriati@vki.ac.be</footnote>

<footnote>†Corresponding author</footnote>
given by providing a detailed study of the competitiveness of the two gases in occupying the free sites on the surface. In a previous study⁵, the nitridation reactions of the ACA model were employed to numerically simulate an experimental campaign carried out at the von Karman Institute (VKI) to study carbon ablation in a nitrogen flow.¹³ The ACA model was calibrated both on molecular beam and on Plasmatron data, but strong assumptions were made on modeling the BL developing in front of the sample in Plasmatron experiments. Furthermore, the model was constructed deterministically, without accounting for experimental and parametric uncertainties. In turn, some observed features (e.g., surface recession rates) were not fully in agreement with the experimental values.

This issue underlined the necessity of including the experimental as well as parametric uncertainties during the calibration of the model parameters. In this context, it was shown that Bayesian methods offer a suitable framework for calibrating both catalytic and ablative phenomenological surface reactions.⁷,⁸

In this work, we propose to calibrate in a stochastic fashion the reactions involving nitrogen of the ACA model.²⁸ The reaction rates, along with their uncertainties, are characterized by means of highly informative molecular beam data. The predictive capability of the calibrated model at higher pressure is assessed by using the model to numerically reproduce Plasmatron experiments. They are simulated with a CFD solver that includes an ablative Boundary Condition (BC) to properly describe the BL in front of the test sample. The molecular beam data are explained through the same analytical model proposed by Prata et al.²⁸ The parameters' posterior distribution is obtained by constructing a Markov Chain Monte Carlo (MCMC) using the molecular beam data as observed quantities. The obtained distributions are then propagated through the CFD to obtain the posterior uncertainty on the blowing mass rate predictions, which are compared to the experimental counterparts.

## 2. Experimental Data

In this section, we first introduce a molecular beam-surface scattering experiment involving a nitrogen beam over a carbon surface. As discussed in the introduction, the nature of the experiment makes it suitable for calibrating the elementary reactions of an FRC model. Nevertheless, the low pressure characterizing the experiment does not guarantee that the resulting model is accurate under re-entry conditions. For this reason, the predictive capability of the model is assessed using plasma wind tunnel data, which consists of a dissociated nitrogen flow impinging on a graphite sample.

### 2.1 Reaction probabilities from Molecular Beam-Surface Scattering

The data used to calibrate the model come from a molecular beam-surface scattering experiment performed by Murray et al.²³ to infer the reaction efficiencies of oxygen and nitrogen on an ablative surface. Being interested in calibrating only reactions involving nitrogen, we will focus on the results obtained using a nitrogen beam.

The nitrogen beam is produced using a 95 mbar mixture (2.5% N₂ in He) discharged by a high-pressure radio frequency source and expanded through a 0.48 mm diameter nozzle. It is then collimated by means of a skimmer with a diameter of 0.8 mm, followed by a 0.4 × 0.4 mm² aperture. The resulting beam has an average velocity of ≈ 2000 m/s and it is composed, in molar fractions, of 18% atomic nitrogen and 82% of molecular nitrogen. The beam is directed over a heated sample, placed in a high-vacuum room, where pumps guarantee pressures of the order of 10 × 10⁻⁷ torr. In addition, the temperature of the sample is measured with a calorimeter.

A rotatable mass spectrometer detects the number density distribution of the scattered products as a function of time, N(t), at different final angles, θf, corresponding to an incident angle, θi. The relative curves are generally referred to as Time-Of-Flight (TOF).

The data were collected at increasing temperatures (from 800 to 1873 K). Before each test, the sample was left to anneal at 1873K for about 1 h. Results showed that the surface undergoes only recombination and nitridation reactions. Their reaction probabilities are computed as:

$$ \gamma^r = \frac{\text{flux of } r \text{ product}}{\text{flux of N atoms onto the surface}} = \frac{n f(\text{product}_r)}{f(N + 2\text{N}_2 + \text{CN})}, \quad \text{with } r = [\text{N}_2, \text{CN}]. \quad (1) $$

Symbol n refers to the number of nitrogen atoms in the products.

### 2.2 Surface recession from high-speed imaging in the VKI Plasmatron

Data from plasma wind tunnel experiments are compared with the model predictions to show that it can be representative of atmospheric entry conditions. We briefly describe the experimental campaign performed by Helber et al.¹³
to study the temperature-dependency of the nitridation reaction efficiency on carbon surfaces in the VKI Plasmatron facility.⁴

The experiment consisted of a plasma nitrogen flow over a 25 mm radius hemispherical graphite sample. The setup included a digital camera for the determination of the stagnation point recession rate, a 2-color pyrometer for measuring the surface temperature of the sample, and a spectrometer to retrieve both the gas temperature and the CN concentration. Seven runs were performed. The electrical power supplied to the generator was progressively increased to obtain a higher free-stream enthalpy, and, in turn, sample's surface temperature. Three of the seven runs exhibit a strong thermal non-equilibrium and will not be considered in this study. The BL edge conditions were obtained through a coupled numerical procedure, detailed in Helber et al.,¹³ using the measured dynamic pressure and the heat flux. An overview of the conditions of the four runs considered can be found in Table 1.

<table><caption>Table 1: Overview of Plasmatron test conditions for a nitrogen flow at 1500 Pa over a 25 mm radius hemisphere graphite sample: dynamic pressure $p_d$, generator power $P_{el}$, sample mean temperature $T_w$, mass blowing rate $\dot{m}$, BL edge conditions (density $\rho_e$, temperature $T_e$, and velocity $u_e$).</caption><thead><tr><th>CASE</th><th><i>p</i><sub>d</sub><br>[Pa]</th><th><i>P</i><sub>el</sub><br>[kW]</th><th><i>T</i><sub>w</sub><br>[K]</th><th><i>ṁ</i><br>[g/m<sup>2</sup>s]</th><th><i>ρ</i><sub>e</sub><br>[g/m<sup>3</sup>]</th><th><i>T</i><sub>e</sub><br>[K]</th><th><i>u</i><sub>e</sub><br>[m/s]</th></tr></thead><tbody><tr><td>G4</td><td>231</td><td>280</td><td>2225</td><td>2.49 ± 0.91 (95%)</td><td>0.236</td><td>10005</td><td>554</td></tr><tr><td>G5</td><td>268</td><td>330</td><td>2410</td><td>2.89 ± 0.97 (95%)</td><td>0.225</td><td>10280</td><td>562</td></tr><tr><td>G6</td><td>312</td><td>370</td><td>2535</td><td>4.41 ± 0.80 (95%)</td><td>0.184</td><td>11040</td><td>846</td></tr><tr><td>G7</td><td>330</td><td>390</td><td>2575</td><td>4.56 ± 0.70 (95%)</td><td>0.189</td><td>10970</td><td>859</td></tr></tbody></table>

In the work of del Val et al.⁷, these data were used to infer the phenomenological nitridation and recombination reaction efficiencies. The obtained marginals and joint distributions of the reaction efficiencies are compared to the ones obtained in this study by post-processing the calibrated FRC model predictions.

### 3. Theoretical Modeling

In this section, we introduce the reactions involving nitrogen of the ACA FRC model²⁸. We intend to calibrate the model using molecular beam experiments. They are numerically described using a 0D reactor simulator, where the FRC model is embedded. The FRC model is also contained in the ablative BC of the CFD solver employed to simulate the Plasmatron experiments.

#### 3.1 Finite-rate nitridation model

The ACA FRC model describes the elementary reactions occurring on a carbon surface, such as adsorption, desorption, oxidation, nitridation, and recombination.²⁸ It has the advantage of capturing the pressure dependence of the chemical reactions and the surface coverage. The subset of nitrogen reactions we are interested in, as well as the corresponding rates, are reported in Table 2.

<table><caption>Table 2: Subset of the FRC reactions involving nitrogen in the ACA model.</caption><thead><tr><th>Tag</th><th>Reaction</th><th>Reaction Rate</th><th>Rate Coefficient</th></tr></thead><tbody><tr><td>1</td><td>N + s → N<sub>s</sub></td><td>k<sub>1</sub>[N][s]</td><td><sup>F<sub>N</sub></sup>&frasl;<sub>B</sub> exp <sup>-E<sub>1</sub></sup>&frasl;<sub>T<sub>s</sub></sub></td></tr><tr><td>2</td><td>N<sub>s</sub> → N + s</td><td>k<sub>2</sub>[N<sub>s</sub>]</td><td><sup>2πm<sub>N</sub>k<sub>2</sub><sup>2</sup>T<sub>s</sub><sup>2</sup></sup>&frasl;<sub>N<sub>s</sub>Bh<sup>3</sup></sub> exp <sup>-73971</sup>&frasl;<sub>T<sub>s</sub></sub></td></tr><tr><td>3</td><td>N + N<sub>s</sub> + C<sub>b</sub> → CN + N + s</td><td>k<sub>3</sub>[N][N<sub>s</sub>]</td><td><sup>F<sub>N</sub></sup>&frasl;<sub>B</sub> X<sub>3</sub> exp <sup>-E<sub>3</sub></sup>&frasl;<sub>T</sub></td></tr><tr><td>4</td><td>N + N<sub>s</sub> → N<sub>2</sub> + s</td><td>k<sub>4</sub>[N][N<sub>s</sub>]</td><td><sup>F<sub>N</sub></sup>&frasl;<sub>B</sub> X<sub>4</sub> exp <sup>-E<sub>4</sub></sup>&frasl;<sub>T<sub>s</sub></sub></td></tr><tr><td>5</td><td>N<sub>s</sub> + N<sub>s</sub> → N<sub>2</sub> + 2s</td><td>k<sub>5</sub>[N<sub>s</sub>][N<sub>s</sub>]</td><td>√<sup>N<sub>A</sub></sup>&frasl;<sub>B</sub> F<sub>N,2D</sub>X<sub>5</sub> exp <sup>-E<sub>5</sub></sup>&frasl;<sub>T<sub>s</sub></sub></td></tr><tr><td>6</td><td>N<sub>s</sub> + C<sub>b</sub> → CN + s</td><td>k<sub>6</sub>[N<sub>s</sub>]</td><td>X<sub>6</sub> exp <sup>-E<sub>6</sub></sup>&frasl;<sub>T<sub>s</sub></sub></td></tr></tbody></table>

Symbol s denotes a free site on the surface that can be occupied by a gas-phase element, the subscript (·)$_s$ indicates absorbed species, while the subscript (·)$_b$ surface elements, [s] is the density of the free sites, $k_i$ the reaction rate
coefficient, B the total active site density, T_s the surface temperature, h the Planck constant, k_b the Boltzmann constant, N_A the Avogadro's number, m_N the mass of the nitrogen, F_N = $\sqrt{k_bT}/2\pi m_N$ a quarter of the mean thermal speed of atomic nitrogen, and F_{N,2D} = $\sqrt{\pi k_b T}/2m_N$ the mean thermal speed of the mobile adsorbed nitrogen on the surface, where T is the temperature of the gas. Symbols E_i and X_i are introduced to designate the parameters of the model that we intend to calibrate. The values proposed by Prata et al. are reported in Table 3.

The ACA model was deterministically inferred both from low-pressure and high-pressure experiments and captures the main observed trends. It includes both atomic nitrogen adsorption and desorption (reactions 1 and 2). The adsorption process was modeled as an activated process to reproduce the sharp drop of the N₂ efficiency observed at low temperatures in the molecular beam experiments. The rate coefficient of the desorption process was computed based on the transition state theory. A desorption energy equal to the bound energy of a double bound between C and N was chosen to describe the non-decreasing efficiency trend observed at high temperatures in the molecular beam experiments. Two nitridation reactions were included in the model to capture the significantly different experimental efficiencies at low and high pressure: a gas-phase dependent reaction (reaction 3) allows for capturing the high-pressure behavior, while the gas-independent reaction (reaction 6) for describing the behavior at low-pressure. Lastly, both Eley-Rideal (ER, reaction 4) and Langmuir-Hinshelwood (LH, reaction 5) mechanisms were included in the model. Although molecular beam experiments only showed evidence of LH recombination, the ER recombination was necessary to describe the recombination over a wide range of pressures²⁸.

## 3.2 0D model

To describe the molecular beam experiment, we used the same 0D reactor simulator as the one proposed by Prata et al.²⁸. It takes as input the partial pressure of the nitrogen and its temperature. A value of T = 1000 K and p = 0.024 Pa allows us to reproduce a flux of 5 × 10²⁰ atom m⁻²s⁻¹, matching the flux of the molecular beam experiment. The second input is the surface temperature, which was set according to the experiment.

Considering the mechanism in Table 2, the set of equations to solve is:

$$ \frac{d}{dt}[N] = -k_1[N][s] + k_2[N_s] - k_4[N][N_s], \qquad (2) $$

$$ \frac{d}{dt}[N_2] = k_4[N][N_s] + k_5[N_s]^2, \qquad (3) $$

$$ \frac{d}{dt}[CN] = k_3[N][N_s] + k_6[N_s], \qquad (4) $$

$$ \frac{d}{dt}[N_s] = k_1[N][s] - k_2[N_s] - k_3[N][N_s] - k_4[N][N_s] - 2k_5[N_s]^2 - k_6[N_s]. \qquad (5) $$

The system is closed with the conservation of active sites

$$ B = N_s + s, \qquad (6) $$

and the steady-state assumption

$$ \frac{d}{dt}[N_s] = 0. \qquad (7) $$

The efficiencies of the macroscopic recombination (N + N → N₂) and nitridation (N + C_b → CN) reactions are then computed as:

$$ \gamma_{N_2} = \frac{2(d[N_2]/dt)}{N_N}, \qquad (8) $$

$$ \gamma_{CN} = \frac{(d[CN]/dt)}{N_N}, \qquad (9) $$

where $N_N = n_N \sqrt{k_bT/(2\pi m_N)}$ is the number flux of particles of nitrogen impinging the surface, and $n_N$ the number density of the nitrogen. Their values are compared to their experimental counterparts.
### 3.3 Stagnation line model

The BL developing in front of the sample in the VKI Plasmatron is numerically simulated by means of an in-house software that solves the chemically reacting Navier-Stokes equations with an ablative BC. The BL edge conditions, as well as the surface temperature, are imposed as BC. Their values are reported in Table 2. The output of the CFD computations is the mass blowing rate, which is compared to the experimental counterpart for each case.

The governing equations describing the problem of interest are:

$$ \frac{\partial \rho_i}{\partial t} + \nabla \cdot (\rho_i \mathbf{u} + \mathbf{j}_i) = \dot{\omega}_i, \quad \forall i \in [1, n_s], \qquad (10) $$

$$ \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u} + p \bar{\bar{I}} + \bar{\bar{\tau}}) = \mathbf{0}, \qquad (11) $$

$$ \frac{\partial \rho E}{\partial t} + \nabla \cdot (\rho \mathbf{u} H + \bar{\bar{\tau}} \cdot \mathbf{u} + \mathbf{q}) = 0. \qquad (12) $$

Equation 10 represents the conservation of mass of the species i: vector u stands for the mass-averaged mixture velocity, $\rho_i$ the partial density of species i, $\mathbf{j}_i$ its diffusive fluxes, $\dot{\omega}_i$ its chemical production/destruction rate, and $n_s$ the number of species in the mixture. Equation 11 is the conservation of momentum: $\rho$ indicates the mixture density, $p$ the thermodynamic pressure, and $\bar{\bar{\tau}}$ the viscous stress tensor. Finally, equation 12 is the conservation of energy: symbol $E$ denotes the total energy, $H = E + p/\rho$ the total enthalpy, and $\mathbf{q}$ the total heat flux.

Thermodynamics properties are obtained using the NASA polynomials.¹⁹ Diffusion mass fluxes are computed using the generalized Stefan-Maxwell equations, while viscosity and thermal conductivity are obtained based on the Chapman-Enskog perturbative solution of the Boltzmann equation.¹⁵,¹⁶ Chemical production rates $\dot{\omega}_i$ for species i are computed using finite-rate chemistry, with reaction rate coefficients taken from Olynick et al.²⁵

When dealing with axisymmetric problems, a Dimensionally Reduced Navier-Stokes Equations (DRNSE) formulation¹⁴ can be used to simulate the flowfield along the stagnation line. The system is written in spherical coordinates $(r, \theta, \phi)$, then, exploiting the axial symmetry, $u_\phi = \partial/\partial\phi = 0$ and letting $\theta \to 0$, the DRNSE formulation is obtained. Munafò²⁰ developed a software for the DRNSE equations based on a finite volume discretization method, together with an implicit local time-stepping technique to reach the steady state. Roe scheme is used to compute the numerical fluxes³⁰. The code is interfaced with MUTATION⁺⁺³,³¹, which provides the closure to the properties and to the ablative BC.

The ablative response of the surface can be modeled by solving mass/energy balances on the infinitesimal volume containing both the gas and the surface. When the surface temperature is known, one only needs to solve a Surface Mass Balance (SMB): surface chemistry generates a species gradient at the surface, driving diffusion. Assuming the recession velocity to be negligible compared to the blowing one ($u_g$), the SMB for each species i reads:

$$ \dot{\omega}_i = (\rho_i \mathbf{u}_g + \mathbf{j}_i) \cdot \mathbf{n}, \quad \forall i \in [1, n_s]. \qquad (13) $$

Versor $\mathbf{n}$ is the normal to the surface. When reactions involve material removal, the material recedes, and the products are advected out from the surface (blowing effect). The mass blowing rate and the blowing velocity are:

$$ \dot{m} = \sum_{i=1}^{n_s} \dot{\omega}_i, \qquad (14) $$

$$ \mathbf{u}_g = \frac{\dot{m}}{\rho}. \qquad (15) $$

The system 13 is solved by means of the GSI module of MUTATION⁺⁺³ where the chemical source terms are computed using the ACA model.

### 4. Stochastic inverse problem

In this section, we describe the stochastic methodology used to calibrate the reactions involving nitrogen in the ACA model. Subsequently, we present the technique to assess the capability of the model to predict the recession rates at
higher pressure.

The posterior distribution of the parameters to calibrate, $\pi(\mathbf{x}|\mathbf{y})$, given a set of observations, $\mathbf{y}$, can be obtained by solving the Bayesian inverse problem. According to Bayes' theorem:

$$ \pi(\mathbf{x}|\mathbf{y}) = \frac{\pi(\mathbf{x})\pi(\mathbf{y}|\mathbf{x})}{\int \pi(\mathbf{x})\pi(\mathbf{y}|\mathbf{x})d\mathbf{x}}, \qquad (16) $$

where $\pi(\mathbf{x})$ is the prior distribution of $\mathbf{x}$, $\pi(\mathbf{y}|\mathbf{x})$ is the likelihood of an observable quantity $\mathbf{y}$ given $\mathbf{x}$, and the integral at the denominator is the marginal likelihood or evidence, which ensures the distribution integrates to 1. Based on preliminary investigations, we choose to characterize 5 activation energies ($E_1, E_3, E_4, E_5, E_6$), 4 elementary reaction efficiencies ($X_3, X_4, X_5, X_6$), and the total active site density ($B$) of the model in Table 2. The vector of parameters to characterize is:

$$ \mathbf{x} = [E_1, E_3, E_4, E_5, E_6, \log_{10}(X_3), \log_{10}(X_4), \log_{10}(X_5), \log_{10}(X_6), \log_{10}(B)]. \qquad (17) $$

The elementary reaction efficiencies and the total active site density have been expressed as $\log_{10}$ because their estimation can potentially span several orders of magnitude.

The reaction efficiencies obtained from the molecular beam experiment²³ are chosen as observable quantities to calibrate the model:

$$ \mathbf{y} = [\gamma_{\text{CN}}, \gamma_{\text{N2}}]_i, \quad i = 1, \dots, N_{\text{obs}}, \qquad (18) $$

where $N_{\text{obs}}$ is the number of observations. Assuming that the $2N_{\text{obs}}$ observations are independent and that the model error is negligible with respect to the experimental uncertainty, $\sigma_i$, on the measurement $y_i$ (i-th component of $\mathbf{y}$), the likelihood reads:

$$ \pi(\mathbf{y}|\mathbf{x}) = \prod_{i=1}^{2N_{\text{obs}}} N(y_i|M_i(\mathbf{x}), \sigma_i^2). \qquad (19) $$

The symbol $M_i(\mathbf{x})$ stands for a mathematical model which expresses $\mathbf{y}$ as a function of $\mathbf{x}$. The 0D model described in section 3.2 is used to explain the molecular beam data.

Non-informative uniform distributions on $\mathbf{x}$ over a wide range are chosen to encompass several values available in the literature. The prior distributions, $\pi(\mathbf{x})$, are reported in Table 3.

<table><caption>Table 3: Non-informative uniform prior distribution on the parameters to calibrate.</caption><thead><tr><th>Parameter</th><th>min</th><th>max</th><th>Prata et al.</th><th>Parameter</th><th>min</th><th>max</th><th>Prata et al.</th></tr></thead><tbody><tr><td>E<sub>1</sub></td><td>0</td><td>3000</td><td>2500</td><td>log<sub>10</sub>(B)</td><td>-7</td><td>-3</td><td>-5</td></tr><tr><td>E<sub>3</sub></td><td>0</td><td>15000</td><td>7000</td><td>log<sub>10</sub>(X<sub>3</sub>)</td><td>-4</td><td>2</td><td>0.176</td></tr><tr><td>E<sub>4</sub></td><td>0</td><td>15000</td><td>2000</td><td>log<sub>10</sub>(X<sub>4</sub>)</td><td>-4</td><td>0</td><td>-0.301</td></tr><tr><td>E<sub>5</sub></td><td>0</td><td>35000</td><td>21000</td><td>log<sub>10</sub>(X<sub>5</sub>)</td><td>-4</td><td>0</td><td>-1</td></tr><tr><td>E<sub>6</sub></td><td>15000</td><td>35000</td><td>20676</td><td>log<sub>10</sub>(X<sub>6</sub>)</td><td>5</td><td>12</td><td>8</td></tr></tbody></table>

The posterior distribution is then numerically obtained by constructing an MCMC, using the affine invariant ensemble algorithm¹⁰ implemented in the UQLAB software.¹⁷

We then use the Plasmatron data for assessing the predictive capability of the calibrated model at higher pressure. Specifically, we compare the model predictions of the mass blowing rate to the value experimentally measured:

$$ \mathbf{y}^{\text{ICP}} = \dot{m}_i, \quad i = 1, \dots, N_{\text{ICP}}, \qquad (20) $$

where $N_{\text{ICP}}$ is the number of Plasamtron tests that we use for the assessment. To this end, the obtained posterior joint distribution is propagated using a Monte Carlo technique. First, a subset of points is randomly sampled from the posterior distribution. A CFD computation is then run for each experimental condition using the sampled FRC model parameters to describe the surface chemistry. Estimates of the probability density of mass blowing rates for each experimental point are obtained by building histograms on the CFD output.

## 5. Results

In this section, we first present the results of the stochastic calibration in terms of parameters' marginal and joint distributions. Subsequently, we show the posterior predictive results for a span of surface temperatures. Lastly, we show the results of the model assessment at higher pressure.
## 5.1 Marginal and joint distributions

We use two different sets of molecular beam points for the calibration. The first set (hereafter referred to as 'set A') does not include the reaction efficiencies obtained for $T < 1000$ K. This subset is the one used for constructing the ACA model. The second calibration was performed using all the available points (hereafter referred to as 'set B'). For both conditions, 100 MCMC chains, of 10000 sampled points each, are constructed. The first 20% of points are discarded as part of the burn-in process.

The marginal distributions of the calibrated parameters for the two nitridation reactions are plotted in Figures 1a, 1b, 2a, and 2b, while in Figures 1c, and 2c the two joint distributions relating the reaction efficiency and the activation energy of selected reactions are given. The values obtained by Prata et al. in a previous calibration,²⁸ and the theoretical values obtained by Nieman et al.²⁴ are also shown in the same figures. Both the nitridation reactions are well characterized using both sets of experiments, especially when using set B, as it contains more information. When the model is calibrated using set A, the posterior distribution includes both the value proposed by Prata et al. and the theoretical ones computed by Nieman et al., which makes them consistent with the data. However, when low-temperature experimental efficiencies are added in the calibration (set B), all the distributions shift towards lower values to allow the model to predict a non-decreasing trend of the nitridation efficiency at low temperatures, as observed in the molecular beam experiments, see Figures 6 and 7. In this case, the distribution of the calibrated activation energies departs from the theoretical values, suggesting some incompatibility of these values in exempling the molecular beam nitridation efficiencies at low temperatures through the model proposed.

<figure><img src="image_3.png" /><figcaption>Figure 1: Probability densities of the efficiency and the activation energy of reaction 3. Blue histograms refer to the parameters posterior using set A, black histograms using set B, red lines indicate the prior distributions, green lines the values obtained by Prata et al., and yellow lines the theoretical values obtained by Nieman et al.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2: Probability densities of the efficiency and the activation energy of reaction 6. Blue histograms refer to the parameters posterior using set A, black histograms using set B, red lines indicate the prior distributions, green lines the values obtained by Prata et al., and yellow lines the theoretical values obtained by Nieman et al.</figcaption></figure>

The marginal distributions of the calibrated parameters for the two recombination reactions are plotted in Figures 3a, 3b, 4a, and 4b, while in Figures 3c, and 4c the two joint distributions for the reaction efficiency and the activation energy are shown. The LH mechanism ($N_s + N_s \rightarrow N_2 + 2s$) appears to be the dominant recombination mechanism, in agreement with what was observed by Murray et al.²³ In fact, the ER reaction ($N_s + N \rightarrow N_2 + s$) appears to be poorly characterized with the molecular beam experiments considered. The distribution for the activation energy of the LH mechanism has a peak very close to the value proposed by Prata et al. when set A is used. Such peak
shifts toward a lower value using set B, closer to the theoretical value of Nieman.

<figure><img src="image_5.png" /><figcaption>Figure 3: Probability densities of the efficiency and the activation energy of reaction 4. Blue histograms refer to the parameters posterior using set A, black histograms using set B, red lines indicate the prior distributions, green lines the values obtained by Prata et al., and yellow lines the theoretical values obtained by Nieman et al.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Probability densities of the efficiency and the activation energy of reaction 5. Blue histograms refer to the parameters posterior using set A, black histograms using set B, red lines indicate the prior distributions, green lines the values obtained by Prata et al., and yellow lines the theoretical values obtained by Nieman et al.</figcaption></figure>

Lastly, the activation energy of the absorption reaction, Figure 5a, is also well characterized, with a distribution
tending towards small values, close to the lowest value proposed by Nieman. Such low values of activation energy
allow the model to predict the very high values of recombination efficiencies experimentally observed in the molecular
beam experiments. The total active site density, Figure 5b, has a broad distribution with a peak close to the value
proposed by Prata et al. when using set A. It moves towards small values when set B is used. However, its distribution
is not very well characterized by the molecular beam data employed in the calibration, suggesting the need of further
experimental points.

5.2 Posterior predictive result

The joint posterior distribution is propagated through the 0D model describing the molecular beam to obtain the pre-
dictive posterior both at low and high pressure. The model's posterior prediction using the set A is plotted in Figures 6a
and 6b. The 22 molecular beam experimental points,²³ and the marginals obtained by del Val et al.⁷ are shown on
the same figures. The predictions of the calibrated model are in perfect agreement with the molecular beam points
used for calibration. Large model uncertainties are encountered in the prediction at high pressure as no calibration
point was used. This suggests the need for higher-pressure experimental points to reduce the uncertainty and make the
predictions more accurate. The same can be achieved, as we will see next, by using experimental points characterized
by a surface coverage similar to the one at high pressure. In fact, the uncertainty on the nitridation reaction at high
pressure, shown in Figure 7b, is drastically reduced when using the data from set B. Nitridation reaction efficiencies at
such high pressures also exhibit excellent agreement with the marginals from del Val et al. at high pressure, suggesting
that the two experiments can be accurately explained through the use of the same calibrated model.

This improved characterization obtained when using the set B can be explained by looking at Figure 8, where
the surface coverage at the two different pressures is plotted as a function of the temperature. As one can see, the
value of the surface coverage at low-pressure/low-temperature is close to the value at high-pressure. Furthermore,
<figure><img src="image_7.png" /><figcaption>Figure 5: Probability densities of the activation energy of reaction 1 and total active site density. Blue histograms refer to the parameters posterior using set A, black histograms using set B, red lines indicate the prior distributions, green lines the values obtained by et al., and yellow lines the theoretical values obtained by Nieman et al.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Model prediction using set A. The red and blue color indicates, respectively, the recombination and the nitridation efficiency. The solid lines represent the mean of the prediction at molecular beam pressure (low pressure), while the dashed lines are the mean at the Plasmatron pressure (high pressure). The uncertainty corresponding to the 95% Confidence Interval (CI) on the model prediction at low pressure is shown in the shadowed region of Figure 6a for low pressures and Figure 6b for high pressures. The error bars with square as mean represent the 22 molecular beam experimental points,²³ while the error bars with circle as mean are the marginals obtained by del Val et al.⁷</figcaption></figure>

the experimental uncertainty on the nitridation reaction efficiency is small. Thus, the low-temperature points bring information on the model behavior at high surface coverage. Accordingly, also the uncertainty at low temperatures of the prediction at low pressure decrease, as seem in Figure 7a.

Regarding the main trends, a sharp increase in the recombination efficiency is predicted at low pressure, in agreement with the assumptions of Murray et al.²³ After a plateau, it decreases around 2250 K because the surface coverage sharply decreases at this temperature, as shown in Figure 8. The nitridation efficiency follows also the experimental trend: first, it increases, mostly driven by reaction 3 of Table 2. After a brief downward trend due to reaction 3, the nitridation efficiency shows a steeper upward slope for higher surface temperatures, driven by reaction 6. For high temperatures, nitridation reaction efficiencies show a downward trend due to the surface coverage. At higher pressure, an almost flat trend in the nitridation probability is predicted, due to the nearly constant surface coverage. On the other hand, the recombination efficiency, mostly driven by reaction 4 at high pressure, is still characterized by large uncertainty, as the uncertainty on the molecular beam experiments at low temperatures is high. The inclusion of more accurate data at low-temperature in the calibration would help to decrease the model uncertainty at high pressure.
<figure><img src="image_9.png" /></figure>

Figure 7: Model prediction using set B. The red and blue color indicates, respectively, the recombination and the nitridation efficiency. The solid lines represent the mean of the prediction at molecular beam pressure (low pressure), while the dashed lines are the mean at the Plasmatron pressure (high pressure). The uncertainty corresponding to the 95% CI on the model prediction at low pressure is shown in the shadowed region of Figure 6a for low pressures and Figure 6b for high pressures. The error bars with square as mean represent the 22 molecular beam experimental points,²³ while the error bars with circle as mean are the marginals obtained by del Val et al.⁷

<figure><img src="image_10.png" /></figure>

Figure 8: Mean of the model prediction of the surface coverage using the set B for the calibration. Solid line refers to the prediction at molecular beam pressure and the dashed line at Plasmatron pressure.

## 5.3 Predictions for Plasmatron experiments

The parameters' posteriors are propagated using a Monte Carlo method: 200 points were randomly sampled by the posteriors' distribution and a CFD was run for each of the Plasmatron conditions.

The propagated uncertainty on the mass blowing rate, along with the experimental uncertainty, is shown in Figure 9, for the cases G4, G5, G6, and G7. It can be seen that the propagated mass blowing rates share much of the support with their experimental counterparts. However, the calibrated model appears to be less sensitive to the surface temperature than what was experimentally observed. The propagated joint distribution of $\gamma_{CN}$, $\gamma_{N2}$, relative to the G5 experiment, is shown in Figure 10, along with the joint distribution obtained by del Val et al.⁷ Also in this case, the two distributions agree fairly well, but the predictive distribution tends to higher $\gamma_{CN}$ values than the ones obtained by del Val. Such a result is consistent with what was observed for the distribution of the mass blowing rates, which is higher than the experimental counterpart for the G5 test. Using these data in the calibration process would help in identifying whether the model is capable of capturing such an upward trend.

Overall the model exhibits a good agreement with the experiment at high pressure, even if no value at this pressure was
<figure><img src="image_11.png" /><figcaption>Figure 9: Mass blowing rate distributions. The continuous lines refer to the experimental uncertainty while the histograms show the propagated distributions. Blue, Grey, red, and green colors indicate, respectively, the G4, G5, G6, and G7 case.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10: Joint distribution of γCN, γN2 on G5 case condition. Green and black lines refer, respectively, to the values obtained by del Val et al.⁷ and in this study.</figcaption></figure>

used in the calibration.

## 6. Conclusion

In this study, we performed a Bayesian calibration of the reactions involving nitrogen of the ACA model by means of molecular beam-surface scattering data. The model parameters were characterized: the nitridation reactions, as well as the absorption process, and the LH recombination reaction turn out to be well characterized. However, the probability distributions of some of the activation energies do not encompass the values theoretically computed by Nieman, indicating the potential incompatibility of these values in simulating the molecular beam nitridation efficiencies through the calibrated model. The ER mechanism and the total number of active sites were found to be poorly characterized with the available data.

Despite this, the model posteriors predictions exhibit excellent agreement with the calibration points. It was shown that using low-pressure low-temperature molecular beam data, characterized by high surface coverage, improves the model prediction at high pressure, also characterized by high surface coverage. The model predictions at high pressure are also in agreement with the values obtained by del Val et al., suggesting that the two experiments can be accurately explained through the same model. Lastly, it was noted that substantial uncertainty characterizes the prediction of the recombination probability, especially at high pressure. More accurate molecular beam data at low pressure and low temperature would potentially lead to a reduction of this uncertainty, enhancing our understanding of the process.

The calibrated model was also used to predict the mass blowing rate obtained in Plasmatron experiments, whose BL was simulated using CFD simulations with an ablative BC. The model predictions exhibit an overall good agreement with the experimental uncertainty, indicating that the model is capable of explaining both sets of experimental data.

However, we observed that the calibrated model is less sensitive to the temperature than the experimental data. Incorporating these data points in the calibration process could lead to an improvement in the predictions.

Future work includes using the Plasmatron data points in the calibration process and the use of a different set of experiments for the assessment.

## References

[1] A. Balter-Peterson, F. Nichols, B. Mifsud, and W. Love. Arc jet testing in NASA Ames Research Center thermophysics facilities. In *AlAA 4th International Aerospace Planes Conference*. American Institute of Aeronautics and Astronautics, 12 1992.

[2] P. Barbante. *Accurate and Efficient Modelling of High Temperature Nonequilibrium Air Flows*. PhD thesis, Université Libre de Bruxelles, 2001.
[3] G. Bellas Chatzigeorgis, A. Turchi, A. Viladegut, O. Chazot, P. Barbante, and T. Magin. Development of catalytic and ablative gas-surface interaction models for the simulation of reacting gas mixtures. In *23rd AIAA Computational Fluid Dynamics Conference*. American Institute of Aeronautics and Astronautics, 2017.

[4] B. Bottin, O. Chazot, M. Carbonaro, V. Haegen, and S. Paris. The VKI Plasmatron Characteristics and Performance. page 26, 04 2000.

[5] M. Capriati, K. Prata, T. Schwartzentruber, G. Candler, and T. Magin. Development of a nitridation gas-surface boundary condition for high-fidelity hypersonic simulations. In *14th WCCM-ECCOMAS Congress*. CIMNE, 2021.

[6] A. Del Val and O. Chazot. Stochastic determination of gas-phase chemical reaction rates from experiments in air plasmas. In *AIAA SCITECH 2023 Forum*. American Institute of Aeronautics and Astronautics, 01 2023.

[7] A. del Val, O. Le Maitre, P. Congedo, and T. Magin. Stochastic calibration of a carbon nitridation model from plasma wind tunnel experiments using a bayesian formulation. *Carbon*, 200:199-214, 11 2022.

[8] A. del Val, O. Le Maitre, O. Chazot, P. Congedo, and T. Magin. Inference Methods for Gas-Surface Interaction Models: From Deterministic Approaches to Bayesian Techniques. In *Massimiliano Vasile and Domenico Quagliarella, editors, Advances in Uncertainty Quantification and Optimization Under Uncertainty with Aerospace Applications*, volume 8, pages 349-364. Springer International Publishing, 2021. Series Title: Space Technology Proceedings.

[9] U. Duzel, O. Schroeder, H. Zhang, and A. Martin. Numerical simulation of an arc jet test section. *Journal of Thermophysics and Heat Transfer*, 34(2):393-403, 04 2020.

[10] J. Goodman and J. Weare. Ensemble samplers with affine invariance. *Communications in Applied Mathematics and Computational Science*, 5(1):65-80, 01 2010.

[11] A. Gordeev, A. Kolesnikov, and V. Sakharov. Flow and heat transfer in underexpanded nonequilibrium jets of an induction plasmatron. *Fluid Dynamics - FLUID DYN*, 46:623-633, 08 2011.

[12] R. Goulard. On Catalytic Recombination Rates in Hypersonic Stagnation Heat Transfer. *Journal of Jet Propulsion*, 28(11):737-745, 1958.

[13] B. Helber, A. Turchi, and T. Magin. Determination of active nitridation reaction efficiency of graphite in industriely coupled plasma flows. *Carbon*, 125:582-594, 2017.

[14] A. Klomfass and S. Muller. Calculation of stagnation streamline quantities in hypersonic blunt body flows. *Shock Waves*, 7(1):13-23, 02 1997.

[15] T. Magin and G. Degrez. Transport algorithms for partially ionized and unmagnetized plasmas. *Journal of Computational Physics*, 198(2):424-449, 2004.

[16] T. Magin and G. Degrez. Transport properties of partially ionized and unmagnetized plasmas. *Physical review. E, Statistical, nonlinear, and soft matter physics*, 70:046412, Oct 2004.

[17] S. Marelli and B. Sudret. UQLab: A Framework for Uncertainty Quantification in Matlab. In *Vulnerability, Uncertainty, and Risk*, pages 2554-2563. American Society of Civil Engineers, 2014.

[18] J. Marschall and M. MacLean. Finite-Rate Surface Chemistry Model, I: Formulation and Reaction System Examples. In *42nd AIAA Thermophysics Conference*, 06 2011.

[19] B. McBride, M. Zehe, and S. Gordon. Nasa glenn coefficients for calculating thermodynamic properties of individual species. Technical report, NASA, 10 2002.

[20] A. Munafó and T. Magin. Modeling of stagnation-line nonequilibrium flows by means of quantum based collisional models. *Physics of Fluids*, 26(9):097102, 2014.

[21] V. Murray, B. Marshall, P. Woodburn, and T. Minton. Inelastic and Reactive Scattering Dynamics of Hyperthermal O and O2 on Hot Vitreous Carbon Surfaces. *The Journal of Physical Chemistry C*, 119(26):14780-14796, 07 2015.
[22] V. Murray and T. Minton. Gas-surface interactions of atomic nitrogen with vitreous carbon. *Carbon*, 150:85-92, 2019.

[23] V. Murray, P. Recio, A. Caracciolo, C.e Miossec, N. Balucani, P. Casavecchia, and T. Minton. Oxidation and nitridation of vitreous carbon at high temperatures. *Carbon*, 167:388-402, 10 2020.

[24] R. Nieman, M. Sands, Y. Wang, T. Minton, E. Mussoni, J. Engerer, and H. Guo. Informing air-carbon ablation modeling with theoretical calculations of atomic oxygen and nitrogen interacting with carbon surfaces. *Physical Chemistry Chemical Physics*, 2023.

[25] D. Olynick, Y. Chen, and M. Tauber. Aerothermodynamics of the Stardust Sample Return Capsul. *Journal of Spacecraft and Rockets*, 36:442-462, 1998.

[26] C. Park. Effects of atomic oxygen on graphite ablation. *AIAA Journal*, 14(11):1640-1642, 1976.

[27] S. Poovathingal, T. Schwartzentruber, V. Murray, T. Minton, and G. Candler. Finite-rate oxidation model for carbon surfaces from molecular beam experiments. *AIAA Journal*, 55(5):1644-1658, 2017.

[28] K. Prata, T. Schwartzentruber, and T. Minton. Air-Carbon Ablation Model for Hypersonic Flight from Molecular-Beam Data. *AIAA Journal*, 60(2):627-640, 2022.

[29] C. Purpura, F. De Filippis, P. Barrera, and D. Mandanici. Experimental characterisation of the cira plasma wind tunnel scirocco test section. *Acta Astronautica*, 62(6):410-421, 2008.

[30] P. Roe. Approximate riemann solvers, parameter vectors, and difference schemes. *Journal of Computational Physics*, 43(2):357-372, 1981.

[31] J. Scoggins, V. Leroy, G. Bellas-Chatzigeorgis, B. Dias, and T. Magin. Mutation + + : MUlticomponent Thermodynamic And Transport properties for IONized gases in C++. *SoftwareX*, 12:100575, 2020.