# Quantification of model-form uncertainties affecting the calibration of a carbon nitridation model by means of Bayesian Model Averaging

Anabel del Val, Thierry Magin, Pietro Marco Congedo

► To cite this version:

Anabel del Val, Thierry Magin, Pietro Marco Congedo. 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, 2022, 213, pp.124271. 10.1016/j.ijheatmasstransfer.2023.124271. hal-04393192

HAL Id: hal-04393192

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Quantification of model-form uncertainties affecting the calibration of a carbon nitridation model by means of Bayesian Model Averaging

Anabel del Val<sup>a,b,*</sup>, Thierry E. Magin<sup>a</sup>, Pietro M. Congedo<sup>b</sup>

<sup>a</sup>von Karman Institute for Fluid Dynamics, Chaussée de Waterloo 72, 1640 Rhode-St-Genèse, Belgium

<sup>b</sup>Inria, Centre de Mathématiques Appliquées, École Polytechnique, IPP, Route de Saclay, 91128 Palaiseau Cedex, France

## Abstract

Severe epistemic uncertainties not only can affect the prescription of parameters within a given model but also the choice of models we make to interpret and infer from experimental data. In this work, we incorporate experimental, parametric and model-form uncertainties in the calibration of a carbon nitridation model. The model-form uncertainties considered stem from the different modeling choices that are taken as valid representations of a set of plasma wind tunnel experiments. To this end, we define a Bayesian model averaging strategy where the marginal posteriors of the nitridation reaction efficiencies are weighted by the marginalized likelihoods of the experimental data for each proposed model. First, Bayes factors are computed to possibly discard invalid models. The baseline model, a thermal equilibrium stagnation line flow with nitridation as only surface reaction, performs as well as all the alternative models proposed, which range from adding surface recombination reactions to considering thermal non-equilibrium in the gas and gas-surface interface. The presence of nitrogen recombination reactions is shown to broaden the support of the nitridation marginal posteriors considerably, allowing it to take on larger values. Lastly, a Bayesian model averaged Arrhenius law for the nitridation efficiencies is computed for a range of surface temperatures.

**Keywords:** Uncertainty Quantification, Bayesian Inference, Thermal Protection Systems, Ablation

## 1. Introduction

Space vehicles plunging into dense planetary atmospheres release a significant amount of kinetic energy in a very short window of time. This energy is transformed into thermal energy through the formation of a strong bow shock ahead of the vehicle. This harsh environment can generate extreme heat loads (> 10,000 J/cm²) [1] which are to be partially absorbed by the vehicle's surface. Aerospace engineers have largely figured out how to efficiently protect the integrity of spacecraft as well as astronauts and cargo from the damage derived from such intense heating. Thermal Protections Systems (TPS) are designed and used for this purpose. In particular, ablative materials such as rigid carbon fibers or silicon composites, impregnated with an organic resin matrix, are broadly used to deal with surface heating for high speed atmospheric entries (~ 11km/s).

The characterization of ablation phenomena plays an important role in the development of theoretical models that can predict material degradation and recession of TPS, both important aspects in the design of ablative TPS [2]. As such, experiments and models are often used to understand the physics of ablation and improve our predictive capabilities. On the experimental side, the stochastic nature of the data must be accurately described and modeled to produce reliable experimental data on which to base our analyses. On the modeling side, many different sources of uncertainties can affect the predictions considerably. Objectively characterizing and quantifying such uncertainties is important to make comparisons to experimentally observed quantities useful, providing a more consistent/quantifiable way of defining new research directions.

Overall, the process of inferring ablation parameters from experimental data poses many questions that span beyond our experimental capabilities to the assumptions in our models. We still grapple with understanding our own

*Corresponding author: von Karman Institute for Fluid Dynamics, Chaussée de Waterloo 72, 1640 Rhode-St-Genèse, Belgium. Email address: ana.isabel.del.val.benitez@vki.ac.be (A. del Val)

Preprint submitted to Elsevier

October 10, 2022
experimental facilities to the point that we must admit ignorance about the relevant physical processes for different experimental conditions in some cases. While experiments are giving us answers to questions, sometimes we must admit that it is not clear what are the questions to which the answers pertain and this is significantly slowing down our progress. In this regard, severe lack of knowledge greatly affects our capabilities to build predictive models.

In particular, the carbon nitridation reaction $C_s + N \rightarrow CN + 0.34eV$ on a solid (s) carbon surface is still hard to predict accurately at temperatures above 1,000 K. The models found in the aerothermodynamics literature are derived empirically by fitting experimental data [3–7]. Even in the cases where uncertainty is taken into account by performing stochastic calibrations, the results only reflect experimental and parametric uncertainty within the given structure of a model [8, 9]. In the cases where different models were entertained, such as in the work of Upadhyay et al. [8], the overall resulting uncertainty in the final calibrated model was never discussed on the basis of model predictive performance and complexity, which prevented such model-form uncertainties from being characterized and included in the end results.

In this work, we address the fact that there exists severe lack of knowledge in the calibration of carbon nitridation reaction efficiencies from plasma wind tunnel data. Experimental and parametric uncertainty, together with a selection of model scenarios, are taken into account to produce a comprehensive picture of nitrogen ablation for selected plasma wind tunnel testing conditions. We use the experimental dataset produced by Helber et al. [10] for which they subjected a graphite material sample to a subsonic plasma flow. They recorded material recession rates and local CN mass densities in the boundary layer by means of spectroscopy measurements. The fact that the material tested is graphite in a pure nitrogen flow gives us an unprecedented window for testing hypotheses concerning the flow and wall states without having to add important uncertainties stemming from other physical processes such as pyrolysis and oxidation. We consider several modeling scenarios to explain the experimental data and predict the observed recession rates and CN densities. The baseline model is a thermal equilibrium stagnation line flow with only nitridation as surface mechanism. Calibrations performed for this model indicate that its assumptions hold for the experimental conditions considered. We use the baseline model as reference or null hypothesis against which three alternative models of additive physical richness are compared.

The first alternative model adds nitrogen recombination reactions at the gas-surface interface by incorporating a recombination efficiency parameter. As recombination reactions compete with nitridation reactions for the available atomic nitrogen, it is important to evaluate their impact in our calibrations. Moreover, we are also interested in learning about nitrogen recombination from these experiments. A second alternative model expresses the same physics but adopts a Surface Energy Balance (SEB) to resolve the surface temperature instead of imposing it as a given boundary condition. The SEB model reduces the complexity in terms of number of parameters to be inferred, given that the surface temperature is now a prediction of the model instead of another free parameter to be imposed. Furthermore, the new parametrization that the SEB introduces comes with more constrained parameter spaces such as the material emissivity. On the other hand, this model has to be tested against an additional observation, the surface temperature. The last alternative model also retains all the physics of the previous models but relaxes the conditions of thermal equilibrium at the gas and gas-surface interface. It remains agnostic about the thermal state of the gas by assuming two distinct thermal baths. At the gas-surface interface, it solves a SEB for each bath, introducing yet a new parametrization that also weighs in the comparison of models.

First, each model proposed is calibrated separately for each experimental condition. The hypothesis testing study is then made through the computation of Bayes factors which compare the resulting nitridation reaction efficiencies marginal posteriors at the order of magnitude level between the baseline and each alternative model. Two different aspects affect the computation of such metric: the goodness-of-fit to the observations that each model can predict as well as their complexities in terms of number of model parameters to be inferred and the extent of their prior spaces. The final nitridation reaction efficiency posterior distributions obtained under each successful model are weighted according to their marginalized likelihoods, and summed using Bayesian Model Averaging (BMA) to account for the epistemic uncertainties stemming from these different valid modeling assumptions. Further, an additional original contribution of this work is to combine all available measurements for the calibration of stochastic Arrhenius laws under each considered model, extending the computation of nitridation efficiencies to a broader span of surface temperatures for which there are no reliable experimental data. The marginal posteriors of the Arrhenius laws parameters are then weighted by their corresponding marginalized likelihoods to generate a model-averaged Arrhenius law that incorporates model-form uncertainties for a broad span of surface temperatures. The overall aim of this work is then twofold: to incorporate model-form uncertainties in the characterization of a carbon nitridation model and to
invalidate or, conversely, find new research directions that would eventually lead to dropping models, refining our understanding of the physics in the process.

The article is organized as follows. Sec. 2 addresses the experimental data, theoretical models and numerical methods used in this work. Sec. 3 presents the Bayesian model averaging methodology in detail. Sec. 4 shows and discusses the obtained results. Sec. 5 presents the conclusions and outlook.

## 2. Overview of experimental data and modeling scenarios

We introduce a brief note on the experimental data used for this study followed by the different modeling choices we can make to retrieve the carbon nitridation model we are after. The experiments and models can be found in more detail in other published works [10–12] so we focus only on recalling the important aspects.

### 2.1. On the experimental data

The experimental scenario is depicted in much detail in [10]. In this section, we briefly describe the measurements and their associated uncertainties.

The experiments were performed to rebuild nitridation reaction efficiencies $\gamma_N^{CN}$ for graphite at high surface temperatures. For this, Helber et al. [10] experimentally observed the response of graphite samples to a subsonic plasma flow produced by the von Karman Institute's (VKI) Plasmatron facility (see Fig. 1). This process creates a reacting boundary layer around the exposed sample where species resulting from the ablation reactions are injected. Overall, Helber et al. were interested in measuring the rate of recession of the graphite sample due to the ablation reaction $C_s + N \rightarrow CN$, as well as the radiation signature from the CN molecule injected in the boundary layer. Additional measurements such as dynamic pressure $P_d$ within the plasma jet, and static pressure $P_s$ in the Plasmatron chamber are also included in the dataset.

<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic view of the ablation experimental set-up for the Plasmatron facility. Figure adapted from Helber [? ]</figcaption></figure>
All the different probes and test sample holders were exposed sequentially to the same plasma flow and the nitridation measurements were recorded until steady-state was maintained for 300 s. The following tables summarize the different measurements and uncertainties used in this work. In Table 1, the measured quantities are the dynamic pressure $P_d$, the wall temperature $T_w$, the chamber static pressure $P_s$, the recession rate $\dot{s}$ and the local CN densities $\rho_{CN}$ estimated from spectra measurements locally resolved. For the particular case G6, CN density measurements are not taken into account given the inconsistencies identified [9].

<table><caption>Table 1: Experimental data.</caption><thead><tr><th>Experiment ID</th><th>P<sub>d</sub> [Pa]</th><th>T<sub>w</sub> [K]</th><th>P<sub>s</sub> [Pa]</th><th>˙s [µm/s]</th><th>ρ<sub>CN</sub> [µg/m<sup>3</sup>]</th></tr></thead><tbody><tr><td>G4</td><td>231</td><td>2225</td><td>1500</td><td>1.41</td><td>385</td></tr><tr><td>G5</td><td>268</td><td>2410</td><td>1500</td><td>1.64</td><td>800</td></tr><tr><td>G6</td><td>312</td><td>2535</td><td>1500</td><td>2.60</td><td>-</td></tr><tr><td>G7</td><td>330</td><td>2575</td><td>1500</td><td>2.51</td><td>1050</td></tr></tbody></table>

The experimental uncertainties are assumed independently distributed Gaussian distributions $N(\mu, \sigma)$ with $\sigma$ taken considering the bounds given on the measurements as $2\sigma$ level of confidence [13–15].

<table><caption>Table 2: Experimental uncertainties.</caption><thead><tr><th>Experiment ID</th><th>P<sub>d</sub> [Pa]</th><th>T<sub>w</sub> [K]</th><th>P<sub>s</sub> [Pa]</th><th>˙s [µm/s]</th><th>ρ<sub>CN</sub> [µg/m<sup>3</sup>]</th></tr></thead><tbody><tr><td>G4</td><td>N(231, 2.31)</td><td>N(2225, 11.12)</td><td>N(1500, 22.5)</td><td>N(1.41, 0.26)</td><td>N(385, 52)</td></tr><tr><td>G5</td><td>N(268, 2.68)</td><td>N(2410, 12.05)</td><td>N(1500, 22.5)</td><td>N(1.64, 0.27)</td><td>N(800, 100)</td></tr><tr><td>G6</td><td>N(312, 3.12)</td><td>N(2535, 12.67)</td><td>N(1500, 22.5)</td><td>N(2.60, 0.27)</td><td>-</td></tr><tr><td>G7</td><td>N(330, 3.3)</td><td>N(2575, 12.87)</td><td>N(1500, 22.5)</td><td>N(2.51, 0.25)</td><td>N(1050, 162)</td></tr></tbody></table>

## 2.2. Baseline model

We use a one dimensional form of the full Navier-Stokes equations in order to perform efficient simulations of a high temperature, reacting flow with complex physico-chemical models and boundary conditions. Originally proposed in [16], this formulation is called “Dimensionally Reduced Navier-Stokes Equations (DRNSE)” and, in the present study, it is implemented in an in-house software developed by Munafò [11]. This code is able to reproduce the flow along the stagnation-line of spherical and cylindrical bodies [16].

In the solver, the spatial discretization is performed by means of a cell-centered finite volume method. The numerical advective fluxes are computed by means of the AUSM$^{++}$-up2 scheme [17], which splits the flux into a convective and a pressure terms. Second order accuracy in space is achieved by using an upwind reconstruction to obtain the cell interface variables. The diffusive fluxes and source terms are both evaluated in terms of primitive variables, namely, species densities $\rho_i$, velocity components $u, v$, and temperature $T$. The value of the primitive variables and their gradients at the volume interface are computed by a weighted average and a central finite difference approximation, respectively. The cell center value of the primitive variables to be used in the evaluation of the diffusive source terms are computed using a two-point central finite-difference approximation. A fully implicit scheme is used for the present study. The code is coupled with the Mutation$^{++}$ library [12]. Species thermodynamic properties are obtained from Rigid-Rotor Harmonic-Oscillator computations, and relative mixture quantities are derived from pure species quantities. The transport properties are derived from kinetic theory, which provides relationships for macroscopic transport coefficients based on microscopic collision integrals. The chemical production rates for species, based on elementary chemical reactions, including third body, are calculated by taking the forward reaction rate coefficients specified by the user in an Arrhenius law form. The backward rate coefficient is determined by satisfying the equilibrium relation. The library has been designed, implemented and extensively tested to ensure high-fidelity together with low computational costs.
A coupled ablative boundary condition is also included in a generalized formulation [18] in Mutation++ based on surface mass and energy balances. Several other authors have applied this approach to study the ablative and non-ablative gas-surface interaction of thermal protection materials with the surrounding flowfield [19, 20]. The advective and the diffusive fluxes at the gas-surface interface are obtained by solving the surface species mass balance. This balance is obtained limiting the control volume of the mass conservation equation to the thin lamina representing the gas-surface interface [10]

$$[\rho_{w,i}(\mathbf{u}_g - \mathbf{u}_r) + \rho_{w,i}V_{w,i}] \cdot \mathbf{n} = \dot{\omega}_{w,i}, \quad \forall i \in S, \quad (1)$$

where the source term $\dot{\omega}_{w,i}$ represents the sum for species i of the production rates of individual heterogeneous reactions with the surface $\dot{\omega}_{w,i} = \sum_{r=1}^{n_{Rw}} \dot{\omega}_{w,i}^r$, with $n_{Rw}$ being the number of considered reactions. For this model, we do not need to include a surface energy balance given that the temperature at the wall $T_w$ is measured and imposed.

We need to introduce closure terms for the production of species $\dot{\omega}_{w,i}^r$. The approach for the production term considered in this work is a phenomenological one because there is only so much that we can learn from the present data. This type of models are quite popular in the aerothermodynamics literature for their straightforward implementation [21]. We are only concerned with ablation models for surface reactions involving nitrogen on the side of the flow and carbon on the material side. The model for the production rate is defined

$$\dot{\omega}_{w,i}^r = \gamma_i^r m_i \mathcal{F}_i^\downarrow, \quad \forall i \in S, \forall r \in n_{Rw}, \quad (2)$$

with $\gamma_i^r$ representing the reaction efficiency, $m_i$ the particle i mass, and $\mathcal{F}_i^\downarrow$ the impinging flux of species i.

The inlet boundary condition used for the simulations of such reacting flows comprises species densities $\rho_i$, temperature $T$, and velocity components $u, v$. The information needed to prescribe the boundary condition at the free stream pertains to a supersonic inlet even though we want to carry out subsonic simulations. This methodology has been investigated and discussed in [22], where it was verified that the solution of the Riemann problem at the boundary interface automatically screens the information needed from the inlet.

All in all, we intend to find the flow solutions that are compatible with our experimental observations under the considered model. This means that from the observations reported by Helber et al. [10], we want to infer the inlet conditions and the wall conditions in such a way that we are left with a population of possible flow solutions. After taking into account some physical considerations for the definition of the inlet quantities [9], for given Plasmatron operating conditions, we are left with the following expressions for the simulated counterparts of the observations of recession rate and local CN densities

$$\dot{s} = \dot{s}(P_{\delta}, T_{\delta}, P_{d}, T_{w}, \gamma_{N}^{\text{CN}}), \quad (3)$$

$$\rho_{\text{CN}} = \rho_{\text{CN}}(P_{\delta}, T_{\delta}, P_{d}, T_{w}, \gamma_{N}^{\text{CN}}), \quad (4)$$

where the subscript $\delta$ makes reference to the quantities defined at the boundary layer edge. These expressions constitute the dependencies of the forward model.

## 2.3. Alternative modeling scenarios

Several model assumptions stand at the basis of the baseline model. The flow is assumed in thermal equillibrium and the ablative wall conditions amount to only nitridation reactions. Furthermore, measuring the surface temperature also allows us to test ablation energy balance models by using the wall temperature as the quantity to be reproduced. In the available literature regarding experimental studies for nitridation, flow simulations are a needed component to compute atomic nitrogen concentrations at probe locations. The resulting nitridation reaction efficiencies rely heavily on such computations which are full of assumptions about the flow and material behaviors. Suzuki et al. [5, 23] used flow simulations that did not account for nitridation at the surface. Helber et al. [10] showed that the atomic nitrogen concentration profile along the stagnation line of the boundary layer is significantly different (~10%) when nitridation is turned on/off. Suzuki et al. [23] underestimated nitridation reaction efficiencies this way.

The thermal estate of the gas is not expected to impact the previous estimations of nitridation reaction efficiencies given that the recession rates and surface temperatures are accurately measured. Nevertheless, this set of experimental
data presents a great opportunity to study the evidence for thermal non-equilibrium and the adequacy of surface energy balance models. Together with the thermodynamics and energy wall conditions, nitrogen recombination reactions are added to our models. As recombination reactions compete with nitridation reactions for the available atomic nitrogen, it is important to evaluate their impact in our calibrations. Moreover, we are also interested in learning about nitrogen recombination from these experiments.

Overall, three different alternative model scenarios of additive physical richness are presented. They all build on the previous baseline model, $M_0$:

* $M_1$: This model accounts for nitrogen recombination at the surface by introducing an extra parameter $\gamma_N^{N_2}$ in the surface mass balance through the nitrogen recombination reaction $N(s) + N \rightarrow N_2$. Included in the original model of Zhluktov and Abe [24], nitrogen recombination is a competing mechanism to nitridation, consuming available atomic nitrogen which, in turn, has the potential to impact the former reaction in a coupled fashion. This model also considers the flow to be in thermal equilibrium and imposes the surface temperature. The implicit model relationships between inputs and outputs for this model reads

$$ \dot{s} = \dot{s}(P_{\delta}, T_{\delta}, P_{d}, T_{w}, \gamma_{N}^{CN}, \gamma_{N}^{N_2}), \quad (5) $$

$$ \rho_{CN} = \rho_{CN}(P_{\delta}, T_{\delta}, P_{d}, T_{w}, \gamma_{N}^{CN}, \gamma_{N}^{N_2}). \quad (6) $$

* $M_2$: In this model the thermal state is considered in equilibrium, surface recombination is taken into account, and the surface temperature is now a prediction of the model through a Surface Energy Balance (SEB)

$$ \left[ \rho_w(\mathbf{u}_g - \mathbf{u}_r)h_w + \sum_{i \in S} \rho_{w,i}h_{w,i}V_{w,i} + \lambda_w \nabla T_w + \lambda_w^V \nabla T_w^V \right] \cdot \mathbf{n} = \dot{\Omega}_{w,i}, \quad \forall i \in S. \quad (7) $$

where the terms $\sum_{i \in S} \rho_{w,i}h_{w,i}V_{w,i}$, $\lambda_w \nabla T_w$ and $\lambda_w^V \nabla T_w^V$ are the diffusive and conductive heat fluxes for translational and vibrational thermal baths, respectively. Quantities $\rho_w, \mathbf{u}_g, \mathbf{u}_r$ and $h_w$ are the gas density at the surface, the velocity of the gas leaving the interface, the velocity at which the interface is moving, and the enthalpy at the surface, respectively. The term $V_{w,i}$ refers to the diffusion velocity of species $i$ at the gas-surface lamina, and $\lambda_w, \lambda_w^V$ are the thermal conductivities for translational and internal groups of energies [9]. The direction normal to the surface pointing towards the gas is designated by $\mathbf{n}$ while the surface temperature for the vibrational thermal bath is $T_w^V$. As it is recalled, thermal equilibrium is considered in this model both in the gas phase and gas-surface interface for which the condition $T_w = T_w^V$ is further imposed. The only source term taken into account in Eq. (7) as the effect of the surface in the surrounding gas is the radiative flux. This term follows the Stefan-Boltzmann law for a body in thermodynamic equilibrium at temperature $T_w$

$$ \dot{\Omega}_w = \sigma \epsilon (T_w)^4, \quad (8) $$

where $\sigma$ is the Stefan-Boltzmann constant and $\epsilon$ the emissivity coefficient of the surface.

In this work, the ablation modeling is restricted to carbon preforms as it constitutes the material of choice for wind tunnel testing for nitridation. For carbon preforms, the binding resine is non-existent, allowing us to drop the modeling of the physics associated to pyrolysis in the energy balances [20, 25]. Consequently, Eq. 7 is the result of simplifications that neglect the presence of pyrolysis gases and approximate the conductive heat flux following steady-state ablation conditions in the solid phase and integrating over the semi-infinite material [26, 27].

The implicit model relationships between inputs and outputs for this model reads

$$ \dot{s} = \dot{s}(P_{\delta}, T_{\delta}, P_{d}, \gamma_{N}^{CN}, \gamma_{N}^{N_2}), \quad (9) $$

$$ \rho_{CN} = \rho_{CN}(P_{\delta}, T_{\delta}, P_{d}, \gamma_{N}^{CN}, \gamma_{N}^{N_2}), \quad (10) $$
$$T_w = T_w(P_\delta, T_\delta, P_d, \gamma_N^{\text{CN}}, \gamma_N^{\text{N}_2}). \qquad (11)$$

* $M_3$: This model assumes thermal non-equilibrium by considering two different thermal baths for the translational-rotational and vibrational-electronic-electron temperatures. Furthermore, it is also agnostic about the wall surface energy condition by allowing it to be in non-equilibrium by additionally solving the surface energy balance

$$\left[ \rho_w(\mathbf{u}_g - \mathbf{u}_r)h_w^{\mathcal{V}} + \sum_{i \in S} \rho_{w,i}h_{w,i}^{\mathcal{V}}V_{w,i} + \lambda_w^{\mathcal{V}}\nabla T_w^{\mathcal{V}} \right] \cdot \mathbf{n} = \dot{\Omega}_w^{\mathcal{V}T} + \dot{\Omega}_w^{\mathcal{C}V}, \qquad (12)$$

where the convective and conductive terms in Eq. (7) are now only associated to the internal energies $\mathcal{V}$. The first source term $\dot{\Omega}_w^{\mathcal{V}T}$ is the relaxation of internally excited particles by collision with the surface (heterogeneous quenching). This term was derived by adapting the volumetric Landau-Teller gas-gas collisions formulation [28] to surface gas-surface collisions. The expression reads

$$\dot{\Omega}_w^{\mathcal{V}T} = \sum_{i \in S} \alpha_i \mathcal{F}_i^{\downarrow} m_i (h_{w,i}^{\mathcal{V}} (T^{\mathcal{V}}) - h_{w,i}^{\mathcal{V}} (T)), \qquad (13)$$

where $\alpha_i$ is a parameter that considers that only the $\alpha_i$ fraction of the collisional events of species $i$ leads to a relaxation. The term $\mathcal{F}_i^{\downarrow}$ represents the impinging flux of species $i$ and $m_i$ is the particle mass of species $i$. The second source term in Eq. (12) is a chemical source term which takes into account that not all the energy from chemical reactions is accommodated on the surface. It reads

$$\dot{\Omega}_w^{\mathcal{C}V} = \sum_{i \in S} (1 - \beta_i) (h_{w,i}^{\mathcal{V}} - h_{w,i}), \qquad (14)$$

where the $\beta_i$ coefficient is an adaptation for a two-temperature model of the energy accommodation term defined by Halpern and Rosner [29], and that a $1-\beta$ quota of the heat flux is transferred into the internal energy modes of the products. More details on the two-temperature surface energy balance here presented can be found in the work of Capriati et al. [30].

Finally, surface recombination reactions are also included. The implicit model relationships between inputs and outputs for this model reads

$$\dot{s} = \dot{s}(P_\delta, T_\delta, T_\delta^{\mathcal{V}}, P_d, \gamma_N^{\text{CN}}, \gamma_N^{\text{N}_2}, \alpha, \beta, \epsilon), \qquad (15)$$

$$\rho_{\text{CN}} = \rho_{\text{CN}}(P_{\delta}, T_{\delta}, T_{\delta}^{\mathcal{V}}, P_{d}, \gamma_{\text{N}}^{\text{CN}}, \gamma_{\text{N}}^{\text{N}_2}, \alpha, \beta, \epsilon), \qquad (16)$$

$$T_w = T_w(P_\delta, T_\delta, T_\delta^{\mathcal{V}}, P_d, \gamma_N^{\text{CN}}, \gamma_N^{\text{N}_2}, \alpha, \beta, \epsilon). \qquad (17)$$

The main characteristics of each modeling scenario are depicted in the following Table 3

<table><caption>Table 3: Modeling scenarios considered in this work.</caption><thead><tr><th>Modeling assumptions</th><th>M<sub>0</sub></th><th>M<sub>1</sub></th><th>M<sub>2</sub></th><th>M<sub>3</sub></th></tr></thead><tbody><tr><td>Thermal non-equilibrium</td><td>No</td><td>No</td><td>No</td><td>Yes</td></tr><tr><td>Surface energy balance</td><td>No</td><td>No</td><td>Yes</td><td>Yes</td></tr><tr><td>Surface recombination γ<sub>N</sub><sup>N<sub>2</sub></sup></td><td>No</td><td>Yes</td><td>Yes</td><td>Yes</td></tr></tbody></table>
### 3. Bayesian model averaging methodology

In this section, we present the methods used to identify whether or not the proposed alternative models can explain the experimental data comparatively well, and how we quantify the resulting uncertainty stemming from such different modeling choices.

#### 3.1. Bayesian inference

Performing Bayesian analyses requires the specification and computation of the different ingredients of the Bayes' rule. We can formulate such rule in a general way for each model considered as

$$ \mathcal{P}(\boldsymbol{q}^{(i)}|\boldsymbol{y}_{\text{obs}}^{(i)}, \mathcal{M}_i) = \frac{\mathcal{P}(\boldsymbol{y}_{\text{obs}}^{(i)}|\boldsymbol{q}^{(i)}, \mathcal{M}_i) \mathcal{P}(\boldsymbol{q}^{(i)}|\mathcal{M}_i)}{\int_{\Omega} \mathcal{P}(\boldsymbol{y}_{\text{obs}}^{(i)}|\boldsymbol{q}^{(i)}, \mathcal{M}_i) \mathcal{P}(\boldsymbol{q}^{(i)}|\mathcal{M}_i)d\boldsymbol{q}^{(i)}} \quad \forall i \in \{0, \dots, 3\}, \quad (18) $$

where $\boldsymbol{q}^{(i)}$ is the vector of parameters to be inferred which results from the stagnation line modeling and physical considerations presented in Sec. 2 for each model, and $\boldsymbol{y}_{\text{obs}}^{(i)}$ is the vector of observations for model $i$. Furthermore, for compact notation, we split the observations and parameters to be inferred in two distinct subsets each. The subset $\boldsymbol{q}^{(i)'}$ whose parameters have direct measured counterparts, and $\boldsymbol{y}_{\text{obs}}^{(i)'}$ containing the corresponding measurements. The subset $\boldsymbol{q}^{(i)''}$ contains the parameters that appear implicitly in the model outputs dependencies, while $\boldsymbol{y}_{\text{obs}}^{(i)''}$ is the vector containing the measurements that are to be compared against the corresponding model outputs. We also define the set of model outputs $\boldsymbol{Y}^{(i)}$ which contains the model outputs to be compared against the corresponding $\boldsymbol{y}_{\text{obs}}^{(i)''}$ set, such as recession rates, CN densities and wall temperatures. This partition of the parameters, observations and model outputs sets helps setting up the likelihood function with a more compact notation. Table 4 shows an overview of the different parameters, observations and model outputs considered under each model.

<table><caption>Table 4: Vectors of parameters to be calibrated, observations and model outputs for the different modeling scenarios.</caption><thead><tr><th>Vectors</th><th>M<sub>0</sub></th><th>M<sub>1</sub></th><th>M<sub>2</sub></th><th>M<sub>3</sub></th></tr></thead><tbody><tr><td>Parameters q<sup>(i)</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>q<sup>(i)'</sup></td><td>P<sub>δ</sub>, P<sub>d</sub>, T<sub>w</sub></td><td>P<sub>δ</sub>, P<sub>d</sub>, T<sub>w</sub></td><td>P<sub>δ</sub>, P<sub>d</sub></td><td>P<sub>δ</sub>, P<sub>d</sub></td></tr><tr><td>q<sup>(i)''</sup></td><td>T<sub>δ</sub>, γ<sub>N</sub><sup>CN</sup></td><td>T<sub>δ</sub>, γ<sub>N</sub><sup>CN</sup>, γ<sub>N</sub><sup>N<sub>2</sub></sup></td><td>T<sub>δ</sub>, γ<sub>N</sub><sup>CN</sup>, γ<sub>N</sub><sup>N<sub>2</sub></sup></td><td>T<sub>δ</sub>, T<sub>δ</sub><sup>V</sup>, γ<sub>N</sub><sup>CN</sup>, γ<sub>N</sub><sup>N<sub>2</sub></sup>, α, β, ε</td></tr><tr><td>Observations y<sub>obs</sub><sup>(i)</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>y<sub>obs</sub><sup>(i)'</sup></td><td>P<sub>d</sub><sup>meas</sup>, P<sub>s</sub><sup>meas</sup>, T<sub>w</sub><sup>meas</sup></td><td>P<sub>d</sub><sup>meas</sup>, P<sub>s</sub><sup>meas</sup>, T<sub>w</sub><sup>meas</sup></td><td>P<sub>d</sub><sup>meas</sup>, P<sub>s</sub><sup>meas</sup></td><td>P<sub>d</sub><sup>meas</sup>, P<sub>s</sub><sup>meas</sup></td></tr><tr><td>y<sub>obs</sub><sup>(i)''</sup></td><td>σ<sup>meas</sup>, ρ<sub>CN</sub><sup>meas</sup></td><td>σ<sup>meas</sup>, ρ<sub>CN</sub><sup>meas</sup></td><td>σ<sup>meas</sup>, ρ<sub>CN</sub><sup>meas</sup>, T<sub>w</sub><sup>meas</sup></td><td>σ<sup>meas</sup>, ρ<sub>CN</sub><sup>meas</sup>, T<sub>w</sub><sup>meas</sup></td></tr><tr><td>Model outputs Y<sup>(i)</sup></td><td>σ, ρ<sub>CN</sub></td><td>σ, ρ<sub>CN</sub></td><td>σ, ρ<sub>CN</sub>, T<sub>w</sub></td><td>σ, ρ<sub>CN</sub>, T<sub>w</sub></td></tr></tbody></table>

The definition of the likelihood function encapsulates the forward model given that the parameters to be inferred have an implicit relation with the simulated observations. In Sec. 2, we referenced the models used to simulate the reacting boundary layer encountered in the experiments. Specifically, we can predict recession rates, CN densities and wall temperature (only for model $M_2$ and $M_3$) which we can then compare to the observations. Having our different sets of measurements, parameters and outputs defined, for each experimental condition and model we propose a general likelihood of the form
$$ \mathcal{P}(\mathbf{y}_{\text{obs}}^{(i)} | \mathbf{q}^{(i)}) \propto \prod_{j=1}^{N^{(i)\prime}} \exp\left[-\frac{|\mathbf{y}_{\text{obs},j}^{(i)\prime} - \mathbf{q}_j^{(i)\prime}|^2}{2\sigma_{\mathbf{y}_{\text{obs},j}^{(i)\prime}}^2}\right] \times \prod_{k=1}^{N^{(i)\prime\prime}} \exp\left[-\frac{|\mathbf{y}_{\text{obs},k}^{(i)\prime\prime} - \mathbf{Y}_k^{(i)}|^2}{2\sigma_{\mathbf{y}_{\text{obs},k}^{(i)\prime\prime}}^2}\right], \quad (19) $$

where $\mathbf{q}_j^{(i)\prime}$, $\mathbf{y}_{\text{obs},j}^{(i)\prime}$ represent the jth elements of the $\mathbf{q}^{(i)\prime}$, $\mathbf{y}_{\text{obs}}^{(i)\prime}$ sets, while $\mathbf{y}_{\text{obs},k}^{(i)\prime\prime}$ is the kth element of the $\mathbf{y}_{\text{obs}}^{(i)\prime\prime}$ set. The terms $N^{(i)\prime}$ and $N^{(i)\prime\prime}$ refer to the number of elements in the prime $\mathbf{q}^{(i)\prime}$, $\mathbf{y}_{\text{obs}}^{(i)\prime}$ sets, and the observations double prime $\mathbf{Y}_k^{(i)}$ set, respectively. The term $\mathbf{Y}_k^{(i)}$ is the kth element of the $\mathbf{Y}^{(i)}$ set. The standard deviations of the elements in both sets of measurements are expressed as $\sigma_{\mathbf{y}_{\text{obs},j}^{(i)\prime\prime}}$ and $\sigma_{\mathbf{y}_{\text{obs},k}^{(i)\prime\prime}}$. In practice, given that the forward model evaluations that compose the sets $\mathbf{Y}^{(i)}$ are expensive, we approximate them with Gaussian process (GP) surrogates whose L₂-norm thresholds are kept under 1% [9].

Besides the likelihood function, we need to define prior distributions for all the parameters in each $\mathbf{q}^{(i)}$. In the available literature for tests in the VKI Plasmatron and stagnation line plasma flows [13–15, 31], we can only establish some bounds a priori. In a maximum entropy set-up [32] with this a priori information, the resulting prior distributions are uniform over the prescribed range. Table 5 shows the chosen prior distributions.

<table><caption>Table 5: Parameters to be calibrated and their priors for the different modeling scenarios.</caption><thead><tr><th>Inference parameter</th><th>M<sub>0</sub></th><th>M<sub>1</sub></th><th>M<sub>2</sub></th><th>M<sub>3</sub></th><th>Prior</th></tr></thead><tbody><tr><td>Dynamic pressure P<sub>d</sub></td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td><td>U[200, 360]</td></tr><tr><td>Edge translational temperature T<sub>δ</sub></td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td><td>U[9000, 13000]</td></tr><tr><td>Static pressure P<sub>δ</sub></td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td><td>U[1200, 1700]</td></tr><tr><td>Surface temperature T<sub>w</sub></td><td>Yes</td><td>Yes</td><td>No</td><td>No</td><td>U[2000, 4000]</td></tr><tr><td>Nitridation efficiency γ<sub>N</sub><sup>CN</sup></td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td><td>log U[-4, 0]</td></tr><tr><td>Recombination efficiency γ<sub>N</sub><sup>N2</sup></td><td>No</td><td>Yes</td><td>Yes</td><td>Yes</td><td>log U[-4, 0]</td></tr><tr><td>Emissivity ε</td><td>No</td><td>No</td><td>Yes</td><td>Yes</td><td>U[0.5, 1]</td></tr><tr><td>Energy accommodation coefficient β</td><td>No</td><td>No</td><td>No</td><td>Yes</td><td>U[-1, 1]</td></tr><tr><td>Collisional efficiency α</td><td>No</td><td>No</td><td>No</td><td>Yes</td><td>U[0, 1]</td></tr></tbody></table>

Depending on the algorithms used to solve Eq. (18), we generally do not need to compute the normalization constant in the denominator [33]. Its computation is not necessary to achieve the objective of sampling from the posterior distribution for a given model. When different models are to be compared in their predictive capabilities, the normalization constant becomes central and this is at the core of hypothesis testing studies using Bayesian statistics which we address in the following section.

## 3.2. Hypothesis testing

Performing a hypothesis testing study entails the evaluation of the posterior probabilities of each competing model (or hypothesis in our context). This comparison allows us to weigh the relative merit of each model. At first thought, we might consider the most suitable model to be the one that best fits the data. Difficulties quickly arise as we aim at comparing models of different complexity. More complex models involving the definition of more parameters can generally be more expressive and fit very well the data. In this case, the results of hypothesis testing studies have the philosophical principle of *Ockham's razor* built in them. This principle embodies the fact that the simplest models that can explain comparatively well the data are the preferred ones. The comparison is driven through the posterior ratio

$$ \frac{\mathcal{P}(M_i | \mathbf{y}_{\text{obs}})}{\mathcal{P}(M_0 | \mathbf{y}_{\text{obs}})}, \quad \forall i \in \{1, \dots, 3\}, \qquad (20) $$

for which $M_0$ is used as reference in the denominator. This notation omits using superscripts (i) to denote each model's sets of parameters and observations for clarity. The comparison is performed at the order of magnitude level. If the posterior of the baseline model is several orders of magnitude below the competing model, we would prefer the
competing model and viceversa. If both posteriors are of similar orders of magnitude, the data do not provide enough evidence to preferentially support either model.

The posterior ratio can be decomposed according to the Bayes' rule

$$ \frac{\mathcal{P}(M_i|\mathbf{y}_{\text{obs}})}{\mathcal{P}(M_0|\mathbf{y}_{\text{obs}})} = \frac{\mathcal{P}(\mathbf{y}_{\text{obs}}|M_i)}{\mathcal{P}(\mathbf{y}_{\text{obs}}|M_0)} \times \frac{\mathcal{P}(M_i)}{\mathcal{P}(M_0)}, \quad \forall i \in \{1, \dots, 3\}, \qquad (21) $$

where the posterior ratio is expressed as the product of the ratio of the marginalized likelihoods and priors of the different models. The marginalized likelihood of the data under the assumptions of each model is expressed as follows

$$ \mathcal{P}(\mathbf{y}_{\text{obs}}|M_i) = \int_{\Omega} \mathcal{P}(\mathbf{y}_{\text{obs}}|\mathbf{q}, M_i)\mathcal{P}(\mathbf{q})d\mathbf{q}, \quad \forall i \in \{0, \dots, 3\}, \qquad (22) $$

where the term $\mathcal{P}(\mathbf{y}_{\text{obs}}|\mathbf{q}, M_i)$ is the likelihood in Eq. (18) for the different models under consideration. $\mathcal{P}(\mathbf{q})$ is the prior of the parameters $\mathbf{q}$. The marginalized likelihood is also the normalization constant of the posterior $\mathcal{P}(\mathbf{q}|\mathbf{y}_{\text{obs}}, M_i)$ under each model. The prior ratio $\mathcal{P}(M_i)/\mathcal{P}(M_0)$ is taken as unity in this case. We consider all models to be equally likely to explain the experiments a priori. The results from the hypothesis testing study is entirely based on how well each model fits the data (values of the likelihood $\mathcal{P}(\mathbf{y}_{\text{obs}}|\mathbf{q}, M_i)$), and how complex each model is regarding the number and ignorance about each model parameter (dimensionality of the integral in Eq. (22) and overall support of the priors of $\mathbf{q}$).

Looking at Table 4, it is straightforward to see that the null hypothesis $M_0$ is the simplest model in terms of model parameters, while $M_3$ is the most complex one. For $M_0$ to be definitely invalidated by the experimental evidence, its predictions must be considerably worse than the ones produced by any other model. We have to consider this assessment with care. We already know that the experimental data can be reasonably well explained with $M_0$ [9], and the fact that a competing model does much better is reason enough to abandon the hypotheses behind $M_0$, or at least modify them. If $M_0$ was not good enough to begin with, we would not have any reason to consider it in the first place for our hypothesis testing study.

Typically, the evidences are compared through the Bayes factor

$$ BF_{i0} = \frac{\mathcal{P}(\mathbf{y}_{\text{obs}}|M_i)}{\mathcal{P}(\mathbf{y}_{\text{obs}}|M_0)}, \quad \forall i \in \{1, \dots, 3\}, \qquad (23) $$

where $BF$ refers to the Bayes factor which is generally expressed in logarithmic scale as

$$ \log_{10} BF_{i0} = \log_{10}(\mathcal{P}(\mathbf{y}_{\text{obs}}|M_i)) - \log_{10}(\mathcal{P}(\mathbf{y}_{\text{obs}}|M_0)), \quad \forall i \in \{1, \dots, 3\}. \qquad (24) $$

For strong evidences in favour or against a certain hypothesis, indicative thresholds are found in the literature by Jeffreys [34] and Kass and Raftery [35]. In general, the Bayes factor should be above a threshold of $|\log_{10} BF| = 8$, which indicates that as long as the evidences for each model remain in within some orders of magnitude, we cannot conclude anything from the data. In this regard, it is straightforward to see that $BF_{10} = 7$ indicates that the data are 7 times more likely under $M_1$ than under $M_0$, whereas $BF_{10} = 0.2$ indicates that the observations are 5 times more likely under $M_0$ than under $M_1$.

### 3.3. Model averaging

Given the same prior odds for both models, the Bayes factor is a metric that allows us to quantify the overall best performing model. Apart from Bayes factors, other model selection criteria are commonly used in the literature. One criterion that closely approximates the Bayes factor is the Bayesian Information Criterion (BIC). Developed by Schwarz [36], the BIC serves as an asymptotic approximation to a transformation of the Bayesian posterior probability of a candidate model. The model corresponding to the minimum value of BIC is selected. Similarly, the Akaike Information Criterion (AIC) score rewards models that achieve a high goodness-of-fit score and penalizes them if they become overly complex [37]. The model with the lower AIC score is expected to strike a superior balance. A comprehensive overview of AIC and other popular model selection methods is given by Ding et al. [38].

In this work, we focus on computing Bayes factors as the metric used to discard different thermochemical models and ablation boundary conditions (Sec. 2). An issue that arises is the fact that Bayes factor thresholds are not very
sharp and there can be a grey area where the prior spaces for both models are comparable as well as their goodness-of-fit ratio. There can also be a trade-off between prior space, or complexity, and goodness-of-fit in such a way that the Bayes factor falls in the grey area. This would indicate that the more complex model manages to overcome the complexity penalization by fitting the data much better than the simpler model. In this case, discarding a given model is not possible on the basis of the data and prior information we possess.

Already Chamberlin [39] advocated the concept of “multiple working hypotheses”, implying that there is not a single “right” model; instead, there are several well-supported models that are being entertained. If new experiments are well-designed to maximize their informational content, the results tend to support one or more models while providing less support for others. Repetition of this general approach leads to advances, combining experimental design with comparison of theoretical models. New or more elaborate models are consequently postulated, while models with little empirical support are gradually dropped from consideration.

Acknowledging this fact means incorporating model uncertainty stemming from the existence of more than one plausible model. Given these conditions, we may forego selection, retaining the uncertainty across the model space. In this work, we focus on Bayesian model averaging [40]. The model averaged posterior density for the model parameters $\gamma_N^{CN}$ shared between the baseline and each alternative model reads

$$ \mathcal{P}(\gamma_N^{CN}|\mathbf{y}_{\text{obs}}) = \sum_{i \in \{0,...,3\}} \mathcal{P}(\gamma_N^{CN}|\mathbf{y}_{\text{obs}}, \mathcal{M}_i)\mathcal{P}(\mathcal{M}_i|\mathbf{y}_{\text{obs}}). \quad (25) $$

The result is an overall marginal posterior which gathers the information on nitridation efficiencies provided by each valid model, incorporating this way model-form uncertainties in the calibration of $\gamma_N^{CN}$.

BMA is important for expressing model uncertainty in a consistent way, thus, avoiding over-confident inferences. Introduced by Roberts [41], who for the first time combined two expert opinions in the statistical literature, BMA has a challenging implementation due to the complex probability distributions it deals with and their integrations. Moreover, assigning prior distributions to each specified model can be a source of caveats when the choice is not obvious.

# 4. Results and discussion

In this section, we present and discuss the results of the analyses. First, the hypothesis testing results are laid out followed by the resulting model averaging outcomes. Summary statistics are provided for the final averaged models.

## 4.1. Hypothesis testing

The evidences $\mathcal{P}(\mathbf{y}_{\text{obs}}|\mathcal{M}_i)$ are computed for each model and experimental case (see Table 6). As we recall, given the same prior odds to each competing model, the evidences dominate the computation of the Bayes factors for model selection. In the particular case of experiment G6, only recession rates and wall temperatures (when applicable) are used in the computations, producing evidences of several orders of magnitude lower than for the other experimental conditions. This is not an issue given that the comparisons are made among the models for each experimental case independently and on the basis of the same experimental data. The CN density measurement under G6 conditions was deemed untrustworthy in previous analyses [9], therefore, it is not used in the context of this study.

<table><caption>Table 6: Evidences or marginalized likelihoods $\mathcal{P}(\mathbf{y}_{\text{obs}}|\mathcal{M}_i)$ of the different models for the given experimental conditions.</caption><thead><tr><th>Experiment ID</th><th>M<sub>0</sub></th><th>M<sub>1</sub></th><th>M<sub>2</sub></th><th>M<sub>3</sub></th></tr></thead><tbody><tr><td>G4</td><td>1.12 &times; 10<sup>11</sup></td><td>5.6 &times; 10<sup>10</sup></td><td>6.0 &times; 10<sup>7</sup></td><td>6.64 &times; 10<sup>7</sup></td></tr><tr><td>G5</td><td>7.31 &times; 10<sup>10</sup></td><td>6.76 &times; 10<sup>10</sup></td><td>2.85 &times; 10<sup>5</sup></td><td>7.73 &times; 10<sup>6</sup></td></tr><tr><td>G6</td><td>5.1 &times; 10<sup>4</sup></td><td>5.08 &times; 10<sup>4</sup></td><td>3.45 &times; 10<sup>1</sup></td><td>4.43 &times; 10<sup>1</sup></td></tr><tr><td>G7</td><td>3.7 &times; 10<sup>10</sup></td><td>3.76 &times; 10<sup>10</sup></td><td>5.0 &times; 10<sup>5</sup></td><td>3.7 &times; 10<sup>6</sup></td></tr></tbody></table>

Table 7 shows the resulting Bayes factors. There is no definitive evidence against the baseline model which considers thermal equilibrium and nitridation reactions as the only wall reactions. On the other hand, it means that the SEB model ($\mathcal{M}_2$) does not perform too poorly compared to imposing the measured surface temperature directly,
giving us some confidence in such models for ablative surfaces under the considered Plasmatron conditions. The baseline model $M_0$ does pretty well compared against the other models for which the larger number of parameters to be inferred penalizes the resulting evidence. In the case of model $M_2$, the offset of the temperature prediction is what lowers its evidence given that $M_0$ and $M_2$ have the same number of parameters $q$. $M_2$ also replaces $T_w$, which is a well-defined parameter within its experimental uncertainty for $M_0$, with the emissivity $\epsilon$ which does not have a measured counterpart in these cases, also contributing to the lowering of $M_2$'s evidence. Model $M_2$ needs a perfect prediction of $T_w$ to amount to the evidence of $M_0$ which is not achieved, although the deviation is small to still consider the model good enough compared to the baseline.

<table><caption>Table 7: Bayes factors of the different models against the baseline model $M_0$ for each experimental condition.</caption><thead><tr><th>Experiment ID</th><th>M<sub>1</sub></th><th>M<sub>2</sub></th><th>M<sub>3</sub></th></tr></thead><tbody><tr><td>G4</td><td>-0.3</td><td>-3.27</td><td>-3.22</td></tr><tr><td>G5</td><td>-0.034</td><td>-5.4</td><td>-4.0</td></tr><tr><td>G6</td><td>-0.0024</td><td>-3.17</td><td>-3.06</td></tr><tr><td>G7</td><td>0.01</td><td>-4.87</td><td>-4.0</td></tr></tbody></table>

From this study we cannot conclude that the baseline model is the truth and it indicates what is happening in the flow but rather that the experimental data are consistent with the baseline model as well as the other models. If the true flow was in non-equilibrium, other measurements would be needed in order to find the relevant evidence and truly invalidate the baseline model of thermal equilibrium. The same is applicable to the presence of recombination reactions. Concentrations of $N_2$ in the vicinity of the material surface that cannot be explained by the baseline model could definitely give evidence of recombination reactions and their relative importance. So far, the calibrations of the different models with this particular set of experimental data do not provide enough information about recombination reactions. Fig. 2 depicts the marginal posteriors for $\gamma_{N2}^N$ under the different models that include such parameter. Almost flat posteriors with the same support as the priors are recovered except for the 1T SEB model $M_2$ where the posterior support is reduced. In order for $M_2$ to be able to explain the experimental observations, recombination reaction efficiencies cannot surpass a certain upper limit. Going over such limit would probably throw off the prediction of the surface temperature. Given that the other models have more parameters to control the surface temperature prediction, model $M_2$ is able to place more constraints on $\gamma_{N2}^N$. It can also be appreciated that as the wall temperature measured decreases, the posterior support for $\gamma_{N2}^N$ shrinks toward lower values.

For more insight about model performance, we can look at the resulting residuals of the statistical inverse problem for the measured quantities. The inverse problem we are solving every time we perform an inference amounts to solving the equation

$$y_{\text{obs}} = Y(q) + E, \qquad (26)$$

where the stagnation line model evaluations are replaced by accurate GP surrogates $Y$ for the predictions of recession rates $\dot{s}$, CN densities $\rho_{CN}$, and surface temperatures $T_w$. The error model is chosen as unbiased Gaussian distributions for which we prescribe the different standard deviations $\sigma$ for each measurement. Once the model is calibrated and we have posteriors on the parameters $q$, we propagate them through the model to generate the resulting distribution of $y_{\text{obs}} - Y(q)$. This distribution is what we call the residual and it should match the error model distribution $E$ initially prescribed for each measurement. In this case, we compute the residual as the distribution of $Y(q) - y_{\text{obs}}$ as it has a more intuitive interpretation when looking at the results. This distribution also has to amount to the error distribution given its symmetrical qualities.

Fig. 3 shows the residuals for case G5 for recession rate $\dot{s}^{\text{meas}}$, CN density $\rho_{CN}^{\text{meas}}$ and surface temperature $T_w^{\text{meas}}$ under the different modeling scenarios. The rest of the results for the other experimental cases can be found in Appendix A. Overall, all models retrieve the noise of recession rate and CN density. Model $M_2$ shows the largest deviations in the predictions which are still quite small. The 1T SEB model ($M_2$) tends to over predict recession rates and under predict CN densities in the flowfield. The surface temperature is perfectly captured by the 2T SEB model $M_3$ while the 1T SEB model $M_2$ shows some deviation towards higher surface temperatures. The deviation is ~ 25 K which is rather small compared to the range of surface temperatures we deal with in this problem (> 2,000 K).
<figure><img src="image_3.png" /><figcaption>Figure 2: Posterior marginal distributions of $\gamma_N^{N_2}$ under the different modeling scenarios for G4 (upper left), G5 (upper right), G6 (lower left) and G7 (lower right).</figcaption></figure>

It is interesting to reflect on how the 2T SEB model $M_3$ manages to achieve a lower prediction of the recession rate and surface temperature, and yet perfectly recuperate the error of the CN density. The only possibility for this to happen is through available thermal non-equilibrium mechanisms which are not considered in the 1T SEB model $M_2$. The 2T surface energy balance (Sec. 2) contains an additional term for the conductive contribution of the internal energy to the total energy balance. Through this mechanism, the 2T SEB model finds the combination of collisional efficiency $\alpha$ and energy accommodation factor $\beta$ so that the internal energy contribution cools off the surface. Fig. 4 shows the resulting marginal posterior distributions as well as samples from the joint distribution of the $\alpha$ and $\beta$ parameters. The collisional efficiency $\alpha$ is quite constrained by the experimental data to low values, while the energy accommodation coefficient $\beta$ keeps all the prior support. The resulting distributions are shifted towards values that favour thermal non-equilibrium ($\alpha < 1, \beta < 1$).

Fig. 5 shows the propagated posterior solutions of the stagnation line model for both the 1T SEB and 2T SEB models with their corresponding confidence intervals. Promoting strong thermal non-equilibrium, as seen in the temperature profiles in Fig. 5, makes the surface excite the internal energy modes, consequently lowering the resulting translational surface temperature to match perfectly the experimental data. This feature of the 2T SEB model can easily decouple the surface temperature prediction from the rest of the temperature profile in the flowfield. The result is an overall increment of the flowfield temperatures with respect to the 1T SEB model.

The fact that this change in the temperature profiles can also alter the chemical composition in the flowfield to
<figure><img src="image_4.png" /><figcaption>Figure 3: Residuals of recession rate $s^{\text{meas}}$, CN density $\rho_{\text{CN}}^{\text{meas}}$ and surface temperature $T_w^{\text{meas}}$ for the different models under conditions G5.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Left: 100,000 samples from the joint posterior distribution of $\alpha$ and $\beta$. Right: posterior marginal distributions of $\alpha$ and $\beta$ under conditions G5.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Propagated temperature profiles for the 1T SEB ($M_2$) and 2T SEB ($M_3$) models with their respective confidence intervals. The vertical black line represents the position of the spectrometer measurement.</figcaption></figure>

match the observed $\rho_{CN}$ (Fig. 6) brings forth the importance of another set of assumptions: the gas phase chemistry. In particular, the resulting overall depletion and creation of CN differs between the two models in such a way that a higher CN density is achieved at the spectrometer location for the same concentration of CN at the wall, dictated primarily by the recession rates predicted (schematic representation can be found in Fig. 7). Additional experiments for the proper calibration of the gas phase chemistry would be needed to have a realistic assessment on the evidence of thermal non-equilibrium at the wall.

Still, considering the same chemical mechanism and reaction rates, we could also propose new experimental avenues to find strong evidence, or lack thereof, of thermal non-equilibrium at the wall. In particular, the differences in the CN mole fraction profiles could be exploited further by considering the addition of a spectrometer measurement in a different location along the stagnation line. Such measurement could be performed downstream the current spectrometer measurement, although the accuracy and uncertainty estimation of the resulting CN concentration should be carefully addressed to be able to tell the two models apart. Upstream, the radiative signal would lose intensity and the observation would become more noisy. In any case, it represents a challenge to find stronger evidence in favour of thermal equilibrium or non-equilibrium at the wall.

## 4.2. Bayesian model averaging

All the modeling scenarios presented are so far consistent with the experimental data. We want to incorporate the uncertainty of these multiple modeling choices into the resulting calibrated nitridation reaction efficiencies $\gamma_N^{CN}$. We recall the model assumptions for which we remain uncertain. The baseline model $M_0$ considers a flow in thermal equilibrium for which only one temperature is defined. As heterogeneous chemistry mechanism at the wall it only considers nitridation. Model $M_1$ introduces the possibility of having nitrogen recombination reactions at the surface. Models $M_2$ and $M_3$ incorporate a surface energy balance model for the prediction of surface temperatures which are not imposed in these cases. Further, model $M_3$ assumes the flow to be in a state of thermal non-equilibrium by assuming two distinct thermal baths.

Fig. 8 shows the resulting $\gamma_N^{CN}$ marginal posteriors under the different modeling scenarios. The resulting weighted average posterior is labelled as BMA. Overall, all the posteriors share supports and are consistent with each other. This is to be expected as $\gamma_N^{CN}$ dominates recession rate and surface temperature measurements which are quite accurate in all cases and are not substantially affected by the differential features of the models. We notice a broader support for the thermal equilibrium model with recombination $\gamma_N^{N2}$ in it $M_1$. As a competing mechanism, recombination reactions allow nitridation reactions to have higher efficiencies for the same observed recession rates. The atomic nitrogen available at the surface is more scarce for nitridation reactions when competing with recombination reactions. Higher nitridation efficiencies would still produce the same observed recession rates. The wider support is also extended to the averaged posterior given that model $M_1$ has similar weight as all the other models, consequently widening the support for $\gamma_N^{CN}$. The models with the tightest supports are the baseline $M_0$ and the 1T SEB $M_2$ models. They are
<figure><img src="image_7.png" /><figcaption>Figure 6: Propagated mole fraction profiles for the 1T SEB ($M_2$) and 2T SEB ($M_3$) models with their respective confidence intervals. The vertical black lines represent the position of the spectrometer measurement.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Schematic representation of the modulating effect of the gas phase chemistry in the calibration of the different stagnation line models. The CN concentration at the wall is dominated by the recession rate while the CN density at the spectrometer location is a combined result of the recession rate and the gas phase chemistry.</figcaption></figure>

also the models with the least number of parameters to be inferred, producing more defined peaks and less uncertainty for $\gamma_N^{CN}$.

The most relevant summary statistics for the resulting BMA distributions are gathered in Table 8.

We also extend the model averaging methodology to the computation of an Arrhenius law that captures such model-form uncertainties for a broad span of surface temperatures. In this context, the calibration of an Arrhenius law is performed for each of the models proposed. It takes the form
<figure><img src="image_9.png" /><figcaption>Figure 8: Posterior marginal distributions of $\gamma_N^{CN}$ under the different modeling scenarios and the resulting averaged posteriors for G4 (upper left), G5 (upper right), G6 (lower left) and G7 (lower right).</figcaption></figure>

<table><caption>Table 8: $\gamma_N^{CN}$ posterior statistics of the resulting BMA distributions for the different experimental cases.</caption><thead><tr><th>Experiment ID</th><th>Mean</th><th>95% C. I.</th><th>MAP</th></tr></thead><tbody><tr><td>G4</td><td>0.005</td><td>[0.003, 0.007]</td><td>0.005</td></tr><tr><td>G5</td><td>0.007</td><td>[0.005, 0.013]</td><td>0.007</td></tr><tr><td>G6</td><td>0.012</td><td>[0.008, 0.026]</td><td>0.010</td></tr><tr><td>G7</td><td>0.012</td><td>[0.009, 0.024]</td><td>0.010</td></tr></tbody></table>

$$ \gamma_N^{CN} = A \exp\left(\frac{-T_a}{T_w}\right), \qquad (27) $$

with the pre-exponential factor A, and the activation energy (given in temperature units) $T_a$. For this analysis, both parameters are the same across all experimental cases and considered model-dependent. The resulting marginal posteriors are then averaged a posteriori to compute the Bayesian model averaged Arrhenius law.
The inference problem needs to be re-formulated to account for the new parametrization introduced. For ease of notation, we unpack here a general likelihood function applicable to each model proposed in Sec. 2 by changing the $q$, $y_{obs}$ and $Y$ sets to their specific ones depicted in Table 4. Along with the Arrhenius parameters, we need to add all the nuisance parameters to the calibrations. In this case, each nuisance parameter is represented as a different variable for each experimental case. The compact notation is as follows. For each case (i) in $G = (G_4, G_5, G_6, G_7)$, we define the following likelihood

$$ \mathcal{P}(y_{\text{obs}}|q) \doteq \prod_{i \in G} \mathcal{P}(y_{\text{obs}}^{(i)}|q^{(i)}), \qquad (28) $$

with $\mathcal{P}(y_{\text{obs}}^{(i)}|q^{(i)})$ defined for each model following Eq. (19). Fig. 9 depicts a schematic representation of the methodology. All experimental cases are used to jointly inform the Arrhenius law parameters for each model scenario. Model averaging is performed on such parameters to retrieve the final parameters which incorporate model-form uncertainty.

<figure><img src="image_10.png" /><figcaption>Figure 9: Schematic representation of the relation among the experimental data and model scenarios for the inference of the Arrhenius law parameters.</figcaption></figure>

The prior distribution is taken as independent priors on all parameters, already specified in Sec. 3.1. For $A$ we prescribe the same prior as for $\gamma_N^{CN}$ as it is an equivalent parameter, in this case, $A \sim \mathcal{U}[-4, 0]$. The prior on $T_a$ is chosen wide enough with the same lower bound as the prior on $T_w$, given that we expect the activation energy to be at least as low as the a priori possible lowest $T_w$. It is not easy to prescribe a good upper bound as there is not so much literature on this particular reaction parameter. We choose the upper bound conservatively large based on the results of Suzuki et al. [23] and Helber et al. [10], having $T_a \sim \mathcal{U}[2000, 10000]$.

The results of the Arrhenius parameter marginal posteriors can be seen in Fig. 10. Overall there is good agreement in their support and the resulting BMA distributions follow closely the posterior obtained with the baseline model $M_0$. The summary statistics are gathered in Table 9.

<table><caption>Table 9: Posterior statistics of the resulting BMA distributions for the Arrhenius parameters.</caption><thead><tr><th>Arrhenius parameter</th><th>Mean</th><th>95% C. I.</th><th>MAP</th></tr></thead><tbody><tr><td>A</td><td>0.253</td><td>[0.065, 0.502]</td><td>0.42</td></tr><tr><td>T<sub>a</sub> [K]</td><td>8328.7</td><td>[4996.1, 9944.2]</td><td>9669</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>Figure 10: Calibrated Arrhenius parameters and their weighted average for the different models considered.</figcaption></figure>

Propagating the pre-exponential factor and activation energy posteriors we get an estimation of the variation of the nitridation parameter $\gamma_N^{CN}$ with the surface temperature $T_w$. Fig. 11 shows the different Arrhenius laws and the BMA with the 95% confidence intervals. The right plot of Fig. 11 shows the obtained calibrated laws plotted against the available experimental data. Even though the resulting laws are calibrated with the data from Helber et al. [10], we can project accurate confidence intervals to lower surface temperatures, where the experimental data available differ by several orders of magnitude. The data of Helber et al. are consistent with the experiments from Zhang et al. [6] for lower surface temperatures as well as Suzuki et al.'s [23] data. This apparent agreement enhances our confidence in the tools developed and the analyses performed.

<figure><img src="image_12.png" /><figcaption>Figure 11: Calibrated Arrhenius laws and their weighted average for the different models considered (left). The Arrhenius laws are compared against the available experimental data (right).</figcaption></figure>

A question that can be posed on these results is if defining a model averaged brings anything compared to just using the estimation recovered with the baseline model alone, given that the differences seem to be negligible. Indeed, we can already see from the case by case inferences that the SEB models have several orders of magnitude less evidence weight and this causes the BMA to basically be an average between $M_0$ and $M_1$ which very much agree in their Arrhenius parameters. One could argue that we could have saved ourselves some work by just extrapolating this reasoning from the case by case analyses to the Arrhenius parameters but the answer is not that simple. Arrhenius parameters introduce a new parametrization, one where there is coupling between the prediction of the wall temperature for the SEB models and the Arrhenius parameters themselves, leading to a non-linear constraint on the predicted $\gamma_N^{CN}$ ($\log(\gamma_N^{CN}/A) + T_a/T_w = 0$, with $T_w = T_w(..., \gamma_N^{CN})$) for every pair of Arrhenius parameters. This issue, together with the fact that we infer using all data jointly, further influences the inferences on the parameters we are after. All in
all, it is hard to say a priori what this picture of the different models will look like. It is not until the exercise is done that we come to realise the similarities even under such constraints. Given the logarithmic nature of the scale of the nitridation parameter, averaging can be useful even if small deviations from model to model are found, considering that a small change in scale can mean a sizeable change in $\gamma_{\text{N}}^{\text{CN}}$ and, therefore, the recession rate to be predicted using this calibrated law. To illustrate this further, we can compute the recession rates resulting from using the parameters calibrated assuming the baseline model and the BMA model, respectively. The resulting predictions of recession rates in logarithmic scale for the same set of inputs except the Arrhenius parameters are

$$ \log_{10} s^{\text{baseline}} = -8.16839, \qquad (29) $$

$$ \log_{10} s^{\text{BMA}} = -8.24106, \qquad (30) $$

which when transformed to natural scale, we get a relative difference of ~17% in recession rates, that is one standard deviation of the experimental measurements. These differences can certainly impact the predictive capabilities of the model.

Another important question that arises regarding the results of Fig. 10 is the fact that the marginal distributions, especially those of $T_a$, seem to be constrained by their chosen priors. As said, the choice of priors for the Arrhenius parameters is challenging, given the poor literature on the subject. Even though the supports were chosen conservatively large, the posteriors seem to push their supports towards the higher end of the possible values for $T_a$ (Fig. 10). It is interesting to look at the joint marginal posterior of both A and $T_a$ given that they are strongly correlated. Fig. 12 shows a scatter plot of some samples taken from the joint marginal posterior. The current prior placed on $T_a$ leads to a wall where the samples suddenly stop. Conversely, a wider prior relaxes the situation and now the limit of the posterior support is dictated by the prior placed on A. An Arrhenius law is no other thing than a macroscopic model with physically motivated parameters [42] that describe the frequency of collisions and the portion of those with enough energy to actually react based on the statistics of such encounters. As such, we find that the model becomes degenerate as it describes the same chemistry for large frequencies of collisions and large activation energies and viceversa. This issue is seen in Fig. 12 as the joint posterior goes on and on towards larger values and it only stops when specific constraints are placed, such as the ones stemming from the priors. The initially superficial question becomes more intricate, and it remarks the importance of looking at the microscopic nature of the flow to prescribe physically consistent priors that can tackle this kind of problems, therefore removing the possible bias that can be induced by selecting arbitrary priors.

<figure><img src="image_13.png" /><figcaption>Figure 12: 10,000 samples from the joint marginal posterior of A and $T_a$ for two different choices of priors for $T_a$.</figcaption></figure>
## 5. Concluding remarks and outlook

This work presents a calibrated carbon nitridation model that takes into account experimental, parametric and model-form uncertainties. Experimental and parametric uncertainties are incorporated through dedicated Bayesian inference frameworks. Model-form uncertainties are quantified by performing model selection and averaging on the Bayesian posteriors. The final result is a calibrated Arrhenius law for nitridation reaction efficiencies spanning a range of surface temperatures whose parameters also include uncertainties stemming from different modeling assumptions about the physico-chemical nature of the plasma flow and gas-surface interface.

Different modeling scenarios are proposed in this work. We first define a baseline model against which to compare a set of alternative models in their predictive performances and complexities. The baseline model accounts for flow thermal equilibrium assumptions with nitridation as only surface mechanism; it also imposes the surface temperature, treating it as a free parameter. The model is calibrated and its predictive performance tested on the experimental dataset. The modeling assumptions hold for the conditions of the experiments giving us a solid base on which to build our subsequent hypothesis testing study. The alternative models relax several assumptions embedded in the baseline model. The first alternative model $M_1$ adds nitrogen recombination reactions through an efficiency parameter that must be also calibrated. Model $M_2$ builds on model $M_1$ and includes a surface energy balance to predict the surface temperature instead of treating it as a free parameter to be calibrated. Finally, alternative model $M_3$ considers nitrogen recombination reactions at the gas-surface interface and remains agnostic about the thermal state of the gas and gas-surface interface for which additional balance equations must be solved.

All these modeling choices are assessed in a hypothesis testing study. Overall no strong evidence is found against the baseline model and none of the considered models can be discarded. The calibrations of models with surface recombination can define an upper limit for the recombination efficiency parameter at best. The calibrated solutions of the thermal non-equilibrium model point to strong thermal non-equilibrium at the gas-surface interface, although the evidence is weak.

Model averaging is then performed to incorporate such model-form uncertainties in a Bayesian setting. The marginal posteriors of nitridation efficiencies parameters for each experimental condition are weighted according to their marginalized likelihoods assuming equal prior probabilities for all models. The results show how for cases G5-G7 the tail of the nitridation reaction efficiencies distributions towards higher values expands, enlarging the space of possible values for $\gamma_N^{CN}$. This growth of the $\gamma_N^{CN}$ uncertainty is mainly promoted by model $M_1$ which includes recombination reactions while freely imposing the surface temperature. Furthermore, an averaged Arrhenius law is also obtained by weighting each Arrhenius model parameter by their corresponding marginalized likelihoods for each individual model. This work offers complete and accurate characterization of the uncertainties affecting the calibration of a carbon nitridation model, updating the current state-of-the-art databases.

The outcomes of this work direct our research efforts in this field towards two different aspects: generating new experimental data with which to be able to rule out some models, finding definitive evidence of thermal non-equilibrium as well as learning surface recombination; and incorporating priors for kinetic models, such as Arrhenius laws, that consistently reflect the microscopic nature of the flow without necessarily incurring into model degeneracies entirely controlled by the prior distributions chosen.

## Acknowledgements

This work is fully funded by the European Commission H2020 programme, through the UTOPIAE Marie Curie Innovative Training Network, H2020-MSCA-ITN-2016, Grant Agreement number 722734. The authors would also like to thank Michele Capriati for his guidance in understanding and simulating the 2T SEB model. Prof. Serge Prudhomme is acknowledged for providing insightful comments which improved this work.

## Appendix A. Complementary results

This appendix contains the results concerning the cases not included in Section 4.
## Appendix A.1. Prediction residuals

Figs. A.13-A.15 show the residuals for cases G4, G6 and G7 for recession rate $\dot{s}$, CN density $\rho_{CN}$ and surface temperature $T_w$ under the different modeling scenarios. We recall that we only take into account the recession rate measurements for case G6 for which the CN density residual is not shown.

<figure><img src="image_14.png" /><figcaption>Figure A.13: Residuals of recession rate $\dot{s}$, CN density $\rho_{CN}$ and surface temperature $T_w$ for the different models under conditions G4.</figcaption></figure>

## Appendix A.2. Calibrated flowfield solutions

Figs. A.16-A.21 show the propagated posterior solutions of the stagnation line model for both the 1T SEB and 2T SEB models with their corresponding confidence intervals for cases G4, G6 and G7.

## Appendix A.3. Posterior distributions of the surface parameters for the 2T model

Fig. A.22 shows samples of the 2D joint posterior distributions of the parameters $\alpha$ and $\beta$ in the 2T surface balance model, as well as their marginal posteriors for cases G4, G6, and G7.
<figure><img src="image_15.png" /><figcaption>Figure A.14: Residuals of recession rate $\dot{s}$ and surface temperature $T_w$ for the different models under conditions G6.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure A.15: Residuals of recession rate $\dot{s}$, CN density $\rho_{CN}$ and surface temperature $T_w$ for the different models under conditions G7.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure A.16: Propagated temperature profiles for the 1T SEB ($M_2$) and 2T SEB ($M_3$) models with their respective confidence intervals for case G4. The vertical black line represents the position of the spectrometer measurement.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure A.17: Propagated mole fraction profiles for the 1T SEB ($M_2$) and 2T SEB ($M_3$) models with their respective confidence intervals for case G4. The vertical black lines represent the position of the spectrometer measurement.</figcaption></figure>

## References

[1] C. Davies, M. Arcadi, Planetary Mission Entry Vehicles Quick Reference Guide. Version 3.0., Tech. rep., NASA Technical Report 2006-3401 (2006).

[2] B. Helber, A. Turchi, J. B. Scoggins, A. Hubin, T. E. Magin, Experimental investigation of ablation and pyrolysis processes of carbon-phenolic ablators in atmospheric entry plasmas (2016), International Journal of Heat and Mass Transfer 100 (2016) 810–824. doi:https://doi.org/10.1016/j.ijheatmasstransfer.2016.04.072.

[3] H. W. Goldstein, The reaction of active nitrogen with graphite (1964), The Journal of Physical Chemistry 68 (1) (1964) 39–42. doi:10.1021/j100783a007.

[4] C. Park, D. Bogdanoff, Shock-tube measurements of nitridation coefficients of solid carbon (2006), Journal of Thermophysics and Heat Transfer 20 (3).
<figure><img src="image_19.png" /><figcaption>Figure A.18: Propagated temperature profiles for the 1T SEB ($M_2$) and 2T SEB ($M_3$) models with their respective confidence intervals for case G6. The vertical black line represents the position of the spectrometer measurement.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure A.19: Propagated mole fraction profiles for the 1T SEB ($M_2$) and 2T SEB ($M_3$) models with their respective confidence intervals for case G6. The vertical black lines represent the position of the spectrometer measurement.</figcaption></figure>

<footnote>[5] T. Suzuki, K. Fujita, T. Sakai, Graphite nitridation in lower surface temperature regime (2010), Journal of Thermophysics and Heat Transfer 24(1).</footnote>

<footnote>[6] L. Zhang, D. A. Pejaković, J. Marschall, M. Dougherty, D. G. Fletcher, Laboratory investigation of the active nitridation of graphite by atomic nitrogen (2012), Journal of Thermophysics and Heat Transfer 26 (1) (2012) 10–21. doi:10.2514/1.T3612.</footnote>

<footnote>[7] C. Park, D. W. Bogdanoff, Shock-tube measurement of nitridation coefficient of solid carbon (2006), Journal of Thermophysics and Heat Transfer 20 (3) (2006) 487–492. doi:10.2514/1.15743.</footnote>

<footnote>[8] Uncertainty quantification of a graphite nitridation experiment using a bayesian approach (2011), Experimental Thermal and Fluid Science 35 (8) (2011) 1588–1599. doi:https://doi.org/10.1016/j.expthermflusci.2011.07.010.</footnote>

<footnote>[9] A. del Val, Bayesian calibration and assessment of gas-surface interaction models and experiments in atmospheric entry plasmas, Ph.D. thesis, IPP/VKI (2021).</footnote>
<figure><img src="image_21.png" /><figcaption>Figure A.20: Propagated temperature profiles for the 1T SEB ($M_2$) and 2T SEB ($M_3$) models with their respective confidence intervals for case G7. The vertical black line represents the position of the spectrometer measurement.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure A.21: Propagated mole fraction profiles for the 1T SEB ($M_2$) and 2T SEB ($M_3$) models with their respective confidence intervals for case G7. The vertical black lines represent the position of the spectrometer measurement.</figcaption></figure>

<footnote>[10] B. Helber, A. Turchi, T. E. Magin, Determination of active nitridation reaction efficiency of graphite in inductively coupled plasma flows (2017), Carbon 125 (2017) 582–594.</footnote>

<footnote>[11] A. Munafò, Multi-scale models and computational methods for aerothermodynamics, Ph.D. thesis, VKI/ECP (2014).</footnote>

<footnote>[12] J. B. Scoggins, V. Leroy, G. Bellas-Chatzigeorgis, B. Dias, T. E. Magin, Mutation++: Multicomponent thermodynamic and transport properties for ionized gases in c++ (2020), SoftwareX 12 (2020) 100575. doi:https://doi.org/10.1016/j.softx.2020.100575.</footnote>

<footnote>[13] B. Helber, Material response characterization of low-density ablators in atmospheric entry plasmas, Ph.D. thesis, VUB/VKI (2016).</footnote>

<footnote>[14] F. Panerai, Aerothermochemistry characterization of thermal protection systems, Ph.D. thesis, Università degli Studi di Perugia, von Karman Institute for Fluid Dynamics (2012).</footnote>

<footnote>[15] A. Viladegut, Assessment of Gas-Surface Interaction Modelling for Lifting Body Re-Entry Flight Design, Ph.D. thesis, VKI/UPC (2017).</footnote>

<footnote>[16] A. Klomfass, S. Müller, A quasi-onedimensional approach for hypersonic stagnation-point flows (1996), Technical report, Internal Report.</footnote>
<figure><img src="image_23.png" /><figcaption>Figure A.22: Left: 100,000 samples from the joint posterior distributions of α and β. Right: posterior marginal distributions of α and β under conditions (from top to bottom) G4, G6 and G7.</figcaption></figure>

<footnote>[17] K. Kitamura, E. Shima, Towards shock-stable and accurate hypersonic heating computations: A new pressure flux for ausm-family schemes (2013), Journal of Computational Physics 245 (2013) 62–83. doi:https://doi.org/10.1016/j.jcp.2013.02.046.</footnote>

<footnote>[18] G. Bellas-Chatzigeorgis, P. F. Barbante, T. E. Magin, Energy accommodation coefficient calculation methodology using state-to-state catalysis applied to hypersonic flows (2020), AIAA Journal 58 (1) (2020) 278–290. doi:10.2514/1.J058543.</footnote>

<footnote>[19] P. Chen, N. Zabaras, I. Bilionis, Uncertainty propagation using infinite mixture of gaussian processes and variational bayesian inference (2015), Journal of Computational Physics 284 (2015) 291–333.</footnote>
[20] A. Turchi, D. Bianchi, F. Nasuti, R. Marocco, Ablative material behavior in oxygen/methane thruster environment, in: 43rd AIAA Thermo-physics Conference, 2012.

[21] J. Lachaud, J. Scoggins, T. Magin, M. Meyer, N. Mansour, A generic local thermal equilibrium model for porous reactive materials submitted to high temperatures (2017), International Journal of Heat and Mass Transfer 108 (2017) 1406-1417.
doi:https://doi.org/10.1016/j.ijheatmasstransfer.2016.11.067.
URL https://www.sciencedirect.com/science/article/pii/S0017931016312480

[22] B. Dias, Thermal ablation and radiation modeling of meteor phenomena, Ph.D. thesis, Université Catolique de Louvain/VKI (2020).

[23] T. Suzuki, K. Fujita, T. Sakai, Graphite nitridation in lower surface temperature regime (2010), Journal of Thermophysics and Heat Transfer 24 (1) (2010) 212-215. doi:10.2514/1.43265.

[24] S. V. Zhluktov, T. Abe, Viscous shock-layer simulation of airflow past ablating blunt body with carbon surface (1999), Journal of Thermo-physics and Heat Transfer 13(1) (1999) 50-59.

[25] A. Turchi, A gas-surface interaction model for the numerical study of rocket nozzle flow sover pyrolyzing ablative materials, Ph.D. thesis, Sapienza Università di Roma (2013).

[26] D. Bianchi, F. Nasuti, E. Martelli, Coupled analysis of flow and surface ablation in carbon-carbon rocket nozzles (2009), Journal of Spacecraft and Rockets 46(3).

[27] D. Bianchi, F. Nasuti, M. Onofri, E. Martelli, Thermochemical erosion analysis for graphite/carbon-carbon rocket nozzles (2011), Journal of Propulsion and Power 27(1).

[28] L. Landau, E. Teller, Theory of sound dispersion, Physikalische Zeitschrift der Sowjetunion, 1936.

[29] B. Halpern, D. Rosner, Chemical energy accommodation at catalyst surfaces. flow reactor studies of the association of nitrogen atoms on metals at high temperatures (1978), Journal of the Chemical Society, Faraday Transactions 74 (1978) 1883-1912.

[30] M. Capriati, G. Bellas-Chatzigeorgis, A. Turchi, B. Helber, T. E. Magin, Thermal non-equilibrium modeling for ablative gas-surface interaction, in: HiSST: 2nd International Conference on High-Speed Vehicle Science and Technology, 2020.

[31] G. Bellas, Development of advanced gas-surface interaction models for chemically reacting flows for re-entry conditions, Ph.D. thesis, VKI/Politecnico di Milano (2018).

[32] D. S. Sivia, J. Skilling, Data Analysis - A Bayesian Tutorial, 2nd Edition, Oxford Science Publications, Oxford University Press, 2006.

[33] H. Haario, E. Saksman, J. Tamminen, An adaptive metropolis algorithm (1998), Bernoulli 7 (1998) 223-242.

[34] H. Jeffreys, Theory of Probability, 3rd Edition, Oxford, Oxford, England, 1961.

[35] R. E. Kass, A. E. Raftery, Bayes factors (1995), Journal of the American Statistical Association 90 (430) (1995) 773-795.
doi:10.1080/01621459.1995.10476572.

[36] G. Schwarz, Estimating the dimension of a model (1978), The Annals of Statistics 6 (1978) 461-464.

[37] H. Akaike, A new look at the statistical model identification (1974), IEEE Transactions on Automatic Control 19 (6) (1974) 716-723.
doi:10.1109/TAC.1974.1100705.

[38] J. Ding, V. Tarokh, Y. Yang, Model selection techniques: An overview (2018), IEEE Signal Processing Magazine 35 (6) (2018) 16-34.
doi:10.1109/MSP.2018.2867638.

[39] T. C. Chamberlin, The method of multiple working hypotheses (1965), Science 148 (3671) (1965) 754-759.
doi:10.1126/science.148.3671.754.

[40] J. A. Hoeting, D. Madigan, A. E. Raftery, C. T. Volinsky, Bayesian model averaging: A tutorial (1999), Statistical Science 14 (4) (1999) 382-401.

[41] H. V. Roberts, Probabilistic prediction (1965), Journal of the American Statistical Association 60 (309) (1965) 50-62.
doi:10.1080/01621459.1965.10480774.

[42] W. G. Vicenti, C. H. Kruger, Introduction to Physical Gas Dynamics, Krieger Pub Co., 1975.