# Plasma-Assisted Partial Methane Oxidation. Part I: One-Dimensional Statistical Modeling of a Dielectric Barrier Discharge Reactor

Mohamed Saadana, Zeina Al Zayed, Hanene Oueslati, Nicolas Barléon, Vincent Robin, Élodie Fourré, Catherine Batiot-Dupeyrat, Sylvie Rossignol

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

## ► To cite this version:

Mohamed Saadana, Zeina Al Zayed, Hanene Oueslati, Nicolas Barléon, Vincent Robin, et al., Plasma-Assisted Partial Methane Oxidation. Part I: One-Dimensional Statistical Modeling of a Dielectric Barrier Discharge Reactor. Journal of Physical Chemistry C, 2025, 129 (13), pp.6157-6176. 10.1021/acs.jpcc.4c06369. [hal-05016493](https://hal.science/hal-05016493v1)

HAL Id: hal-05016493

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

Submitted on 2 Apr 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_2.png" /></figure>
# Plasma-assisted partial methane oxidation.
## Part I: One-dimensional statistical modeling of a DBD reactor.

Mohamed Saadana,¹ Zeina Al Zayed,² Hanen Oueslati,¹ Nicolas Barléon,³ Vincent Robin,*,¹ Élodie Fourré,² Catherine Batiot-Dupeyrat,² and Sylvie Rossignol⁴

¹Institut Pprime, ISAE-ENSMA, Université de Poitiers, 1 Av. Clément Ader, 86360 Chasseneuil-du-Poitou, France

²IC2MP, Université de Poitiers, 1 rue Marcel Doré, 86073 Poitiers, France

³CERFACS, 42 Avenue Gaspard Coriolis, 31057 Toulouse, France

⁴IRCER, Université de Limoges, 12 Rue Atlantis, 87068 Limoges, France

* E-mail: vincent.robin@isae-ensma.fr

### Abstract

This paper presents an original one-dimensional statistical model designed to complement experimental data from a practical DBD reactor. The experimental setup consists of a plasma-assisted reactor with gas injection composed solely of CH₄ and O₂. The numerical procedure uses electrical measurements to provide a realistic description of the power consumption and the current flowing through the gas, in a complex scenario where approximately one hundred current peaks are measured per electrical half-period. A plasma kinetic model is then used to analyse the chemistry and characteristic times of both isolated discharges and a representative train of discharges occurring within the reactor. These times are used to optimise the computational costs and to
select the most appropriate kinetic schemes for the gas phase, whether in a plasma or quasi-thermodynamic equilibrium state. This approach also allows the separation of fast transformations (plasma) occurring at constant volume, from slower transformations occurring at constant pressure. The statistical approach, based on a Monte Carlo method, clearly identifies the assumptions required to reduce the real complexity of the DBD reactor to a 1D flow model. The combination of chromatographic measurements at the reactor outlet and numerical simulations provides the heterogeneity factor of the discharges, which is identified as a key parameter in the model. Although the flow can be considered stationary on average, the obtained value reveals a highly heterogeneous spatial distribution of the discharges within the reactor. Thus, the numerical results suggest that the gases passing through the reactor are rarely in a plasma state.

Introduction

The synthesis of high-value-added chemical compounds from methane frequently necessitates an intermediate syngas production stage.<sup>1,2</sup> Among the established processes mentioned in the literature,<sup>3</sup> such as steam reforming,<sup>4</sup> dry reforming,<sup>5</sup> and auto-thermal reforming,<sup>6</sup> the partial oxidation of methane<sup>7-9</sup> is particularly noteworthy due to its stability, high conversion efficiency, and low energy consumption. A highly effective approach for methane partial oxidation is the utilisation of non-thermal plasma (NTP), generated through various discharges<sup>10</sup> including microwave,<sup>11</sup> gliding arc, nanosecond repetitive pulse (NRP),<sup>12</sup> and dielectric barrier discharges (DBD).<sup>13,14</sup>

One key advantage of NTP is its ability to transfer most of its energy to electrons, which attain temperatures in the range of 10<sup>4</sup>-10<sup>5</sup> K (1-10 eV), while the heavy species temperature remains near ambient levels. Among the discharge generation systems, DBD stands out for its simple design (easy setup and maintenance) and by the use of dielectric materials that prevent the formation of an electric arc and enable better spatial distribution of the active species (ions, radicals, electrons, and excited species) within the reactor volume.
Despite recent significant progress in the field of plasma,¹⁵,¹⁶ a deeper comprehension of the active species generated in the gas phase can only be achieved by combining experimental and simulation results. The most commonly employed simulation tools for the analysis of these flows are homogeneous reactors¹⁷⁻²⁰ (0D), which provide an accurate description of the chemical kinetics occurring within the plasma and offer insights into the most probable reaction pathways for a given composition and thermodynamic state. The associated numerical results can be employed to reinforce experimental analysis; however, a direct comparison with practical measurements is seldom feasible.

Advances in computing capabilities have facilitated the development of multidimensional simulation tools that consider species and charges transport,²¹⁻²⁴ allowing the analysis of discharge dynamics and their spatial structure. However, these simulations are limited by small computational domains and short simulation times, making validation with experiments challenging. For more information on the numerical tools used in plasma simulations, see the review article by Maitre et al.²⁵

In contrast, there have been considerably fewer attempts²⁶,²⁷ to develop a one-dimensional model capable of simulating a practical DBD reactor as a whole system, which is yet crucial for integrating numerical data with experimental measurements. The development of such a numerical simulation tool poses several challenges. The large number of chemical reactions and species involved in the plasma results in prohibitively long computation times, and the corresponding detailed kinetic schemes are validated for specific conditions, but their use in practical systems with a wide range of operating conditions can introduce significant bias. Furthermore, the fastest chemical reactions in the plasma occur on timescales of the order of nanoseconds or even less, while the residence time in conventional reactors is several seconds. This huge discrepancy in timescales requires a computational time step that is several tens of orders of magnitude smaller than the total simulation time, to ensure that both the local chemical evolution and the chemical composition at the reactor outlet are accurately computed. Finally, the multidimensional topology of the discharges, whether filamentary,
streamer or localized micro-discharges, results in spatial heterogeneities, necessitating so-
phisticated modeling strategies.

In light of the aforementioned challenges, the development of a fully predictive simulation
tool for an entire DBD reactor is unreachable. However, integrating experimental data into
numerical simulations for which the most physically relevant modeling parameters have been
identified represents a technique under development for many different type of flows.²⁸⁻³⁰ One
challenge associated with this approach is how experimental information can be integrated
into the model to adjust numerical parameters.

Accordingly, this work employs a twofold strategy. Firstly, the experimental current
measurements and overall electric power consumption, which are commonly available exper-
imental data, are processed in order to facilitate their use in a one-dimensional numerical
simulation. Secondly, the reduction to a one-dimensional DBD reactor requires a statistical
approach to model the spatial heterogeneity of the discharges within the volume.

The experimental setup considered for simulations is a DBD reactor where the gas in-
jected consisted solely of CH₄ and O₂. This composition was selected for its potential to
capitalize on recent progress in the field of plasma-assisted combustion,²⁰,²²,³¹⁻³⁵ thereby fa-
cilitating the model validation process. These previous studies focused on methane oxidation
chemistry through a two-way plasma/combustion kinetic coupling applied to homogeneous
reactors as well as multi-dimensional simulations of NRP discharges. The model developed
by the Cerfacs group,²⁰,³⁵ which employs a coupling of a kinetic solver Cantera³⁶ embedding
a Boltzmann solver Bolos,¹⁷ was selected for this study.

The first section of this paper provides a concise overview of the DBD reactor facil-
ity, the available measurements and the electrical setup. Subsequently, the plasma-kinetic
model is presented, its relevance verified against literature results³¹ and an isolated discharge
is analyzed. Finally, the statistical modeling strategy specifically developed for this joint ex-
perimental/numerical study of the DBD reactor is presented, with a particular focus on the
reactor's heterogeneity behavior.
# Description of the DBD reactor setup

## Geometry, operating conditions and measurements

The reactor consists of a co-axial cylindrical tube constructed from quartz, featuring a central high-voltage electrode made of stainless steel (refer to Figure 1). Surrounding the outer diameter of the reactor is a copper sheet, functioning as the external electrode. The reactor's wall serves as a dielectric to prevent arc formation. Within the reactor, the plasma reactive zone is 7 cm in length, with a 2 mm gap between the two electrodes. The diameter of the central electrode is 1.2 cm, resulting in a discharge volume of 6.16 cm³.

The experiments were conducted at atmospheric pressure and a temperature of 400 K by placing the reactor in a temperature-controlled oven. The injected gas mixture, consisting of CH₄ and O₂ with mole fractions of 30% and 70%, respectively, was also preheated to 400 K. The temperature at the reactor outlet was not directly measured. However, the low conversion rate observed at the outlet, combined with the heat exchanges occurring within the oven, suggests that the rise in gas temperature is likely negligible.

Four distinct experiments were conducted at various overall flow rates of 70, 60, 50, and 40 ml min⁻¹, corresponding to residence times *t*ᵣ within the reactor of 5.27, 6.17, 7.37, and 9.22 s, respectively.

The gaseous products were analyzed using a four-way gas chromatograph (micro-GC Varian Quad CP 4900). The analysis employed Poraplot Q (10 m, 0.15 mm), CP-Sil 5 CB (8 m, 0.15 mm), CP-Wax 52 CB (4 m, 0.25 mm), and molecular sieve 5 A (10 m, 0.32 mm) columns with thermal conductivity detectors (TCD) for the online analyze of various gases (O₂, H₂, CO, CO₂, CH₄, C₂H₆, C₂H₄, C₃H₈, C₃H₆, CH₃OH, CH₃OCH₃, and C₂H₅OH) within a 3 min timeframe.
<figure><img src="image_3.png" /><figcaption>Figure 1: Experimental Set-up, Plasma Reactor.</figcaption></figure>

## Electrical setup for Lissajous curve acquisition

Plasma generation employed a sinusoidal voltage $V_s(t) = V_{s0} \cos(2\pi ft)$, with a frequency $f$ of 1 kHz and an amplitude $V_{s0}$ of 9 kV. This voltage and frequency were selected empirically to achieve a power consumption of 5 W, ensuring the generation of a stable cold plasma while minimizing gas heating.

A low-frequency generator (GBF) from TTi®, model TG1010A, was coupled to a signal amplifier from TREK®, model 20/20A, to produce the required voltage $V_s(t)$.

The reactor charge $Q(t)$ is measured by the capacitor connected to the copper counter-electrode, which has a sufficiently large capacitance ($C_c=4.5$ nF) to collect all the charges coming from the plasma:³⁷

$$Q(t) = C_c V_c(t) = C_r V_r(t) \qquad (1)$$

where $C_r$ and $V_r(t)$ are respectively the reactor's capacity and voltage and $V_c(t)$ the voltage across the capacitor.

Two high-voltage probes (Lecroy PPE 20 kV) were used to measure $V_s(t)$ and $V_c(t)$, see Figure 1. These measurements were recorded on a Lecroy oscilloscope (Waversurfer 3054, 500 MHz), as shown in Figure 2-left where the reactor's voltage is obtained as follows:
$$V_r = V_s - V_c.$$

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

<figure><img src="image_5.png" /><figcaption>Figure 2: Left: The measured voltage across the capacitor ($V_C$), the voltage across the reactor ($V_r$), and the corresponding experimental current. Right: Lissajous diagram of the DBD reactor.</figcaption></figure>

The power is calculated from the integration of the Lissajous curve, see Figure 2-right, issued from the Manley method:³⁸

$$P = \frac{1}{T} \int_{T} V_{r} I dt = f \oint_{T} V_{r} dQ = f C_{c} \oint_{T} V_{c} dV_{r} \quad (2)$$

where $T$ is the electrical period and $I$ the current. This curve plots the voltage across the capacitor $V_c$ as a function of the voltage applied to the reactor $V_r$. The reactor capacitance $C_r$, which takes two distinct values ($C_{r1}$, $C_{r2}$) depending on the presence or absence of plasma discharges, was determined from the slopes of the cycle lines in the Lissajous curve, as illustrated in Figure 2 right. A current probe (Stangene, model 0.1-0.5) was also used to obtain the precise shape of the current peaks during the plasma phase, see Figure 2 bottom left.

The current $I(t)$ and power $P$ represent the electrical measurements that were incorporated into the one-dimensional numerical model of the reactor. To simplify both the
measurement process and the numerical processing, only the first half-period of the current was used and repeated for each half-period.

# Plasma chemical kinetic computational model

## Solver description

The modeling of this reactor relies on a recently developed solver³⁴ that describes the temporal evolution of chemical species in a homogeneous (0D) reactor subjected to plasma discharges. Since the objective of DBD reactors is to achieve higher selectivity of specific chemical species, relying solely on non-thermal plasma chemistry is not sufficient. Thus, the set of chemical reactions involved is divided into two distinct subsets: one associated with molecules reactions, subset C (superscript c), describing classical chemical kinetic; and another one associated with electron collision reactions, subset B (superscript b), describing the plasma. The solver is based on the Cantera solver³⁶ in which the EEDF solver bolos¹⁷ has been implemented²⁰,³⁵ as done in previous plasma-assisted combustion studies.

For the zero dimensional simulations, the set of equations for the gas flow reduces to:

$$ \frac{\partial \rho Y_k}{\partial t} = \dot{\omega}_k, \quad \frac{\partial \rho \mathcal{E}}{\partial t} = \dot{\omega}_T + \dot{Q} \qquad (3) $$

With $\rho$ is the density, $Y_k$ the mass fraction, $\mathcal{E}$ the energy, $\dot{Q}$ the heat source term. The mass reaction rate of species k ($\dot{\omega}_k$) and the heat release rate $\dot{\omega}_T$ are separated into two distinct contributions, related to each reactions subset:

$$ \dot{\omega}_k = \dot{\omega}_k^c + \dot{\omega}_k^b = \sum_{j=C} \dot{\omega}_{kj}^c + \sum_{j=B} \dot{\omega}_{kj}^b \qquad (4) $$

$$ \dot{\omega}_k = \sum_{j=1}^{M} \nu_{kj} W_k \dot{Q}_j \quad \text{with} \quad \dot{Q}_j = K_j \prod_{k=1}^{N} [X_k]^{\nu_{kj}'} \qquad (5) $$

$$ \dot{\omega}_T = - \sum_{j \in P} \dot{Q}_j (\Delta u_j - \epsilon_{th}^j) - \sum_{j \in C} \dot{Q}_j \Delta u_j \qquad (6) $$
Where $\nu_{kj}$ is the net stoechiometric coefficient ($\nu_{kj} = \nu''_{kj} - \nu'_{kj}$), $\nu'_{kj}$ the stoechiometric coefficient for reactants, $\nu''_{kj}$ the stoechiometric coefficient for products, $W_k$ the molar mass of species k, $\dot{Q}_j$ the molar production rate of reaction j, $X_k$ the molar fraction, $\epsilon_{th}^j$ the energy lost during process j, $\Delta u_j$ the reaction energy change per mole and $K_j$ is the reaction rate. The definition of $K_j$ will depend on the nature of the reaction (B or C).

The discharges are characterized by an intense electric field that imparts a notably higher energy to electrons compared to other heavier species. Thus, a two-temperature fluid model is considered, treating the energy of molecules and electrons separately. The energy of electron population within a fluid element initially conforms to a Maxwellian distribution. However, frequent collisions with neutrals can cause changes in the distribution and generate ions, radicals, and excited species. Consequently, the Electron Energy Distribution Function (EEDF) is determined by solving the Boltzmann equation¹⁷ which accounts for the effects of the electric field and collisions from a set of collisional cross-sections³⁹:

$$ \frac{\partial}{\partial \epsilon} \left( \tilde{W} F_0 - \tilde{D} \frac{\partial F_0}{\partial \epsilon} \right) = \tilde{S} \qquad (7) $$

Where $F_0$ is the EEDF, $\epsilon$ the electron energy, $\tilde{W}$ the convection part representing the cooling by elastic collisions with less energetic particles (neutrals or electrons), $\tilde{D}$ the diffusive part representing heating by the field and by elastic collisions with more energetic particles and $\tilde{S}$ the source term depending on collisions. These parameters strongly depend on the electric field, the electron density ($n_e$) and the set of cross-sections.

The EEDF $F_0$, obtained by solving Eq. (7), is then used to compute the reaction rate $K_j$ and the effective electron temperature $T_e$ by integrating over electron energy $\epsilon$.

Eventually, for the set of reactions B and C the reaction rates are given by:

$$ K_j^b = \gamma \int_0^\infty \epsilon \sigma_j F_0 d\epsilon \qquad (8) $$

$$ K_j^c = B_j T_g^{\beta_j} \exp \left( -\frac{E_{aj}}{RT_g} \right) T_e^{\beta'_j} \exp \left( -\frac{E'_{aj}}{RT_e} \right) \qquad (9) $$
The electron temperature is obtained from the EEDF by:

$$T_e = \frac{2}{3} \int_0^\infty \epsilon^{3/2} F_0 d\epsilon \qquad (10)$$

Where $\gamma = \sqrt{2e/m_e}$ is a constant, $\epsilon = (v/\gamma)^2$ the electron energy, $\sigma_j$ the cross-section, $B_j$ the pre-exponential constant, $T_g$ the gas temperature, $\beta_j$ and $\beta'_j$ the temperature exponents and $E_{aj}$ and $E'_{aj}$ the activation energies.

This solver, referred to as the ‘Cantera-Boltzmann’ solver, and its associated chemical mechanism require validation before use. After validation, the Cantera-Boltzmann solver is used to analyze the chemical and temporal characteristics of the discharges produced in the DBD.

## Validation of the chemical kinetic mechanism

The plasma chemistry mechanism was established by using a set of cross-sections data for electron-neutral collisions taken from different data-sets of the LXCat database,⁴⁰ the Phelps and Hayashi databases⁴¹,⁴² have been used for $O_2$ and $CH_4$ respectively. The detailed plasma mechanism²⁰,³⁵ developed for methane-air plasma assisted combustion has been considered. It contains 71 species and 680 reactions. The plasma mechanism has been then combined with the GRI 3.0 oxidation mechanism⁴³ yielding to a combined plasma-combustion mechanism of 100 species and 964 reactions. This combined mechanism have been validated for a wide range of temperature [300-1500 K] and pressure [0.1-1 bar] by Cheng et al.²⁰. Given that the chemical scheme employed in this study is applied under conditions significantly different from those of plasma-assisted combustion, our initial step was to assess its relevance in the context of a DBD reactor.

This verification was conducted by comparing the results obtained with this solver to numerical data on methane oxidation in a DBD single-pulse reactor, as presented in the literature.³¹ The simulation was conducted using a constant volume homogeneous reactor
<figure><img src="image_6.png" /><figcaption>Figure 3: Time evolution of excited species and radicals generated by a single pulse in CH₄/O₂/Ar mixture, predicted by the model (solid lines) and by Z.Ecket et al.³¹ (dashed lines).</figcaption></figure>

evolving in time. The initial conditions consists of a mixture of methane and oxygen diluted in argon with mole fractions of CH₄: 0.0016, O₂: 0.003, and Ar: 0.995, at a temperature of 750 K and a pressure of 101 325 Pa. The electric field and current profiles used in the simulation were taken from the previous study.³¹ Figure 3 illustrates the predicted mole fractions of excited species and radicals generated by a single nanosecond discharge pulse in methane.

The minor discrepancies observed between the results of Eckert et al.³¹ and our results can be attributed to slight variations in the reaction scheme, as for example metastable argon self-quenching reaction, see Equation (11), which are absent from the kinetic scheme used here.³⁴

$$ \begin{gathered} Ar^* + Ar \rightarrow Ar + Ar \qquad (11) \\ Ar^* + Ar + Ar \rightarrow Ar + Ar + Ar \end{gathered} $$

These reactions by consuming slightly more Ar*, limit the production of H and O by colli-
sional quenching of metastable argon, see Equation (12):

$$ 
\begin{align*}
& \mathrm{Ar}^* + \mathrm{O}_2 \rightarrow \mathrm{Ar} + \mathrm{O} + \mathrm{O} \\
& \mathrm{Ar}^* + \mathrm{CH}_4 \rightarrow \mathrm{Ar} + \mathrm{CH}_2 + \mathrm{H} + \mathrm{H} \tag{12} \\
& \mathrm{Ar}^* + \mathrm{CH}_4 \rightarrow \mathrm{Ar} + \mathrm{CH} + \mathrm{H}_2 + \mathrm{H} \\
& \mathrm{Ar}^* + \mathrm{H}_2 \rightarrow \mathrm{Ar} + \mathrm{H} + \mathrm{H}
\end{align*}
$$

However, considering the striking similarity in the overall behavior of these two chemical mechanisms,³¹,³⁴ we confidently adopted the one developed by the Cerfacs group for the numerical simulations considered in this study.

## Characteristic times of a plasma discharge

In the field of combustion, the thickness of a flamelet as well as self-ignition delay time – the concepts behind the vast majority of turbulent flame models⁴⁴,⁴⁵ – must be precisely defined in order to construct a relevant statistical model of such reactive flows. Similarly, to model DBD using a statistical approach, defining the duration of a plasma discharge is essential. Furthermore, as in the case of flamelets or homogeneous gas self-ignition, the 0D plasma discharge serves as a fundamental element in this study, maintaining a quasi-stable temporal structure throughout the flow. This subsection therefore aims to describe the different characteristic times involved in a single discharge.

Naturally, various physical parameters, such as composition, temperature and the applied electric field, influence the shape of the discharge time structure. Some of the parameters in question have been determined through experimental means, such as the initial composition of the gas and the temperature. However, other quantities, such as the electric field experienced by the gas, cannot be directly measured and are difficult to determine from experiments. Consequently, in order to define the characteristic times of a discharge and evaluate their value, it is first necessary to conduct a parametric study using constant volume homogeneous reactors where the reduced electric field $E/N$ applied is varied. Here, $E$
is the electric field and N is the number of particles per unit volume, N = 1.84 × 10²⁵ m⁻³.

Similarly, the precise duration of each current peak, as illustrated in Figure 2 bottom left, is challenging to ascertain through experimental means. However, we assume that the maximum current of each peak p can be obtained from the current measurements. This value can then be used in numerical simulations by imposing the maximum electron density (EDN), denoted as $n_e|_{p,max}$, in accordance with the generalized formulation of Ohm's law:

$$n_e|_{p,max} = \frac{j|_{p,max}}{\mu_e q_e E} \qquad (13)$$

where $n_e|_{p,max}$ is the maximum value of EDN of each peak p, $\mu_e$ is the electron mobility, $q_e$ is the charge of the electron and $j|_{p,max}$ is the peak current density. This current density can be approximated by the measured peak current intensity, $j|_{p,max} = I|_{p,max}/S_c$, where $S_c = v/d$ is a characteristic surface evaluated from the volume of the reactor (v) and the inter-electrode distance (d). For this parametric study, the maximum EDN is maintained at a constant value throughout the parametric analysis ($n_e|_{p,max} = 3.5 \times 10^{15} m^{-3}$), that corresponds to a representative current peak observed in the experiment, calculated using Equation 13.

Within the DBD reactor, when the electric field applied to the gas exceeds a specific threshold (the breakdown voltage), the Townsend avalanche is initiated, resulting in a sudden increase in the EDN until current passes through the gas, and the electric field declines.

In order to represent this mechanism, each discharge in the parametric study is simulated in a 0D homogeneous reactor using the Cantera-Boltzmann solver and the combined plasma-combustion mechanism previously described. The simulation proceeds as follows: a constant reduced electric field (122 Td in the example shown in the graph on the left of Figure 4) is applied to the gas until the EDN reaches the critical value $n_e|_{p,max}$ ($3.5 \times 10^{15} m^{-3}$). At this point, referred to as the breakdown time $\tau_b$, the electric field is set to zero. Once the voltage dropped, the electrons recombined to ions, and the production rates of excited species ceased, as illustrated in Figure 4-left. Consequently, two additional characteristic times were
<figure><img src="image_7.png" /><figcaption>Figure 4: Definition of breakdown $\tau_b$, current $\tau_c$ and discharge $\tau_d$ times characterising the plasma discharge. Right: values of these three characteristic times as function of the reduced electric field. Left: Definition of the breakdown time ($\tau_b$), current time ($\tau_c$), and discharge time ($\tau_d$), characterizing the plasma discharge based on the time evolution of the reduced electric field, the EDN, and the reaction rates (as defined in Eq. 4) of the excited species $O_2(b_1)$ and $O_2(A_3)$.</figcaption></figure>

introduced: the current flow time, denoted by $\tau_c$, and the plasma discharge delay, designated by $\tau_d$. The first is defined as the time interval when the number of electrons exceeds a value of $n_e|_c$ which is taken as 3% of $n_e|_{p,max}$. The second is the total duration of the plasma discharge, encompassing the period of electric field application along with subsequent electronic and molecular recombinations and de-excitation. The criterion considered for the discharge end is when the chemical reaction rate of all plasma species (ions and excited species) falls below 5% of its maximum value. The numerical parameters employed to define the characteristic times of a single plasma discharge are summarized in Table 1.

<table><caption>Table 1: Numerical parameters used to define the characteristic times of a single plasma discharge</caption><thead><tr><td>Parameter</td><td>Value</td><td>Description</td></tr></thead><tbody><tr><td>$n_e|_{p,max}$</td><td>$3.5 \times 10^{15} \text{ m}^{-3}$</td><td>Representative EDN peak</td></tr><tr><td>$n_e|_c$</td><td>$3\% n_e|_{p,max}$</td><td>Discharge triggering condition on EDN</td></tr><tr><td>$\dot{\omega}_k^p / \dot{\omega}_{k(max)}^p$</td><td>$5\%$</td><td>End of plasma phase (at $\tau_d$)</td></tr></tbody></table>

The right graph of Figure 4 shows the three characteristic times resulting from the parametric study as a function of the applied reduced electric field. It can be observed that these
characteristic times, which are highly dependent on the reduced electric field, exhibit a range
of values, spanning from 10 ns to 1 µs. Consequently, they exhibit a period that is much less
than that of the sinusoidal signal applied to the electrodes. Additionally, simulations are
not feasible for values of the reduced electric field that are insufficient. This threshold value,
estimated to be approximately 115 Td, is of a comparable magnitude to the breakdown volt-
age (≈124 Td, see next subsection), although slightly inferior. In order to corroborate this
outcome with greater precision and to complete this investigation of characteristic times, a
more detailed chemical analysis of the discharge has been conducted.

Chemical analysis of a single discharge

The Phelps and Hayashi databases$^{41,42}$ have been used in conjunction with BOLSIG+ to estimate the breakdown voltage and the energy fraction used by each type of chemical process. The elastic and inelastic processes are described in detail for CH$_4$ and O$_2$ in Tables 2) and 3), respectively.

<figure><img src="image_8.png" /><figcaption>Figure 5: Reduced Townsend coefficients (left) and discharge energy branching (right) as a function of the reduced electric field. Results obtained with the Phelps and Hayashi databases$^{41,42}$ and BOLSIG+.</figcaption></figure>

Figure 5-left depicts the reduced Townsend coefficients as a function of the reduced elec-
tric field with $\alpha$ and $\eta$ representing the ionization coefficient and the attachment coefficient, respectively. First, the ionization coefficient sharply increases with $E/N$, whereas in contrast, the attachment coefficient shows a slight increase at lower $E/N$ values [0-50 Td], and decreases slightly at $E/N > 50$ Td. Breakdown occurs when ionization becomes the dominant process over attachment, which corresponds approximately to 124 Td, see Figure 5-left. This result, obtained with BOLSIG+ (which accounts only for plasma reactions), is slightly higher than the 115 Td obtained with the Cantera-Boltzmann solver (see Figure 4), which also includes neutral and ionic reactions.

Additionally, the resulting discharge energy branching is displayed as a function of the reduced electric field in Figure 5-right. For a reduced electric field strength of approximately 124 Td, the figure shows that the majority of the electron energy is expended on the dissociation of $O_2$ to O molecules (70%), which corresponds to reactions R3, R13, R14 and R15 in table 3. In contrast, a smaller proportion is allocated to the dissociation of $CH_4$ to $CH_3$ and H (11%) corresponding to reaction R4 in table 2, and another 11% on the vibrational excitation of $O_2$, (R4-12 in table 3).

<table><caption>Table 2: List of processes in the Hayashi database$^{42}$ for $e^-$-$CH_4$ processes and the associated reactions considered in this work. $\Delta\epsilon$ is the energy threshold of a process</caption><thead><tr><th>#</th><th>Process</th><th>Reaction</th><th>Δε [eV]</th></tr></thead><tbody><tr><td>R1</td><td>CH<sub>4</sub> → CH<sub>4</sub><sup>-</sup></td><td>e<sup>-</sup> + CH<sub>4</sub> → CH<sub>4</sub><sup>-</sup></td><td>0.0</td></tr><tr><td>R2</td><td>CH<sub>4</sub> → CH<sub>4</sub>·(0·159eV)</td><td>e<sup>-</sup> + CH<sub>4</sub> → e<sup>-</sup> + CH<sub>4</sub>(v<sub>24</sub>)</td><td>0.162</td></tr><tr><td>R3</td><td>CH<sub>4</sub> → CH<sub>4</sub>·(0·37eV)</td><td>e<sup>-</sup> + CH<sub>4</sub> → e<sup>-</sup> + CH<sub>4</sub>(v<sub>13</sub>)</td><td>0.361</td></tr><tr><td>R4</td><td>CH<sub>4</sub> → CH<sub>4</sub>·(7·9eV)</td><td>e<sup>-</sup> + CH<sub>4</sub> → e<sup>-</sup> + CH<sub>3</sub> + H + ΔE</td><td>7.9</td></tr><tr><td>R5</td><td>CH<sub>4</sub> → CH<sub>4</sub><sup>+</sup></td><td>e<sup>-</sup> + CH<sub>4</sub> → 2e<sup>-</sup> + CH<sub>4</sub><sup>+</sup> + ΔE</td><td>12.9</td></tr></tbody></table>
<table><caption>Table 3: List of processes in the Phelps database$^{41}$ for e⁻-O₂ processes and the associated reactions considered in this work.</caption><thead><tr><th>#</th><th>Process</th><th>Reaction</th><th>Δε [eV]</th></tr></thead><tbody><tr><td>R1</td><td>O<sub>2</sub> → O<sub>2</sub>(rot)</td><td>e<sup>-</sup> + O<sub>2</sub> → e<sup>-</sup> + O<sub>2</sub> + ΔE</td><td>0.02</td></tr><tr><td>R2</td><td>O<sub>2</sub> (ATT)</td><td>e<sup>-</sup> + O<sub>2</sub> + M → O<sub>2</sub><sup>-</sup> + M</td><td>0.0</td></tr><tr><td>R3</td><td>O<sub>2</sub> → O<sup>-</sup> + O</td><td>e<sup>-</sup> + O<sub>2</sub> → O<sup>-</sup> + O</td><td>0.0</td></tr><tr><td>R4</td><td>O<sub>2</sub> → O<sub>2</sub>(v<sub>1</sub>)</td><td>e<sup>-</sup> + O<sub>2</sub> → e<sup>-</sup> + O<sub>2</sub>(v<sub>1</sub>)</td><td>0.19</td></tr><tr><td>R5</td><td>O<sub>2</sub> → O<sub>2</sub>(v<sub>1</sub>res)</td><td>e<sup>-</sup> + O<sub>2</sub> → e<sup>-</sup> + O<sub>2</sub>(v<sub>1</sub>)</td><td>0.19</td></tr><tr><td>R6</td><td>O<sub>2</sub> → O<sub>2</sub>(v<sub>2</sub>)</td><td>e<sup>-</sup> + O<sub>2</sub> → e<sup>-</sup> + O<sub>2</sub>(v<sub>2</sub>)</td><td>0.38</td></tr><tr><td>R7</td><td>O<sub>2</sub> → O<sub>2</sub>(v<sub>2</sub>res)</td><td>e<sup>-</sup> + O<sub>2</sub> → e<sup>-</sup> + O<sub>2</sub>(v<sub>2</sub>)</td><td>0.38</td></tr><tr><td>R8</td><td>O<sub>2</sub> → O<sub>2</sub>(v<sub>3</sub>)</td><td>e<sup>-</sup> + O<sub>2</sub> → e<sup>-</sup> + O<sub>2</sub>(v<sub>3</sub>)</td><td>0.57</td></tr><tr><td>R9</td><td>O<sub>2</sub> → O<sub>2</sub>(v<sub>4</sub>)</td><td>e<sup>-</sup> + O<sub>2</sub> → e<sup>-</sup> + O<sub>2</sub>(v<sub>4</sub>)</td><td>0.75</td></tr><tr><td>R10</td><td>O<sub>2</sub> → O<sub>2</sub>(a<sub>1</sub>)</td><td>e<sup>-</sup> + O<sub>2</sub> → e<sup>-</sup> + O<sub>2</sub>(a<sub>1</sub>)</td><td>0.977</td></tr><tr><td>R11</td><td>O<sub>2</sub> → O<sub>2</sub>(b<sub>1</sub>)</td><td>e<sup>-</sup> + O<sub>2</sub> → e<sup>-</sup> + O<sub>2</sub>(b<sub>1</sub>)</td><td>1.627</td></tr><tr><td>R12</td><td>O<sub>2</sub> → O<sub>2</sub>(4·5eV)</td><td>e<sup>-</sup> + O<sub>2</sub> → e<sup>-</sup> + O<sub>2</sub>(A<sub>3</sub>)</td><td>4.5</td></tr><tr><td>R13</td><td>O<sub>2</sub> → O<sub>2</sub>(6·0eV)</td><td>e<sup>-</sup> + O<sub>2</sub> → e<sup>-</sup> + 2O + ΔE</td><td>6.0</td></tr><tr><td>R14</td><td>O<sub>2</sub> → O<sub>2</sub>(8·4eV)</td><td>e<sup>-</sup> + O<sub>2</sub> → e<sup>-</sup> + O + O(<sup>1</sup>D) + ΔE</td><td>8.4</td></tr><tr><td>R15</td><td>O<sub>2</sub> → O<sub>2</sub>(9·97eV)</td><td>e<sup>-</sup> + O<sub>2</sub> → e<sup>-</sup> + O + O(<sup>1</sup>S) + ΔE</td><td>9.97</td></tr><tr><td>R16</td><td>O<sub>2</sub> → O<sub>2</sub><sup>+</sup></td><td>e<sup>-</sup> + O<sub>2</sub> → 2e<sup>-</sup> + O<sub>2</sub><sup>+</sup></td><td>12.06</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>Figure 6: Molar production rate at four times $t_1$: the maximum intensity of the discharge ($t_1 = \tau_b$), $t_2$: a few nanoseconds after the pulse ($t_2 = \tau_b + 5$ ns), $t_3$: at the end of the plasma phase ($t_3 = \tau_d$), $t_4$: further away from the discharge ($t_4 = \tau_d + 3$ µs).</figcaption></figure>

To further analyze the chemical processes, the rates of progress ($Q_j$) of the 15 most important reactions obtained with Cantera-Boltzmann at four different times during and after the discharge are shown in Figure 6. The four times analysed, $t_1, t_2, t_3, t_4$, correspond
respectively to (i) the maximum intensity of the discharge ($t_1 = \tau_b$), (ii) a few nanoseconds after the pulse ($t_2 = \tau_b + 5$ ns), (iii) at the end of the plasma phase as defined in the previous section ($t_3 = \tau_d$) and (iv) further away from the discharge ($t_4 = \tau_d + 3$ ns). Within the current pulse, at $t_1$, as expected, the dissociation reactions of $O_2$ and $CH_4$ and the excitation of $O_2$ are the most significant (see figure 6 top left), as obtained with BOLSIG+. Immediately after the pulse, at $t_2$, the dominant processes are the quenching of excited states of atomic oxygen, facilitating $CH_4$ dissociation and $O_2$ excitation. Chemical reactions of neutral species become more important, while electron impact reactions hardly occur. At the end of the plasma phase $t_3$, as expected, neutral reactions become more important, in particular the dissociation of $CH_4$ by the OH radical and the formation of $O_3$, but there is still an important reaction of de-excitation of $O_2$. Finally, further away from the discharge, at $t_4$, common gaseous chemical reactions with neutral species become dominant, especially those involving radicals (O, OH, $CH_3$, ...).

# Statistical modeling of the DBD reactor

## Averaged equations and basic assumptions

The currently available experimental data does not permit a detailed understanding of the behaviour of the discharges within the reactor. To simulate this device using a statistical approach, it is thus necessary to introduce modelling parameters. The aim of this section is therefore to properly set out the theoretical analysis, based on the characteristic times and lengths of the physical mechanisms, which will enable these parameters to be introduced. Given that the practical goal of research on these reactors is to enhance the selectivity of the target product, our primary focus lies on analyzing the total and partial mass balance
equations:

$$ \frac{\partial \rho}{\partial t} + \frac{\partial \rho u_i}{\partial x_i} = 0 \qquad (14) $$

$$ \frac{\partial \rho Y_k}{\partial t} + \frac{\partial \rho u_i Y_k}{\partial x_i} = \frac{\partial \mathcal{F}_{ik}}{\partial x_i} + \dot{\omega}_k \qquad (15) $$

where $u_i$ the velocity vector ($i \in [1, 3]$), and $\mathcal{F}_{ik}$ the diffusion flux of species $k$ (molecular and drift). The convection and molecular diffusion terms in equations (14-15) are small compared to the unsteady term because the characteristic times of the plasma discharges are negligible compared to the residence time of a fluid particle in the reactor. However, they cannot be neglected a priori as they are essential to characterise the amount of target species leaving the reactor. Considering the absence of evidence suggesting spatial homogeneity during discharges, all three dimensions are crucial. However, given the shape of the reactor, it makes more sense here to use cylindrical coordinates $x$ (longitudinal), $r$ (radial), $\theta$ (angular).

Given the constant flow rate, the stability of chemical measurements at the reactor outlet, and the consistent light emission observed, the reactor's behavior can be considered statistically stationary. Consequently, the thermodynamic quantities of interest are the ensemble averages, akin to a RANS (Reynolds Averaged Navier-Stokes) approach for turbulent flow. So the averaging process encompasses both spatial and temporal dimensions. From a spatial perspective, we assume statistical homogeneity in both the angular $\theta$ and radial $r$ directions.

The first assumption implies that discharges are either uniform over the surface of the electrodes or heterogeneous but randomly distributed in a manner where a large number of discharges cover the entire electrode surface, thus achieving statistical homogeneity.

The second assumption, namely statistical homogeneity in the radial direction, implies that the chemical effect of the plasma discharge, which propagates radially from one electrode to the other, remains similar regardless of the radial position $r$. To evaluate the validity of this hypothesis, 2D simulations were conducted to determine the distribution of the electric field across a reactor slice, i.e., in the radial ($r$) direction, in a non-reactive case. The
SOMAFOAM²³ code was used to compute the electric field ($\mathbf{E} = -\nabla\phi$) between the two electrodes by solving Poisson's equation for the electric potential:

$$ \Delta\phi = -\frac{\rho}{\varepsilon} \qquad (16) $$

where $\rho$ is the charge density and $\varepsilon$ the electric permittivity. This numerical tool accounts for different media, each with its own permittivity value (The relative permittivity for gas and dielectric are respectively 1.006 and 4.3), using a multi-zone approach in which the zones are connected through appropriate boundary conditions. This feature enables an accurate representation of the influence of the dielectric material and the asymmetric configuration on the electric field distribution in the gas. The typical electric field obtained is presented in Figure 7-left, while Figure 7-right shows the reduced radial electric field profiles for different electrode voltages.

<figure><img src="image_10.png" /><figcaption>Figure 7: Reduced electric field within a section of the reactor (left) obtained from non-reactive 2D simulation with $V_{r0} = 9$ kV and corresponding profile along diameter with three different values of $V_{r0}$ (right).</figcaption></figure>

As evidenced by these results, the assumption of statistical homogeneity in the radial direction is clearly inconsistent with reality, as each discharge propagates from one electrode to the other within a variable electric field.

Furthermore, this assumption also leads to neglect interactions between the discharge and the reactor walls, as well as the potential presence of a boundary layer in the composi-
tion resulting from catalytic reactions. Despite these limitations, this assumption has been retained for the sake of simplicity, but it is evident that further research is required in order to propose a more accurate model.

Eventually, only the longitudinal direction $x$ is considered from an average standpoint, thus the balance equations become, on average:

$$ \frac{d \bar{\rho} \tilde{u}}{dx} = 0 \qquad (17) $$

$$ \frac{d \bar{\rho} \tilde{u} \tilde{Y}_k}{dx} = \frac{d \tilde{\mathcal{F}}_k}{dx} - \frac{d \bar{\rho} u'' \tilde{Y}_k''}{dx} + \bar{\omega}_k \qquad (18) $$

Given the potential density variations resulting from changes in chemical composition, the Favre averaged operator $\tilde{\mathcal{F}}$ have been introduced. Although the flow is not turbulent, fluctuations in velocity $u''$ and composition $Y_k''$ induced by the plasma (ionic winds) are likely to be non-negligible. As a result, the first and second terms of the right hand side of Eq. (18), representing the unresolved fluxes due to diffusion and turbulence-induced-plasma, respectively, require specific models. However, in this study, our focus is on modelling the average chemical source term, $\tilde{\omega}_k$. Consequently, for the sake of simplicity, the unresolved fluxes are neglected, thereby reducing the complexity of the problem to a classical plug flow reactor simulation. This can be solved either in the longitudinal direction, $x$, or in the residence time space, $t_r$ ($dx = \tilde{u}dt_r$):

$$ \frac{d\tilde{Y}_k}{dx} = \frac{\tilde{\Omega}_k}{\tilde{u}} \quad \text{or} \quad \frac{d\tilde{Y}_k}{dt_r} = \tilde{\Omega}_k \qquad (19) $$

where the Favre averaged chemical rate is $\tilde{\Omega}_k = \bar{\omega}_k / \bar{\rho}$.

Finally, the main hypotheses are summarized in Table 4.
<table><caption>Table 4: Summary of Hypotheses Adopted</caption><thead><tr><td>Hypothesis</td><td>Justification</td></tr></thead><tbody><tr><td>statistically stationary</td><td>constant measurements</td></tr><tr><td>angular θ spatial homogeneity</td><td>uniform spatial distribution of discharges</td></tr><tr><td>radial r spatial homogeneity</td><td>simplicity</td></tr><tr><td>negligible unresolved fluxes</td><td>simplicity</td></tr></tbody></table>

## Numerical implementation of measured electric signals

The simulation of the reactor, in conjunction with that of its electrical circuit, is a topic that extends beyond the scope of this study, which focuses on the one-dimensional modelling of the reactor based on a statistical approach. Consequently, it is not possible to directly exploit the voltage signals measured, see Figure 1. Instead, the following section outlines the methodology employed to derive the voltage experienced by the gas and the current in simulations from the measured intensity signal.

Experimental measurements yield the sinusoidal voltage applied to the electrodes, denoted here $V_a(t)$, and the corresponding current $I(t)$ flowing through the gas.³⁷ A single period of this voltage (red line) is presented in Figure 8. The current (orange line) displays a multitude of peaks that are characteristic of plasma discharges in DBD reactors. Nevertheless, the electric field experienced by the gas remains unknown due to the presence of dielectric material between the gas and electrodes, the current flowing during plasma discharges which modifies the electric field, and the accumulation of charges on walls.

For reasons of evident simplicity, we postulated that at each measured current peak (78 for a half-period in this case) there is a corresponding drop to zero of the voltage experienced by the gas within the discharge. Once the electrical effects of the discharge have dissipated, the voltage exhibits a rapid increase until the commencement of the subsequent discharge cycle. During this phase of rising voltage, its value reaches a level comparable to the breakdown voltage.

Consequently, the voltage experienced by the gas is limited by the breakdown voltage. The voltage signal is therefore in the form of a square wave, as illustrated in Figure 4-left,
<figure><img src="image_11.png" /><figcaption>Figure 8: Left: measured applied voltage $V_s$ (red), measured current peaks (orange), and breakdown voltage limit (dashed black). Right: current peaks (orange) and reduced electric fields experienced by the gas $(E/N)_g$ for the first five peaks of the half-period.</figcaption></figure>

which depicts a single discharge. The Figure 8-right illustrates the reduced electric field $(E/N)_g$ applied to the gas that was employed in numerical simulations of several plasma discharges. The duration of each peak is not determined by the measured signal but by the numerical simulation itself, in order to obtain the same value of the maximum current as that measured for each peak (orange dots in Figure 8-right), as given by Equation 13.

The simulation enables the determination of the growth of the EDN, $n_e$, in response to the applied electric field of each peak. Once $n_e$ reaches its local maximum $n_e|_{p,\text{max}}$, the electric field is set to zero, and it increases again just before the occurrence of the next peak. This process is repeated for each subsequent peak. Consequently, only the number of peaks, their instant of appearance and their respective maximum value are derived from the measured signal and utilised in the simulations.

## Simulation of the discharge train over one half period

Having defined the characteristic times of a discharge and the maximum current of each discharge, it is now possible to simulate a homogeneous reactor which would be subjected to the same electric current as that measured by applying to the gas a maximum reduced electric field value slightly higher than the breakdown value. Following the procedure described
previously for determining the duration of each discharge, the maximum reduced electric field value retained, $E/N|_{g,max} \approx 122$ Td, was determined iteratively in order to obtain a total power consumption of 5 W, as measured experimentally. It can be observed that this value is in very close proximity to the breakdown value determined with BOLSIG+, thereby providing evidence that the methodology employed to integrate the experimental current signal into the simulations is relevant and effective.

Figure 9 presents a comparison of the temporal evolution of the current obtained in simulation with those obtained in experiments.

Since the maximum value of each measured current peak is a target for the simulation (see Equation 13), any discrepancies between the simulation and experiments are attributable solely to numerical errors. As shown in Figure 10, a parametric study on the time step demonstrates that, for a time step of $10^{-9}$s, the current peak is reproduced with sufficient accuracy. Accordingly, the results in Figure 9 were obtained using this sufficiently small time step to ensure the target was achieved with the desired precision. The simulation thus exhibits excellent quantitative agreement with the experimental data.

<figure><img src="image_12.png" /><figcaption>Figure 9: Current signal from measurements (left) and obtained from numerical simulations (right) for the first half period.</figcaption></figure>

Figure 11 illustrates the instantaneous specific power ($p = n_e q_e \mu_e E^2$) obtained over half a period (left) and a zoom for two different discharges (right). As anticipated, the power maximums associated with each discharge have significant magnitude and are correlated to
<figure><img src="image_13.png" /><figcaption>Figure 10: Numerical current ($I_{num}$) for different time steps ($\Delta t$) compared with the experimental target current ($I_{exp}$). Current of discharge 6 (left) and current of discharge 15 (right).</figcaption></figure>

the current, see Figure 11 right. We can also observe that the shape of the curves remains the same, regardless of the amplitude, which is a consequence of the shape chosen for the applied voltage.

<figure><img src="image_14.png" /><figcaption>Figure 11: Specific power ($kW/m^3$) obtained from numerical simulations for the first half period (left) and zoom for discharges 6 and 15 (right).</figcaption></figure>

Figure 12 illustrates the temporal evolution of the mole fractions of the principal chemical species ($CH_4$, $O_2$, $CO$, $CO_2$) during the initial half-period. The concentrations of the reactants, $CH_4$ and $O_2$, show a sudden drop with each discharge, due to electron collisions, and continue to decrease gradually between discharges. In contrast, the concentrations of
the two main products, CO and CO₂, do not show sudden variations with each discharge, but only increase gradually between discharges. Over the brief interval corresponding to

<figure><img src="image_15.png" /><figcaption>Figure 12: Molar progress variable ($X_k - X_{k|0}$) of O₂ and CH₄, molar fraction of CO and CO₂ for the first-half period with all applied current peaks. $X_{O2|0} = 0.7$ and $X_{CH4|0} = 0.3$.</figcaption></figure>

this half-period (less than one 1 ms), the quantities consumed (respectively produced) are markedly low. While slight variations in the behavior of plasma discharges may occur over time, the reactor's low conversion suggests that this behavior remains relatively consistent throughout. It is therefore the accumulation of these very small quantities over the periods that the model must be able to reproduce, which makes simulations of the entire reactor and comparison with experimental data particularly challenging. Even the slightest error can accumulate and result in significant drift over time. The procedure that consists of calibrating the model parameters on the experimental results is therefore essential. Although the model cannot be considered predictable, it is nonetheless a powerful tool for analysing the reactor's behavior and complementing experimental data.
# Modeling of spatial heterogeneities

## Monte Carlo methodology and assumptions

If the discharges were homogeneous in volume, then the homogeneous reactor approach that led to the results of the previous section, see Figure 9, obtained by simulating half a period, could be applied to the entire residence time of the reactor.

Furthermore, given that plasma discharges are highly unlikely to be homogeneous, this approach is no longer valid. In order to account for the aforementioned heterogeneities, a Monte Carlo method was employed. This approach involves seeding a sufficient number of particles to describe all the discharges, with each particle only experiencing a part of these discharges. Thus the mean state can be regarded as the average of an ensemble of *M* homogeneous reactor states.

The Figure 13 depicts the sum of the current signals emanating from the seeded particles (six in this instance). The summed signal is identical to the measured signal, see Figure 9-left. However, in order to limit the number of particle required for the Monte-Carlo method, the distribution of the current peaks across the particles is not random. In this illustrative example, comprising six particles, the first particle experiences the first peak, the second particle the second peak, and so on until the seventh peak, which is assigned to the first particle. This assignment sequence is repeated until all the peaks have been assigned to one and only one particle. This strategy allows for a reduction in the number of particles required while maintaining the total power consumed. Consequently, in this instance, each particle underwent 13 peaks ($N_p = 13$), see Figure 13, in order to obtain the total number of discharges occurring over half a period, which was found to be $N_{tot} = 78$.

Figure 14(a) illustrates the time evolution of the molar fraction of $O_2$ over one period. The results indicate that the majority of particles exhibit a similar trend in $X_{O_2}$ values, with the curves almost overlapping. However, there is one exception, namely particle 3, for which significantly less oxygen is consumed. This behavior is attributed to the lower values of the
<figure><img src="image_16.png" /><figcaption>Figure 13: Overall current corresponding to a Monte-Carlo simulation using 6 particles (Left) and example of current corresponding to three different particles (right).</figcaption></figure>

maximum current peaks, see Figure 13 (right). The power consumed by this particle over one period is $p \approx 0.07 \text{ W cm}^{-3}$, whereas the others exhibit a significantly higher consumption rate of $p \approx 0.14 \text{ W cm}^{-3}$.

Figure 14(b) shows the same graph as figure 14(a), but where the simulation is continued over six periods (12 half-periods), always applying the same set of discharges to each particle. The small deviations observed in one period accumulate in each period and lead to a significant divergence in the composition of each particle. For particle 3, the trajectory is even completely different from that of the other particles.

Since there is no reason why a particle should always experience the same set of discharges, a rotation of the sets of discharges is applied to each half-period. Thus the set of discharges experienced by particle 2 in the first half-period will be experienced by particle 1 in the second half-period, and so on. Figure 14(c) shows the evolution of the $O_2$ mole fraction of each particle, accounting for this rotation, over the six periods. In this example (13 peaks per particle) a complete rotation occurs after 3 periods (6 half-periods), so the 6 curves intersect exactly after 3 ms. This means that each particle will have experienced all 78 peaks after 3
<figure><img src="image_17.png" /><figcaption>Figure 14: Evolution of the O2 molar progress variable ($X_k - X_{k|0}$) with $X_{O2|0} = 0.7$, for a case with 6 particles (each with 13 peaks, M = 6) over one period (a) and over 6 periods: without rotation of peaks (b), with rotation (c) and comparison of different M values with rotation (d).</figcaption></figure>
periods.

All particles will therefore behave similarly, and if the number of periods is sufficient, the composition of a single particle will be sufficient to describe the evolution of the composition in the reactor. The essential modelling parameter is therefore the fraction of the set of measured discharges that is actually experienced by a particle. This parameter is thus proportional to the spatial heterogeneity of the discharges occurring within the reactor. Consequently, it is possible to introduce a discharge heterogeneity factor, which is an essential parameter for modelling DBD reactors:

$$ \beta = 1 - \frac{N_p}{N_{tot}} = 1 - \frac{1}{M} \qquad (20) $$

Where $N_p$ is the number of discharges experienced by each particle, $N_{tot}$ the total number of discharges experienced by the reactor and $M$ the number of particles seeded in simulations. The total number of discharges $N_{tot}$ is assumed to be sufficiently large to preserve statistical homogeneity. Therefore, the knowledge of this heterogeneity factor allows the determination of the minimal number of particle required for Monte-Carlo simulations. In the case of spatial homogeneity, where $\beta = 0$, $M = 1$, $N_p = N_{tot}$, only one particle that has experienced all discharges is required. Conversely, in the case of maximum heterogeneity, each particle experiences one discharge per period: $\beta = (N_{tot} - 1)/N_{tot}$, $M = N_{tot}$, $N_p = 1$ (or even less if the particle has not experienced a discharge for each period). As this factor is not directly measurable, it will be determined iteratively in order to ensure that the simulation provides an accurate representation of the chemical composition at the reactor outlet, in accordance with the experimental results. Figure 14(d) shows results for four simulations with different values of the heterogeneity factor: $\beta \in \{0.83, 0.89, 0.92, 0.96\}$ which corresponds respectively to a number of seeded particles of $M \in \{6, 9, 13, 26\}$. For each case, the residence time of the simulation is 12 ms which corresponds to ten periods.

The results presented in Figure 14(d) indicate that when each discharge occupies a significant
icant portion of the reactor volume—i.e., when the reactor approaches a spatially homogeneous state (low values of the heterogeneity factor $\beta$)—the number of discharges experienced by each injected particle is high and tends to align closely with the measured peak count, thereby leading to a high chemical conversion. Conversely, when the discharges are distributed in a highly heterogeneous manner within the reactor, the probability of an injected particle being traversed by a discharge is low, resulting in a low level of chemical conversion.

## Reactor type transitions: from plasma to standard chemistry

As mentioned earlier, one of the challenges in simulating this reactor is the need to use a very small time step, which must be less than $10^{-9}$ s in order to accurately account for the plasma chemistry, while the total residence time in the reactor is about ten seconds, corresponding to ten thousand periods and more than ten billion time steps. Consequently, it is necessary to implement various techniques for reducing the computational time.

As illustrated in Figure 4, the characteristic time of a discharge is less than a microsecond. Consequently, for one hundred discharges (approximately a half-period), it corresponds to less than one hundred microseconds. Thus, there remains more than 400 µs per half-period (80% of the time) during which the thermodynamic state is in quasi-equilibrium, i.e. out of the plasma phase, and does not require such small time steps. This ratio is further accentuated when heterogeneity in discharges is considered (see preceding section), with each particle undergoing only a subset of the total number of discharges. In order to optimise the calculation time, an adaptive time step was therefore employed. During the course of a discharge, the time step is on the order of $10^{-9}$ s, while out of the plasma phase, the time step is on the order of $10^{-6}$ s. This variation is considerable and is carried out numerically in a progressive manner.

In cases where the reactor behaviour is highly heterogeneous and where particles undergo only a few discharges per period, the majority of the gas states are out of the plasma phase. It is therefore worth considering whether the use of a chemical solver and dedicated
kinetic scheme designed for plasma is justified. Indeed, kinetic schemes are suitable for use
in the specific conditions for which they were developed and it is known that they can lead
to significant biases when used in different conditions. Consequently, in light of the char-
acteristic timescales of the discharges described above, we have also opted to implement a
procedure that enables the type of reactor, the solver and kinetic scheme to be switched
throughout the calculation. Consequently, the reactor using the Cantera-Boltzmann solver
is employed exclusively during the discharge phase and is refered as R1, while the second
reactor, which utilizes a standard Cantera solver (a homogeneous reactor coupled with the
standard GRI30 scheme<sup>43</sup>), is applied out of the plasma phase and is referred to as R2.
Furthermore, this sequence of reactors enables the separation of the constant-volume trans-
formations, whose characteristic times are very short, typically for plasma discharges, from
the constant-pressure transformations, which correspond to the operation of the DBD reac-
tor. Accordingly, R1 is a constant volume reactor while R2 is a constant pressure reactor.
This reactor change also accounts for heat losses, as it is assumed that the reactor operates
closer to isothermal conditions than adiabatic ones. Accordingly, the temperature is reset to
exactly 400 K at each switch from R1 to R2.

However, precautions must be taken when using this technique. Since some species are
present in one chemical scheme and absent in the other, perfect continuity of reaction rates is
not guaranteed. To check the effect of this reactor switch, a chemical analysis was performed
at two different times: t₃ = τ₉ (less than 1 µs), corresponding to the end of the plasma phase
as defined in Figure 4, and t₅ = 0.25 ms, corresponding to an infinitely long time compared
to t₃. Two simulations were carried out: one using only the Cantera-Boltzmann solver with
the GRI30/plasma mechanism until t = t₅ (R1) and another one using the coupled Cantera-
Boltzmann solver during the plasma phase only (until t = t₃) and the standard Cantera
solver (with standard GRI30 mechanism) outside the plasma phase, from t₃ to t₅ (R1 + R2).

Figure 15 illustrates the results of the chemical analysis conducted in the two simulations.
In the R1+R2 simulation, the two times under consideration, t3 and t5, are located in the
<figure><img src="image_18.png" /><figcaption>Figure 15: Molar production rate at two different times, $t_3$: at the end of the plasma phase ($t_3 = \tau_d$) and $t_5$: long time after the discharge ($t_5 = 0.25$ ms) and for single reactor case (R1) and switching reactor case (R1 + R2).</figcaption></figure>
<figure><img src="image_19.png" /></figure>

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

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

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

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

<figure><img src="image_24.png" /><figcaption>Figure 16: Molar progress variable ($X_k - X_{k|0}$) of O₂ (a) and CH₄ (b), molar fraction of CO (c) and CO₂ (d), temperature (e) and pressure (f) for single reactor case (R1) and switching reactor case (R1 + R2) for a residence time of 6 ms. $X_{O2|0} = 0.7$ and $X_{CH4|0} = 0.3$.</figcaption></figure>
R2 reactor. The time, t3, corresponds to the first time-iteration following the switch from R1 to R2. The objective is to assess the influence of the reactor change on the chemical reaction rates. The precise values of the reaction rates of the most important chemical reactions common to the chemical schemes of R1 and R2, as shown in this figure, are summarized in Tables 5-6. The results show that, for the two post-discharge times studied, although the de-excitation reactions of O₂ are not considered when using the standard solver and GRI30 (R2), the conversion rates of the other reactions present in the two reactors are perfectly equivalent, see Tables 5-6.

<table><caption>Table 5: Molar production rates at $t_3$ of the 5 most important neutral species reactions, for single reactor case (R1), and switching reactor case (R1 + R2).</caption><thead><tr><th rowspan="3">Reaction</th><th colspan="2">Q&#775;<sub>j</sub> [mol.m<sup>3</sup>.s<sup>-1</sup>]</th></tr><tr><th colspan="2">t<sub>3</sub></th></tr><tr><th>R1</th><th>R1+R2</th></tr></thead><tbody><tr><td>CH<sub>4</sub> + OH &rarr; CH<sub>3</sub> + H<sub>2</sub>O</td><td>0.69</td><td>0.69</td></tr><tr><td>O + 2O<sub>2</sub> &rarr; O<sub>2</sub> + O<sub>3</sub></td><td>0.55</td><td>0.55</td></tr><tr><td>H + 2O<sub>2</sub> &rarr; HO<sub>2</sub> + O<sub>2</sub></td><td>0.37</td><td>0.37</td></tr><tr><td>H + O<sub>2</sub> + M &rarr; HO<sub>2</sub> + M</td><td>0.21</td><td>0.21</td></tr><tr><td>CH<sub>3</sub>O + O<sub>2</sub> &rarr; CH<sub>2</sub>O + HO<sub>2</sub></td><td>0.022</td><td>0.022</td></tr></tbody></table>

<table><caption>Table 6: Molar production rates at $t_5$ of the 5 most important neutral species reactions, for single reactor case (R1), and switching reactor case (R1 + R2).</caption><thead><tr><th rowspan="3">Reaction</th><th colspan="2">Q&#775;<sub>j</sub> &times; 10<sup>-3</sup> [mol.m<sup>3</sup>.s<sup>-1</sup>]</th></tr><tr><th colspan="2">t<sub>5</sub></th></tr><tr><th>R1</th><th>R1+R2</th></tr></thead><tbody><tr><td>2CH<sub>3</sub>(+M) &rarr; C<sub>2</sub>H<sub>6</sub>(+M)</td><td>1.2</td><td>1.2</td></tr><tr><td>CH<sub>3</sub>O + O<sub>2</sub> &rarr; CH<sub>2</sub>O + HO<sub>2</sub></td><td>0.6</td><td>0.6</td></tr><tr><td>CH<sub>4</sub> + OH &rarr; CH<sub>3</sub> + H<sub>2</sub>O</td><td>0.6</td><td>0.59</td></tr><tr><td>CH<sub>3</sub> + HO<sub>2</sub> &rarr; CH<sub>3</sub>O + OH</td><td>0.59</td><td>0.59</td></tr><tr><td>CH<sub>3</sub> + HO<sub>2</sub> &rarr; CH<sub>4</sub> + O<sub>2</sub></td><td>0.0156</td><td>0.0155</td></tr></tbody></table>

This tends to show that neglecting reactions which only involve different excited states of O₂, has a negligible effect on the chemical composition of the mixture. This method of switching reactors is therefore interesting because it avoids the drift of minority species caused by the use of a plasma solver and associated plasma chemical scheme that is not
adapted to the actual state of the gas between discharges. Although negligible for a single discharge, this drift can be significant for thousands of discharges.

Figure 16 shows the evolution of mole fractions, pressure and temperature over the first 6 periods using the R1 reactor alone (solid line) and the switch reactor method R1/R2 (dashed line). There is a drift in the results between the two methods, but it can also be seen that the use of a constant volume reactor (R1) results in a significant increase in pressure and temperature. The observed drift between the two methods is therefore mainly due to this increase in temperature. In practice, however, the reactor is operated at constant pressure and temperature variations are considered negligible. Constant volume transformations only make sense for very short times when the acoustic waves do not have time to leave the reactor volume. It can therefore be concluded that the most realistic results are obtained using the reactor switch method.

## Evaluation of the heterogeneity factor from measurements

The heterogeneity factor $\beta$ is a modeling parameter that must be determined from both experimental and numerical results. Its value is determined by comparing the mass fractions obtained from the simulations with the gas chromatography measurements available at the reactor outlet. However, this process is not straightforward, as it requires parametric simulations over the entire reactor, i.e. over tens of thousands of periods, with multiple discharges during each period. Therefore, in order to limit the computational cost, it was decided to simulate only the first hundred periods and then to use an extrapolation of the mass fractions at the reactor outlet for comparison with the measurements.

In the absence of a known value for the heterogeneity factor, an iterative process was employed. This process involves starting with an arbitrary initial value for $\beta$ (first guess), running a simulation, extrapolating the results to allow comparison with experimental data, and then correcting the value of $\beta$. The process is repeated until a satisfactory agreement is achieved. Finally, it was found that the heterogeneity factor of $\beta = 0.96 (M=26)$ provides
satisfactory results for O₂, CH₄, CO₂ and CO. Only the results from the final iteration of this process are presented here.

The extrapolation of this last simulation results ($\beta = 0.96$) is based on the assumption that the averaged progress variable $\tilde{C}_k$ of the main species $k$ ($k \in \{\text{CH}_4, \text{O}_2, \text{CO}_2\}$), whose definition is based on the mass fractions, follow a power law of order $n$. Accordingly, the equation (19) writes:

$$ \frac{d\tilde{C}_k}{dt_r} = A_k \tilde{C}_k^{n_k}, \quad k \in \{\text{CH}_4, \text{O}_2, \text{CO}_2\} \quad \text{and} \quad \tilde{C}_k = |\tilde{Y}_k - \tilde{Y}_k|_{(t_r=0)} \quad (21) $$

Integration of Equation (21) gives:

$$ \tilde{C}_k = \left( (1 - n_k) A_k t_r \right)^{\frac{1}{1-n_k}} \quad (22) $$

The values of the extrapolation parameters $A_k$ and $n_k$, obtained by fitting the extrapolation curves to the simulation results for the three species (CH₄, O₂, CO₂), are collected in Table 7.

<table><caption>Table 7: Extrapolation parameters $n_k$ and $A_k$</caption><thead><tr><th colspan="3">M = 26</th></tr><tr><th></th><th>A<sub>k</sub></th><th>n<sub>k</sub></th></tr></thead><tbody><tr><td>O<sub>2</sub></td><td>0.0067</td><td>0.155</td></tr><tr><td>CH<sub>4</sub></td><td>0.0023</td><td>0.139</td></tr><tr><td>CO<sub>2</sub></td><td>0.0001</td><td>0.675</td></tr></tbody></table>

There is excellent agreement over the hundred periods between the simulation results and the extrapolation function given by equation (22), see Figure 17.
<figure><img src="image_25.png" /><figcaption>Figure 17: Extrapolation of O₂, CH₄ and CO₂ using equation (22).</figcaption></figure>

The mass fraction of CO was obtained from the carbon balance at the reactor outlet, neglecting the effect of minor carbon species on the overall balance. The numerical results underestimate the conversion of species when compared to experimental results, see Figure 18, but the shape of the evolution as a function of residence time is relevant. Studies are still needed to refine the value of this heterogeneity factor, which is very high, i.e. close to its maximum of unity. The aim of the future studies will be to carry out complete simulations of the reactor without extrapolation.
<figure><img src="image_26.png" /><figcaption>Figure 18: Extrapolation of O₂ (a), CH₄ (b), CO (c) and CO₂ (d) mass fractions for the 26 particles case. Comparison of the extrapolation with the experimental result at different residence times (5.27, 6.14, 7.37 and 9.22 ms).</figcaption></figure>

Finally, a summary of the complete simulation procedure is provided in Figure 19.

Conclusions

This study presents a comprehensive methodology for the one-dimensional statistical modeling of a practical DBD reactor, from the processing of the measured electrical signals used in simulations to the Monte Carlo method that accounts for the spatial heterogeneities of the plasma discharges. The proposed statistical approach clearly identifies the assumptions needed to simplify the real complexity of the DBD reactor to a 1D flow, laying the ground-
<figure><img src="image_27.png" /><figcaption>Figure 19: Summary of the modeling strategy.</figcaption></figure>

work for more refined future models capable of incorporating phenomena such as catalysis.

The numerical results provide a realistic description of the electrical power consumption and the current flowing through the gas in a complex case where approximately one hundred plasma discharges occur per electrical half-cycle. The combination of chromatographic measurements at the reactor outlet with numerical simulations yields the heterogeneity factor of the discharges, which is identified as a key parameter in the model. Although the flow appears stationary on average, the results indicate a highly heterogeneous spatial distribution of the discharges within the reactor.

A parametric analysis of an isolated discharge was also carried out to accurately define the characteristic times. These times are then used to optimize the computational cost and to select the most appropriate kinetic schemes for the gas phase, whether in a plasma or quasi-thermodynamic equilibrium state. This method also distinguishes between fast trans-
formations occurring at constant volume and slower transformations occurring at constant pressure.

The numerical results suggest that the gases passing through the reactor are rarely in a plasma state, highlighting the need for further research to develop kinetic schemes that are suitable not only for the chemical behavior of discharges, but also to the quasi-equilibrium states that follow. The accuracy of these post-discharge kinetic schemes is crucial to ensure accurate simulations of yield and selectivity for high-value chemical compounds.

Finally, additional analysis is required to assess the numerical biases that accumulate with each discharge, which can lead to significant drifts when attempting to simulate the entire reactor. Despite its current limitations, the numerical tool developed in this study provides valuable insights to complement experimental data.

## Acknowledgement

This work was supported by the France 2030 program (LABEX INTERACTIFS, reference ANR-11-LABX-0017-01), the Région Nouvelle-Aquitaine, the ANR (PRC grant ANR-21-CE43-0001-02) and PRCI grand ANR-23-CE05-0025. This work pertains to the French government programme "Investissements d'Avenir" (EUR INTREE, reference ANR-18-EURE-001).

## References

(1) Enger, B. C.; Lødeng, R.; Holmen, A. A review of catalytic partial oxidation of methane to synthesis gas with emphasis on reaction mechanisms over transition metal catalysts. *Applied Catalysis A: General* **2008**, *346*, 1–27.

(2) Blankenship, A.; Artsiusheuski, M.; Sushkevich, V.; van Bokhoven, J. A. Recent trends,
current challenges and future prospects for syngas-free methane partial oxidation. *Nature Catalysis* **2023**, *6*, 748–762.

(3) Rostrup-Nielsen, J. R.; Sehested, J.; Nørskov, J. K. Hydrogen and synthesis gas by steam-and C02 reforming. *Advances in Catalysis* **2002**, *65–139*.

(4) Zhang, H.; Sun, Z.; Hu, Y. H. Steam reforming of methane: Current states of catalyst design and process upgrading. *Renewable and Sustainable Energy Reviews* **2021**, *149*, 111330.

(5) Jang, W.-J.; Shim, J.-O.; Kim, H.-M.; Yoo, S.-Y.; Roh, H.-S. A review on dry reforming of methane in aspect of catalytic properties. *Catalysis Today* **2019**, *324*, 15–26.

(6) Aasberg-Petersen, K.; Christensen, T. S.; Nielsen, C. S.; Dybkjær, I. Recent developments in autothermal reforming and pre-reforming for synthesis gas production in GTL applications. *Fuel Processing Technology* **2003**, *83*, 253–261.

(7) Hickman, D.; Schmidt, L. D. Steps in CH4 oxidation on Pt and Rh surfaces: High-temperature reactor simulations. *AIChE Journal* **1993**, *39*, 1164–1177.

(8) Elbadawi, A. H.; Ge, L.; Li, Z.; Liu, S.; Wang, S.; Zhu, Z. Catalytic partial oxidation of methane to syngas: review of perovskite catalysts and membrane reactors. *Catalysis Reviews* **2021**, *63*, 1–67.

(9) Loenders, B.; Engelmann, Y.; Bogaerts, A. Plasma-catalytic partial oxidation of methane on Pt (111): a microkinetic study on the role of different plasma species. *The Journal of Physical Chemistry C* **2021**, *125*, 2966–2983.

(10) Heijkers, S.; Aghaei, M.; Bogaerts, A. Plasma-based CH4 conversion into higher hydrocarbons and H2: modeling to reveal the reaction mechanisms of different plasma sources. *The Journal of Physical Chemistry C* **2020**, *124*, 7016–7030.
(11) Pham, T. P.; Ro, K. S.; Chen, L.; Mahajan, D.; Siang, T. J.; Ashik, U.; Hayashi, J.-i.; Pham Minh, D.; Vo, D.-V. N. Microwave-assisted dry reforming of methane for syngas production: a review. *Environmental chemistry letters* **2020**, *18*, 1987–2019.

(12) Montesano, C.; Faedda, M.; Martini, L. M.; Dilecce, G.; Tosi, P. CH4 reforming with CO2 in a nanosecond pulsed discharge. The importance of the pulse sequence. *Journal of CO2 Utilization* **2021**, *49*, 101556.

(13) De Bie, C.; Van Dijk, J.; Bogaerts, A. The dominant pathways for the conversion of methane into oxygenates and syngas in an atmospheric pressure dielectric barrier discharge. *The Journal of Physical Chemistry C* **2015**, *119*, 22331–22350.

(14) Bouchoul, N.; Fourre, E.; Duarte, A.; Tanchoux, N.; Louste, C.; Batiot-Dupeyrat, C. Plasma-metal oxides coupling for CH4-CO2 transformation into syngas and/or hydrocarbons, oxygenates. *Catalysis Today* **2021**, *369*, 62–68.

(15) Bogaerts, A.; Tu, X.; Whitehead, J. C.; Centi, G.; Lefferts, L.; Guaitella, O.; Azzolina-Jury, F.; Kim, H.-H.; Murphy, A. B.; Schneider, W. F., et al. The 2020 plasma catalysis roadmap. *Journal of physics D: applied physics* **2020**, *53*, 443001.

(16) Maitre, P.-A.; Bieniek, M. S.; Kechagiopoulos, P. N. Plasma-enhanced catalysis for the upgrading of methane: A review of modelling and simulation methods. *Reaction Chemistry & Engineering* **2020**, *5*, 814–837.

(17) Hagelaar, G.; Pitchford, L. C. Solving the Boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models. *Plasma sources science and technology* **2005**, *14*, 722.

(18) Luque, A. BOLOS—BOLTZMANN EQUATION SOLVER Open Source library. 2025-02-10; https://github.com/aluque/bolos.
(19) Pancheshnyi, S.; Eismann, B.; Hagelaar, G. J. M.; Pitchford, L. C. computer code ZDPlasKin (University of Toulouse, LAPLACE, CNRS-UPS-INP, Toulouse, France, 2008). 2025-02-10; https://http://www.zdplaskin.laplace.univ-tlse.fr/.

(20) Cheng, L.; Barleon, N.; Cuenot, B.; Vermorel, O.; Bourdon, A. Plasma assisted combustion of methane-air mixtures: Validation and reduction. *Combustion and Flame* **2022**, *240*, 111990.

(21) Abdollahzadeh, M.; Pascoa, J.; Oliveira, P. Implementation of the classical plasma-fluid model for simulation of dielectric barrier discharge (DBD) actuators in OpenFOAM. *Computers & Fluids* **2016**, *128*, 77–90.

(22) Bechane, Y.; Fiorina, B. Numerical investigations of turbulent premixed flame ignition by a series of Nanosecond Repetitively Pulsed discharges. *Proceedings of the Combustion Institute* **2021**, *38*, 6575–6582.

(23) Verma, A. K.; Venkattraman, A. SOMAFOAM: An OpenFOAM based solver for continuum simulations of low-temperature plasmas. *Computer Physics Communications* **2021**, *263*, 107855.

(24) Barléon, N.; Cheng, L.; Cuenot, B.; Vermorel, O. A phenomenological model for plasma-assisted combustion with NRP discharges in methane-air mixtures: PACMIND. *Combustion and Flame* **2023**, *253*, 112794.

(25) Maitre, P.-A.; Bieniek, M. S.; Kechagiopoulos, P. N. Plasma-enhanced catalysis for the upgrading of methane: A review of modelling and simulation methods. *Reaction Chemistry & Engineering* **2020**, *5*, 814–837.

(26) De Bie, C.; Verheyde, B.; Martens, T.; van Dijk, J.; Paulussen, S.; Bogaerts, A. Fluid modeling of the conversion of methane into higher hydrocarbons in an atmospheric pressure dielectric barrier discharge. *Plasma Processes and Polymers* **2011**, *8*, 1033–1058.
(27) Snoeckx, R.; Aerts, R.; Tu, X.; Bogaerts, A. Plasma-based dry reforming: a computational study ranging from the nanoseconds to seconds time scale. *The Journal of Physical Chemistry C* **2013**, *117*, 4957–4970.

(28) Rochoux, M. C.; Cuenot, B.; Ricci, S.; Trouvé, A.; Delmotte, B.; Massart, S.; Paoli, R.; Paugam, R. Data assimilation applied to combustion. *Comptes rendus. Mécanique* **2013**, *341*, 266–276.

(29) Meldi, M.; Poux, A. A reduced order model based on Kalman filtering for sequential data assimilation of turbulent flows. *Journal of Computational Physics* **2017**, *347*, 207–234.

(30) Labahn, J. W.; Wu, H.; Coriton, B.; Frank, J. H.; Ihme, M. Data assimilation using high-speed measurements and LES to examine local extinction events in turbulent flames. *Proceedings of the Combustion Institute* **2019**, *37*, 2259–2266.

(31) Eckert, Z.; Tsolas, N.; Togai, K.; Chernukho, A.; Yetter, R. A.; Adamovich, I. V. Kinetics of plasma-assisted oxidation of highly diluted hydrocarbon mixtures excited by a repetitive nanosecond pulse discharge. *Journal of Physics D: Applied Physics* **2018**, *51*, 374002.

(32) Casey, T. A.; Han, J.; Belhi, M.; Arias, P. G.; Bisetti, F.; Im, H. G.; Chen, J.-Y. Simulations of planar non-thermal plasma assisted ignition at atmospheric pressure. *Proceedings of the Combustion Institute* **2017**, *36*, 4155–4163.

(33) Deak, N.; Bellemans, A.; Bisetti, F. Plasma-assisted ignition of methane/air and ethylene/air mixtures: Efficiency at low and high pressures. *Proceedings of the Combustion Institute* **2021**, *38*, 6551–6558.

(34) Barleon, N.; Cheng, L.; Cuenot, B.; Vermorel, O.; Bourdon, A. Investigation of the impact of NRP discharge frequency on the ignition of a lean methane-air mixture using
fully coupled plasma-combustion numerical simulations. *Proceedings of the Combustion Institute* **2022**, 5521–5530.

(35) Barléon, N.; Cheng, L.; Cuenot, B.; Vermorel, O. A phenomenological model for plasma-assisted combustion with NRP discharges in methane-air mixtures: PACMIND. *Combustion and Flame* **2023**, 253, 112794.

(36) Goodwin, D. G.; Speth, R. L.; Moffat, H. K.; Weber, B. W. Cantera: An Object-oriented Software Toolkit for Chemical Kinetics, Thermodynamics, and Transport Processes. 2021; Version 2.5.1.

(37) Allegraud, K.; Guaitella, O.; Rousseau, A. Spatio-temporal breakdown in surface DBDs: evidence of collective effect. *Journal of Physics D: Applied Physics* **2007**, 40, 7698.

(38) Manley, T. The electric characteristics of the ozonator discharge. *Transactions of the electrochemical society* **1943**, 84, 83.

(39) Frost, L.; Phelps, A. Rotational excitation and momentum transfer cross sections for electrons in H 2 and N 2 from transport coefficients. *Physical Review* **1962**, 127, 1621.

(40) Pitchford, L. C.; Alves, L. L.; Bartschat, K.; Biagi, S. F.; Bordage, M.-C.; Bray, I.; Brion, C. E.; Brunger, M. J.; Campbell, L.; Chachereau, A., et al. lxcat: An open-access, web-based platform for data needed for modeling low temperature plasmas. *Plasma Processes Polymers* **2017**, 14, 1600098.

(41) Phelps, A. Phelps database. 2025-02-10; www.lxcat.net/Phelps.

(42) Pitchford, L. C. Hayashi database. 2025-02-10; www.lxcat.net/Hayashi.

(43) Smith, G.; Golden, D.; Frenklach, M.; Moriarty, N.; Eiteneer, B.; Goldenberg, M.; Bowman, C.; Hanson, R.; Song, S.; Gardiner, W.; Lissianski, J.; Qin, Z. GRI-Mech. 2025-02-10; http://www.me.berkeley.edu/gri_mech/.
(44) Van Oijen, J.; Donini, A.; Bastiaans, R.; ten Thije Boonkkamp, J.; De Goey, L. State-of-the-art in premixed combustion modeling using flamelet generated manifolds. *Progress in Energy and Combustion Science* **2016**, 57, 30–74.

(45) Kim, Q.; Robin, V.; Mura, A.; Champion, M. Implications of laminar flame finite thickness on the structure of turbulent premixed flames. *Journal of Fluid Mechanics* **2016**, 787, 116–147.
<figure><img src="image_28.png" /><figcaption>Figure 20: TOC Graphic</figcaption></figure>