# A semi-analytical model for Plasma-Assisted Combustion: Application to NH3-H2-air mixtures

Harry Gillingham, Thomas Lesaffre, Nicolas Barléon

► To cite this version:

Harry Gillingham, Thomas Lesaffre, Nicolas Barléon. A semi-analytical model for Plasma-Assisted Combustion: Application to NH3-H2-air mixtures. Combustion and Flame, 2025, 283, pp.114528. 10.1016/j.combustflame.2025.114528. hal-05320527

HAL Id: hal-05320527

https://hal.science/hal-05320527v1

Submitted on 17 Oct 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
A semi-analytical model for Plasma-Assisted
Combustion: Application to NH₃-H₂-air mixtures

Harry **Gillingham**<sup>a,b</sup>, Thomas **Lesaffre**<sup>a</sup>, Nicolas **Barléon**<sup>a,*</sup>

<sup>a</sup>CERFACS, 42 Avenue Gaspard Coriolis, 31057 Toulouse, Cedex 1, France

<sup>b</sup>KTH Royal Institute of Technology, Stockholm 114 28, Sweden

Abstract

The transition to low-carbon energy systems results in strong interest for ammonia and hydrogen as alternative carbon-free fuels. Plasma-Assisted Combustion (PAC) offers a promising pathway to overcome their limitations in terms of combustion properties and pollutant emissions. However, incorporating non-equilibrium plasma effects into turbulent flame simulation remains computationally expensive: efficient and accurate modelling approaches are therefore essential. A Semi-Analytical (SA) model for PAC, named **PACMIND-SA**, enabling cost-effective prediction of discharge energy partitioning in ammonia–hydrogen–air mixtures, has been developed. Unlike some existing phenomenological approaches, PACMIND-SA does not require detailed time-resolved kinetic simulations for discharge energy branching evaluation, while detailed kinetics is only required once for the initial processes identification. Instead, it estimates the energy branching into chemical, thermal, and vibrational channels using electron-impact reaction data and a reaction rate-based analysis, informed by the electron energy distribution function. The model is validated against detailed plasma-kinetic simulations for various conditions, including single-pulse discharges in air, ignition of NH₃-H₂ mixtures, NOₓ formation in various mixtures, MILD combustion, oxy-combustion and plasma-assisted ammonia reforming. PACMIND-SA accurately reproduces key thermo-chemical trends such as ignition delays, radical formation, and pollutant emissions (NOₓ). In particular, the error on ignition delay time for a single pulse against the detailed mechanism is found between -9% and +6% over a wide range of E/N ∈ [100 – 500] Td and NH₃-H₂ blending ratios (from pure ammonia to pure hydrogen). Tests performed on multiple pulses, 100 pulses at 20 kHz, showed a very good agreement with less than 1% error on the ignition delay time. Despite the assumption of constant thermochemical states during the pulse, the model shows excellent agreement across a wide range of discharge energies and fuel compositions. Additionally, the overall methodology can be applied to other type of fuel, and thus is not limited to the applied NH₃-H₂ cases considered in this work, as shown in this work with CH₄-air mixture. The low cost associated to the model parameter evaluation (O(ms)), three order of magnitudes lower than with the original PACMIND methodology (O(s)), could allow their on-the-fly evaluation in CFD solver, leveraging the need for lookup tables within a limited state space. Tools to generate the model coefficients are available at https://chemistry.cerfacs.fr/en/plasma/.

Keywords: Reduced-Order Model, Nanosecond Repetitively Pulsed Discharges, Chemistry, Ammonia-Hydrogen Blends, Methane

* Nicolas Barléon: CERFACS, 42 Avenue Gaspard Coriolis, 31057 Toulouse, Cedex 1, France

Novelty and significance statement

1

2

3

4

5

6

7

8

9

This work introduces a novel semi-analytic formu-
lation for a phenomenological plasma-assisted com-
bustion model that does not require detailed kinetic
simulations any more once the main processes are
identified. Instead, it incorporates the information
from detailed chemical kinetic mechanisms directly,
enabling efficient and cost-effective discharge energy
branching evaluation. This approach significantly re-

duces computational overhead, facilitating integration
into CFD solvers. Its generalizable framework al-
lows for adaptation to different fuel/mixture composi-
tions (ammonia, hydrogen and methane tested in this
work), enhancing its practical use, provided that de-
tailed kinetic mechanisms are available. This new
model marks a step toward scalable and reliable mod-
eling of plasma-assisted combustion.

1. Introduction

Combustion technologies are increasingly shifting
toward carbon-free fuels such as ammonia. However,
the use of ammonia for large-scale energy production

10

11

12

13

14

15

16

17

18

19

20

21

e-mail: barleon@cerfacs.fr (Nicolas Barleon)
through its oxidation is still at the research level, still needing to overcome the high reactive nitrogen oxide emissions and low reactivity [1]. A common strategy to enhance the reactivity of ammonia-air mixtures involves blending ammonia with small amounts of hydrogen, either via ammonia cracking or pure hydrogen addition [2, 3, 4], even though this can lead to increased pollutant emissions. Plasma-Assisted Combustion (PAC) is intensively studied for the unique properties of non-equilibrium plasma to enhance combustion process [5, 6], including ignition [7], flammability limit extension, flame stabilization [8], mitigation of thermo-acoustic instabilities [9, 10] and pollutant reduction [11]. Recent works highlighted the plasma potential to accelerate ammonia oxidation [12, 13, 14], extend its flammability limits and reduce its lean and rich Minimum Ignition Energy [15], increase the ammonia lean blowout limit [16, 17] and reduce NOx emissions [17, 11]. Despite these encouraging findings, the current strategies for Plasma-Assisted Ammonia Combustion (PAAC) do not reach the expected performance as obtained with hydrocarbon fuels [18]. To improve PAAC performances, a deeper understanding of the physico-chemical processes is required, for which numerical simulations can provide critical insights.

Recent numerical studies on PAAC have mainly employed detailed kinetic models in zero-dimensional (0D) reactors [19, 20, 21]. While these simulations offer valuable insights into plasma-induced N₂O formation in ammonia-containing mixtures [22], they fail to capture multidimensional effects that influence flame and plasma behaviours and its interaction with the surrounding flow. To address these limitations, reduced chemical mechanisms have been developed [23, 24], enabling fully coupled plasma-combustion simulations that investigate plasma-flame-flow interactions [25, 26]. However, such detailed approaches remain computationally prohibitive for complex three-dimensional turbulent configurations, despite recent progress [27].

To bridge this gap and enable modeling of practical burner systems, phenomenological approaches have emerged over the past decade. Initiated by the work of Castela et al. [28], these methods are based on the concept that plasma energy can be partitioned into three primary channels: i) fast gas heating, ii) slow gas heating and iii) rapid radical production, which can be incorporated as source terms into the reactive compressible Navier-Stokes equations. Originally formulated in air [28], this approach has since been extended to the combustion of fuel-air mixtures [29, 30, 31] or plasma-assisted fuel reforming [32], drawing on analyses of detailed plasma-assisted combustion kinetics [33]. In particular in [34], a streamer code has been coupled to a phenomenological discharge energy branching to improve the spatio-temporal accuracy of the energy deposition. However, this approach for computing the plasma power density is computationally expensive and restricted to small-scale canonical configurations, making it impractical for large-scale burner applications and thereby limiting the overall benefits of the phenomenological model. Alternatively, plasma-affected reactive mixture composition could be modified based on zero-dimensional plasma-combustion simulation to perform DNS of turbulent flame influenced by plasma [35]. Whereas these approaches have been widely used to model plasma-actuated turbulent flames [10, 36, 37, 38], they required a priori knowledge of the thermochemical state of the plasma region in order to perform the associated detailed simulations that can be computationally expensive (i.e., a few minutes).

In this work, a semi-analytical approach is proposed to alleviate the need for pre-computing the model parameters, building upon the previously developed PACMIND model [29]. The approach is based on the concept of discharge power branching, which evaluates how electron energy is distributed among various pathways, including the excitation and ionization of heavy molecules. Subsequently, an analytical evaluation of reaction rates, obtained from a detailed kinetic mechanism, is used to estimate the final ground-state species produced following the plasma discharge as well as the fast and slow gas heating. The methodology is applied to NH₃-H₂ mixtures, with optional water vapour addition corresponding to burnt gases, without loss of generality, and is particularly suited for non-equilibrium plasma such as Nanosecond Repetitively Pulsed (NRP) discharges, Micro-Wave (MW) discharges or Dielectric Barrier Discharges (DBD) operating in the range of E/N ∈ [100–500] Td. The methodology is also evaluated for CH₄ and a different mechanism to validate its wide fuel applicability and mechanism flexibility.

The paper is organized as follows. Section 2 presents the numerical framework and the detailed mechanism used to develop and validate the proposed model. Section 3 provides a detailed description of the new methodology, named PACMIND-SA, for computing the model parameters, along with an analysis of discharge energy partitioning into chemical and thermal channels. Section 4 validates the model across various mixtures by comparing its predictions with detailed kinetic simulations. Finally, Section 5 examined the model's cost and the remaining challenges associated with it.

## 2. Numerical framework

### 2.1. Governing equations for detailed plasma simulations

The mass fractions $Y_k$ and the gas temperature T of a constant volume reactor are governed by:

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112
$$m \frac{dY_k}{dt} = V W_k \dot{\omega}_k \qquad (1)$$

$$mc_v \frac{dT}{dt} = V \sum_k W_k \dot{\omega}_k u_k + \dot{R}_{VT}^k \qquad (2)$$

where m is the mass of the mixture, t the time, $\dot{\omega}_k$ the kth species molar production rate considering both combustion and plasma reactions, $W_k$ its molar mass, $u_k$ its internal energy, V the reactor volume and $c_v$ the mass heat capacity at constant volume and $\dot{R}_{VT}^k$ is the vibrational-translational relaxation resulting in a slow gas heating. To account for non-equilibrium plasma kinetics, the Boltzmann equation solver named Bolos (based on BOLSIG+ algorithm [39]) is employed to obtain the Electron Energy Distribution Function (EEDF) from which electron temperature, electron impact collision rates and other swarm parameters are computed based on a set of complete and consistent cross-sections from the LxCat database [40]. These plasma-combustion equations are solved in a self-consistent way using the unified Cantera-Bolos solver presented in [24] with a dedicated plasma reactor. To account for vibrational excitation, only the mean vibration energy density $e_{vib}^k$ of the kth species is considered following

$$\frac{de_{vib}^k}{dt} = \dot{e}_{vib,k}^p - \dot{R}_{VT}^k \qquad (3)$$

where $\dot{e}_{vib,k}^p$ is the vibrational excitation from electron impact reactions or vibrational energy residuals. This approach has been considered for N₂, O₂, H₂ and NH₃ instead of accounting for each vibrational states showing good agreement with literature reference as shown in Section 2.2.

## 2.2. Detailed mechanism validation

This work is based upon the validated NH₃-air plasma-assisted combustion mechanism from [19] extended to H₂ in [41]. The mechanism originally validated on ZDPlasKin-ChemKin simulation has been converted to Cantera format using a mean vibrational energy approach for the vibrational excitation. In order to validate this step, the simulation of a single pulse is performed and compared with the results presented in [20]. A stoichiometric NH₃-air mixture is considered with an initial charge species concentration of 9.8 × 10¹³ m⁻³. A square voltage pulse is applied with a constant reduced electric field of 250 Td, which is switched off when the discharge energy density $e^p$ reaches a value of 400 J · m⁻³. The discharge energy density is computed from

$$e^p = \int_t e \mu_e n_e E^2 dt \qquad (4)$$

where e is the elementary charge, $\mu_e$ the electron mobility, $n_e$ the electron density and E the electric field. The results of the major excited species and radicals

are shown in Fig. 1. The time evolution of the radicals and excited species shows a very good agreement, including the electrons, which shows that the pulse duration is also retrieved. For the temperature, the initial fast gas heating phase is well captured, and there is an overall very good agreement for its temporal evolution.

<figure><img src="image_2.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(b)</figcaption></figure>

Fig. 1: Time evolution of a) the main radicals and the gas temperature and b) the main excited species for the detailed kinetic simulation of (—) this work and (---) the reference case from [20].

# 3. Model development

## 3.1. Model description and objectives

The model developed in this work, PACMIND-SA, aims to model the main thermo-chemical impact of non-equilibrium plasma discharge in a NH₃/H₂/N₂/O₂/H₂O mixture. To do so, the overall plasma energy density $e^p$ is divided in three energy channels following the ideas presented in [28, 29]: namely the fast chemical energy effects $e_{chem}^p$, the fast gas heating $e_{heat}^p$ and slow gas heating $e_{vib}^p$:

$$e^p = e_{chem}^p + e_{heat}^p + e_{vib}^p \qquad (5)$$
As done for methane-air combustion in [29], the chemical effect is divided in several processes $j \in [1, P]$ from which the chemical energy density reads

$$e_{\text{chem}}^{p} = \sum_{j} e_{\text{chem},j}^{p} \qquad (6)$$

For the vibrational energy, we considered the excitation of N₂, O₂, H₂ and NH₃. The vibrational energy of NH₃ is relaxed in 10 to 100 ns [19], while for N₂, O₂ and H₂, the relaxation time is of the order of the ms [42]. Hence, the vibrational energy of NH₃ is modeled as fast gas heating. In the studied conditions (E/N ≥ 100 Td), vibrational excitations of oxygen and hydrogen molecules are marginal so that only the vibrational excitation of N₂ is retained in the final model. These assumptions could be easily relaxed and do not alter the generality of the model development.

The model has been implemented in a zero-dimensional reactor in order to compare its accuracy with the detailed kinetic mechanism. The governing equations for species mass fractions, gas temperature and mean non-equilibrium vibrational energy density read:

$$m \frac{dY_k}{dt} = V W_k (\dot{\omega}_k + \dot{\omega}_k^{\text{mod}}) \qquad (7)$$

$$mc_v \frac{dT}{dt} = V \sum_k W_k \dot{\omega}_k u_k + \dot{R}_{VT}^k + \dot{\epsilon}_{\text{heat}}^p \qquad (8)$$

$$\frac{\mathrm{d}e_{\text{vib}}^{\mathrm{N}_2}}{\mathrm{d}t} = \dot{\epsilon}_{\text{vib,N}_2}^{\mathrm{p}} - \dot{R}_{\text{VT}}^{\mathrm{N}_2} \qquad (9)$$

where $\dot{\omega}_k^{\text{mod}} = \sum_{j \in P} \nu_{kj} \dot{Q}_j$ is the chemical source terms resulting from the plasma discharge with $\dot{Q}_j$ the molar production rate from the plasma processes included in the model, $\dot{\epsilon}_{\text{heat}}^p$ corresponds to the energy source term from the fast gas heating and the treatment of vibrational excitation is similar to the detailed model presented in Sec. 2.1 except that only N₂(v) is considered in the present model. For consistency in the comparison of the model and the detailed PAC simulations, the same NH₃/H₂ oxidation mechanism from [43] was used.

## 3.2. Chemical pathway analysis

As discussed in the previous section, the model rely on a set of global processes that have to be defined. A chemical pathway analysis has been conducted on the detailed plasma-combustion mechanism in NH₃ and NH₃-H₂ blends to identify the dominant chemical processes resulting from the plasma discharge. Thus, the model only requires a detailed mechanism to perform this pre-analysis task. The main chemical effects of a plasma discharge in air are the ultra fast O₂ and N₂ dissociation and NO formation, as previously shown in [29]. The O₂ dissociation process primarily occurs due to electron impact reaction or reactive collisional quenching of N₂* leading to

the overall process (P1).

$$O_2 \xrightarrow{\alpha_1} 2O \qquad (P1)$$

The N₂ dissociation mainly results from electron impact reactions forming highly excited pre-dissociating singlet electronic states

$$e^{-} + N_{2} \longrightarrow e^{-} + N + N(^{2}D) \qquad (10)$$

which is summarized by the global process (P2).

$$N_2 \xrightarrow{\alpha_2} 2N \qquad (P2)$$

The excited nitrogen N(²D) formed in reaction (10) is prone to react with the oxygen molecule to form NO molecule modeled through the global process (P3). Nitrous oxide N₂O can also be formed during plasma process via collisional quenching of N₂(A₃) with O₂ [44] leading to the global process (P4). Regarding the NO₂, the plasma seems to have a negligible effect on its production with the selected mechanism. Thus, it is not taken into account in the present model.

$$O_2 + N_2 \xrightarrow{\alpha_3} 2NO \qquad (P3)$$

$$O_2 + 2N_2 \xrightarrow{\alpha_4} 2N_2O \qquad (P4)$$

During the pulse, ammonia is mainly dissociated via electron impact dissociation, either leading to NH₂ + H or to NH + H₂. Additionally, the ammonia dissociation can also come from the quenching of excited states of N₂, electron recombination with NH₃⁺ or via the production of NH₄⁺ followed by its ion-ion recombination with O⁻, O₂⁻ or OH⁻. This is summarized by two global processes (P5) and (P6).

$$NH_3 \xrightarrow{\alpha_5} NH_2 + H \qquad (P5)$$

$$NH_3 \xrightarrow{\alpha_6} NH + H_2 \qquad (P6)$$

Whereas N₂* is the dominant excited species produced in air-dominated mixture in the range of conditions of interest for NRP discharge, other excited species such as O(¹D, ¹S) or N(²D) are produced during the discharge. These excited atoms quickly react with ammonia following two main reactions:

$$NH_3 + N(^2D) \longrightarrow NH + NH_2 \qquad (11)$$

$$NH_3 + O(^1D) \longrightarrow NH_2 + OH \qquad (12)$$

which can be globally modeled through two overall processes:

$$NH_3 + \frac{1}{2}N_2 \xrightarrow{\alpha_7} NH + NH_2 \qquad (P7)$$

$$NH_3 + \frac{1}{2}O_2 \xrightarrow{\alpha_8} NH_2 + OH \qquad (P8)$$

The model is extended to H₂-containing mixtures
1 to allow its use in NH₃–H₂ blends. As for ammo-
nia, H₂-dissociation process (P9) is considered in the
2 model (either through electron impact reaction or col-
3 lisional quenching of N₂*) as well as OH formation
4 (P10) following reactions with excited atomic oxygen.
5

$$
H_2 \xrightarrow{\alpha_g} 2 H \qquad (P9)
$$

$$
6 \qquad H_2 + 0.5 O_2 \xrightarrow{\alpha_{10}} OH + H \quad (\text{P10})
$$

7 In burnt gases, H₂O is the only major species added
8 to the mixture and the plasma effect can be modeled
9 as in [29] following the global process (P11).

$$
H_2O \xrightarrow{\alpha_{11}} OH + H \quad (\text{P11})
$$

10 Finally, for lean mixtures, O₂ is still present in the
11 burnt gases. This allows reactions between excited
12 atomic oxygen and water, leading to the global pro-
13 cess (P12).

$$
H_2O + \frac{1}{2} O_2 \xrightarrow{\alpha_{12}} 2 OH \quad (\text{P12})
$$

14 The final model ends with 13 species present in con-
15 ventional NH₃–H₂ oxidation mechanisms along with
16 12 processes summarized in Tab. 1.

3.3. *Production rate of the global processes*

18 A method for computing the energy fraction de-
19 posited in the jth chemical process αj has been pre-
20 viously presented in [29] and is referred to as the
21 PACMIND model. This method is based on the post-
22 processing of a detailed simulation to calculate the
23 mean rate of progress of the global process j from
24 which αj can be computed as

$$
\alpha_j = \frac{\overline{\dot{Q}_j} T_p \Delta u_j}{e^p} \qquad (13)
$$

25 where $\overline{\dot{Q}_j}$ is the mean production rate of the process
26 over the pulse duration $T_p$. One of the drawbacks
27 of this method is that it requires to solve for Eqs. 1
28 and 2 over time, as depicted in Fig. 2.a, which usually
29 requires a few minutes depending on the mechanism
30 size. Consequently, the model parameters cannot be
31 computed at run time in a multi-species Navier-Stokes
32 solver and require the knowledge of all the thermody-
33 namic states a priori.

34 In this work, an alternative and faster method de-
35 picted in Fig. 2.b is presented to compute the α<sub>j</sub> pa-
36 rameters. Rather than solving time-dependent species
37 and energy equations with a zero-dimensional solver,
38 the new model relies exclusively on reaction rate coef-
39 ficients without time integration. Consequently, a sig-
40 nificant reduction of the computational cost associated
41 to the model is obtained as detailed in Section 5. The
42 methodology is based on the idea that the discharge
43 energy branching (i.e., the electron energy transferred

to the surrounding molecules) can be separated in five
44 major channels, namely the chemical $e_{chem}^p$, the heat-
45 ing $e_{heat}^p$, the vibrational $e_{vib}^p$, the electronic excitation
46 $e_{elec}^p$ and the ionization $e_{ion}^p$. The fractional powers dis-
47 sipated by electrons can be computed using an EEDF
48 solver knowing the initial thermodynamic state (Y<sub>k</sub>, T)
49 and the reduced electric field E/N [42].
50

The contributions to $e_{chem}^p$, $e_{heat}^p$, and $e_{vib}^p$ from electron impact reaction are directly linked to the model parameters. However, it is necessary to evaluate the conversion pathways of $e_{elec}^p$ (mainly $N_2^+$ for air-dominated mixtures) and $e_{ion}^p$ into the three considered channels knowing the reactions and their associating rates of progress that mainly depend on (Y<sub>k</sub>, T), in order to avoid including excited and ionized species in the combustion chemical mechanism. The proposed method, PACMIND-SA, relies on two main assumptions: 1) the thermodynamic state of the gas does not vary during a pulse and 2) the radical concentrations remain negligible compared to the main reactants and products (i.e. NH<sub>3</sub>/H<sub>2</sub>/O<sub>2</sub>/N<sub>2</sub>/H<sub>2</sub>O). Thus, quenching of excited and ionized species occurs only with these dominant species in addition to radiative quenching. These assumptions may no longer hold at high pulse energy densities, as will be discussed in Sections 4.1 and 4.2.

As shown in Sec. 3.2, the global processes are de-
71 rived from reactions between neutral and electron or
72 excited species. At any instant, the molar production
73 rate ω<sub>M*</sub> of an excited species (or an ion) denoted M*
74 reads
75

$$
\dot{\omega}_{M^*} = \sum_j (v''_{M^*j} - v'_{M^*j}) \dot{Q}_j \quad (14)
$$

$$
\equiv \dot{\omega}_{p,M^*} - \dot{\omega}_{c,M^*} \quad (15)
$$

where $v'_{M^*j}$ and $v''_{M^*j}$ are respectively the reactant and product stoichiometric coefficients of $M^*$ in the $j^{th}$ reaction, and $\dot{\omega}_{p,M^*}$, $\dot{\omega}_{c,M^*}$ are respectively the production and consumption rates of the species $M^*$. It is assumed that, for $t=0$ and $t \to \infty$, the gas is only constituted of neutral species in their ground state. Consequently, the overall production $\omega_{p,M^*}$ and consumption $\omega_{c,M^*}$ of $M^*$ are equal:

$$
\omega_{p,M^*} \equiv \int_t \sum_j v''_{M^*j} \dot{Q}_j dt = \int_t \sum_j v'_{M^*j} \dot{Q}_j dt \equiv \omega_{c,M^*} \quad (16)
$$

The species M* is assumed to be consumed through
a set of reactions $j \in C_{M^*}$ defined by

$$
M^* + A_j \rightarrow \text{Products}
$$

whose progress rate reads $\dot{Q}_j = k_j[M^*][A_j]$. Thus,
the fraction of M* consumed through the reaction j,
84
85
<table><caption>Table 1: Description of the global processes considered in the PACMIND-SA phenomenological model for ammonia-hydrogen-air flames. Δu<sub>j</sub> corresponds to the internal energy change of the process.</caption><thead><tr><th></th><th>Params</th><th>Global process</th><th>Δu<sub>j</sub> [kJ/mol]</th></tr></thead><tbody><tr><td>P1</td><td>α<sub>1</sub></td><td>O<sub>2</sub> ——> 2 O</td><td>498.3</td></tr><tr><td>P2</td><td>α<sub>2</sub></td><td>N<sub>2</sub> ——> 2 N</td><td>945.4</td></tr><tr><td>P3</td><td>α<sub>3</sub></td><td>O<sub>2</sub> + N<sub>2</sub> ——> 2 NO</td><td>182.5</td></tr><tr><td>P4</td><td>α<sub>4</sub></td><td>O<sub>2</sub> + 2 N<sub>2</sub> ——> 2 N<sub>2</sub>O</td><td>163.2</td></tr><tr><td>P5</td><td>α<sub>5</sub></td><td>NH<sub>3</sub> ——> NH<sub>2</sub> + H</td><td>455.9</td></tr><tr><td>P6</td><td>α<sub>6</sub></td><td>NH<sub>3</sub> ——> NH + H<sub>2</sub></td><td>402.8</td></tr><tr><td>P7</td><td>α<sub>7</sub></td><td>NH<sub>3</sub> + 0.5 N<sub>2</sub> ——> NH + NH<sub>2</sub></td><td>594.9</td></tr><tr><td>P8</td><td>α<sub>8</sub></td><td>NH<sub>3</sub> + 0.5 O<sub>2</sub> ——> NH<sub>2</sub> + OH</td><td>277.3</td></tr><tr><td>P9</td><td>α<sub>9</sub></td><td>H<sub>2</sub> ——> 2 H</td><td>436.0</td></tr><tr><td>P10</td><td>α<sub>10</sub></td><td>H<sub>2</sub> + 0.5 O<sub>2</sub> ——> OH + H</td><td>255.5</td></tr><tr><td>P11</td><td>α<sub>11</sub></td><td>H<sub>2</sub>O ——> OH + H</td><td>499.2</td></tr><tr><td>P12</td><td>α<sub>12</sub></td><td>H<sub>2</sub>O + 0.5 O<sub>2</sub> ——> 2 OH</td><td>360.8</td></tr></tbody></table>

1 denoted β<sub>M<sup>*</sup>j</sub>, can be computed as<fnref n="1" />:

$$ \beta_{M^*j} = \frac{v'_{M^*j} \dot{Q}_j}{\sum_{i \in C_{M^*}} v'_{M^*i} \dot{Q}_i} = \frac{v'_{M^*j} k_j x_{A_j}}{\sum_{i \in C_{M^*}} v'_{M^*i} k_i x_{A_i}} \quad (17) $$

2 where x<sub>A<sub>j</sub></sub> is the molar fraction of the species A<sub>j</sub>. With

3 the assumption that the mixture composition and the

4 temperature do not vary much, β<sub>M<sup>*</sup>j</sub> is constant over

5 time. Since $\dot{Q}_j = (\beta_{M^*j}/v'_{M^*j}) \dot{\omega}_{c,M^*}$, it yields:

$$ Q_j = \int_t \dot{Q}_j dt = \frac{\beta_{M^*j}}{v'_{M^*j}} \omega_{c,M^*} = \frac{\beta_{M^*j}}{v'_{M^*j}} \omega_{p,M^*} \quad (18) $$

6 which can be used to compute the overall production

7 Q<sub>j</sub> of the reaction j as a fraction of the overall production

8 of M<sup>*</sup>.

Initially, all the excited species (or ions) M<sup>*</sup> are produced via electron impact reactions. Then, their quenching can produce other excited species with lower energies, higher level re-population pathways being neglected. For ions, Penning and associative ionization are also neglected. If M<sup>*</sup> is only produced by electron impact reaction, its overall production would read

$$ \omega_{p,M^*} = \int_{t_1}^{t_2} \sum_{j} v''_{M^*j} k_j(\sigma_j)[e^{-}][X_{A_j}] dt \quad (19) $$

$$ = I_e N \sum_{j} v''_{M^*j} k_j(\sigma_j) x_{A_j} \quad (20) $$

where k<sub>j</sub>(σ<sub>j</sub>) is the reaction rate associated to the electron impact cross section σ<sub>j</sub>, N the particle number density and I<sub>e</sub> ≡ ∫<sub>t<sub>1</sub></sub><sub>t<sub>2</sub></sub> [e<sup>-</sup>] dt. In general, in addition to the electron impact reactions, the quenching of other excited species A<sub>k</sub><sup>*</sup> (not necessarily the same ground state species as M<sup>*</sup>) can contribute to the production of M<sup>*</sup>. This contribution to the production of M<sup>*</sup> corresponds to a fraction of the production of the A<sub>k</sub><sup>*</sup>

<figure><img src="image_4.png" /><figcaption>Fig. 2: Principle used to compute the energy branching following a) the original methodology PACMIND [29] based on a posteriori analysis of detailed kinetic simulation and b) the new methodology PACMIND-SA based on the electron power density from the EEDF solver and reaction rates analysis.</figcaption></figure>

<footnote n="1">¹First order reaction such as optical transition are discussed in Supplementary Materials</footnote>
whose only higher excited energy levels are consid- ered:

$$
\omega_{p,M^*} \equiv I_e N \sum_j v''_{M^*_j} k_j(\sigma_j) x_j \\
+ \sum_{A_k^* \in \mathcal{E}_{M^*}} \sum_{j \in R_{A_k^*} \to M^*} \frac{\beta_{A_k^*,j}}{v'_{A_k^* j}} \omega_{p,A_k^*} \quad (21)
$$

where $\mathcal{E}_{M^*} \equiv \{A_k^* | E_{A_k^*} > E_{M^*}\}$ is the set of species
$A_k^*$ which have a higher excitation energy than $M^*$ and
$R_{A_k^*} \to M^*$ is the set of reactions producing $M^*$ from $A_k^*$.

As the production of the higher excited species
is proportional to $I_e N$, the same proportionality re-
cursively applies to all lower excited levels. Thus,
by defining the reduced overall production $\omega'_{p,M^*} \equiv$
$\omega_{p,M^*}/I_e N$:

$$
\omega'_{p,M^*} = \sum_j v''_{M^*_j} k_j(\sigma_i) x_j \\
+ \sum_{A_k^* \in \mathcal{E}_{M^*}} \sum_{j \in R_{A_k^*} \to M^*} \frac{\beta_{A_k^*,j}}{v'_{A_k^* j}} \omega'_{p,A_k^*} \quad (22)
$$

Following a decaying excitation energy level order, it
yields by induction:

$$
Q_j = \frac{\beta_j}{v'_{M_j^*}} \omega'_{p,M_j^*} I_e N \quad (23)
$$

The global process parametrized by $\alpha_p$ introduced
in Section 3.2 actually corresponds to a set of reac-
tions between electrons (subscript $e$) and ground state
species initially present in the mixture (subscript $G$)
($R_{e,G}$) or an excited species (or ions) (subscript $M^*$)
with a ground state species $G$ ($R_{M^*,G}$) leading to the
global process $p \in [1, P]$:

$$
Q_{\alpha_p} = \sum_{i \in R_{e,G}} w_{pi} Q_i + \sum_{i \in R_{M^*,G}} w_{pi} Q_i \quad (24)
$$

where $w_{pi}$ is the weight associated to a reaction $i \in$
$[1, N_r]$ from the detailed mechanism relative to the
global process $p \in [1, P]$ defined by:

$$
w_{pi} = \min_k \left\{ \frac{v''_{ki}}{v''_{kp}} \mid \forall k \in N_s, v''_{kp} > 0 \right\} \quad (25)
$$

where $N_s$ and $N_r$ are respectively the number of
species in the mixture and the number of reactions in
the detailed mechanisms.

As an example for the first global processes
$O_2 \rightarrow 2O$ (i.e., $v''_{Op} = 2$), the reactions

$$
e^{-} + O_{2} \rightarrow e^{-} + 2O \qquad (26)
$$

$$
e^{-} + O_{2} \rightarrow e^{-} + O + O(^{1}D) \quad (27)
$$

have a weight of 1 and 0.5, respectively, the excited
atomic oxygen O(1D) being treated separately from
O as it does not necessary relax back to O without

any chemical effect. Additionally if a reaction *i* perfectly fits with a global process *p* (i.e, $w_{pi} = 1$), then the weights of this reaction for all the other partially matching global processes are set to zero following the modified weight $w'_{pi}$

$$
\begin{align*}
w'_{pi} &= w_{pi} && \text{if } m < 1 \\
&= w_{pi} \delta_{w_{pi},m} && \text{else}
\end{align*}
$$

where $m = \max_i (w_{pi})$ and $\delta_{w_{pi},m}$ is the Kronecker symbol. This solves for instance reactions producing $\text{NH}_2 + \text{OH}$ that could otherwise be incorrectly considered for the process P11.

The first sum in Eq. 24 can be computed in the
same manner as in Eq. 20 while Eq. 23 is used for
the second sum, resulting in:

$$
\begin{align}
Q_{\alpha_p} = NI_e \sum_{i \in R_{e,G}} w_{pi} v_i^{\prime\prime} k_i(\sigma_i) x_i & \nonumber \\
& + NI_e \sum_{i \in R_{M^*,G}} w_{pi} \frac{\beta_{M^*,i}}{v^{\prime}_{M^*,i}} \omega'_{p,M^*} \tag{28} \\
& \equiv I_e N Q'_{\alpha_p} \nonumber
\end{align}
$$

where $Q'_{\alpha_p}$ is the reduced overall production of the global process $j$. The model parameters computation is detailed Section 3.4.

3.4. Energy balance and model closure

In general, the solution of the EEDF solver for low
ionization level verifies a power density balance [39]:

$$
\dot{e}^p = \dot{e}_{\text{elas}}^p + \dot{e}_{\text{inel}}^p + \dot{e}_{\text{growth}}^p \quad (29)
$$

where $\dot{e}^p$ is the total discharge power, $\dot{e}_{\text{elas}}^p$ the elastic power loss, $\dot{e}_{\text{inel}}^p$ the inelastic power loss and $\dot{e}_{\text{growth}}^p$ the growth power loss. For the reduced electric fields relevant to NRP discharge applications (100-500 Td), most of the electron energy is transferred to vibrational and electronic state of the gaseous molecules. In that case the discharge energy can be approximated by the inelastic energy loss that reads:

$$
\begin{align}
e^p &= \int \sum_i k_i(\sigma_i) [e^{-}][X_i] \Delta u_i dt \nonumber \\
&= NI_e \sum_i k_i(\sigma_i) x_i \Delta u_i \tag{30} \\
&\equiv NI_e e^{p'} \nonumber
\end{align}
$$

where $\Delta u_i$ is the threshold energy of the $i^{th}$ process and $e^{p'}$ is the normalized discharge energy.

The energy fractions $\alpha_p$, $\alpha_{\text{heat}}$ and $\alpha_{\text{vib}}$ can now be computed by comparing their energy with the overall discharge energy. The normalized vibrational energy $e_{\text{vib}}$ is simply retrieved from electron impact reaction and excited species quenching following:
$$e'_{\text{vib}} = \sum_{i \in \mathcal{V}} k_i(\sigma_i)x_i\Delta u_i + \sum_{j \in R_{M^* \to \Lambda(v)}} \epsilon_{res,j}^{\text{vib}} \frac{\beta_{M^*,j}}{\nu'_{M^*j}} \omega'_{p,M^*} \quad (31)$$

where $\mathcal{V}$ is the set of electron impact vibrational excitation yielding $\alpha_{\text{vib}} = e'_{\text{vib}}/e^{\prime\prime}$ and $R_{M^* \to \Lambda(v)}$ is the set of reaction leading to a vibrational energy residual in the products (e.g., $N_2(C_3) + O_2 \longrightarrow N_2(v) + 2O$) as considered in the detailed mechanism of [24] with $\epsilon_{res,j}^{\text{vib}}$ the vibrationnal energy residual. For the global chemical processing, the coefficients can be computed as:

$$\alpha_j = \frac{Q'_j \Delta u_j}{e^{\prime\prime}} \quad (32)$$

The computation of $Q'_j$ (see Sec. 3.3), $e'_p$ and thus $\alpha_j$ can be done only with the knowledge of the mixture composition and thermodynamic state and a chemical mechanism with all the rate coefficients. It does not require simulation through time, leading to an improved efficiency in the computation of the model coefficients. Finally, knowing the part of the discharge power stored in chemical energy and vibrational modes, it is assumed that the remaining energy entirely results in fast gas heating to close the model:

$$\alpha_{\text{heat}} = 1 - \alpha_{\text{vib}} - \sum_{j \in \mathcal{J}, P} \alpha_j \quad (33)$$

$$= 1 - \alpha_{\text{vib}} - \alpha_{\text{chem}} \quad (34)$$

Once the $\alpha(E/N, P, T, X_k)$ coefficients are known, they can be introduced in a CFD solver to compute the plasma source terms in the compressible reactive Navier-Stokes equations following the description given in [29].

## 3.5. Model parameters

### 3.5.1. Influence of $E/N$

For a given mixture, the most influential discharge parameter affecting the energy branching is the reduced electric field $E/N$. The variations in the model's coefficients $\alpha$ are shown in Fig. 3 for a stoichiometric ammonia mixture at 300 K and atmospheric pressure. As vibrational excitation energy is primarily stored in $N_2(v)$, its contribution sharply decreases with increasing $E/N$, consistent with observations in air and hydrocarbon-air mixtures. The fast gas heating contribution increases with $E/N$, rising from approximately 20% at 100 Td to 42% at 500 Td. Ammonia vibrational excitation, which is treated as fast gas heating due to its rapid relaxation under the present conditions, accounts for a significant portion of this effect at low reduced electric fields. Regarding chemical effects, a significant increase is observed in the range $E/N \in [100, 200]$ Td, followed by a plateau at around 57% of the discharge

energy. This fast chemical effect is primarily driven by global processes such as ammonia dissociation ($\alpha_5$ and $\alpha_6$) and oxygen dissociation ($\alpha_1$). The dedicated NO and $N_2$O processes, represents an energy fraction around $10^{-4}$ and $10^{-6}$ respectively. Unlike the previous PACMIND method, the updated PACMIND-SA approach can decompose the global process fraction into distinct physical reactions, including electron impact, quenching, and ion recombination (see Fig. 4). This could help to highlight the most important reactions within the mechanism that require careful definition of their reaction rates. Other chemical pathways exhibit a general increase with $E/N$, but each contributes only to a small fraction of the total discharge energy.

The case of a stoichiometric hydrogen-air mixture is also considered and presented in Fig. 5. For this mixture, the fast gas heating fraction increases with the reduced electric field, reaching approximately 50% of the discharge energy at 250 Td. The chemical energy fraction peaks at 45% around 250 Td and then decreases to 40% at 500 Td. Among the chemical pathways, oxygen dissociation represents the most significant contribution to chemical energy transfer. It is followed by hydrogen dissociation, which increases between 100 and 220 Td, then saturates at higher reduced electric fields reaching about 15% of the discharge energy. As previously with ammonia, the NOx processes represent a low energy fraction, around $10^{-3}$ and $10^{-4}$ for $\alpha_3$ and $\alpha_4$ respectively.

The model parameters variation with $E/N$ for various $\text{NH}_3/\text{H}_2$ blends is provided in the Supplementary Materials for conciseness. Notice that these coefficients rely on the specific kinetic mechanism considered in this work. Comparisons have been made with the previous methodology to compute the $\alpha$ coefficients with an overall good agreement.

<figure><img src="image_5.png" /><figcaption>Fig. 3: Model parameters against the reduced electric field for a stoichiometric $\text{NH}_3$-air mixture: ● $\alpha_{\text{heat}}$, ★ $\alpha_{\text{vib}}$, ■ $\alpha_{\text{chem}}$, ▼ $\alpha_1$, ◁ $\alpha_2$, ▲ $\alpha_5$, ▲ $\alpha_6$, ♦ $\alpha_7$, ● $\alpha_8$</figcaption></figure>

### 3.5.2. Burnt gases

The model parameters have been computed for a mixture of 0.69 $N_2$ / 0.30 $\text{H}_2\text{O}$ / 0.01 $\text{O}_2$ (in volume),
<figure><img src="image_6.png" /><figcaption>Fig. 4: Contribution of chemical process in the global ammonia dissociation process ● α₅ for a stoichiometric ammonia-air mixture: ▲ Electron impact, ▲ Quenching of N₂*, ♦ Via NH₃⁺/NH₄⁺</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5: Model parameters against the reduced electric field for a stoichiometric H₂-air mixture: ● αheat, ★ αvib, ■ αchem, ▼ α₁, ◁ α₂, + α₉, × α₁₀</figcaption></figure>

representative of NH₃/H₂ burnt gases under slightly lean conditions, and are shown in Fig. 6. The chemical energy behave like in the H₂ case peaking around 200 Td, whereas the fast gas heating contribution continuously increases within [100-500] Td. It is noteworthy that, although O₂ accounts for only 1% of the molar composition, it contributes to approximately 20% of the discharge energy for E/N > 200 Td. In practice, this cannot hold for high energy discharge as the possible oxygen dissociation energy will be lower than what the model predicts. This is primarily due to the fact that quenching of electronically excited N₂ states (N₂(B), N₂(C), and N₂(a)) by H₂O is not included in the kinetic mechanism from [19]. As a result, an important part of the electronic excitation energy is consumed via highly efficient collisions with O₂. In other words, the branching ratio βN₂,j in Eq. 17 is significantly higher than 1% for the reactions j ≡ N₂(B, C, a) + O₂ → N₂ + 2 O (12% for N₂(B), 28% for N₂(C) and 68% for N₂(a)). This observation may change with updated kinetic mechanisms that account for reactive collisional quenching of electroni-

<figure><img src="image_8.png" /><figcaption>Fig. 6: Model parameters against the reduced electric field for a mixture of 0.69 N₂ / 0.30 H₂O / 0.01 O₂: ● αheat, ★ αvib, ■ αchem, ▼ α₁, ◁ α₂, ◁ α₃, ● α₁₁, × α₁₂</figcaption></figure>

cally excited N₂ species with H₂O. Alternatively, the model coefficients could be regularly updated, at least if the oxygen molecule concentration falls to zero.

# 4. Validation and discussion

## 4.1. Discharges in air

The new methodology to compute the PACMIND-SA parameters is first assessed in a single pulse occurring in air by comparison with the detailed kinetic simulation and the previous PACMIND methodology [29]. The reduced electric field is set to 200 Td and the pulse energy has been set to 10⁴ and 10⁶ J · m⁻³ by adjusting the pulse duration in the detailed simulation. The discharge energy is deposited in τpulse = 50 ns in the phenomenological approach. The model parameters are shown in Tab. 2, evaluating the impact of the discharge energy in the PACMIND methodology whereas a single set of parameters is defined for the PACMIND-SA methodology as it assumes a constant thermodynamic states and does not depend on the discharge energy density. For the low discharge energy density, both methodologies lead to very similar discharge energy branching. The resulting radicals and gas temperature behaviors are shown in Fig. 7a. The species production is well reproduced by both modeling approaches, with a very good match of decay rate during the afterglow. Differences observed in the fast gas heating is primarily due to the model closure which consider all the missing processes as gas heating in PACMIND-SA whereas in PACMIND 10% of the energy is still stored in electronic modes at 100 ns when the parameters are evaluated. For the high energy density, slight parameter differences are observed between PACMIND and PACMIND-SA that can be observed in Fig. 7b. For the PACMIND simulation, which account for the discharge energy density, a very good match is observed for all the radicals and the gas temperature. For PACMIND-SA, atomic oxygen production is slightly higher than in
the detailed simulation. The resulting delay of O decay results in small differences in the gas heating around $10^{-4}$ s due to an earlier relaxation of the vibrational energy. The differences in NO and N productions are mainly due to the ineffectiveness of the reaction $N_2^* + O \rightarrow NO + N$ in the PACMIND-SA approach when considering an initial mixture without atomic oxygen (i.e., $\dot{Q}_j = 0$).

In general, the new procedure PACMIND-SA offers a good accuracy compared to detailed simulation even at high energy density with a very low computational cost for the model parameter evaluation compared to PACMIND. In Section 4.2, the PACMIND-SA model accuracy will be verified in ignition cases of various NH$_3$/H$_2$/N$_2$/O$_2$ mixtures.

<table><caption>Table 2: Model parameters in air using the previous method (PACMIND) and the new one (PACMIND-SA).</caption><thead><tr><th>Parameter</th><th>PACMIND</th><th>PACMIND-SA</th></tr></thead><tbody><tr><td>$e^p$ [J·m<sup>-3</sup>]</td><td>1 × 10<sup>4</sup></td><td>1 × 10<sup>6</sup></td><td>-</td></tr><tr><td>α<sub>heat</sub></td><td>0.39</td><td>0.43</td><td>0.38</td></tr><tr><td>α<sub>vib</sub></td><td>0.16</td><td>0.16</td><td>0.16</td></tr><tr><td>α<sub>1</sub></td><td>0.41</td><td>0.35</td><td>0.43</td></tr><tr><td>α<sub>2</sub></td><td>0.036</td><td>0.044</td><td>0.024</td></tr><tr><td>α<sub>3</sub></td><td>0.006</td><td>0.018</td><td>0.005</td></tr></tbody></table>

## 4.2. Single-pulse ignition

In this section, different ignition cases are considered as detailed in Tab. 3. Pure ammonia, pure hydrogen and blends are considered by defining the fuel blend composition

$$ \gamma H_2 + (1 - \gamma) NH_3 $$

where $\gamma$ is the molar blending factor. In all cases, a high energy density of $e^p = 10^6$ J·m$^{-3}$ is used to ignite the mixture in a single pulse.

<table><caption>Table 3: NH$_3$-H$_2$ molar blending factor $\gamma$ and equivalence ratio $\phi$ for Cases A, B, and C sharing common $e^p = 10^6$ J·m$^{-3}$, $T_0 = 300$ K, $P_0 = 1$ atm and $E/N = 250$ Td.</caption><thead><tr><th></th><th>Case A</th><th>Case B</th><th>Case C</th></tr></thead><tbody><tr><td>γ</td><td>0</td><td>1</td><td>0.3</td></tr><tr><td>φ</td><td>0.8</td><td>0.5</td><td>0.8</td></tr></tbody></table>

The PACMIND-SA model is compared to the detailed simulations of the three cases in Fig. 8. In all cases, the phenomenological approach is able to well reproduce the ignition phase with an ignition delay error of -5 %, -6 % and -2 % respectively for the A, B and C cases. The time evolution of the main radicals of the case C is shown in Fig. 9. After the transient period [0-100] ns corresponding to the energy

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

$$ (a) e^p = 10^4 \text{ J} \cdot \text{m}^{-3} $$

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

$$ (b) e^p = 10^6 \text{ J} \cdot \text{m}^{-3} $$

<caption>Fig. 7: Time evolution of the main radicals and the gas temperature for a single discharge in air with a) $e^p = 10^4$ J·m$^{-3}$ and b) $e^p = 10^6$ J·m$^{-3}$ for (---) the detailed kinetics, (- - -) PACMIND and (--- -) PACMIND-SA.</caption>

deposition, a very good agreement is observed for the consumption and subsequent production of the radicals. Similar conclusion can be drawn for the Cases A and B (see Supplementary Materials). Interestingly, the error introduced at high energy densities by the constant thermochemical state assumption in the PACMIND-SA model has a limited impact on the ignition process.

In order to validate the model on a broad range of conditions, the ignition delay for both PACMIND-SA and the detailed kinetics is computed for a range of reduced electric field $E/N \in [100, 500]$ Td and blending ratio $\gamma \in [0, 1]$ and a discharge energy of $10^6$ J·m$^{-3}$. The relative error between the detailed and PACMIND-SA ignition delay time, respectively $\tau_{\text{ig}}^{\text{detailed}}$ and $\tau_{\text{ig}}^{\text{PACMIND-SA}}$, is defined as:

$$ \text{Err} = \frac{\tau_{\text{ig}}^{\text{PACMIND-SA}} - \tau_{\text{ig}}^{\text{detailed}}}{\tau_{\text{ig}}^{\text{detailed}}} \quad (35) $$

The resulting errors are shown in Fig. 10. The relative error globally increases with the reduced electric
<figure><img src="image_11.png" /><figcaption>Fig. 8: Time evolution of the gas temperature during the ignition of the cases A, B and C for (—) the detailed kinetics and (---) PACMIND-SA.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 9: Time evolution of the main radicals during the ignition of the Case C for (—) the detailed kinetics and (---) PACMIND-SA.</figcaption></figure>

field for ammonia-containing mixtures. This could be due greater importance of ionization at high reduced electric fields not considered in detail with PACMIND-SA, especially with complex $NH_x^+$ processes. Overall, the ignition delay error is bounded by -9.4% and +6.4% which shows the good capability of PACMIND-SA to properly reproduce the detailed kinetics in an ignition scenario on a broad range of blending ratios and reduced electric fields.

### 4.3. NOₓ and N₂O productions

The model is also validated with respect to its ability to properly capture the right amount of NOₓ and N₂O emissions. Two scenarios are considered corresponding to the cases A (pure NH₃) and B (pure H₂) with only the pulse energy changed from $e^p = 10^6$ J·m⁻³ to $e^p = 10^5$ J·m⁻³ to avoid ignition before 1 s. For each scenario, a detailed simulation is performed as well as two simulations with the PACMIND-SA model: i) the full 12 processes model and ii) a PACMIND-SA-woNOx model performed with $α_3$ and $α_4$ set to 0, the processes 3 and

<figure><img src="image_13.png" /><figcaption>Fig. 10: Relative error (%) on a single pulse ignition delay between PACMIND-SA and detailed kinetics for a range of reduced electric field and blending ratios. The equivalence ratio is set to 0.8</figcaption></figure>

4 being dedicated to the NO and N₂O. The corresponding energy fractions, although small, are considered as fast gas heating (i.e., $α_{heat} = α_3 + α_4$) to have the same energy density in both models. The additional gas temperature increase in the PACMIND-SA-woNOx model is significantly lower than 1 K and does not increase the thermal NOₓ pathways. The results are presented in Figures 11a and 11b.

For the ammonia case (see Figure 11a), a very good match is obtained for NO, N₂O and NO₂. It is interesting to note that in the case of pure NH₃, the deactivation of the dedicated NOₓ processes does not influence the time evolution of NOₓ as the PACMIND-woNOx model correctly reproduces the trends. Indeed in this case, a predominant part of NO is produced during the oxidation processes via NH that is produced by the discharge. This NO acts as a precursor in the formation pathways of NO₂ and N₂O. Thus to properly capture the NOₓ emissions in plasma-assisted ammonia oxidation process, the correct prediction of NH radicals resulting from NH₃ plasma-dissociation is sufficient. However, in the case of H₂ (see Figure 11b), the plasma-NOₓ global processes are fundamental to retrieve the right amount of NOₓ. This is an indication that for H₂, the plasma process and not the resulting heavy molecule reactions are responsible for production of NOₓ. This shows that dedicated NOₓ processes are necessary to capture in all scenarios the right amount of NOₓ.

### 4.4. Multi-pulse ignition

Finally, the model is validated for multiple pulses. The same three cases as in Tab. 3 are considered except that the energy is deposited equally over 100 pulses (i.e., $e^p = 10^4$ J·m⁻³) at a frequency of 20 kHz. For the simulations using the PACMIND-SA model, the coefficients were updated on the fly at the beginning of each pulse, taking advantage of
<figure><img src="image_14.png" /><figcaption>(a) NH₃-air</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) H₂-air</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 12: Time evolution of the gas temperature during the ignition of the cases A, B and C via 100 discharges at 20 kHz with (—) the detailed kinetics and (---) PACMIND-SA.</figcaption></figure>

detailed chemical mechanisms, the approach has been applied to CH₄-air ignition using the mechanism presented in [24]. The ignition of a lean CH₄-air mixture at T₀ = 300 K, P₀ = 1 atm, φ = 0.67, and subjected to a single pulse of $e^p = 10^6 \text{ J} \cdot \text{m}^{-3}$, was simulated using both the detailed mechanism and the PACMIND-SA model, considering the same global processes as in [29] for methane-air mixtures. The resulting ignition delay times and relative errors are shown in Fig.13. Overall, excellent agreement is obtained between the two simulations, with relative errors below 10% across a wide range of reduced electric fields. These results further reinforce the capability of PACMIND-SA to handle a wide variety of fuels and non-equilibrium plasma kinetic mechanisms.

the low computational cost for the model parameters. The time evolution of the temperature is presented in Fig. 12 and shows a remarkable agreement between PACMIND-SA and the detailed kinetics with an ignition delay error of -1%, -0.35% and ≤ 0.1 % for the cases A, B and C, respectively. Dynamics of the updated model coefficients can be found in the Supplementary Materials.

1

2

3

4

5

6

7

8

9

## 4.5. Other applications

10

11

12

13

14

15

16

17

18

19

The model is eventually tested in a series of limit-case scenarios, taking advantage of the ability of PACMIND-SA to work for any NH₃/H₂/N₂/O₂/H₂O mixture. These scenarios correspond to a pulse in the context of MILD combustion, oxycombustion, plasma assisted ammonia reforming as well as combustion with water injection. The results presented in the Supplementary Materials shows an overall very good agreement on the temperature evolution as well as radicals production.

20

## 4.6. Methane-air ignition

21

22

To further validate the methodology and demonstrate its flexibility with respect to both fuel type and

<figure><img src="image_17.png" /><figcaption>Fig. 13: Ignition delay times for a detailed simulation and the PACMIND-SA model of a methane-air mixture at φ = 0.67, $e^p = 10^6 \text{ J} \cdot \text{m}^{-3}$ with various E/N.</figcaption></figure>

# 5. Computational cost and remaining challenges

38

39

40

41

42

43

44

The current version of the PACMIND-SA model is able to reduce drastically the computation cost of the model parameters compared to the original PACMIND model. Indeed, it takes around 0.05 s to PACMIND-SA to compute the model parameters instead of ≥ 100 sec for PACMIND. In more detail,
the EEDF solver takes approximately 99.5% of the overall PACMIND-SA cost. Thus, optimization of the EEDF solver (i.e., currently a python version of bolos) could be a way to boost the performance of the model. With the aim of integrating the model in a CFD solver, the parameters would have to be computed at the cell level, which could be an expensive and challenging task. To control this cost, EEDF computations could be constrained by changes in the main parameters of the solver such as the target species mass fractions or the reduced electric field based on what is done in the ZDPlasKin solver or more recently in a streamer code [45]. Moreover, as discharges are only applied in a very limited region of the domain, Dynamic Load Balancing strategy could be consider to balance the cost of the EEDF solver over all the processors as suggested in [45]. In-Situ Adaptive Tabulation (ISAT) [46] could also be employed to store the coefficients, incrementally populating the table as the discharge region explores the entire physical space.

Perspectives in the model development have also been identified. First, the global processes identified a priori based on a path-flux analysis could be generalized following approaches such as presented in [47]. Second, the radicals produced by the discharge could be introduced as main species, thus helping to consider them as quenchers that could improve the model at high energy density.

## 6. Conclusion

This work introduced PACMIND-SA, a semi-analytical model designed to include the primary thermo-chemical effects of non-equilibrium plasma discharges in ammonia–hydrogen–air mixtures. Building upon previous phenomenological approaches of PACMIND [29], the new methodology named PACMIND-SA eliminates the need for time-consuming detailed simulations, by precomputing the discharge energy distribution over the fast gas heating, the fast chemical effect and the slow gas heating using reaction rates derived from detailed a kinetic mechanism and the analysis of the Electron Energy Distribution function. It captures the key plasma–combustion interactions while maintaining computational efficiency, making it suitable for integration into large-scale CFD solvers. While mainly applied to NH₃–H₂ blends, the model is intended to be generic, and the methodology behind it can be applied to a wide variety of mixture as long as a detailed kinetic mechanism exist for it, as shown by its application to methane–air ignition. Variation of the N₂/O₂ ratio, already supported within the model, could be relevant for applications involving oxygen-enriched combustion. In addition, dilution with other species such as Ar or He could be considered within the model.

Validation against detailed kinetic simulations across a range of conditions, including single-pulse discharge in air, ignition via a single or multiple discharges in different NH₃–H₂ blends, and NOx emissions

predictions, demonstrated the accuracy and robustness of the proposed approach. Despite relying on simplifications such as constant thermochemical states during the pulse, PACMIND-SA yields results in excellent agreement with detailed models, even at high energy densities. The model also captures the distinct plasma-driven NOx formation mechanisms in hydrogen and ammonia mixtures, underlying its flexibility and accuracy.

Finally, PACMIND-SA represents a significant step toward the reliable, scalable, physics-informed, and computationally affordable modeling of plasma-assisted combustion, with direct applicability to the development of cleaner and more efficient energy systems based on low-carbon fuels.

## Acknowledgments

The research leading to these results has received funding from the ANR PRCI grant JETHPAC (ANR-23-CE05-0025).

## Model parameters availability

The model parameters can be computed for a given thermochemical state and reduced electric field via an executable file that will be made available at https://chemistry.cerfacs.fr/en/plasma/ and described in the Supplementary Materials.

## CrediT authorship contribution statement

**Harry Gillingham:** Writing – original draft; Conceptualization; Data curation; Formal analysis; Investigation; Methodology; Software; Validation; Visualization – **Thomas Lesaffre:** Writing – review & editing; Supervision; Validation – **Nicolas Barléon:** Writing – original draft, review & editing; Conceptualization; Investigation; Methodology; Software; Supervision; Validation; Visualization; Funding acquisition.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Supplementary material

Supplementary material associated with this article can be found in the online version.

## References

[1] H. Kobayashi, A. Hayakawa, K. K. A. Somaratne, E. C. Okafor, Science and technology of ammonia combustion, Proc. Combust. Inst. 37 (1) (2019) 109–133.

[2] C. Lhuillier, P. Brequigny, F. Contino, C. Mounaïm-Rousselle, Experimental study on ammonia/hydrogen/air combustion in spark ignition engine conditions, Fuel 269 (2020) 117448.
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

[3] B. Mei, J. Zhang, X. Shi, Z. Xi, Y. Li, Enhancement of ammonia combustion with partial fuel cracking strategy: Laminar flame propagation and kinetic modeling investigation of nh3/h2/n2/air mixtures up to 10 atm, Combust. Flame 231 (2021) 111472.

[4] N. N. Shohdy, M. Alicherif, D. A. Lacoste, Transfer functions of ammonia and partly cracked ammonia swirl flames, Energies 16 (3) (2023) 1323.

[5] Y. Ju, W. Sun, Plasma assisted combustion: Dynamics and chemistry, Prog. Energ. Combust 48 (2015) 21–83.

[6] D. A. Lacoste, Flames with plasmas, Proc. Combust. Inst. 39 (4) (2023) 5405–5428.

[7] S. Starikovskaia, D. A. Lacoste, G. Colonna, Non-equilibrium plasma for ignition and combustion enhancement, Eur. Phys. J. D 75 (2021) 231.

[8] J.-B. Perrin-Terrin, N. Vaysse, D. Durox, R. Vicquelin, S. Candel, C. O. Laux, A. Renaud, Plasma-assisted combustion of hydrogen swirling flames: Extension of lean blowout limit and nox emissions, Proc. Combust. Inst. 40 (1–4) (2024) 105546.

[9] S. A. Shcherbanev, Q. Malé, B. Dharmaputra, R. Solana-Pérez, N. Noiray, Effect of plasma-flow coupling on the ignition enhancement with non-equilibrium plasma in a sequential combustor, J. Phys. D Appl. Phys. 55 (2022) 425202.

[10] Q. Malé, S. Shcherbanev, M. Impagnatiello, N. Noiray, Stabilization of a thermoacoustically unstable sequential combustor using non-equilibrium plasma: Large eddy simulation and experiments, Proc. Combust. Inst. 40 (1–4) (2024) 105277.

[11] J. Choe, W. Sun, T. Ombrello, C. Carter, Plasma assisted ammonia combustion: Simultaneous NOx reduction and flame enhancement, Combust. Flame 228 (2021) 430–432.

[12] A. Shioyoke, J. Hayashi, R. Murai, N. Nakatsuka, F. Akamatsu, Numerical investigation on effects of nonequilibrium plasma on laminar burning velocity of ammonia flame, Energy amp; Fuels 32 (3) (2018) 3824–3832.

[13] M. Shahsavari, A. A. Konnov, A. Valera-Medina, M. Jangi, On nanosecond plasma-assisted ammonia combustion: Effects of pulse and mixture properties, Combust. Flame 245 (2022) 112368.

[14] G. Faingold, O. Kalitzky, J. K. Lefkowitz, Plasma reforming for enhanced ammonia-air ignition: A numerical study, Fuel Communications 12 (2022) 100070.

[15] S. Shy, V. Mai, Y. Chen, H. Hsieh, Nanosecond repetitively pulsed discharges and conventional sparks of ammonia-air mixtures in a fan-stirred cruciform burner: Flammability limits and ignition transition, Appl. Energy Combust. Sci. 15 (2023) 100164.

[16] Y. Tang, D. Xie, B. Shi, N. Wang, S. Li, Flammability enhancement of swirling ammonia/air combustion using ac powered gliding arc discharges, Fuel 313 (2022) 122674.

[17] N. N. Shohdy, M. Alicherif, T. F. Guiberti, E.-t. Essebbar, D. A. Lacoste, Quantitative comparison of the benefits of cracking and nanosecond repetitive pulsed discharges on the lean blow-off, emissions, and topology of ammonia premixed swirl flames, Combust. Flame 269 (2024) 113678.

[18] J. Sun, Y. Bao, K. Zhang, A. A. Konnov, M. Richter, E. Kristensson, A. Ehn, A comprehensive study on dynamics of flames in a nanosecond pulsed discharge. part ii: Plasma-assisted ammonia and methane combustion, Combust. Flame 275 (2025) 114076.

[19] H. Zhong, X. Mao, N. Liu, Z. Wang, T. Ombrello, Y. Ju, Understanding non-equilibrium n2o/nox chemistry in plasma-assisted low-temperature nh3 oxidation, Combust. Flame 256 (2023) 112948.

[20] X. Mao, H. Zhong, N. Liu, Z. Wang, Y. Ju, Ignition enhancement and no formation of nh3/air mixtures by non-equilibrium plasma discharge, Combust. Flame 259 (2024) 113140.

[21] G. Rekkas-Ventiris, P. Sabia, G. Sorrentino, A. Bellemans, Nanosecond pulsed plasma-assisted MILD combustion of ammonia, Proc. Combust. Inst. 40 (1–4) (2024) 105384.

[22] N. Liu, Q. Chen, X. Lu, J. Pan, H. Wei, X. Mao, N2o reduction with selective excitation of species by non-equilibrium plasma in an nh3/air mixture, Combust. Flame 276 (2025) 114164.

[23] A. Bellemans, N. Kincaid, N. Deak, P. Pepiot, F. Bisetti, P-DRGEP: a novel methodology for the reduction of kinetics mechanisms for plasma-assisted combustion applications, Proc. Combust. Inst. 38 (2021) 6631–6639.

[24] L. Cheng, N. Barleon, B. Cuenot, O. Vermorel, A. Bourdon, Plasma assisted combustion of methane-air mixtures: Validation and reduction, Combust. Flame 240 (2022) 111990.

[25] N. Barleon, L. Cheng, B. Cuenot, O. Vermorel, A. Bourdon, Investigation of the impact of nrp discharge frequency on the ignition of a lean methane-air mixture using fully coupled plasma-combustion numerical simulations, Proc. Combust. Inst. 39 (4) (2023) 5521–5530.

[26] Z. Shi, X. Mao, Z. Wang, Y. Ju, Plasma assisted nh3/h2/air ignition in nanosecond discharges with non-equilibrium energy transfer, Combust. Flame 275 (2025) 114053.

[27] R. Shetty, L. Massa, Plasma-fluid energy transfers in plasma assisted ignition, Combust. Sci. Technol. (2024) 1–38.

[28] M. Castela, B. Fiorina, A. Coussement, O. Gicquel, N. Darabiha, C. O. Laux, Modelling the impact of non-equilibrium discharges on reactive mixtures for simulations of plasma-assisted ignition in turbulent flows, Combust. Flame 166 (2016) 133–147.

[29] N. Barléon, L. Cheng, B. Cuenot, O. Vermorel, A phenomenological model for plasma-assisted combustion with nrp discharges in methane-air mixtures: Pacmind, Combust. Flame 253 (2023) 112794.

[30] X. Shao, D. A. Lacoste, H. G. Im, A computational study of plasma-assisted ammonia/hydrogen combustion using a phenomenological model, in: AIAA SCITECH 2024 Forum, 2024.

[31] P. N. Johnson, T. S. Taneja, Z. Sun, S. Yang, K. Kim, C.-B. Kweon, Phenomenological modeling of plasma-assisted methane ignition using plasma energy fraction manifolds, in: AIAA SCITECH 2025 Forum, 2025.

[32] P. N. Johnson, T. S. Taneja, S. Yang, 3d phenomenological modeling of plasma-assisted methane reforming, in: AIAA SCITECH 2024 Forum, 2024.

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

[33] R. J. Dijoud, N. Laws, C. Guerra-Garcia, Mapping the performance envelope and energy pathways of plasma-assisted ignition across combustion environments, Combust. Flame 271 (2025) 113793.

[34] P. N. Johnson, T. S. Taneja, S. Yang, Effects of spatiotemporal plasma power distribution on the modeling of ignition kernel evolution in quiescent and turbulent methane/air mixtures, J. Phys. D Appl. Phys 57 (45) (2024) 455201.

[35] M. Shahsavari, N. Chakraborty, S. Zhong, A. Valera-Medina, M. Jangi, Comparative effects of plasma and preheating in assisting premixed ammonia/air flames: A dns study, Fuel 381 (2025) 133645.

[36] Y. Bechane, B. Fiorina, Numerical analysis of turbulent flame enhancement by nanosecond repetitively pulsed plasma discharges, Proceedings of the Combustion Institute 39 (4) (2023) 5465-5476.

[37] N. Barléon, B. Cuenot, O. Vermorel, Large-eddy simulation of swirled flame stabilisation using nrp discharges at atmospheric pressure, Appl. Energy Combust. Sci. 15 (2023) 100163.

[38] Q. Malé, N. Barléon, S. Shcherbanev, B. Dharmaputra, N. Noiray, Numerical study of nitrogen oxides chemistry during plasma assisted combustion in a sequential combustor, Combust. Flame 260 (2024) 113206.

[39] G. J. M. Hagelaar, L. C. Pitchford, Solving the Boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models, Plasma Sources Sci. T. 14 (2005) 722-733.

[40] L. C. Pitchford, L. L. Alves, K. Bartschat, S. F. Bigagi, M. Bordage, I. Bray, C. E. Brion, M. J. Brunger, L. Campbell, A. Chachereau, B. Chaudhury, L. G. Christophorou, E. Carbone, N. A. Dyatko, C. M. Franck, D. V. Fursa, R. K. Gangwar, V. Guerra, P. Haefliger, G. J. M. Hagelaar, A. Hoesl, Y. Itikawa, I. V. Kochetov, R. P. McEachran, W. L. Morgan, A. P. Napartovich, V. Puech, M. Rabie, L. Sharma, R. Srivastava, A. D. Stauffer, J. Tennyson, J. de Urquijo, J. van Dijk, L. A. Viehland, M. C. Zammit, O. Zatsarinny, S. Pancheshnyi, Lxcat: an open-access, web-based platform for data needed for modeling low temperature plasmas, Plasma Process. Polym. 14 (1-2) (2016).

[41] N. Liu, B. Mei, X. Mao, Z. Wang, Z. Sun, Y. Xu, Z. Shi, Y. Ju, Kinetics of low temperature plasma assisted NH3/H2 oxidation in a nanosecond-pulsed discharge, Proc. Combust. Inst. 40 (1-4) (2024) 105353.

[42] A. Starikovskiy, N. Aleksandrov, Plasma-assisted ignition and combustion, Prog. Energ. Combust. 39 (2013) 61-110.

[43] L. S. Thorsen, M. S. Jensen, M. S. Pullich, J. M. Christensen, H. Hashemi, P. Glarborg, V. A. Alekseev, E. J. Nilsson, Z. Wang, B. Mei, N. Liu, Y. Ju, High pressure oxidation of nh3-heptane mixtures, Combust. Flame 254 (2023) 112785.

[44] M. Capitelli, F. C.M., B. Gordiets, A. Osipov, Plasma Kinetics in Atmospheric Gases, Springer, Berlin, 2000.

[45] X. Shao, D. A. Lacoste, H. G. Im, A unified fluid model for nonthermal plasmas and reacting flows, Plasma Sources Sci. T. 34 (8) (2025) 085016.

[46] S. Pope, Computationally efficient implementation of combustion chemistry using in situ adaptive tabulation, Combust. Theory Model 1 (1) (1997) 41-63.

[47] P. N. Johnson, T. S. Taneja, S. Yang, Plasma-based global pathway analysis to understand the chemical kinetics of plasma-assisted combustion and fuel reforming, Combust. Flame 255 (2023) 112927.

63

64

65

66