# Combining Molecular Beam and Plasmatron Data for Robust Calibration of Nitridation Mechanism

Michele Capriati, Anabel del Val, Thomas Schwartzentruber, Timothy Minton, Pietro Marco Congedo, Thierry Magin

► To cite this version:

Michele Capriati, Anabel del Val, Thomas Schwartzentruber, Timothy Minton, Pietro Marco Congedo, et al.. Combining Molecular Beam and Plasmatron Data for Robust Calibration of Nitridation Mechanism. AIAA Journal, 2025, pp.1-15. 10.2514/1.J065221. hal-05249808

HAL Id: hal-05249808

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

Submitted on 11 Sep 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Combining molecular beam and Plasmatron data for robust
calibration of nitridation mechanism

Michele Capriati *

von Karman Institute for Fluid Dynamics, Waterloosesteenweg 72, Sint-Genesius-Rode, 1640, BE
École polytechnique, IPP, Route de Saclay, 91120 Palaiseau, France

Anabel del Val † and Thomas E. Schwartzentruber ‡

University of Minnesota, 110 Union Street SE, Minneapolis, 55455, Minnesota, USA

Timothy K. Minton §

University of Colorado, 3775 Discovery Dr, Boulder, 80303, Colorado, USA

Pietro M. Congedo ¶

École polytechnique, IPP, Route de Saclay, 91120 Palaiseau, France

Thierry E. Magin ‖

von Karman Institute for Fluid Dynamics, Waterloosesteenweg 72, Sint-Genesius-Rode, 1640, BE
Université Libre de Bruxelles, Avenue F.D. Roosevelt 50, 1050, Brussels, Belgium

The modeling of gas-surface interaction phenomena is crucial for accurately predicting the heat flux and the mass loss experienced by hypersonic vehicles. Gas-surface interactions refer to the phenomena occurring between the reacting gas and the thermal protection material. An important part of the modeling concerns the description of the surface chemical reactions. In this regard, we propose a novel methodology to infer the parameters underlying such surface chemistry models. It combines uncertainty quantification techniques with state-of-the-art modeling and different types of experiments. The methodology is used to calibrate, in a Bayesian sense, the rates of the elementary reactions between a nitrogen gas and a carbon surface. We rely both on molecular beam and plasma wind tunnel observations. The former provides detailed data on the chemical mechanisms, but is characterized by pressures non-representative of atmospheric entries. By contrast, plasma wind tunnel experiments are conducted at representative pressures, but contain only macroscopic information. The parameters' posterior distributions are then propagated through the models representing the two experiments. The calibrated model is able to satisfactorily explain both experiments, highlighting the robustness of the proposed methodology.

<footnote>*Research Engineer, Aeronautics and Aerospace Department. Corresponding author: michele.capriati@vki.ac.be</footnote>

<footnote>†Assistant Professor, Department of Aerospace Engineering and Mechanics</footnote>

<footnote>‡Professor, Department of Aerospace Engineering and Mechanics</footnote>

<footnote>§Professor, Ann and H.J. Smead Department of Aerospace Engineering Sciences</footnote>

<footnote>¶Research Director, Inria, Centre de Mathématiques Appliquées</footnote>

<footnote>‖Professor, Aeronautics and Aerospace Department</footnote>
Nomenclature

$B$ = total active sites

$e$ = thermodynamic energy

$E$ = total energy

$H$ = total enthalpy

$h$ = Planck constant

$j$ = diffusive mass flux

$K_H$ = Homman's correction

$k_B$ = Boltzmann constant

$m$ = particle mass

$\dot{m}$ = mass blowing rate

$N$ = number flux

$NA$ = Avogadro number

$NRMSE$ = normalized root mean square error

$p$ = pressure

$P_{el}$ = generator power

$q$ = heat flux

$t$ = time

$\tau$ = viscous stress tensor

$T$ = temperature

$u$ = velocity

$V_T$ = quarter of mean thermal velocity

$V_{T,2D}$ = mean thermal velocity of mobile absorbed particle

$X$ = elementary reaction efficiency

$Z$ = stationary Gaussian process

$\gamma$ = reaction efficiency

$\mu$ = mean

$\pi$ = distribution

$\rho$ = density

$\sigma$ = standard deviation

$\dot{\omega}$ = chemical production rate

Subscripts
d = dynamic

e = boundary layer edge

g = gas

v = validation

w = wall

# I. Introduction

Spacecraft entering a planet's atmosphere are exposed to extreme conditions, including 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) which interacts with the surface of the vehicle. Vehicles are thus protected with Thermal Protection Materials (TPM) for a safe landing. Two different classes of TPM 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 spacecraft. On the other hand, ablative TPMs dissipate energy mostly by degradation and they experience mass loss.

Gas-Surface Interaction (GSI) describes the set of phenomena occurring between the reacting gas and the solid surface. The GSI should be accurately modeled to predict the in-flight heat flux and mass loss. One key aspect of GSI modeling is the description of the chemical reactions happening at the surface. Among others, reactions involving atomic nitrogen are important in entry applications, for example in Earth and in Titan atmospheres. In fact, the nitridation reaction ablates the TPM, while the exothermic recombination of atomic nitrogen into molecular nitrogen increases the heat flux experienced by the spacecraft.

Several experimental works addressed the characterization of the efficiencies of the two reactions. For instance, the nitridation efficiency was estimated, over a wide range of temperatures and pressures, by means of experiments in shock tube [1], in furnace-heated quartz flow tube coupled to microwave discharge [2], in Inductively Coupled Plasma (ICP) facilities [3-7], and using a molecular beam-surface scattering apparatus [8]. The studies [3, 5, 8] also included the determination of nitrogen recombination efficiency. The values proposed in these works largely differ from each other, suggesting a poor understanding of both the nitridation and the recombination reactions.

It is important to remark that, in most of the above-mentioned studies, coupled numerical simulations supported the experiments in determining those Quantities of Interest (QoIs) that were not experimentally accessible, adding an extra layer of modeling uncertainty. For example, the concentration of the atomic nitrogen in ICP facilities, needed to compute the reaction efficiencies, was determined through numerical simulations in [4, 6, 7]. In such computations, the surface recombination was neglected, potentially biasing the determination of the gas concentrations. Furthermore,
Computational Fluid Dynamics (CFD) simulations contain assumptions also on the gas-phase modeling, for instance on the rates of the homogeneous chemistry and on the transport properties. Consequently, the assumptions contained in the modeling of these processes directly affect the characterization of the surface reactions. However, it is virtually impossible to characterize all these sources of uncertainties, along with the surface reaction probabilities, given the limited amount of observations that can be acquired in these facilities. Contrary to ICP, in the molecular beam-surface scattering experiments, the ablative sample is placed in a vacuum chamber and bombarded with a beam of molecules or atoms. Thus, few gas collisions are expected to happen and phenomena such as gas chemistry or transport are negligible. As a consequence, simulating these experiments requires no assumptions in the modeling of the gas. However, while these experiments manage to isolate the GSI phenomenon, we cannot fully trust that they are representative of hyper-sonic flight conditions because the very low pressures involved may result in a surface coverage state that differs from flight.

From a modeling standpoint, phenomenological Arrhenius laws for the nitridation reaction were proposed by Suzuki et al. [6], and, later, by Helber et al. [4]. Such Arrhenius laws capture only the temperature dependence of the reaction and are expected to be valid at pressures close to the ones involved in the calibration experiments. Recently, Prata et al. [9] derived the Air-Carbon Ablation (ACA) model based both on molecular beam and ICP experiments. To capture the different probabilities observed at different temperatures and pressures, both these dependencies were accounted for in the model. Furthermore, the model was calibrated to capture the main trends of the reaction efficiencies experimentally observed. Among these reactions, the ones involving atomic nitrogen were calibrated using the molecular beam experiments performed by Murray et al. [10], and the ICP experiments from Lutz [5] and Helber et al. [4]. It is important to note that the ACA model fitting and calibration only used nominal quantities inferred from the experiments and did not include any experimental error or quantified uncertainty.

Few Uncertainty Quantification (UQ) works have addressed a robust characterization of the nitridation reaction. Upadhyay et al. [11] employed a Bayesian method to infer the reaction probability accounting for the experimental and modeling uncertainties in the data from Zhang et al. [2]. Four different models, differing in the underlying modeling assumptions, were used for the calibration. The temperature dependence of the nitridation reaction was captured, for each model, by calibrating a power law. A Bayesian framework was also used by del Val et al. [12] to infer the nitridation efficiency from the ICP data from Helber et al. [4]. First, a sensitivity analysis was performed. It highlighted that the computation of the experimental observations was highly dependent on the value of the nitridation efficiency, while the Boundary Layer (BL) edge conditions had little influence. After that, the nitridation probabilities were individually inferred from each surface temperature-dependent measurement, then an Arrhenius law was calibrated considering them jointly. No heterogeneous recombination was considered. In a subsequent work by the same authors [13], this aspect was tackled, along with other model-form uncertainties. Although the analysis did not highlight evidence of the
recombination reaction, it showed that its uncertainty broadens the posterior distribution of the nitridation probability. All the above-mentioned inferences relied on experimental observations obtained in a single facility, involving similar operating pressures. Hence, a phenomenological model was adequate to describe the task at hand. To date, no stochastic calibration has been performed attempting to capture the reaction pressure dependence by informing a Finite-Rate Chemistry (FRC) model combining observations from very different experiments.

We propose a new methodology to infer the reaction rates of FRC models. It relies on UQ techniques to combine I) molecular beam data, which are highly informative, and II) plasma wind tunnel observations, as they are representative of atmospheric entry conditions.

The proposed methodology is employed to calibrate the elementary reactions involving nitrogen of the ACA model [9] by means of both the molecular beam experiments performed by Murray et al. [10], and the Plasmatron experiments from and Helber et al. [4]. First, the Plasmatron data were used to assess the predictive capability of the model calibrated only on molecular beam data. Then, we considered both sets of data in a joint inference framework. Unlike the work from Prata et al. [9], we simulated the Plasmatron experiments with a CFD solver, where an ablative Boundary Condition (BC) is embedded to properly describe the BL developing in front of the test sample. The molecular beam data were reproduced with the same analytical model proposed by Prata et al. [9]. The parameters' posterior distributions were obtained by drawing samples by means of Markov Chain Monte Carlo (MCMC). The posterior distributions were then propagated through both the CFD and the analytical model by Prata et al. [9] to obtain the posterior predictions. These are compared to the experimental counterparts.

The article is structured as follows. In Section II, we present two different experimental campaigns, whose data served as input for solving the inverse problem. In Section III, we introduce the models describing such experiments. The UQ tools and methodologies that we employed to solve the forward and inverse problem are described in Section IV, while Section V is devoted to showing the UQ results. Finally, conclusions and future perspectives are drawn in Section VI.

## II. Experiments

In this section, we describe the facilities and experiments whose data were employed to inform the FRC model. Although other data are available [1, 5–7], we decided to focus our attention on the molecular beam and ICP experiments performed, respectively, by Murray et al. [8] and Helber et al. [4].
## A. Reaction probabilities from molecular beam-surface scattering

The first class of data employed in the solution of the inverse problem comes from a molecular beam experimental campaign by Murray et al. [8].

The main components of the experimental apparatus are I) a beam source, II) a heated sample, over which the beam is directed, III) a rotatable mass spectrometer to detect the number density distribution of the scattered products, and IV) a chopper wheel to modulate the scattered products. The sample is placed in a vacuum chamber, where pumps guarantee pressure of the order of $10 \times 10^{-7}$ torr, essential to avoid gas-phase collisions. The temperature of the sample is measured through a calorimeter.

In the experiment, continuous beams were employed to study the reactions between a nitrogen gas and a vitreous carbon surface. The nitrogen beam was produced using a 95 mbar mixture (2.5% N₂ in He) discharged by a high-pressure radio frequency discharge 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 and after by a 0.4 × 0.4 mm² aperture. The resulting beam is composed by $x_N = 0.18$ and $x_{N2} = 0.82$ with an average velocity of ≈ 2000 m/s. The data were collected at increasing temperatures (from 800 to 1873 K), while annealing the sample at 1873 K for around 1 hour before each test.

Experimental results showed that the surface undergoes only recombination and nitridation reactions. Their reaction probabilities were determined as:

$$ \gamma^r = \frac{\text{flux of } r \text{ product}}{\text{flux of N atoms onto the surface}}. \qquad (1) $$

Since the experiment was conducted in vacuum, one can neglect the gas-phase collisions and employ the 0D model proposed by Prata et al. [9] to describe the phenomenon (presented in Section III.A).

## B. Surface recession from high-speed imaging in ICP facility

The FRC model was also informed by employing the data from a nitridation experimental campaign performed by Helber et al. [4] in the von Karman Institute (VKI) Plasmatron facility. Such a facility can duplicate the stagnation-point heat flux and recession experienced by an object traveling at hypersonic speed by means of a subsonic flow [14–16], ensuring that a model, learned on these data, is representative of atmospheric entry conditions.

The facility generates high-purity/high-enthalpy flows by means of a 160 mm diameter ICP torch powered by a high-frequency, high-power, high-voltage (400 kHz, 1.2 MW, 2 kV) generator [17]. The facility is started with Argon, as it is easy to ionize. Once the plasma is stabilized, the gas is switched to the operating one (air, N₂ or CO₂). The gas exits in the test chamber, where vacuum pumps maintain the pressure at the operating value. A holder places the test sample into the centerline of the plasma jet.

The experiment consisted of a nitrogen plasma flow over a superfine-grain/high-density graphite ($\rho_s = 1760 \text{ kg/m}^3$),
machined in-house to have a 25 mm radius hemispherical shape with a 25 mm cylindrical after-body. It was placed at a distance of 445 mm from the exit of the torch. The experimental set-up included a digital camera for the determination of the stagnation-point recession, a two-color pyrometer for measuring the surface temperature of the sample, and a spectrometer to record the violet spectral band of the CN molecule in the BL. It resolved the radial distribution of the emission versus the wavelength of the signal. Such a 2D matrix can be Abel-inverted to obtain the locally-resolved CN-emission. By post-processing it, it was possible to retrieve both the gas temperature and the CN concentration at a point belonging to the BL.

Seven runs were performed. From one to the other, the electrical power supplied to the generator was progressively increased to obtain a higher free-stream enthalpy, and, in turn, a higher surface temperature of the sample. No mass loss of the sample was observed in three of the seven runs. Such runs will not be considered for the calibration proposed within this work. An overview of the conditions of the four considered runs is given in Table 1, where, beside the experimental uncertainties, the values of the BL edge temperature and velocity rebuilt by Helber et al. [4] are given.

<table><caption>**Table 1** Overview of Plasmatron test conditions for a nitrogen flow at 1500 ± 50 Pa over a 25 mm radius hemisphere carbonaceous sample: dynamic pressure $p_d$, generator power $P_{el}$, sample mean temperature $T_w$, mass blowing rate $\dot{m}$, and reconstructed BL edge conditions (temperature, $T_e$, and velocity, $u_e$). Conditions from Helber et al. [4].</caption><thead><tr><th rowspan="2">CASE</th><th colspan="4">Experimental</th><th colspan="2">Rebuilt</th></tr><tr><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>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 ± 4.6</td><td>280</td><td>2225 ± 22</td><td>2.49 ± 0.91</td><td>10005</td><td>554</td></tr><tr><td>G5</td><td>268 ± 5.4</td><td>330</td><td>2410 ± 24</td><td>2.89 ± 0.97</td><td>10280</td><td>562</td></tr><tr><td>G6</td><td>312 ± 6.2</td><td>370</td><td>2535 ± 26</td><td>4.41 ± 0.80</td><td>11040</td><td>846</td></tr><tr><td>G7</td><td>330 ± 6.6</td><td>390</td><td>2575 ± 26</td><td>4.56 ± 0.70</td><td>10970</td><td>859</td></tr></tbody></table>

### III. Numerical modeling

In this section, we first introduce the ACA ablation model. Then, we present the 0D model employed to describe the molecular beam experiments, and the stagnation line code used to simulate the BL in the Plasmatron experiments.

#### A. Air-carbon ablation model

The ACA model proposed by Prata et al. [9] comprises 20 elementary reactions occurring between a carbon surface and an air gas. We are interested in a subset of reactions involving atomic nitrogen. Such a subset, as well as the corresponding rates, is reported in Table 2.
<table><caption>Table 2 Subset of the FRC reactions involving nitrogen in the ACA model.</caption><thead><tr><th>Index</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>V<sub>T,N</sub></sup>&frasl;<sub>B</sub> exp <sup>-2500</sup>&frasl;<sub>T<sub>w</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>b</sub><sup>2</sup>T<sub>w</sub><sup>2</sup></sup>&frasl;<sub>N<sub>A</sub>Bh<sup>3</sup></sub> exp <sup>-73971</sup>&frasl;<sub>T<sub>w</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>V<sub>T,N</sub></sup>&frasl;<sub>B</sub> 1.5 exp <sup>-7000</sup>&frasl;<sub>T<sub>w</sub></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>V<sub>T,N</sub></sup>&frasl;<sub>B</sub> 0.5 exp <sup>-2000</sup>&frasl;<sub>T<sub>w</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>w</sub>]</td><td>√<sup>N<sub>A</sub></sup>&frasl;<sub>B</sub> V<sub>T,2D,N</sub> 0.1 exp <sup>-21000</sup>&frasl;<sub>T<sub>w</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>10<sup>8</sup> exp <sup>-20676</sup>&frasl;<sub>T<sub>w</sub></sub></td></tr></tbody></table>

In the table, the symbol s denotes a free site on the surface that can be occupied by a gas-phase element, the subscript (·)s indicates an absorbed species, while the subscript (·)b indicates a carbon surface element. B is the total active site density, kᵢ the rate coefficient of the reaction i, T<sub>w</sub> the surface temperature, m<sub>N</sub> the mass of the atomic nitrogen, V<sub>T,N</sub> a quarter of the mean thermal speed of atomic nitrogen, and V<sub>T,2D,N</sub> the mean thermal speed of the mobile adsorbed nitrogen on the surface. Symbols k<sub>B</sub>, h, and N<sub>A</sub> represent the Boltzmann and Planck constants, and Avogadro's number. The model predictions were compared both to high-pressure [4, 5] and to low-pressure [10] data and it was shown to capture their main trends.

The adsorption process (reaction 1) was modeled as an activated process to reproduce the drop of the recombination efficiency observed at low temperatures in the molecular beam experiments. The rate coefficient of the desorption process (reaction 2) was computed based on transition state theory. A desorption energy equal to the bound energy of a double bond 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 gas-phase dependent (reaction 4) and gas-phase independent (reaction 5) mechanisms were included in the model to describe the nitrogen recombination.

## B. 0D model

Molecular beam experiments can be described using the 0D model proposed by Prata et al. [9]. According to the law of mass action, the set of surface production rates read:

$$ \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}[\text{CN}] = k_3[\text{N}][\text{N}_s] + k_6[\text{N}_s], \qquad (4) $$

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

One can close such a system with the conservation of active sites, $B = [\text{N}_s] + [\text{s}]$, and the steady-state assumption, $\partial[\text{s}]/\partial t = 0$. Further inputs to the model are the surface temperature, and the gas pressure and temperature, used to compute the concentration of the atomic nitrogen. These two values were set to produce a flux close to the one observed in the molecular beam experiments ($p_{\text{beam}} = 0.024$ Pa and $T_{\text{beam}} = 1000$ K). Finally, the efficiencies of the macroscopic recombination ($N + N \rightarrow N_2$) and nitridation ($N + C_b \rightarrow \text{CN}$) reactions can be computed as:

$$ \gamma_{\text{N}_2} = \frac{2d[\text{N}_2]/dt}{\mathcal{N}_{\text{N}}}, \qquad (6) $$

$$ \gamma_{\text{CN}} = \frac{d[\text{CN}]/dt}{\mathcal{N}_{\text{N}}}, \qquad (7) $$

where $\mathcal{N}_{\text{N}}$ is the number flux of particles of nitrogen impinging on the surface.

Such a model can be also employed to estimate the behavior at high pressure by imposing the surface temperature and the gas pressure. Full details can be found in Prata et al. [9].

### C. Stagnation line model

Contrary to the molecular beam experiments, the ones in the VKI Plasmatron facility consist of a gas flow impinging on a carbon surface. Thus, to describe the BL developing in front of the sample, one has to solve the chemically reacting Navier-Stokes equations:

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

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

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

Equation (8) represents the species mass conservation: symbol $\mathbf{u}$ is the mass-averaged mixture velocity, $\rho_i$ the partial density of species $i$, $\mathbf{j}_i$ its diffusive mass flux, $\dot{\omega}_i$ its chemical production/destruction rate, and $n_s$ the number of species in the mixture. Equation (9) represents the conservation of momentum: symbol $p$ stands for the thermodynamic pressure of the mixture, $\bar{I}$ the unit tensor, and $\bar{\bar{\tau}}$ the viscous stress tensor. Finally, Equation (10) represents the conservation of energy: symbol $E = e + \frac{1}{2}\mathbf{u} \cdot \mathbf{u}$ is the total energy, sum of the thermodynamic contribution ($e$) and the kinetic one ($\frac{1}{2}\mathbf{u} \cdot \mathbf{u}$), $H = E + p/\rho$ is the total enthalpy, and $\mathbf{q}$ the total heat flux.

When dealing with axisymmetric problems, a Dimensionally Reduced Navier-Stokes Equations (DRNSE) formula-
tion [18] can be used to simulate the flowfield only along the stagnation line, which is computationally efficient. The system is first written in spherical coordinates (r, θ, φ), then, exploiting the axial symmetry, $u_φ = ∂/∂φ = 0$ and letting $θ → 0$, the DRNSE formulation is obtained. Based on such a formulation, Munafò [19, 20] developed the stagnation line code. It is based on a finite volume discretization method. An implicit local time-stepping technique ensures rapid convergence to the steady state. In this study, numerical fluxes are evaluated based on the Roe scheme [21], thermodynamics properties using the NASA polynomials [22], diffusion mass fluxes through the generalized Stefan-Maxwell equations, while viscosity and thermal conductivity are obtained based on the Chapman-Enskog perturbative solution of the Boltzmann equation [23, 24]. Chemical production rates $\dot{\omega}_i$ for species i are computed using finite-rate gas-phase chemistry, with reaction rate coefficients taken from Olynick et al. [25].

Two BCs were imposed. The inlet was characterized by imposing the BL edge conditions, while the reacting carbon surface was modeled using a surface mass balance. It reads:

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

where $\dot{\omega}_i$ and $\mathbf{j}_i$ respectively represent the chemical surface production, computed employing the ACA model, and the diffusive mass flux, $\mathbf{u}_g = \dot{\mathbf{m}}/\rho$ is the gas blowing velocity at the interface, and $\mathbf{n}$ indicates the outward normal. Symbol $\dot{\mathbf{m}}$ stands for the mass blowing rate, which can be computed as:

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

The stagnation line code is fully coupled to Mutation$^{++}$, which is called for the computation of thermodynamics, chemistry, and transport properties, as well as for solving the ablative BC.

Output of the CFD computations is the mass blowing rate, which is compared to the experimental counterpart for each case.

# IV. Uncertainty quantification methodology

In this section, we describe the UQ tools and methodologies employed to solve both the forward and the inverse problem, and build a surrogate model of the CFD computations.

## A. Uncertainty propagation

Uncertainty propagation refers to the process of characterizing the uncertainty on a QoI, $y = M(x)$, by propagating the uncertainties on the input $x$ through the model $M$. Among others, Monte Carlo (MC) [26, 27] is the most
straightforward and popular method. It requires quasi-randomly sampling $N_s$ independent realizations of the input, $\mathbf{x}^{(i)}$, according to its joint Probabilistic Density Function (PDF). Then, the PDF of $\mathbf{y}$ can be inferred through the collection of the $N_s$ independent model realizations, $\mathcal{M}(\mathbf{x}^{(i)})$. Specifically, the numerical approximations of the first two moments (mean and variance) can be computed as:

$$ \bar{\mu}_y = \frac{1}{N_s} \sum_{i=1}^{N_s} \mathcal{M}(\mathbf{x}^{(i)}), \qquad (13) $$

$$ \bar{\sigma}_y^2 = \frac{1}{N_s} \sum_{i=1}^{N_s} \left( \mathcal{M}(\mathbf{x}^{(i)}) - \bar{\mu}_y \right)^2 . \qquad (14) $$

## B. Stochastic inverse problem

The stochastic inverse problem refers to characterizing the uncertainty on some model parameters $\mathbf{x}$ given the uncertainty on some observations $\mathbf{y}$. In this work, we employed a Bayesian framework. According to Bayes' theorem, the posterior distribution of the parameters $\mathbf{x}$ that we intend to calibrate reads:

$$ \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 (15) $$

The quantity $\pi(\mathbf{x})$ is the prior distribution of $\mathbf{x}$ and represents what we know about $\mathbf{x}$ before performing the inversion, it can come from previous studies or expert judgment. The quantity $\pi(\mathbf{y}|\mathbf{x})$ is the likelihood function and expresses how likely the realization $\mathbf{y}$ is obtained by the input $\mathbf{x}$ through a numerical model $\mathbf{y} = \mathcal{M}(\mathbf{x})$. It relates the QoIs to the experimental observations. Finally, the denominator is the marginal likelihood or evidence, which ensures the distribution integrates to 1.

In our study, we assumed that I) the $N_{obs}$ experimental observations are independent, II) the model error is negligible compared to the experimental uncertainty, $\sigma_i$, on the measurement $y_i$, and III) the experimental uncertainty is normally distributed around the mean ($\mathcal{N}(0, \sigma_i^2)$). The likelihood reads:

$$ \pi(\mathbf{y}|\mathbf{x}) = \prod_{i=1}^{N_{obs}} \mathcal{N}(y_i|\mathcal{M}_i(\mathbf{x}), \sigma_i^2). \qquad (16) $$

In practical applications, analytically computing the evidence of Equation (15) is not trivial, and the posterior distribution is obtained by employing a sampling method, such as the MCMC [28]. The MCMC algorithm requires several calls to the forward model. For efficiency reasons, one generally approximates expensive numerical responses with surrogate models.
### C. Surrogate model

Although efficient, the stagnation line computations, needed to describe the Plasmatron experiments, are too costly to be directly used in a sampling strategy. For this reason, we construct a Kriging surrogate model for approximating its response.

We have a model $M_i(\mathbf{x})$, which provides the response $y_i$ (e.g. mass blowing rates) as a function of the input quantity $\mathbf{x}$ (CFD BCs). We can sample $n_t$ points ($\mathcal{X} = \{\mathbf{x}^{(1)}, ..., \mathbf{x}^{(n_t)}\}$) and compute the corresponding model realizations ($\mathcal{Y}_i = \{y_i^{(1)}, ..., y_i^{(n_t)}\}$). We can use these pairs to train a surrogate model, which can be then called to provide the response at an unseen $\mathbf{x}$ point.

The Kriging model [29–31] assumes $M_i$ to be the realization of a Gaussian process:

$$M_i(\mathbf{x}) \approx M_i^{\mathrm{K}}(\mathbf{x}) = \beta_i^{\mathrm{T}} \mathbf{f}_i(\mathbf{x}) + Z_i(\mathbf{x}). \qquad (17)$$

The first term of the RHS is the process mean, also called trend, whose unknown $\beta_i$ parameters need to be determined. The second term $Z_i$ is a stationary Gaussian process characterized by a variance $\sigma_i^2$. It is defined through a correlation function $(R(x, x'; \theta_i))$ between the points $x$ and $x'$. Parameters $\sigma_i^2$ and $\theta_i$ also need to be determined. Once the process is conditioned to the training points, one can compute the posterior mean of the process, $\hat{M}^{\mathrm{K}}(\mathbf{x})$, at any point, with confidence given by the process’ standard deviation, $\hat{\sigma}^{\mathrm{K}}(\mathbf{x})$. They read:

$$\hat{M}^{\mathrm{K}}(\mathbf{x}^{(*)}) = \beta f^{\mathrm{T}}(\mathbf{x}^{(*)}) + r^{\mathrm{T}}(\mathbf{x}^{(*)}) R^{-1} (\mathcal{Y} - \beta F), \qquad (18)$$

$$(\hat{\sigma}^{\mathrm{K}})^2(\mathbf{x}^{(*)}) = \sigma^2 \left( 1 - r^{\mathrm{T}}(\mathbf{x}^{(*)}) R^{-1} r(\mathbf{x}^{(*)}) + u^{\mathrm{T}}(\mathbf{x}^{(*)}) (F^{\mathrm{T}} R^{-1} F)^{-1} u(\mathbf{x}^{(*)}) \right) \qquad (19)$$

where $F$ is the matrix of the trends computed at the training points, $r(\mathbf{x}^{(*)})$ is the vector of cross-correlations between $\mathbf{x}^{(*)}$ and each of the training points, $R$ is the $n_t \times n_t$ correlation matrix between all the training points. Their components respectively are $F_{ij} = f_j(\mathbf{x}^{(i)}), r_i = R(\mathbf{x}^{(*)}, \mathbf{x}^{(i)}), R_{ij} = R(\mathbf{x}^{(i)}, \mathbf{x}^{(j)}).$ Finally, $u(\mathbf{x}^{(*)}) = F^{\mathrm{T}} R^{-1} r(\mathbf{x}^{(*)}) - f(\mathbf{x}^{(*)})$.
We employed an ordinary Kriging. In this case, the functions reduce to $f = 1$, and $F = 1$.

As we will see in Section V.C, one of the bottlenecks of the Kriging is the so-called ‘curse of dimensionality’: the number of training points needed to produce an accurate surrogate model increases exponentially with the dimension of the input space.

## V. Calibration procedure and results

In this section, we first present the assessment of the ACA model predictions against Plasmatron data. Then, we outline the procedure that we employed to robustly calibrate the rate parameters of the ACA model and we present the
results.

A. Assessment of ACA model predictions against Plasmatron experiments

We assessed the predictive capability of the ACA model against the four Plasmatron experiments presented in
Section II.B.

Stagnation line computations were performed to describe the experiments, allowing us to obtain the numerical mass blowing rate. Input of the simulations are the surface temperature, $T_w$, and the BL edge velocity, $u_e$, pressure $p_e$, and temperature, $T_e$. Both the surface temperature and BL edge pressure, assumed equal to the chamber one, are experimentally accessible. The BL edge velocity can be obtained from the experimental dynamic pressure ($u_e = \sqrt{2p_d/K_H\rho_e}$, where $K_H$ is the Homman's correction), provided that the BL edge density, $\rho_e$, is known. Assuming thermo-chemical equilibrium, the latter can be computed through the ideal gas law, imposing the known BL edge pressure and the BL edge temperature, which is not accessible from an experimental standpoint.

The numerical uncertainties on *ṁ* for the four different experimental conditions were obtained as follows:

1) We propagated the experimental uncertainties on $T_w$, $p_d$, and $p_e$ reported in Table 1.

2) We prescribed a large uniform uncertainty (± 1000 K) to the BL edge temperature around the value rebuilt by Helber et al. [4], as it was considered less trustworthy. In fact, this value was obtained through a coupled numerical procedure that involved several model assumptions. This uncertainty was also propagated.

3) We described the chemical behavior of the carbon surface employing the ACA model. We used the reaction rate parameters proposed by Prata et al. [9], without prescribing any model uncertainty.

A total of 200 simulations were performed for each condition by randomly sampling the input from the stochastic space. The propagated uncertainties are plotted in Figure 1, along with the experimental counterparts, against the surface temperature. As one can see, the predicted mass blowing rates are lower than what was experimentally observed. Furthermore, the numerical predictions fail to capture the upward trend of $\dot{m}$ with respect to the surface temperature.

We note that the ACA model was originally calibrated using these experimental data. However, the model predictions were not directly compared to experimental mass blowing rates, but with the values of the nitridation probabilities derived by Helber et al. [4]. The latter were computed through a coupled numerical-experimental methodology, assuming no surface recombination of the nitrogen atoms. By contrast, the ACA model predicts a high recombination probability (see Figure 2) in agreement with what was observed in the experiments performed by Murray et al. [10] and Lutz [5]. Furthermore, the analysis perfomed by Prata et al. [9] to compare ACA results to the Plasmatron experiments made strong assumptions on the modeling of the BL developing in front of the sample. Specifically, the analysis did not account for the reduction of atomic nitrogen at the surface due to surface reactions. As a result, although the
<figure><img src="image_2.png" /><figcaption>Fig. 1 Uncertainties on the stagnation-point mass blowing rates for the conditions in Table 1. Black and red error bars: experimental and numerical uncertainties.</figcaption></figure>

equivalent nitridation probabilities are in fair agreement with the values proposed by Helber et al. [4] (see Figure 2) the high recombination probability reduces the atomic nitrogen available at the surface for nitridation, and hence the mass blowing rate ($\dot{m} \propto \gamma_{CN}n_N$).

<figure><img src="image_3.png" /><figcaption>Fig. 2 Uncertainties on the nitridation and recombination probabilities for the conditions in Table 1. Numerical-experimental nitridation values rebuilt by Helber et al. [4] with black error bars. Numerical uncertainties obtained in this work with error bars: nitridation probabilities in red, recombination in blue.</figcaption></figure>
## B. Calibration procedure

The results obtained in the previous section highlighted the necessity of I) directly exploiting pure experimental observations to calibrate the model (for example, the blowing mass rates), and II) employing an accurate model to reproduce the near-surface state in the Plasmatron experiments. Furthermore, the ACA model was constructed deterministically, without accounting for experimental and parametric uncertainties. To improve the robustness of the predictions, we employed a Bayesian formalism to compute the posterior distributions of the ACA reaction rate parameters.

We chose to learn four elementary reaction efficiencies ($X_3, X_4, X_5, X_6$), and the total active site density ($B$) of the ACA model; the subscript of the efficiencies denotes the index of the reaction in Table 2. Additionally, although it is common practice to consider uncertain only the pre-exponential coefficients of the reaction rates when dealing with homogeneous FRC [32, 33], we decided to infer also five activation energies ($E_1, E_3, E_4, E_5, E_6$) as different values can be found in literature (see Table 3). Hence, the vector of parameters that we intend to learn is:

$$ \mathbf{x} = [E_1, E_3, E_4, E_5, E_6, X_3, X_4, X_5, X_6, B]. $$

The prior distributions on $\mathbf{x}$ are reported in Table 3: non-informative uniform distributions over a wide range of values are chosen to encompass the data available in the literature. Note that we prescribed log-uniform distributions to the elementary reaction efficiencies and the total active site density because their priors span several orders of magnitude.

<table><caption>**Table 3** Non-informative uniform prior distributions on the parameters to calibrate, along with the values proposed by Prata et al. [9] and Nieman et al. [34].</caption><thead><tr><th>Parameter</th><th>min</th><th>max</th><th>Prata</th><th>Nieman</th></tr></thead><tbody><tr><td>E<sub>1</sub>, K</td><td>0</td><td>3000</td><td>2500</td><td>0, 2436, and 4525*</td></tr><tr><td>E<sub>3</sub>, K</td><td>0</td><td>15000</td><td>7000</td><td>10676</td></tr><tr><td>E<sub>4</sub>, K</td><td>0</td><td>15000</td><td>2000</td><td>0</td></tr><tr><td>E<sub>5</sub>, K</td><td>0</td><td>35000</td><td>21000</td><td>6150</td></tr><tr><td>E<sub>6</sub>, K</td><td>15000</td><td>35000</td><td>20676</td><td>30286</td></tr><tr><td>log<sub>10</sub>(X<sub>3</sub>)</td><td>-4</td><td>2</td><td>0.176</td><td>-</td></tr><tr><td>log<sub>10</sub>(X<sub>4</sub>)</td><td>-4</td><td>0</td><td>-0.301</td><td>-</td></tr><tr><td>log<sub>10</sub>(X<sub>5</sub>)</td><td>-4</td><td>0</td><td>-1</td><td>-</td></tr><tr><td>log<sub>10</sub>(X<sub>6</sub>)</td><td>5</td><td>12</td><td>8</td><td>-</td></tr><tr><td>log<sub>10</sub>(B)</td><td>-7</td><td>-3</td><td>-5</td><td>-</td></tr></tbody></table>

The first data employed for the calibration comprises the efficiencies estimated by Murray et al. [10] in the molecular

<footnote>*The three values correspond to three different adsorption sites.</footnote>
beam-surface scattering experiments:

$$ \mathbf{y}^{\text{MB}} = [\gamma_{\text{CN}}, \gamma_{\text{N2}}]_i \quad i = 1, \dots, 11. $$

We divided this group into two sets. The first one does not include the data obtained for $T_w < 1000$ K (hereafter referred to as 'set A'). It is the same subset used by Prata et al. [9] to construct the ACA model. The second set includes all the available points (hereafter referred to as 'set B'). For the set A, the likelihood reads:

$$ \pi^{\text{MB,A}}(\mathbf{y}^{\text{MB}}|\mathbf{x}) = \prod_{i=5}^{22} \mathcal{N}(y_i^{\text{MB}}|\mathcal{M}_i^{\text{MB}}(\mathbf{x}), (\sigma_i^{\text{MB}})^2), \quad (20) $$

while, for the set B, it is:

$$ \pi^{\text{MB,B}}(\mathbf{y}^{\text{MB}}|\mathbf{x}) = \prod_{i=1}^{22} \mathcal{N}(y_i^{\text{MB}}|\mathcal{M}_i^{\text{MB}}(\mathbf{x}), (\sigma_i^{\text{MB}})^2). \quad (21) $$

The symbol $\mathcal{M}_i^{\text{MB}}(\mathbf{x})$ stands for the 0D model, described in Section III.B, which expresses $y_i^{\text{MB}}$ as a function of $\mathbf{x}$. The symbol $\sigma_i^{\text{MB}}$ represents the experimental uncertainty on each molecular beam measurement. Since the 0D model is analytical, no model error was considered. Along with the parameters $\mathbf{x}$, the surface temperature and the impinging nitrogen flux are inputs of the model. No uncertainty was considered on these quantities.

The second group comprises the blowing mass rates observed in the Plasmatron experiments:

$$ \mathbf{y}^{\text{ICP}} = \dot{m}_i \quad i = 1, \dots, 4. $$

In this case, the likelihood reads:

$$ \pi^{\text{ICP}}(\mathbf{y}^{\text{ICP}}|\mathbf{x}) = \prod_{i=1}^{4} \mathcal{N}(y_i^{\text{ICP}}|\mathcal{M}_i^{\text{ICP}}(\mathbf{x}), (\sigma_i^{\text{ICP}})^2). \quad (22) $$

The symbol $\mathcal{M}_i^{\text{ICP}}(\mathbf{x})$ stands for the surrogate model, trained on stagnation-line CFD computations, which expresses $y_i^{\text{ICP}}$ as a function of $\mathbf{x}$ (see Section V.C). The symbol $\sigma_i^{\text{ICP}}$ represents the experimental uncertainty on each Plasmatron measurement. The surrogate model is characterized by an error around the 1.6% (see Section V.C), much below the experimental uncertainty (ranging from 20% to 35%). For this reason, the model uncertainty was not accounted for during the inference. The settings of the simulations were discussed in Section V.A. Along with the ACA model rate coefficients, the surface temperature and the conditions at the edge of the BL are inputs of the model.

A total of four different calibrations were performed. They are summarized in Table 4, along with the likelihoods
used to solve the Bayesian inverse problem. In the first two tests, we used only molecular beam-surface scattering data to calibrate the model. The Plasmatron data were employed to assess the predictive capability of the model at high pressure. In the third test, we calibrated the model only on Plasmatron data to verify that such data can be explained through the ACA model. Finally, in the last test, the calibration was performed imposing both sets of observations.

<table><caption>**Table 4 Matrix of the inversions performed.**</caption><thead><tr><th>Test</th><th>Molecular beam</th><th>Plasmatron</th><th>Likelihood</th></tr></thead><tbody><tr><td>A</td><td>Calibration</td><td>Assessment</td><td>π<sup>MB,A</sup></td></tr><tr><td>B</td><td>Calibration</td><td>Assessment</td><td>π<sup>MB,B</sup></td></tr><tr><td>C</td><td>-</td><td>Calibration</td><td>π<sup>ICP</sup></td></tr><tr><td>D</td><td>Calibration</td><td>Calibration</td><td>π<sup>MB,B</sup> · π<sup>ICP</sup></td></tr></tbody></table>

The quantities to calibrate were tightly correlated in a high-dimensional space. For this reason, the posterior distributions were obtained by constructing an MCMC using the Affine Invariant Ensemble Sample (AIES) algorithm [35]. The UQLab software [36] was employed both to construct the MCMC chains and to train the surrogate model.

### C. Surrogate model construction

A Kriging model was constructed to accelerate the evaluation of the mass blowing rate as a function of both the ACA model and the input conditions. First, we attempted to create a global surrogate model able to approximate the CFD response for the four experiments. Thus, in addition to the 10-dimensional $x$ input, the model was parameterized also according to the surface temperature, $T_w$, and the temperature, $T_e$, velocity, $u_e$ and pressure, $p_e$, at the edge of the BL. The resulting model

$$ \mathcal{M}_i^{\text{ICP}}(\mathbf{x}) \approx \mathcal{M}^{\text{K, ICP}}(\mathbf{x}, T_{w,i}, T_{e,i}, u_{e,i}, p_{e,i}) $$

has a total input dimension of 14. The training space was chosen large enough to accommodate both the prior on $x$ and the four experimental conditions in Table 1, along with the relative uncertainties. Such a space was normalized to simplify the construction of the surrogate model. Its accuracy was evaluated by computing a normalized error based on 150 Latin Hypercube Sampling (LHS) validation points. Such an error reads:

$$ NRMSE = \sqrt{\frac{1}{N_v} \sum_{i=0}^{N_v} \frac{(Y_{v,i} - \hat{Y}_i)^2}{Y_{v,i}^2}}, \quad (23) $$

where $\hat{Y}_i = \mathcal{M}^{\text{K, ICP}}(\mathbf{x}, T_{w,i}, T_{e,i}, u_{e,i}, p_{e,i})$ is the Kriging model evaluation at the validation point $i$, while $Y_{v,i}$ is the respective CFD response. The large input dimension, along with the large priors prescribed, deteriorates the rate of convergence of the model. In fact, as one can see in Figure 3, the error drops from around 7%, when 200 LHS training points are employed, to around the 2.3% with 6200 LHS points. Such a surrogate model was not very efficient to be
used coupled to the sequential AIES algorithm.

<figure><img src="image_4.png" /><figcaption>**Fig. 3** NRMSE of the global surrogate model built on 14 dimensions with dashed red line. NRMSE of the 6-dimensional G6 surrogate model with solid black line.</figcaption></figure>

For this reason, we employed it to perform a Sobol analysis, aiming at identifying the most sensitive variables and reducing the input dimension. The resulting total Sobol indices are shown in Figure 4. The quantities at the edge of the BL show poor influence on the blowing mass rate. Similarly, the surface temperature exhibits a small influence. It is important to note that such an index was computed on the entire experimental range (≈ 600 K), while the experimental uncertainty on each test is only around 20 K. Hence, the influence of the uncertainty on the surface temperature for each single experiment can be considered negligible. A similar behavior was also observed by del Val et al. [12].

<figure><img src="image_5.png" /><figcaption>**Fig. 4** Total Sobol indices for the mass blowing rates.</figcaption></figure>
Based on this analysis, we decided to build a surrogate model for each test, neglecting the uncertainties on the surface temperature and on the parameters at the edge of the BL. This reduced the input dimension to 10. Further reduction in dimension was achieved by performing a change of variables. In fact, the reaction rates are in the form:

$$k_i \propto X_i \exp\left(-\frac{E_i}{T_w}\right), \qquad (24)$$

allowing us to introduce the variable:

$$P_i = X_i \exp\left(-\frac{E_i}{T_{w,j}}\right) \quad i = 3, 4, 5, 6, \qquad (25)$$

where $T_{w,j}$ is the surface temperature for the four different conditions ($T_w = [T_{w,G4}, T_{w,G5}, T_{w,G6}, T_{w,G7}$]). Such surface temperature is considered constant in each surrogate model. The resulting surrogate model:

$$M_j^{\text{ICP}}(\mathbf{x}) \approx M_j^{\text{K, ICP}}(P_3, P_4, P_5, P_6, E_1, B)$$

has only an input dimension of 6. The convergence of the model built on the G6 experiment was assessed by means of 150 LHS points. It is shown in Figure 3 where the error drops from around 5% when 200 LHS points are used, to around 1.2% with 5200 LHS points. For efficiency reasons, we decided to use the surrogate model built on 3000 points with a NRMSE = 1.6%.

A bridging function was used during the inference process to transform the points sampled by the MCMC algorithm ($X_3, E_3, X_4, E_4, X_5, E_5, X_6, E_6$) to the corresponding reduced variables ($P_3, P_4, P_5, P_6$) employing the surface temperature of each experiment.

## D. Results of the calibration

In this section, we describe the results obtained by calibrating the model employing the four sets of observations presented in Table 4. For all the conditions, 100 MCMC chains were constructed with an AIES algorithm. For each of them, we sampled enough points to have, after the burn-in process, a low multivariate potential scale reduction factor from the Gelmen-Rubin diagnostic ($R_A = 1.1, R_B = 1.1, R_C = 2.5, R_D = 1.2$). For the case C, which shows a slightly high value of $R$, we assessed the chain convergence by visually inspecting the posterior distributions. They exhibited steady behavior. In the following, we first present the marginal distributions of the calibrated parameters, then the posterior predictions of the model.
1. *Posterior marginal distributions*

The marginal distributions of the calibrated parameters for the two nitridation reactions are plotted in Figures 5(a), 5(b), 6(a), and 6(b), while Figures 5(c), and 6(c) show the two joint distributions relating the reaction efficiency and the activation energy of the selected reactions. The values reported by Prata et al. [9] and by Nieman et al. [34] are also shown in the same figures.

<figure><img src="image_6.png" /><figcaption>**Fig. 5** Probability densities of efficiency and activation energy of reaction 3 (N + N<sub>s</sub> + C<sub>b</sub> → CN + N + s). Blue histograms refer to the parameters posterior using set A, black histograms using set B, gray histograms using set C, purple histograms using set D, red lines indicate the prior distributions, green lines the values obtained by Prata et al. [9], and orange lines the values by Nieman et al. [34].</figcaption></figure>

As one can see in Figures 5 and 6, both nitridation reactions are well characterized when the model is calibrated on
molecular beam data (set A and B). This is especially true when using set B, as it contains information also at low
temperatures. On the other hand, a poor characterization of the parameters is observed when the model is calibrated
using only the Plasmatron data (set C). However, when the Plasmatron data are combined with the molecular beam ones
(set D) an improved characterization is achieved, especially for reaction 3 (see Figure 5), whose parameters' distribution
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Probability densities of efficiency and activation energy of reaction 6 (N<sub>s</sub> + C<sub>b</sub> → CN + s). Same legend convention as in Figure 5.</figcaption></figure>

presents more pronounced peaks. It is consistent as reaction 3 was included in the model to capture the ablative behavior
at high pressures.

In Figure 5(b) we see that when the model is calibrated by means of set A, the posterior distribution of the activation energy of reaction 3 includes both the values proposed by Prata et al. [9] and by Nieman et al. [34], which makes them consistent with the data employed in the calibration. However, when low-temperature experimental efficiencies are included (set B), all the distributions shift towards lower values. As we will discuss later, this allows the model to predict the increasing trend of the nitridation efficiency observed at low temperatures in the molecular beam experiments (see Figures 10 and 11). In this case, the distribution of the calibrated activation energies departs from the values proposed by Prata et al. [9] and by Nieman et al. [34], suggesting some incompatibility of the latter in explaining the molecular beam nitridation efficiencies at low temperatures through the proposed model. The same holds true when the model is calibrated with set D.
The marginal distributions of the calibrated parameters for the two recombination reactions are shown in Figures 7(a), 7(b), 8(a), and 8(b), while in Figures 7(c), and 8(c) the two joint distributions for the reaction efficiency and the activation energy are given.

As one can see in Figure 7, the gas-phase dependent recombination reaction (reaction 4, $N_s + N \rightarrow N_2 + s$) appears to be poorly characterized when the molecular beam and the Plasmatron data are used separately. However, when these data are used jointly (set D), the reaction 4 is predicted to be unlikely. In fact, the joint distribution moves toward low values of the efficiency and high values of the activation energy (see Figure 7(c)).

<figure><img src="image_8.png" /><figcaption>**Fig. 7** Probability densities of the efficiency and the activation energy of reaction 4 ($N + N_s \rightarrow N_2 + s$). Same legend convention as in Figure 5.</figcaption></figure>

On the other hand, the gas-phase independent recombination reaction (reaction 5, $N_s + N_s \rightarrow N_2 + 2s$) presents more characterized distributions (see in Figure 8). Thus, it appears to be the dominant recombination mechanism, in agreement with what was observed by Murray et al. [8]. The distribution for the activation energy of the reaction 5 has a peak very close to the value proposed by Prata et al. [9] when set A is used. The peak shifts toward a lower value
employing set B and set D, closer to the value computed by Nieman et al. [34]. Also in this case, the parameters show a
poor characterization when using set C.

<figure><img src="image_9.png" /><figcaption>**Fig. 8 Probability densities of efficiency and activation energy of reaction 5 (N<sub>s</sub> + N<sub>s</sub> → N<sub>2</sub> + 2s). Same legend convention as in Figure 5.**</figcaption></figure>

The activation energy of the absorption reaction, see Figure 9(a), is also well characterized, with a distribution
tending towards small values, close to the lowest value proposed by Nieman et al. [34]. Interestingly, the more the
model is refined the more this value approaches zero. Such a low value of activation energy allows the model to predict
the very high values of recombination efficiencies experimentally observed in the molecular beam experiments.

Finally, the total active site density, see Figure 9(b), has a broad distribution with a peak close to the value employed by Prata when employing the set A. It moves towards small values using sets B and D, while is poorly characterized with set C, although tending to high values. It is important to note that Prata et al. [9] chose the value of the total active site density as a reference constant and proceeded to fit all other ACA parameters with this constant fixed. It is not clear what effects changing the value of the total active site density would have on the physical accuracy of the ACA model.
Further modeling research would be required to understand this effect.

<figure><img src="image_10.png" /><figcaption>**Fig. 9** Probability densities of activation energy of reaction 1 (N + s → Ns) and total active site density. Same legend convention as in Figure 5. The three different values by Nieman et al. [34] corresponding to three distinct adsorption sites.</figcaption></figure>

## 2. Model posterior predictions

### I) From molecular beam data

The posterior distributions obtained by calibrating the model using set A and set B were propagated through the 0D model to obtain the efficiencies posterior predictions both at low and high pressure.

The model posterior predictions from set A are plotted in Figures 10(a) and 10(b). The molecular beam experimental uncertainties [8], and the marginals on the nitridation reaction efficiency obtained by del Val et al. [13] are shown in the figures. The latter were obtained by directly calibrating these efficiencies employing the same Plasmatron experiments. The predictions of the calibrated model are in perfect agreement with the molecular beam points used for the calibration. Large model uncertainties are encountered in the predictions at high pressure as no high-pressure calibration point was used. This suggests the necessity for higher-pressure experimental points to reduce this uncertainty. The same can be achieved 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 11(b), is drastically reduced when the model is calibrated with set B. Furthermore, the nitridation efficiencies at such high pressure exhibit excellent agreement with the marginals from del Val et al. [13], suggesting that the two experiments can be accurately described through the use of the same calibrated model.

This improved characterization can be explained by looking at Figure 12, where the surface coverage at the two different pressures is plotted as a function of the surface temperature. As one can see, the value of the surface coverage at

<footnote>†The surface coverage is defined as the ratio between the number of adsorbed atoms/molecules on a surface and the total number of active sites.</footnote>
<figure><img src="image_11.png" /><figcaption>Fig. 10 Predictions of the model calibrated with set A at pressures representative of molecular beam (left) and Plasmatron (right) experiments. The red and blue color indicates the recombination and the nitridation efficiencies. The solid lines represent the mean of the prediction. The shadow areas correspond to the prediction uncertainty (95% Confidence Interval (CI)). The error bars with square as mean represent the molecular beam from Murray et al. [8] (last 9 points used for the calibration), while the error bars with circle as mean are the marginals obtained by del Val et al. [13] (not used for the calibration).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11 Predictions of the model calibrated with set B at pressures representative of molecular beam (left) and Plasmatron (right) experiments. Same legend convention as in Figure 10. Molecular beam data from Murray et al. [8] used for the calibration.</figcaption></figure>

low-pressure/low-temperature is close to the almost flat value predicted at high pressure. Furthermore, the experimental uncertainty on the nitridation reaction efficiency at low temperatures is small. Thus, the low-temperature points contain information on the model behavior at high surface coverage.

On the other hand, the recombination efficiency at high pressure is still characterized by large variation as the relative experimental uncertainty at low-pressure/low-temperatures is high. The inclusion in the calibration of more accurate measurements at low temperatures would certainly help in decreasing the model uncertainty at high pressure.
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. [8]. After a plateau, it decreases around 2250 K because the surface coverage sharply decreases at this temperature, as shown in Figure 12. 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 decreasing surface coverage. At higher pressure, an almost flat trend in the nitridation probability is predicted, due to the nearly constant surface coverage.

<figure><img src="image_13.png" /><figcaption>Fig. 12 Model predictions of the surface coverage using the set B for the calibration. Solid line refers to the mean predictions at molecular beam pressure and the dashed line at Plasmatron pressure. The uncertainty corresponding to the 95% CI is shown with shadow area.</figcaption></figure>

The parameters' posteriors were also propagated through the surrogate models to obtain the uncertainties on the stagnation-point mass blowing rate for the cases G4, G5, G6, and G7. Their values, along with the experimental uncertainties, are shown in Figures 13(a) and 13(b). It can be seen that the propagated uncertainties are large when the model is calibrated using set A. They dramatically decrease when the calibration is performed by means of the set B. Furthermore, such distributions share much of the support with their experimental counterparts, which is interesting as no Plasmatron data was used in the inference. However, the calibrated model appears to be less sensitive to the surface temperature (which increases from G4 to G7) than what was experimentally observed.

## II) From Plasmatron data

Although Plasmatron data poorly inform the parameters, such an analysis allowed us to assess whether the model is able to correctly capture the experimental data. The posterior distribution of the model parameters calibrated with set
<figure><img src="image_14.png" /><figcaption>**Fig. 13** Mass blowing rate distributions for the conditions in Table 1. The continuous lines refer to the experimental uncertainty, while the histograms show the propagated distributions. Blue, gray, red, and green colors indicate the G4, G5, G6, and G7 case.</figcaption></figure>

C was propagated through the four surrogate models to obtain the posterior predictions of the stagnation-point mass blowing rate for the cases G4, G5, G6, and G7. These are shown in Figure 14, along with the experimental uncertainties. Overall, the propagated uncertainties are in good agreement with their experimental counterparts, and the upward trend with respect to the surface temperature (which increases from G4 to G7) is observed. Only for the G5 case, the mean prediction is slightly higher than the experimental one. This analysis verified that the model is capable of well explaining the Plasmatron experimental data.

### III) *From molecular beam and Plasmatron data*

The parameters posteriors obtained calibrating the model with set D were propagated through the 0D model. The posterior predictions are shown in Figure 15. While the main behavior of the model is similar to the one obtained when the calibration was performed with set B, the uncertainty at low pressures is slightly reduced at low and high temperatures. The uncertainties in both the nitridation and the recombination reactions at high pressure exhibit a larger reduction, as a natural consequence of the inclusion of the high-pressure Plasmatron data.

The parameters' posterior was also propagated through the four surrogate models representing the case G4, G5, G6, and G7. The results are shown in Figure 16. As one can see, the posterior predictions of the stagnation-point mass blowing rate are contained in the relative experimental uncertainties. However, the predicted upward trend with respect to the surface temperature (which increases from G4 to G7) is negligible.

This behavior can be explained by looking at Figures 17(a) and 17(b), where the specific contributions of reaction 3 and reaction 6 to the overall nitridation reaction are shown along with their uncertainties. They were computed using
<figure><img src="image_15.png" /><figcaption>Fig. 14 Mass blowing rate distributions for the conditions in Table 1 obtained propagating the parameters posteriors calibrated with set C. Same legend convention as in Figure 13.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15 Predictions of the model calibrated with set D at pressures representative of molecular beam (left) and Plasmatron (right) experiments. Same legend convention as in Figure 10. Molecular beam data from Murray et al. [8] used for the calibration. Mass blowing rates from the same experimental campaign employed by del Val et al. [13] are used for the calibration.</figcaption></figure>

the 0D analytical model, imposing a pressure of 1600 Pa.

The results obtained by propagating the PDFs relative to the calibration performed using set C (only Plasmatron data) are shown in Figure 17(a). At high surface temperatures reaction 6 is dominant over reaction 3. As the magnitude of reaction 6 increases with the surface temperature, the upward trend of $\dot{m}$ is also observed. On the other hand, such
<figure><img src="image_17.png" /><figcaption>Fig. 16 Mass blowing rate distributions for conditions reported in Table 1 obtained propagating the parameters posteriors calibrated with set D. Same legend convention as in Figurw 13.</figcaption></figure>

dominant behavior of reaction 6 appears to be not compatible with molecular beam data. It is evident in Figure 17(b), where the posterior predictions obtained by propagating the PDFs from set D (Plasmatron + molecular beam data) are plotted. In this case, the largest contribution comes from reaction 3. Such a reaction is roughly dependent on $\exp(-E_3/T_w)/\sqrt{T_w}$, with $E_3 \approx 1400$. It exhibits a flat behavior at surface temperatures close to the ones measured in the Plasmatron, explaining the flat behavior observed in the mass blowing rate prediction.

Since the model cannot predict such an upward behavior when learning also from the molecular beam data, the posterior predictions converge in a reduced part of the experimental uncertainty, where all four measurements agree with each other.

### 3. Joint distributions of the reaction probabilities

Finally, 200 points were sampled from the posterior distributions relative to the calibrations performed using set B, C, and D. They were propagated directly through the CFD model, for the conditions G4, G5, G6, and G7. The numerical results were post-processed, using the Equation 6 and 7, to obtain the equivalent probability of the nitridation and the recombination reactions. The posterior joint distributions ($\gamma_{CN}, \gamma_{N2}$) for the three different calibrations and for the four different cases are shown in Figure 18. On the same plots, the joint distributions obtained by del Val et al. [13] are shown. These were obtained by performing the calibration directly on these efficiencies.

First of all, we can notice the excellent agreement between the joint distributions from del Val et al. [13] and the ones
<figure><img src="image_18.png" /><figcaption>**Fig. 17** Uncertainties on the contributions of reactions 3 (solid blue line) and 6 (dashed red line) to the nitridation reaction at high pressures.</figcaption></figure>

corresponding to the calibration using the set C. It is consistent because, although the two calibrated models are different, they attempt to predict the same behavior on the same calibration points. The joint distributions from the calibration on set B share part of the support with the distribution from del Val et al. [13]. It is an expected result, from what was discussed about the posterior predictions on the mass blowing rates, that are almost in agreement with the experimental values. Lastly, the joint distributions from the calibration on set D are included in the joint distribution from del Val et al. [13].

Interestingly, as a result of calibrating the model using jointly the Plasmatron and the molecular beam data, we were able to reduce the posterior uncertainties on the two efficiencies with respect to the results by del Val et al. [13]. The largest reduction was observed in the recombination reaction, whose value is predicted to not exceed 0.1.

## VI. Conclusion

This paper demonstrates a new paradigm in multi-physics ablation research where uncertainty quantification techniques are applied to multiple types of experiments in order to both inform model parameters and assess the uncertainty in model predictions. The proposed methodology is applied to infer the rates of the reactions involving atomic nitrogen of the state-of-the-art ACA model. First, we performed numerical simulations of a Plasmatron experiment. The ablative behavior of the sample was described by means of the ACA model, which turned out to predict lower stagnation-point mass blowing rates than what was experimentally observed.

For this reason, we calibrated the model using a Bayesian framework. Based on preliminary investigations, we decided to characterize the pre-exponential factors and the activation energies of selected reactions, and the total active site density. Four different calibrations were performed. The first two employed the molecular beam data for the inference and the Plasmatron data to assess the model predictive capability at high pressure. The third inference relied only on
<figure><img src="image_19.png" /><figcaption>Fig. 18 Joint distribution of the reaction probabilities.</figcaption></figure>

the Plasmatron data, while the last relied on all the observations jointly. Regarding the numerical models, we used the 0D reactor proposed by Prata et al. [9] to explain the molecular beam data, while a surrogate model, built on CFD simulations equipped with the specific ablative BC, was employed to describe the Plasmatron experiment. To further increase the efficiency of the method, the surrogate model was built for each experimental point individually, after having performed an appropriate dimension reduction of the input space.

Overall, molecular beam data were found to be highly informative for calibrating the FRC model. In contrast, Plasmatron data appeared to be inadequate to infer such a complex model. However, the inclusion of these data in the calibration process was important to provide information at pressures typical of ablation applications. Specifically, both the nitridation reactions turned out to be well characterized with the data at hand, especially when exploiting together information from molecular beam and Plasmatron experiments. Regarding the recombination reactions, the gas-phase independent mechanism appeared to be the dominant mechanism, especially when the Plasmatron data were combined with the molecular beam ones. In fact, the analysis suggested the gas-phase dependent mechanism to be unlikely. Also
the activation energy of the adsorption process was found to be well characterized. Its distribution moved towards the lowest value suggested by Nieman et al. [34] as the observation data set is enriched. Finally, the probability of the total number of active sites also tended towards low values, especially when including low-temperature points from the molecular beam experiments.

The parameters posteriors obtained for each set of calibration were then propagated through both the 0D model and the surrogate models representing the four conditions in the Plasmatron experiment. Interestingly, when the model is calibrated using all the molecular beam points, the propagated mass blowing rates share much of the support with their experimental counterparts. This improved calibration was explained by the fact that high-pressure and low-pressure/low-temperature conditions are characterized by a very similar surface coverage. Thus, the molecular beam points at low pressure and low temperature are highly informative for the high-pressure behavior. When the model was calibrated only on Plasmatron data, the propagated uncertainty on the mass blowing rates compared almost perfectly with the experimental uncertainties. This suggested that the FRC can well explain the Plasmatron data. When the model was calibrated using both molecular beam and Plasmatron data a good agreement was found with all the measurements used for the calibration. However, the upward trend of the mass blowing rate with respect to the surface temperature appeared to be much less pronounced than what was experimentally observed.

Finally, we compared the joint distributions of the post-processed phenomenological efficiencies with the values obtained by del Val et al. [13]. Excellent agreement was found when the model was calibrated only on Plasmatron data. When using all the measurements, the joint distributions are a subset of the ones obtained by del Val et al. [13]. Interestingly, the use of a higher-fidelity ablation model and the consequent possibility to calibrate it on different experiments, allowed us to reduce the posterior uncertainties obtained by del Val et al. at high pressures.

We remind that in the last analysis, we used all the data to infer the model, leaving no experimental point for validation. Future perspectives include the assessment of the calibrated model including new Plasmatron data obtained at different pressures. Another interesting direction would be to calibrate the whole ACA model, including those reactions involving atomic oxygen.

## Acknowledgements

The research of Capriati and Magin was supported by the Air Force Office of Scientific Research (AFOSR) under the grant FA8655-23-1-7021. Schwartzentruber was supported by a Space Technology Research Institutes grant from NASA's Space Technology Research Grants Program. The authors also thank George Bellas for his assistance in implementing the ACA model in Mutation++.
References

[1] Park, C., and Bogdanoff, D., "Shock-Tube Measurement of Nitridation Coefficient of Solid Carbon," *Journal of Thermophysics and Heat Transfer*, Vol. 20, No. 3, 2006, pp. 487-492. https://doi.org/10.2514/1.15743.

[2] Zhang, L., Pejaković, D., Marschall, J., Dougherty, M., and Fletcher, D., "Laboratory Investigation of the Active Nitridation of Graphite by Atomic Nitrogen," *Journal of Thermophysics and Heat Transfer*, Vol. 26, No. 1, 2012, pp. 10-21.
https://doi.org/10.2514/1.T3612.

[3] Vignoles, G., Turchi, A., Bianchi, D., Blaineau, P., Lamboley, X., Le Quang Huy, D., Levet, C., Caty, O., and Chazot, O., "Ablative and catalytic behavior of carbon-based porous thermal protection materials in nitrogen plasmas," *Carbon*, Vol. 134, 2018, pp. 376-390. https://doi.org/10.1016/j.carbon.2018.03.087.

[4] Helber, B., Turchi, A., and Magin, T., "Determination of active nitridation reaction efficiency of graphite in inductively coupled plasma flows," *Carbon*, Vol. 125, 2017, pp. 582-594. https://doi.org/10.1016/j.carbon.2017.09.081.

[5] Lutz, A., "Experimental investigation and analysis of high-enthalpy nitrogen flow over graphite," Ph.D. thesis, University of Vermont, 2015.

[6] Suzuki, T., Fujita, K., and Sakai, T., "Graphite Nitridation in Lower Surface Temperature Regime," *Journal of Thermophysics and Heat Transfer*, Vol. 24, No. 1, 2010, pp. 212-215. https://doi.org/10.2514/1.43265.

[7] Suzuki, T., Fujita, K., Ando, K., and Sakai, T., "Experimental Study of Graphite Ablation in Nitrogen Flow," *Journal of Thermophysics and Heat Transfer*, Vol. 22, No. 3, 2008, pp. 382-389. https://doi.org/10.2514/1.35082.

[8] Murray, V., Recio, P., Caracciolo, A., Miossec, C., Balucani, N., Casavecchia, P., and Minton, T., "Oxidation and nitridation of vitreous carbon at high temperatures," *Carbon*, Vol. 167, 2020, pp. 388-402. https://doi.org/10.1016/j.carbon.2020.05.076.

[9] Prata, K., Schwartzentruber, T., and Minton, T., "Air-Carbon Ablation Model for Hypersonic Flight from Molecular-Beam Data," *AIAA Journal*, Vol. 60, No. 2, 2022, pp. 627-640. https://doi.org/10.2514/1.J060516.

[10] Murray, V., and Minton, T., "Gas-surface interactions of atomic nitrogen with vitreous carbon," *Carbon*, Vol. 150, 2019, pp. 85-92. https://doi.org/10.1016/j.carbon.2019.04.117.

[11] Upadhyay, R., Miki, K., Ezekoye, O., and Marschall, J., "Uncertainty quantification of a graphite nitridation experiment using a Bayesian approach," *Experimental Thermal and Fluid Science*, Vol. 35, No. 8, 2011, pp. 1588-1599. https://doi.org/10.1016/j.expthermflusci.2011.07.010.

[12] del Val, A., Le Maître, O., Congedo, P., and Magin, T., "Stochastic calibration of a carbon nitridation model from plasma wind tunnel experiments using a Bayesian formulation," *Carbon*, Vol. 200, 2022, pp. 199-214. https://doi.org/10.1016/j.carbon.2022.07.069.
[13] del Val, A., Magin, T., and Congedo, P., "Quantification of model-form uncertainties affecting the calibration of a carbon nitridation model by means of Bayesian Model Averaging," *International Journal of Heat and Mass Transfer*, Vol. 213, 2023, p. 124271. https://doi.org/10.1016/j.ijheatmasstransfer.2023.124271.

[14] Kolesnikov, A., "Conditions of simulation of stagnation point heat transfer from a high-enthalpy flow," *Fluid Dynamics*, Vol. 28, 1993, pp. 131-137. https://doi.org/10.1007/BF01055676.

[15] Barbante, P., and Chazot, O., "Flight Extrapolation of Plasma Wind Tunnel Stagnation Region Flowfield," *Journal of Thermophysics and Heat Transfer*, Vol. 20, No. 3, 2006, pp. 493-499. https://doi.org/10.2514/1.17185.

[16] Turchi, A., Matesanz Saiz, J., Magin, T., and Chazot, O., "Duplication of hypersonic stagnation-region aerothermochemistry and gas-surface interaction in high-enthalpy ground testing," *Experiments in Fluids*, Vol. 62, 2021. https://doi.org/10.1007/s00348-021-03320-6.

[17] Bottin, B., Chazot, O., Carbonaro, M., Haegen, V., and Paris, S., "The VKI Plasmatron Characteristics and Performance," 2000, p. 26.

[18] Klomfass, A., and Müller, S., "Calculation of stagnation streamline quantities in hypersonic blunt body flows," *Shock Waves*, Vol. 7, No. 1, 1997, pp. 13-23. https://doi.org/10.1007/s001930050057.

[19] Munafò, A., and Magin, T., "Modeling of stagnation-line nonequilibrium flows by means of quantum based collisional models," *Physics of Fluids*, Vol. 26, No. 9, 2014, p. 097102. https://doi.org/10.1063/1.4894842.

[20] Munafo, A., "Multi-Scale models and computational methods for aerothermodynamics," Ph.D. thesis, Ecole Centrale Paris, 2014.

[21] Roe, P., "Approximate Riemann solvers, parameter vectors, and difference schemes," *Journal of Computational Physics*, Vol. 43, No. 2, 1981, pp. 357-372. https://doi.org/10.1016/0021-9991(81)90128-5.

[22] McBride, B., Zehe, M., and Gordon, S., "NASA Glenn coefficients for calculating thermodynamic properties of individual species," Tech. rep., NASA, 10 2002.

[23] Magin, T., and Degrez, G., "Transport algorithms for partially ionized and unmagnetized plasmas," *Journal of Computational Physics*, Vol. 198, No. 2, 2004, pp. 424-449. https://doi.org/10.1016/j.jcp.2004.01.012.

[24] Magin, T., and Degrez, G., "Transport properties of partially ionized and unmagnetized plasmas," *Physical review. E, Statistical, nonlinear, and soft matter physics*, Vol. 70, 2004, p. 046412. https://doi.org/10.1103/PhysRevE.70.046412.

[25] Olynick, D., Chen, Y., and Tauber, M., "Aerothermodynamics of the Stardust Sample Return Capsule," *Journal of Spacecraft and Rockets*, Vol. 36, 1998, pp. 442-462. https://doi.org/10.2514/2.3466.

[26] Caflisch, R., "Monte Carlo and quasi-Monte Carlo methods," *Acta Numerica*, Vol. 7, 1998, p. 1-49. https://doi.org/10.1017/S0962492900002804.
[27] Metropolis, N., and Ulam, S., "The Monte Carlo Method," *Journal of the American Statistical Association*, Vol. 44, No. 247, 1949, pp. 335-341. https://doi.org/10.1080/01621459.1949.10483310.

[28] Hastings, W., ""Monte Carlo sampling methods using Markov chains and their applications"", *Biometrika*, Vol. 57, No. 1, 1970, pp. 97–109. https://doi.org/10.1093/biomet/57.1.97.

[29] Sacks, J., Welch, W., Mitchell, T., and Wynn, H., "Design and Analysis of Computer Experiments," *Statistical Science*, Vol. 4, No. 4, 1989, pp. 409–423. https://doi.org/10.1214/ss/1177012413.

[30] Rasmussen, C., and Williams, C., *Gaussian processes for machine learning*, The MIT Press, 2006.

[31] Forrester, A., Sobester, A., and Keane, A., *Engineering Design via Surrogate Modelling: A Practical Guide*, Wiley, 2008.

[32] Tryoen, J., Congedo, P. M., Abgrall, R., Villedieu, N., and Magin, T., "Bayesian-based method with metamodels for rebuilding freestream conditions in atmospheric entry flows," *AIAA Journal*, Vol. **52**, No. 10, 2014, pp. 2190 – 2197. https://doi.org/https://doi.org/10.2514/1.J052831.

[33] Bose, D., Wright, M., and Gökçen, T., "Uncertainty and Sensitivity Analysis of Thermochemical Modeling for Titan Atmospheric Entry," *37th AIAA Thermophysics Conference*, 2004. https://doi.org/10.2514/6.2004-2455.

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

[35] Goodman, J., and Weare, J., "Ensemble samplers with affine invariance," *Communications in Applied Mathematics and Computational Science*, Vol. 5, No. 1, 2010, pp. 65–80. https://doi.org/10.2140/camcos.2010.5.65.

[36] Marelli, S., and Sudret, B., "UQLab: A Framework for Uncertainty Quantification in Matlab," *Vulnerability, Uncertainty, and Risk*, American Society of Civil Engineers, 2014, pp. 2554–2563. https://doi.org/10.1061/9780784413609.257.