A novel two-stage kinetic model for surface DBD simulations in air

In this work, a novel 0D model for the evaluation of O3 and NO2 produced by a surface dielectric barrier discharge (SDBD) in a closed environment is presented. The model is composed by two coupled sub-models, a discharge sub-model and an afterglow one. The first one, simulating the discharge regime and consequently including electron impact reactions, aims to calculate the production rates of a set of key species (atomic oxygen, excited states of molecular oxygen and molecular nitrogen). These latter are the input of the afterglow sub-model, that simulates the afterglow regime. We introduce a methodology to relate the production rates of the above mentioned species to the input power of the SDBD reactor. The simulation results are validated by a comparison with experimental data from absorption spectroscopy. The experimental measurements are carried out as follows. First, the discharge is turned on until the NO2 number density reaches steady state. Then, the discharge is turned off for several minutes. Finally, the discharge is turned on again to observe the effects of the NO2 concentration on ozone dynamics. The entire process is done without opening the box. The system operating in all the above-listed conditions is simulated for three different levels of input power.


Introduction
Non-thermal plasmas are used in many different fields, ranging from industrial applications to aerospace and biomedical ones [1,2]. The description of chemical and physical phenomena that characterize low temperature plasmas * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
in molecular gasses such as air is not trivial. Performing experimental activities on plasma-based devices is often both costly and challenging. Furthermore, it is not infrequent that an important physical quantity cannot be measured at all. In this context, modeling activities can be very helpful to achieve a better comprehension of the physical mechanisms being studied. Historically, different methodologies have been adopted to model non-thermal plasmas [3]. These can be broadly grouped in three groups, i.e. kinetic (which the popular particle-in-cell approach) [4][5][6], fluid [7][8][9][10] and-if the transport phenomena are negligible for the given applicationsglobal models [11]. These latter, popularized by Liebermann and colleagues during the mid-90s, allow focusing on the chemical characterization of the plasma. Indeed, thanks to the 0D hypothesis, these models allow considering large numbers of chemical reactions while retaining reasonable computational times. During the years, such kind of models were adopted to investigate a lot of different phenomena and problems. Specifically, in [12] Lee and Liebermann studied the chemistry of several gasses and mixtures for different values of key parameters such as absorbed power and pressure. Global models have also extensively applied to study CO 2 . For example, Koelman et al studied the CO 2 splitting in plasma for the production of CH 4 using renewable energy sources [13]. Later on, CO 2 conversion has also been studied using a quasi-1D plug flow model [14]. In a slightly different field, Van Gaens and Bogaerts adopted the same methodology to simulate a plasma jet in an air-argon mixture, investigating the effect of gas temperature, flow speed, power density and air humidity on the chemistry [15]. Dorai and Kushner used a global model to investigate the effects of plasma treatment on wetting and adhesion properties of polymers [16]. Murakami et al [17] analyzed the effect of humid air impurity in a plasma jet with helium-oxygen mixture.
In recent years, global models have also been adopted to study multiphase domains comprising gases and liquids. Lietz and Kushner used a global model to study the interaction between a plasma and a water-covered tissue [18]. Heirman et al also combined a global model with a 2D fluid model to assess the active species production in water solution treated with a plasma jet [19].
In this work, we present a study on the kinetic mechanisms that take place in a surface dielectric barrier discharge (SDBD) in air. In particular, we investigate the production of reactive species of oxygen and nitrogen (RONS), since these species play a central role in determining the antimicrobial and antitumor properties of a non-thermal air plasma [20]. Since the reaction paths for the generation of such biologically active species are very complex and can considerably vary depending on the given physical conditions (e.g. presence of humidity), different modeling strategies can be found in literature. Indeed, one can either develop a detailed modelusually characterized by extensive reaction sets-or capture the main physical properties of the given studied system with reduced models. While the advantages of a detailed set are obvious, they come at the expense of an increased complexity of implementation and computational burden. However, advocating for reduced models is the inherent uncertainty with which reaction rates are measured. Turner has investigated this issue in [21].
As an example of a detailed model, Sakiyama et al [22] assessed the concentration of chemical species created in humid air by a SDBD. The model-comprising 624 reactions-describes both the discharge and the afterglow region by coupling two dedicated 0D-models and adopting a multiple time step procedure. The authors show the steady state value of the computed number densities of several neutral species after 1000 s of discharge application.
Conversely, Shimizu et al used a simple parametric model in [23] to assess the ozone production in a SDBD. Despite the limited number of reactions considered, the experimental validation performed by the authors showed that the ozone number density production was correctly simulated.
A parametric model was also adopted by Park et al [24] to model a similar physical configuration. Using the proposed reaction scheme, the authors of [24] were able to correctly assess the O 3 and NO 2 number density time-evolution using only 36 reaction paths. The model will be discussed in more detail in section 4.1.
In this work, we propose a novel two-stage model to study the main plasma chemistry processes of a SDBD reactor operating with atmospheric pressure air in a closed environment. The simulations performed with the developed model are validated against experimental measurements obtained by absorption spectroscopy. The discharge sub-model simulates the discharge and include electron impact reactions and ionized species. This sub-model features 103 reactions, including several kinetic channels for electronically excited molecular nitrogen states. These latter are crucial in correctly asses NO, and consequently NO 2 , number densities. The afterglow sub-model-comprising 45 reactions-represents the afterglow regime, focusing on the neutral species kinetics. The model simulates the transient behavior of several key neutral species-such as O 3 and NO 2 -over hundreds of seconds and is built around two main ideas. First, we assume that the output of the discharge sub-model is a production rate, rather than a number density. This rate is used as an input for the afterglow sub-model, and periodically updated depending on the simulated discharge conditions. Second, the measured power fed to the discharge by the plasma generator is considered as an input for the computational model. We develop an original numerical technique to accomplish this task. We show what the model reasonably captures a number of key features of the number density time-evolution of O 3 and NO 2 , measured for three different values of the electrical input power levels. The study shows that including kinetic processes involving excited states of atomic and molecular nitrogen is crucial for the assessment of the O 3 and NO 2 number densities, at least for the studied configurations.
The SDBD experimental setup is described in section 2, while the experimental results obtained for different values of the input power are reported in section 3. Subsequently, the device described in section 2 is modeled using the approach proposed by Park et al. Then, following a discussion of the obtained results and a comparison with the experimental measurements, the characteristics and the implementation details of the novel two-stage model proposed by the authors are presented in section 4.2. Then, in section 5, the device is simulated using the model proposed by the authors, for three different values of the input power. The results are again compared to the measurements and discussed in detail. Finally, conclusions are presented in section 6.

Experimental setup
To benchmark the proposed two-stage model, absorption spectroscopy has been employed to measure the O 3 and NO 2 densities, i.e. two of the key neutral species of the system, generated by a SDBD. Absorption spectroscopy is a well-established method which makes use of the Lambert-Beer's law to determine the number density of a specific gas species as follows: Lσ abs (1) where n is the number density of the absorbing species expressed in mol cm −3 , I and I 0 are the light intensities measured with and without gas-light interaction, L is the interaction length, and σ abs is the absorption coefficient. The latter can be found in literature and, in this work, the absorption coefficients for O 3 and NO 2 are taken from [25] and [26] respectively. The SDBD used in this work is a 9 × 5 cm 2 alumina dielectric surface, 1 mm thickness, with 9 metal fingers 7 cm long, 1 mm large and 4 mm apart made of an alloy of ruthenium, nickel, silver and palladium (figure 1(a)). The plasma is lit up by a sinusoidal signal, between 8 and 10 kVpp high voltage power supply working at power levels up to 150 W [27] and at a frequency of (8.5 ± 0.5) kHz [28,29].
The experimental setup consists of a SDBD positioned inside a 4-way cross vacuum chamber of about 3 l, sealed at two opposite sides by quartz windows and at the other ends by open-close valves, as shown in figure 1(b). The light sources, i.e. the UV-A (UV Philips Lighting) and UV-C (UV Lawtronics) lamps, and the respective radiometric detectors used for the light intensity measurements are installed in front of the opposite quartz windows. The absorption length was measured to be L = 20.0 ± 0.1 cm, which refers to the linear distance between the quartz windows as determined with a caliber. As a radiometric device, we employed the HD 31 datalogger from DeltaOhm coupled with two radiometric probes that collect the irradiated power per surface unit (µW cm −2 ) in two wavelength ranges: the LP471-P-A probe measures in the 315-400 nm range, and the LP471-UVC in the 220-280 nm range. The details of the lamp emissions and the radiometric probes can be found in [30]. In figure 2, we compare the absorption coefficients of O 3 [25] (shown in blue) and NO 2 [26] (in red) with the normalized sensitivities of LP471-UVC (in blue) and LP471-P-A (in red) probes, respectively. The sensitivities are represented by dotted lines, and scaled by 10 20 . Additionally, the lamp emissions of UV-C (in blue) and UV-A (in red) are depicted by dashed lines, also in normalized and scaled. Since the UV-A lamp emission is significantly stronger than that of the UV-C lamp, the LP471-UVC radiometer in combination with the UV-C lamp is sensitive only to the absorption of O 3 . On the other hand, the LP471-P-A probe in combination with the UV-A lamp is sensitive only to the absorption of NO 2 . Convoluting the wavelength dependence of both the lamp emission and the radiometer sensitivity, we can derive the following effective absorption coefficients for O 3 and NO 2 : σ abs , O 3 = 1.127 × 10 −17 cm 2 mol −1 and σ abs , NO 2 = 5.66 × 10 −19 cm 2 mol −1 . To account for potential calibration drifts in the radiometer calibration, an error of 10% was assumed in the estimation of σ abs , O 3 and σ abs , NO 2 . However, this uncertainty is only introduced via σ abs as n solely depends on intensity ratios (as shown in equation (1)). Nevertheless, a 1% error due to non-linearities, as well as a ±1 digit reading error, were included in accordance with DeltaOhm specifications. These uncertainties are then propagated to determine error bars of n O3 and n NO2 as shown in the following.
The gas-valves are used to fill the chamber with pure air (79% of N 2 and 21% of O 2 and other negligible pollutants) and evacuate it respectively before and after each experiment. Moreover, two temperature probes are installed: one in contact with the ground side of the SDBD, and the other in the middle of the chamber volume to measure the gas temperature. Herein, these measurements will be referred to as contact and gas temperature, respectively. Finally, to ensure uniformity of the gas, a vent has been placed inside the chamber.

Experimental results
For these experiments we filled the cross chamber box with pure air, and then we closed the box. While registering the light of the UV-A and the UV-C lamps, referring to figure 3, we light on the plasma at a fixed input power for 600 s (a), 240 s (b) and 230 s (c), then we switch it off, and we power on the plasma again after 750 s ((a) and (b)) and 850 s (c). In that way, we can investigate the differences occurring in the O 3 and NO 2 production during the two cycles.
Different O 3 and NO 2 regimes are identified at different input powers (figure 3).
Looking at the first plasma cycle, shown in figure 4, the O 3 slope at the ignition is steeper for increasing power levels. The maximum O 3 concentration decreases with the input power. Also, the decrease happens at lower treatment times for higher power levels. This trend is similar to the results obtained by Park et al [24]. As we will show with the simulations, the nitrogen oxidation chain is responsible for the O 3 depletion and the NO 2 fast increase.
In the second plasma cycle, the O 3 mimics the first cycle for P in = 35 W and stays under the instrument sensitivity for P in = 100 W and 120 W, where NO 2 concentration is above 0.5 × 10 16 cm −3 . The NO 2 continues to slowly increase if O 3 is high, while it decreases if the O 3 is null.
The contact temperature, that is the temperature of the dielectric of the SDBD, it is strictly dependent on the power, up to 60 • C for P in = 12 W, 90 • C for P in = 35 W, 190 • C for P in = 100 W and 220 • C for P in = 120 W. The gas temperature slightly increases from the initial value and stabilizes in a short time, reaching a value always below 40 • C.

Simulations
First, we performed a first assessment of the described configuration using the parametric model proposed in [24] by Park et al. The model will be briefly discussed in section 4.1. After that, to achieve a better comprehension of the phenomenon, we develop and implement a novel two-stage model, containing an extended set of reactions involving a larger number of chemical species. A description is provided in section 4.2.

Preliminary assessment
Park and colleagues developed a parametric model (Park model from here onward) to study a SDBD configuration similar to the one described in the previous sections. The Park model, which involves only neutral species, has been described and validated against experiments by the original authors in [24]. Despite the Park model cannot be used to predict the temporal behavior of the number densities of the considered chemical species, due to its inherent simplicity and the limited number of reactions (36), it can be a very useful tool to individuate the key mechanisms governing the phenomenon under investigation. These reactions-which are included as a subset of the two-stage model developed by the authors that will be discussed later on-are reported in table 2, marked by a 'P' letter. The Park model uses 5 fitting parameters, i.e. the number densities of the species O, N 2 (A) and O 2 (a), the steady state vibrational temperature of the vibrationally excited nitrogen molecules T v0 and the time constant of the vibrational temperature τ v . The time constant τ v is used to calculate the evolution of T v , as in equation (2): where the gas temperature T g is considered constant at 300 K. The time-evolution of the nitrogen vibrational temperature governs the number density of the vibrationally excited nitrogen molecules according to the following relations: where ∆ϵ v is the vibrational energy for an harmonic oscillator (0.29 eV), k B is the Boltzmann constant, e is the elementary charge. The Park model assumes the number densities of O, N 2 (A) and O 2 (a) species, n O , n N2(A) , n O2(a) , respectively, to be constant in time and-as aforementioned-does not account for electron impact reactions. Moreover, the power supplied to the reactor is not taken into account, and only the lower energy excited states of molecular nitrogen and oxygen (O 2 (a) and N 2 (A)) are considered.
We implemented the Park model using the ZDPlaskin Fortran module [31], and used it to simulate several chemical species during the first ignition of the discharge as a preliminary assessment of the physical phenomenon under investigation. The initial number densities are the ones typical of dry air, i.e. n N2 = 1.96 × 10 19 cm −3 , n O2 = 5.2 × 10 18 cm −3 . The other number densities, not considering of course the fitting parameters, are set to 1 × 10 4 cm −3 . Note that, changing the latter value between 1 × 10 2 cm −3 and 1 × 10 6 cm −3 does not change the simulation results. The obtained results are shown in figure 5, together with the experimental measurements performed for an input power of 120 W. The numerical values of the model parameters used to fit the experimental results are reported in the figure caption. As it can be observed, the Park We observed that the value of τ v strongly affects the time it takes n O3 to reach its maximum value, while T v0 (and the steady state value of the number density of the vibrationally excited nitrogen) mainly affects the slope of the O 3 number density curve in the descending phase and the steady state value of the NO 2 number density.  The production rates evaluated as the algebraic sum of all the reaction rates for a given species of O, N 2 (A) and O 2 (a) yielded by the Park model are shown in figure 6. Note thatas discussed-their number densities are kept constant by hypothesis. Practically, this implies that there is an availability of these species allowing the number densities to remain constant, even when they are consumed by chemical reactions. The rates obtained with this kind of hypothesis will be compared to the ones obtained using the proposed two-stage model in section 5.

Two-stage model description
A novel two-stage model has been developed and implemented. The first stage, hereinafter referred to as the discharge sub-model, has the purpose of simulating the kinetics induced by hot electrons in the SDBD. The second stage, called the afterglow sub-model, represents the chemical kinetics of the species produced by the SDBD occurring in the test chamber. The reaction set adopted for the discharge model is reported in table 1 and is constituted by an extension of the set used in [39], with the addition of the excitation and de-excitation mechanisms of N 2 (a'), N 2 (a), N 2 (w), N 2 (W), N 2 (B') taken from [45]. All the mentioned excited species of molecular nitrogen play a central role, as reported in [52], to correctly assess the NO number density and, consequently, the ones of O 3 and NO 2 . The cross-sections used in the simulations are retrieved from the LXCat database [58][59][60].
It is worth noticing that we did not include the spontaneous de-excitation of N 2 (W), N 2 (B') and N 2 (w). Indeed, the excited state N 2 (W) spontaneously de-excitates through reactions N 2 (W) -> N 2 , characterized by a reaction rate of n N2(W) 1.5 × 10 −1 , that is 20 times lower than the reaction rate of process 13b in the considered conditions (with n O ≈ 1 × 10 11 cm −3 ). Including such reaction would not affect significantly the results. Regarding the excited species N 2 (w) and N 2 (B'), they spontaneously de-excitate through reactions N 2 (w) -> N 2 (a') and N 2 (B')-> N 2 (A'). The first one causes the creation of the excited state N 2 (a'), that still feed the formation of NO via reaction 16b, with the same reaction rate of channel 17b, producing NO from N 2 (w). As a result, accounting for the intermediate state N 2 (a') does not change significantly NO formation rate. A similar discussion can be done regarding the second cited reaction, that forms the excited state N 2 (A). Nevertheless, the situation is more complex in this second case since N 2 (A) contributes in more processes compared to N 2 (a'). Not including the spontaneous de-excitation of the N 2 (B') can bring to a over-estimation of this excited state.
However, this can compensate (regarding NO production) the under estimation of the vibrational state of N 2 , explained in section 5.5. The fact that we manage to fit the experimental data for different power inputs can be a proof that the made assumption are not inconsistent.
The aim of the discharge sub-model is to estimate production rates of a set of key species (O, O 2 (a), O 2 (b), N 2 (A), N 2 (B), N 2 (C), N 2 (a'), N 2 (a), N 2 (W), N 2 (w), N 2 (B')), produced by the discharge. As will be further detailed, these will then be used as an input for the afterglow sub-model.
The afterglow sub-model is constituted by the reaction set reported in table 2. The number density of the vibrationally excited nitrogen molecules is still calculated using equation (3) as in [23]. For what concerns O 2 (v) number density, reaction 1 is only considered to correctly calculated the electron energy distribution function (EEDF), as done in [39]. Considering different studies regarding NO formation, e.g. [61], and other concerning the kinetic modeling of an SDBD [22,24], it seems reasonable to assume that O 2 (v) does not play an important role in the mechanisms leading to the formation of nitrogen oxides.
When one of the number densities n s in the afterglow satisfies the relation ns−ns0 ns+ns0 > 10%, where n s0 is the initial number density in the afterglow, the discharge model is run again starting from the number densities yielded by the afterglow model. Then, the production rates obtained from the discharge model are used to update the afterglow input. This procedure is iteratively repeated.
The rate coefficients for the electron-collision processes, which depend on the EEDF, are calculated using the Boltzmann solver BOLSIG+ [62]. The local field   [40] 56 [45] 80 N 2 (a ′ ) + M −→ N 2 (a) + M    [51] [45] (Continued.)  [50] 45b ) 0.5 [48] approximation (LFA) is used, in that rate coefficients are expressed as a function of the local electric field. In a 0D plasma simulation, the approach used to model the physical effects of electron impact processes is a critical point to correctly assess the generation of chemical species. Several approaches are reported in literature, developed to suit the operational conditions of different plasma devices. As an example, in [63] the plasma generated by a radio-frequency source is simulated by imposing an electron temperature such that the electron number density is constant over time. A different approach is used by Vagapov and Schweigert in [64], where a periodic behavior of the electron temperature T e is assumed to account for the effects of streamers nearby a surface in a plasma jet device. Each period of T e is divided into four parts: prior to the streamer onset, the electron temperature is kept at 300 K, corresponding to the gas temperature. Then, T e is rapidly raised to 3 eV. This value is maintained constant for a short time (3-5 µs) after the streamer onset, before returning again at gas temperature. Adopting a different strategy, Heirman and colleagues defined a relation between the electrical input power and the electric field [19], which, in turn, is used to compute the rate coefficients of electron impact reactions under the LFA. The electric field-power relation used in [19] reads as: where σ is the electrical conductivity of the plasma. All the listed approaches have been tested within the two-stage model that we introduced in this work. However, none of them yielded results in agreement with the measurements performed on the SDBD considered in this study. Consequently, we implemented a novel method to describe the physical behavior of the electrons with respect to the discharge, that is based on considering the discharge itself as a collection of streamers. The effects of a single streamer are simulated by applying a reduced electric field with an exponentially decaying behavior: where E 0 and τ E are the maximum value of the reduced electric field and the time-constant of the field decrease over time, respectively. The electric field E affects both the peak values of produced number densities over time and the ratios between the different species. The peak electric field E 0 and the time constant τ E have been chosen to both fit the experimental measurements and to obtain a streamer duration and a maximum electron number density that are typical of a real streamer. In our simulations t str ≈ 100 ns and n e,max ≈ 1 × 10 12 cm −3 . E 0 results to be slightly higher than the value allowing the ignition of the streamer for the conditions considered in this study. As a consequence, n e starts decreasing as soon as the electric field value reaches the 90% of E 0 . The streamer is considered extinguished when the electron number density returns to its initial value n e0 . In order to evaluate the number densities of the species which are delivered per unit time (production rates) to the afterglow sub-model by the discharge phase, we assume that the discharge is constituted by a number of streamers, responsible for the electron impact reactions in table 1. The number of particles of a given species s created per unit time within the discharge region due to a streamer (Ṅ s ) can be expressed as: where n s and n s0 are the number densities of the species s after and before the streamer, respectively. V str is the volume of a single streamer, andṄ str is the number of streamers taking place per unit time.
We then proceed to express the number of streamers per unit time as a function of the input electric power (P in ) of the SDBD reactor and the energy density of a streamer (w str ). The resulting expression is:Ṅ In (7), the efficiency η has been introduced as an additional parameter. The efficiency has to be intended as a global efficiency of the system in transferring the input electric power to the afterglow volume in the test chamber. The efficiency is considered as a fitting parameter accounting for all the loss mechanisms affecting the process, including losses in electric supply system, dielectric losses, fluxes due to the three-dimensional nature of the system under investigation. Combining equations (6) and (7), and dividing for the volume of the afterglow region V aft , we obtain the production rate of the given species s: V aft has been chosen according to the dimension of the test chamber employed in the experimental procedure, i.e. 3 dm 3 . Note that, the energy density of a streamer ω str is related to the chosen electric field E(t) through the relation where t str ≈ 100 ns in our simulations. Equation (9) bring to calculated value of ω str ≈ 1.4 × 10 3 J m −3 . During the discharge phase, an output time-step of 1 × 10 −11 s was used to perform the time-integration of the number densities (note that the employed implicit stiff-ODE solver DVODE, integrated in ZDPlaskin, employs further sub-steps to grant the userrequired accuracy [65]). During the afterglow simulation, the time step is adapted according to the relative variation of the (neutral) species number densities. Overall, the described procedure allows simulating hundreds of seconds of device operation in few minutes of CPU time. Summing up, the described two-stage algorithm-as well as the mentioned time-step criteria-can be visualized in the flow chart in figure 7.

Model validation and discussion
In this section, we describe how the developed two-stage model is applied to the assessment of the SDBD configuration in section 2. We then proceed to discuss the results in comparison with the experimental measurements.
We simulated the device for an input power of P in = 120 W. The values of E 0 = 223.7 Td and τ E = 1E − 6 s-selected following the criteria described in section 4.2-have been employed in the simulation. The values of the other parameters are T v0 = 4340 K, τ v = 17 s and η = 9.1%.
The developed two-stage model yields results that are in reasonable agreement with the experimental data from absorption spectroscopy. The results referring to the first ignition of the discharge are reported in figure 8. A discussion on the obtained results is provided in the following section.

Initial discharge ignition
For the sake of clarity, we subdivide the simulation of the SDBD device operation in three distinct phases. In this paragraph we describe phase one, i.e. the first discharge ignition taking place during the first ∼230 s, highlighted in yellow in figure 3(c).
With regard to figure 8, as soon as the discharge is turned on, the O 3 number density starts to rapidly increase. As well known [23,24], this behavior is due to reaction 2b . As one can see, the O 3 growth carries for the first 35 s of phase one, when a maximum is reached. Then, a rapid decrease of O 3 number density can be observed. As a result, n O3 falls below 1 × 10 15 cm −3 after 110 s from the discharge onset.
The described O 3 decrease is mainly due to the intense N 2 (v) production driven by equations (2) and (3). Indeed, N 2 (v) produces NO through reaction 18b . Note that reaction 18b is already active from the discharge ignition, but its rate becomes relevant (and dominant) only after 30 s due to the exponential increase of T v and N 2 (v) in equations (2) and (3).
Concerning the temporal behavior of O 3 , a key role is played by the excited states of atomic and molecular nitrogen. Indeed, these contribute to inhibit the O 3 production in two different ways. First, reaction 45b (N(D) + O 2 −→ NO + O(D)) and 10b-17b produces NO which, as discussed, reacts with O 3 . Second, reactions 10b-17b reduce the number density of O which, in turn, is a reactant of reaction 2b, constituting the main kinetic channel for the O 3 production.

Comparison with the Park model
First of all, as already anticipated, the number densities of O, N 2 (A) and O 2 (a) are not constant in the proposed twostage model. This constitutes one of the major differences with respect to the Park model, where these three quantities are held constant throughout the simulation. The time-evolution of n O , n N2(A) and n O2(a) is shown in figure 9. Note that these values refers to the previously described results in figure 8, performed for P in = 120 W. As it can be verified, the three obtained number densities are not constant over time, and exhibit variations of several orders of magnitude over the course of the simulation.
In particular, the N 2 (A) number density is much lower than the (constant) one enforced in the Park model, see section 4.1, fitting the same experimental data set. Consequently, additional pathways other than reaction 10b (namely 11b-17b, 45b) must be considered to correctly describe the NO and NO 2 number density production mechanisms.
Finally, coherently with the observations on the number densities, we also note that the obtained rates for O, O 2 (a) and N 2 (A)-reported in figure 10-are markedly different with respect to the ones yielded by the Park model in figure 6.
For what concerns the Park model, the O 2 (a) rate is practically fixed by the O 3 concentration through reaction 31b The same applies to the production rates of N 2 (A) and O.
Differently, in the proposed two-stage model the number densities of the considered species decrease when these react with other species. This leads to the observed smaller changes over time of the rates.
Finally, the magnitude of the N 2 (A) rate yielded by our model is significantly lower than the one of the Park model. This is due to the considerably lower values of N 2 (A) number density (see the above discussion on the inclusion of N 2 excited states).

Second discharge ignition
In the following, we describe the second and third phases of the simulation. The second phase starts at the end of the first ignition, described in the previous paragraph (∼230 s). Then, the  third phase starts after ∼1100 s, corresponding to the second ignition of the discharge (highlighted in yellow in figure 3(c)).
At the beginning of the second phase (discharge off), the ozone number density has already dropped below 1 × 10 15 cm −3 . Also, turning off the discharge implies the cessation of the electron impact reactions leading to the production of atomic O within the discharge volume. As a result, the NO 2 number density remains nearly constant over time.
When the discharge is turned on again (third phase), the NO 2 number density starts to decrease. In fact, residual NO 2 molecules prevent O 3 formation by both quenching O 3 through reaction 29b and reacting with O through reactions 6b-8b. After some seconds, however, reaction 18b is activated due to the N 2 (v) raise described by equations (2) and (3). Once again, N 2 (v) triggers the O 3 depletion observed after 35 s from the new discharge onset, with the same mechanism discussed for phase one. Note that this mechanism explains the different behaviors of O 3 number density during the third phase at P in = 120 W and P in = 35 W in figure 3. In fact, we observed that, in the experimental conditions considered in this study  conditions, a number density of NO 2 ≈ 8 × 10 15 cm −3 can prevent the ozone formation through the mechanism that has been explained. The described time-evolution of the number densities of O 3 , NO 2 and N 2 (v) is reported in figure 11.

Simulation of two different input power levels
As a further validation of the model, the developed reaction scheme has been used to simulate two other experiments, performed with lower input power levels, i.e. 100 W and 35 W. The adopted parameters are reported in table 3, below the corresponding values for 120 W. Note that only the maximum electric field value E 0 , the time constant τ v and the steady state value of the vibrational temperature T v0 have been changed (following the same criteria described in section 4.2).
The obtained number densities during the first discharge ignition (phases one) are compared to the corresponding measurements in figures 12 and 13, for P in =100 W and P in = 35 W, respectively. The results show a good quantitative agreement for P in =100 W. The comparison is also satisfactory for O 3 at 35 W.     The proposed model manages to correctly represents crucial parts of the O 3 and NO 2 dynamic, such as the peak of ozone number density, the decreasing phase of O 3 and the steady state value of the NO 2 number density. This fact demonstrates that the proposed approach is not tied to a specific operational condition of the modeled device, and can provide physical insight in a variety of different conditions without major adjustments of computational parameters, which values are reported in table 3.

Limitations of the model
Uncertainties related to the measurements, the chosen reaction set and the corresponding rate coefficient are common in all the papers that deal with plasma modeling, plasma simulations and experimental validation. However, in this section we are going to report the limitations that specifically affect the presented two-stage model, i.e. the estimation of the vibrationally excited states of nitrogen molecules. During the second ignition (phase three), the largest difference between the predicted and measured O 3 number density is obtained for the largest input power, i.e. 120 W. This is probably due to the fact that equations (2) and (3) cannot properly describe the evolution of N 2 v, particularly for what concerns the first few seconds of simulation. Indeed, in reality, some of the vibrational levels that contribute to NO formation-in particular the levels between v = 13 and v = 15-are populated due to electron impact collisions [33]. As a result, they are rapidly created in the discharge region as soon as the SDBD reactor is turned on. This is disregarded when using equations (2) and (3), yielding values of the vibrational nitrogen number density below 1 × 10 6 cm −3 for t < 2 s after the discharge onset. This effect is more relevant when larger input power levels are considered, since-according to equation (6)-greater power levels correspond to larger production rates. The same mechanism explains the systematic offset between the predicted and measured NO 2 number density observed at the beginning of the first discharge ignition (phase one), see figures 8, 12 and 13 for reference. Moreover, according to data in literature [23], the N 2v distribution function does not follow a equilibrium distribution (contrary to what is postulated in (2) and (3)). As a consequence, the exponential behavior adopted to describe such number densities (equation (3)) could not be fully correct.
It is worth noticing that the vibrational temperatures obtained from the Park approach and the model presented by the authors in this paper result to be quite high. This feature has been observed also in [23]. Indeed, in the considered operative conditions, the vibrational excited states do not follow an equilibrium distribution. We adopted the equation reported in [66] to estimate the vibrational temperature of the lower levels of N 2 (v): where E 1 = 0.29 eV is the value of the energy of the first vibrational level, T low,v is the temperature of the lowest levels, ∆E = 14.32 cm −1 is the anharmonism value for nitrogen molecules [67] and T gas is the gas temperature. The obtained temperature is T low,v = 2791 K and, as shown in table 3, it is considerably lower than the ones utilized in our calculation. However, as observed in [66], the higher vibrational levels, which participate in reaction 18b, tend to acquire a vibrational temperature which is much higher than the one of the lower levels. Therefore, the vibrational temperature results showed in table 3, can be interpreted as an indication of a nonequilibrium regime in the vibrational level distribution. One way to overcome these limitations could be to implement a detailed set of reaction, including vibrational excitation and de-excitation processes. However, this would also probably reduce the practicality of a relatively simple kinetic model as the one proposed in the work; in addition, the low accuracy with which the reaction rates of such kinetic processes are often known experimentally would compromise the accuracy of the model.

Conclusions
In this work, we presented a novel 0D two-stage model to simulate the plasma chemistry of an SDBD reactor operating at atmospheric pressure with dry air. The afterglow submodel is developed as an extension of a validated model previously presented by Parket et al in [24]. Both models have been applied to the simulation of a sinusoidal SDBD reactor operating with dry air at atmospheric pressure in a closed environment. Absorption spectroscopy measurements have been performed on the SDBD reactor to measure the evolution of O 3 and NO 2 over ∼1 × 10 3 s. The discharge is extinguished after several hundreds of seconds from the ignition (phase one). Then, a second ignition is performed (phase three) after a cool down time (phase two).Both models have been used to simulate the first ignition of the SDBD reactor. The results and the underlying assumptions of the two models have been compared. The proposed two-stage model was also used to simulate the device during phases two and three (cool down and second ignition). The obtained results are in qualitative good agreement with the measurements for both considered phases.
The novel developed two-stage model, contrary to the Park one, allows to account self-consistently of the SDBD reactor input power and implements a physical description of the discharge. These additions consent to represent the SDBD reactor and its kinetics accurately for different supply conditions. In particular, in this work are reported simulations referring to three different levels of the input power i.e. 35 W, 100 W and 120 W.
The results show that the N 2 (v) number density is likely not fully represented by an exponential function (made under equilibrium hypothesis). Nevertheless, the obtained good agreement with the experimental data shows that such a simplified description of the phenomenon may still provide useful results. The study also highlights that the excited states of atomic and molecular nitrogen play a central role in correctly assessing the O 3 and NO 2 number densities in the studied configurations, and should be included in models of this kind.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).