Two-temperature balance equations implementation, numerical validation and application to H2O–He microwave induced plasmas

Global Models are widely used to study reaction kinetics in low-temperature plasma discharges. The governing conservative equations are simplified into a system of ordinary differential equations in order to provide computationally feasible conditions to study complex chemistries with hundreds of species and thousands of reactions. This paper presents a detailed two-temperature global model for a H2O–He mixture. The model developed in this work uses a statistical thermodynamics approach to solve the heavy particles energy equation self-consistently together with the electron energy and particles balance equations in order to improve the description of reactive plasma environments. Three analytical test cases are presented to validate and demonstrate the capability of this newly developed functionality embedded in PLASIMO software’s [1] global model module. The developed H2O–He models are compared with the reported results for a radio-frequency plasma [2] and then with experimental measured electron densities and gas temperature for a microwave induced plasma. In addition, conversion and energy efficiencies of hydrogen and hydrogen peroxide productions are compared with experimental values (only for hydrogen) for a pure H2O microwave induced plasma and with available literature results. This comparison underlines the challenges toward finding an optimal plasma configuration and conditions for production of hydrogen from water. The three analytical test cases for validation of the gas-temperature balance implementation in the PLASIMO global model and the detailed developed H2O–He model can be used as benchmarks for other global plasma models. The PLASIMO input files for the presented H2O–He model are available as supplementary materials (https://stacks.iop.org/PSST/30/075007/mmedia); for any future update, please consult the PLASIMO website, https://plasimo.phys.tue.nl.


Introduction
Earth global warming due to emission of anthropogenic green house gases becomes a threatening problem for mankind. Producing hydrocarbon fuels like CH 4 from reduction of CO 2 is a promising path to solve at least one part of this critical problem by helping to achieve CO 2 neutrality. In order to produce hydrocarbons, hydrogen and carbon sources are needed. The conversion of CO 2 to CH 4 can be done with the help of microwave plasma and catalytic processes. The main motivation to make a global model for H 2 O and He in this paper is to make a step in understanding underlying mechanisms in H 2 O-CO 2 microwave induced plasmas and production of H 2 from H 2 O, which was previously reported with quite high energy efficiencies in the order of 35%-40% [3].
A lot of research has been already carried out to understand the governing chemistry of CO 2 in both microwave and dielectric barrier discharges [4][5][6][7][8], but to the best of our knowledge, plasma modeling research on H 2 O in microwave induced discharges is still lacking. H 2 O is mainly considered in plasma discharges that have applications in biotechnology at atmospheric pressure. Liu et al [2] studied the He-H 2 O discharge in atmospheric pressure by investigation into the effect of H 2 O concentration on a He discharge with a global model that contains 46 species and 577 reactions. Vasko et al [9] compared their global model and 1D fluid model of H 2 O-He with experimental results and the global model of [2]. They used the chemistry from Liu et al [2] and improved it by adding some vibrational states of H 2 O. Ding et al [10] focused on atmospheric He-H 2 O discharges by using a hybrid two-temperature global model, and in another paper [11], they determined the main reaction pathways in H 2 O-He for atmospheric pressure capacitive discharges. van Gaens et al [12] studied the kinetics of an argon plasma jet in atmospheric pressure humid air. They considered a subset of water chemistry that is provided in [2] while adding some of the vibrational states of H 2 O molecule to their species list. Liu et al [13] did research on an Ar-H 2 O mixture in cold atmospheric radiofrequency plasma as a cheaper alternative to He-H 2 O plasma for production of reactive oxygen species due to the scarcity of He as a resource. They used a global model with 57 species and 1228 reactions. The main chemical pathways were recognized for some important species like O, OH, OH(A), and H 2 O 2 . Recently, Luo et al [14] investigated the plasma kinetics in a nanosecond pulsed filamentary discharge in atmospheric Ar-H 2 O and pure H 2 O plasmas. They studied the main mechanisms in the production of H 2 O 2 , which is one of the most important species in biomedical applications of the plasmas.
Production of hydrogen and hydrogen peroxide from water by microwave induced plasma with energy efficiencies around 35%-40% and 45%, respectively, were reported in several works in the eighties [3,15,16] but no attempt was made to reproduce such results yet. Fast thermal quenching is counted as one of the most important factors for improving energy efficiency and conversion to reach optimal operational conditions of a microwave reactor.
In this paper, the main goal is to provide a predictive, proper (two-temperature medium), and qualitative (zerodimensional) model to study microwave induced plasmas in H 2 O. For the development of such a complex model, a systematic approach is undertaken by starting with a mixture of water and helium. Therefore, for the first steps in verification and validation of the model, the available and previously validated kinetic models from the literature for cold atmospheric pressure He/H 2 O plasmas can be used. The global model for He-H 2 O mixture from Liu et al [2] is used for the implementation and verification of the model presented here. Then, the He-H 2 O plasma chemistry is refined and the global model is extended to the two-temperature media of microwave plasmas where the gas temperature is self-consistently calculated in each step of the model instead of using the classical constant gas temperature value assumption. The newly developed functionality for self-consistent calculation of the heavy particles gas temperature is explained and validated with three analytical test cases. This extension is particularly important not only for microwave plasmas but also for any application of low-temperature plasmas that are sensitive to variation of gas temperature, as the latter can strongly affect the plasma chemistry and more particularly neutral species kinetics. In a later stage, the predictability of the model is assessed by comparing the electron density and gas temperature values in a He + 11% H 2 O microwave induced discharge at 500 mbar over a range of power 300-1000 W with measured values in a microwave plasma reactor. In addition, conversion of water to hydrogen and hydrogen peroxide in a pure microwave induced water plasma discharge and their energy efficiencies are calculated by modeling the plasma afterglow and compared with measured experimental values (only for production of hydrogen).

Global model description
A global model is also known as a zero-dimensional model because only temporal variations of variables are taken into account without considering any spatial resolution. The plasma properties are either solved locally or spatially integrated over the plasma volume. The modeling of species densities variations versus time (or their steady state values) is the main goal of these models. Therefore, the continuity equation in the particle density form is solved for each species involved in the model. Transport terms are modeled as frequencies. The heavy particles temperature (gas temperature) and electron temperature can be set as constant values, or energy equations can be solved to have better approximation of these two temperatures. The gas temperature is often kept constant in models of low temperature plasmas, but in microwave plasmas, gas heating is not negligible. Therefore, in order to have a better representation of the plasma, the heavy particles energy can be added to the system of ordinary differential equations (ODEs) such that the gas temperature is updated at each iteration. The momentum equation is usually not directly resolved in global models. The density of each species and their temperatures are assumed to be volume-averaged (or local) variables; it means they are homogeneous in the considered volume of the plasma in global models. In the following sections, the balance equations that are solved by the PLASIMO global model are discussed in detail.

Species density balance
The simplified particle balance for a species i in a global model is derived from the particle density balance equation where n i is the density of species i, and ω i is the source term that describes the net production of that species in chemical reactions. The vector u i is the velocity of the species, and it is defined as the sum of the bulk mass-average velocity u b and the species' diffusive velocity u d b i . After integrating equation (1) over the volume, for a cylindrical shape plasma with volume V, side area surface S side , and cross area A plasma , the particle balance equation can be written as where u b , u b s , u d b i , and u d b is are the bulk particle-average velocity normal to the cross section of the plasma, the bulk particleaverage velocity normal to the side area of the plasma, the diffusive velocity of species i normal to the cross section of the plasma, and the diffusive velocity of species i normal to the side area of the plasma, respectively. After dividing by the volume we get

Electron energy balance
The electron energy balance is added to the system of ODEs if the electron temperature is not assumed constant. The complete form of electron energy equation and detailed derivation can be found in [17, p 21-24]. There are a couple of assumptions in the global model. First it is assumed that the kinetic energy due to the bulk velocity and diffusive velocity are negligible compared to the thermal energy of electrons. The work done by shear stress tensor, hydrostatic pressure, and gravitation also are negligible. Usually in low temperature plasmas and in the case of absence of an external magnetic field, the terms related to the magnetic field can be set to zero.
Other terms (heat conduction, collision terms, and ohmic heating n e q e E · u e where n e , q e , E, u e represent electron particle density, electron charge, electric field, and electron velocity, respectively) should be modeled as frequencies, if they are important. Therefore, the electron energy balance is written as where P in represents the power density that electrons absorb from the electric field (n e q e E · u e ). For a volume averaged model, it can be replaced with P/V, the ratio between the total absorbed power and the plasma volume. The collision term is divided into two different contributions, elastic and inelastic collision processes. The averaged transferred energy per collision, ε el , can be calculated for any arbitrary elastic collision between electron and a species X, e + X → e + X, as [18, pp 47-51] ε el = 2m eX m e + m X with m eX the reduced mass. Because the electron mass is so small compared to the heavy particles, the above equation can be approximated as Therefore, for an elastic collision e + X → e + X, the elastic energy source of the electron energy equation is where k el stands for the rate coefficient of the elastic collision. The rate coefficient can be calculated with the help of the electron energy distribution function (EEDF) and cross section of the corresponding elastic collision.
In electron inelastic collisions, the internal energy state of the species X is changed or the latter is converted into other species. For example, the species X is vibrationally and/or electronically excited, ionized or dissociated into other species. The source (that can be positive or negative) of the electron energy for these types of collisions depends on the process. If the reaction is endothermic and the electron is on the left-hand side of the reaction, it means that the electron collides with the heavy particle/s and looses energy: this collision is a sink for the electron energy balance. This loss term can be calculated for any arbitrary reaction as follows where ΔE is defined as In equations (7) and (8) • The summation is done on each side of the reaction over all species involved on that side, X i , • α e j and β e j are the stoichiometric coefficients for the electrons on the left and right hand sides of reaction j, respectively, • α X i j and β X i j are the stoichiometric coefficients for the species X on the left and right hand sides of the reaction j, respectively, • E i is the formation energy of the species X with respect to the zero energy level reference of the system, • R inel is the rate of the inelastic collision.
Based on the above definition, the stoichiometry of species X or electron for reaction j can be defined as the index i can be used both for electrons and heavy particles. Note that ΔE is positive in (8), but because the reaction is endothermic, the negative sign appears in the calculation of the source term for the electron energy (equation (7)). There can be an exception in calculation of endothermic source term for electron dissociation processes. It is possible that the necessary energy which electron should have to overcome the energy barrier of the reaction is higher than the ΔE. In that case, this excess of energy should be added to ΔE in order to have correct estimation of the sink. For example in electron dissociation of CO 2 , an electron needs about 10 eV to dissociate CO 2 to CO and O. However, the ΔE is 5.5 eV, in order to take into account correctly the electron energy sink of the process another 4.5 eV should be added to this ΔE.
For exothermic reactions where an electron is on the righthand side, the inelastic source can be calculated as equation (7). The negative sign should be kept because ΔE is negative in these cases, and electrons gain energy in these types of processes (the source term for electron energy). The explained approach for the inelastic processes' source term calculation is the one used in the PLASIMO chemical source term calculation module.
In general, the rate of electron impact reactions can be obtained either from experiments or be calculated based on cross sections (σ) of reactions and an EEDF. In the case that cross sections for given electron impact processes are available, the electron collision rate with species X can be calculated as R| e−X = n e n X σv , with n X the particle density of species X and where σ(v) and f (v) are the cross section and an electron distribution function in the velocity domain. In other words σv is the rate coefficient of a collision process of electrons and species X. The EEDF is calculated using the two-term approximation [19, pp 76-92] by BOLSIG+ [20] in this paper. After providing cross sections for electron impact collisions and specifying a range for E/N for BOLSIG+, it returns the calculated rate coefficients for each process in a tabular form as a function of the mean electron energy. For all presented results in this paper, the ionization degree, the gas temperature, and the chemical composition for a chosen range of E/N in BOLSIG+ input are set in such a way that steady state conditions of the under studying plasma are represented properly.

Heavy particles energy balance
If the gas temperature is treated as a variable parameter, the heavy particles energy balance is added to the system of ODEs to update the gas temperature in each iteration of the model. The heavy particles energy is derived from the conservative form of the total heavy particles energy (bulk mean kinetic energy, thermal energy, and chemical energy) equation, for more details [17, pp 22-24]. Couple of assumptions are made to reach the simplified global model version of the heavy particles energy. It is assumed the work due to shear stress, electromagnetic, gravitation are negligible in comparison with source terms of elastic and heavy particles collisions. The hydrostatic pressure work is combined with the advection term. In addition, the bulk kinetic energy is negligible in comparison with U m , which represents the contributions of thermal and chemical energies for heavy particles. All electron related terms that appear in the conservative heavy particles equation for twotemperature media are assumed negligible [17, p 24]. Transport and heat transfer are modeled as frequencies similar to the particle balance equation. The statistical thermodynamic approach is chosen in our model to calculate U m for heavy particles. Species' formation energies are not considered on the lhs of (12), so we have a source term due to heavy particle collisions on the rhs Therefore, the heavy particles energy balance can be written as where U m and H m stand for energy and enthalpy of heavy particles, respectively. U m and H m are in J m −3 , and enthalpy can be calculated as U m + p, where p stands for the hydrostatic pressure of a mixture. In order to take into account the heat flux due to the inter-diffusion phenomena [17, p 23], [21, p 320], [22], and [23, p 589 & p 785], H i is defined to represent the enthalpy of the species i. Please consider that equation (12) is an oversimplified version of conservative form of total heavy particles energy equation [17, p 24], and various physical phenomena are captured only approximately in a global model. The elastic source term has the same value as its corresponding term in the electron energy equation (3) with the opposite sign. S ch stands for the energy source (that can have either a positive or a negative value) from heavy particles collisions. Because species' formation energies are not considered on the lhs of equation (12), S ch should be taken into account properly in order to have a correct representation of chemical energy transformations in a mixture. If a reaction is endothermic and no electron is on the lhs of the reaction, the reaction energy difference ΔE becomes a sink term for the heavy particles energy. If a reaction is exothermic and no electron is on the rhs of the reaction, −ΔE becomes a source term for the heavy particles energy. For both cases the chemical source term due to an arbitrary reaction j with the above specifications can be calculated as with R j the rate of reaction j. There is an exception again due to electron dissociative reactions when an electron needs more energy than the ΔE to overcome the reaction barrier. This extra energy is multiplied by the corresponding rate of the electron dissociative reaction and added as a source to the heavy particles energy equation. Finally, S h represents additional heat transfer mechanisms, such as conductive heat losses.
In each iteration of the model, the right-hand side of the heavy particles energy balance (12) is passed to an ODE solver and as a result new value for U m is obtained. Then, the heavy particles temperature T g is obtained from the heavy particles energy U m . The heavy particles internal energy U m , which represents the contributions of thermal and chemical energies, is derived from statistical thermodynamics with k B Boltzmann constant and N s the total number of included species in the global model. q T i and q int i stand for the translational and internal partition sum of species i, respectively. Subsequently, ε T i and ε int i represent the mean translational and mean internal energy of species i. Detailed derivation of equation (14) can be found in appendix A.
2.3.1. Gas temperature calculation from the heavy particles energy. In order to update the heavy particles (gas) temperature from updated heavy particles energy U m , The Brent method [24]-a root finding algorithm that combines bisection, secant and inverse quadratic interpolation-is implemented in PLASIMO to find the root of the following equation where U n+1 m is an updated heavy particles energy in each iteration of the global model, U T m t (T g ) is the total translation energy-first term on the rhs of (14)-and U int m t (T g ) is the total internal energy-contribution of electronic and ro-vibrational, second term on the rhs of (14). A temperature range is defined for the Brent algorithm in such a way that the above equation gives the opposite sign for the lowest temperature and highest temperature in the defined range. Then, in the form of a C++ structure, the above equation is passed to the Brent algorithm, and in each iteration, the temperature is updated until the root (T g ) is found, which satisfies U n+1 The calculated value of the temperature is used as the updated temperature for the next iteration of the global model.

Gas temperature modeling validation
In order to assess the performance of the implemented algorithm for the calculation of heavy particles temperature, three test cases have been designed. In this section, all of them are explained, and analytical results are compared with PLASIMO simulations.
3.1. One-temperature mixture with ro-vibrational contribution to the mixture energy Imagine we have a mixture composes of two species H and H 2 . At time zero, we only have H, and by marching in time, all the atomic hydrogen recombine to molecular hydrogen with the help of reaction 2H → H 2 . This reaction is exothermic, so we expect that the gas temperature increases in time. In addition, if the initial temperature of the mixture is higher than characteristic vibrational (θ V = hcν/k B (T g θ V )) and rotational temperatures (θ R = hcB/k B (T g θ R ))-see appendix A-we can have an approximation for the vibrational and rotational mean energy for H 2 as 2k B T g , one k B T g for the vibrational and the another one for the rotational motion.
Because in practice the above approximations cannot achieve except for real high temperatures, in the developed plugins in PLASIMO, manually the equations for partition sums of rotational and vibrational modes are changed to Therefore, the analytical internal energy should be equal to k B T g at any temperature for both rotational and vibrational modes.
If the rate coefficient and the energy difference between left and right hand-side of reaction 2H → H 2 are represented by k f and ΔE = 2E f H − E f H 2 (this reaction is exothermic and ΔE becomes the source term for the mean energy equation), respectively, we can write the system of ODEs for this mixture as In addition, the mixture mean internal energy U m is The analytical solution for the H particle balance is where n H 0 represents the density of H at initial time. For H 2 , we have this solution where n H 2 0 stands for the initial value of H 2 . The analytical solution of gas temperature has the following form and we can write The derivation of gas temperature solution can be found in appendix B. The explained example is modeled by PLASIMO with k f = 6.04 × 10 −18 m 3 s −1 , n H 0 = 1 × 10 25 m −3 , n H 2 0 = 1 × 10 −20 m −3 , T g 0 = 300 K, B = 60.853 cm −1 , andν = 4401.21 cm −1 . Figure 1 presents the results of the comparison between PLASIMO and the analytical result.

Two-temperature model with only elastic collision
Imagine that there is a mixture of H 2 and e. The only interaction between these two species is due to the elastic collision e + H 2 → e + H 2 . There are no reactions and no external power source. Therefore, the densities remain constants. The system of ODEs composes of the below equations where k e is the rate coefficient for the elastic process, n e and n H 2 are densities of electron and hydrogen, respectively. The reduced mass is also defined as The vibrational and rotation partition sums for H 2 are set similar to the previous example. Therefore, we can write In order to find an analytical solution for the above system of ODEs, we can define some parameters to simplify the task The final result for gas temperature has the following form where T g 0 and T e 0 represent the gas and electron temperature at the initial step. The analytical solution of T e has the following form The D is defined based on the initial condition The analytical derivations of electron and gas temperatures can be found in appendix C. The designed model is run by PLASIMO to check the validity of the gas temperature implementation in two-temperature medium. The following parameters are set for the input: k e = 1 × 10 −20 m 3 s −1 , T e 0 = 1300 K, and T g 0 = 300 K. The model, as expected, reaches the equilibrium point where the gas temperature and the electron temperature have the same value. Figure 2 presents the results of the PLASIMO calculations in comparison with the analytical solutions.

One-temperature mixture with power density
In order to validate the calculation of mixture energy based on partition sum in PLASIMO implementation, a single-species gas H 2 is designed with the rotational and vibrational partition sums as follows There is no reaction, and a constant power density is imposed to the mixture energy. Therefore, the only equation that should be solved is where P is power density with unit W m −3 . The analytical solution of the above equation is where U m 0 is the mixture energy at initial condition. Because there is no reaction considered in the model, the particle density of H 2 remains constant during the simulation. The energy for H 2 based on above partition sums can be calculated as The value for rotational and vibrational frequencies are set as B = 60.853 cm −1 andν = 4401.21 cm −1 . The purpose of this example is to check the numerical calculation of partition sum and its first moment (T ∂q ∂T ) that is needed for the calculation of the energy. The first moment is numerically calculated with the central finite difference with dT = 1 K in PLASIMO. By knowing the analytical value of U m (30), in order to have an analytical solution for T g , we need only to solve a quadratic equation. Therefore, the analytical gas temperature has the following form where

H 2 O-He global model in radio-frequency atmospheric plasma
The presented chemistry set in Liu et al study [2] is implemented in PLASIMO global model as a first step. It contains 46 species and 577 reactions. In order to verify the implementation, similar conditions to the atmospheric radio-frequency plasma of [2] are set to benchmark the developed model. The heavy particles temperature is assumed constant and has a value of 300 K. The plasma is formed between two circular electrodes with radius 1 cm and gap 500 μm.  figure 15 of Liu et al study [2]. In the presented global model in that study, no electron energy is solved [2]. They set a constant electron density for the whole simulation, and by solving the electron particle balance in such a way that electron density remains constant, the electron temperature is updated in each iteration of the model. The algorithm used in Liu et al research [2] is completely different than the PLASIMO global model. In the developed model, the particle balance is solved for each species except the electron. The density of the electrons is updated in each iteration by quasi-neutrality. It means where N + s and N − s represent all included positive and negative ions, respectively, in the model. For the electron temperature, the electron energy balance is solved as explained in the section 2.2. The wall processes for neutral, excited, and positive ions are considered in Liu et al study [2] in table 3. The same processes are implemented in the global model. The rate of the wall processes for neutral and excited species is modeled as where α is the probability of the reaction occurrence; m i and n i are the mass and density of the neutral particle. The flux is calculated based on a Maxwellian distribution [25, pp 176-180].
No wall processes for negative ions are considered, and the rate of the wall processes for positive ions are modeled as where u B is the Bohm velocity for positive ion [26, p 169]; n + i and m + i are the particle density and mass of the positive ions. The wall processes rates are multiplied by S electrod /V in order to convert to the proper sink term for the particle balance of relevant species. S electrode is the area of the electrode, and V is the volume of the plasma. In Liu et al work [2], probability for wall processes for positive ions are set to 1. However, when clusters of H 2 O collides with the wall, it is not clear in Liu et al study [2], what the products of such collisions are. In the developed model in this study, the products are considered as neutrals H 2 where F k is the volume flow rate that is set to 1 × 10 −4 m 3 min −1 [2]. The side diffusion loss is also considered for the neutral particles. The rate is calculated as (34). However, for the sink term of particle balance instead of the electrode area, the side area of the cylinder, plasma volume, is used. BOLSIG+ is used to calculate the rate coefficient of the electron impact reactions based on available referred cross sections in Liu et al work [2]. The EEDF is calculated in the beginning of the simulation with a mixture composed of only H 2 O and He and a fixed ionization degree 1 × 10 −8 . Because in the PLASIMO global model, the electron energy is self-consistently solved, electron energy sink terms for the positive ions wall processes are considered. The energy loss is equal to 2k B T e + 50 eV for each ion which is lost in the wall processes. This value is similar and consistent with the assumed energy loss in Liu et al study [2]. Because quasi-neutrality is used to calculate the electron density in the PLASIMO global model, in addition to 50 eV that is the energy penalty for each positive ion loss, 2k B T e is also considered for the corresponding electron that is lost in a wall process due to the quasi-neutrality. Note that, these amount of energy losses are based on the assumptions that Liu et al [2] made in their study.
The steady state results are reported in Liu et al study [2]. The particle balance equations are solved until the difference between updated values and old particle densities become less than a defined threshold. In order to reach steady state, either the input or source for the feed gas should be considered, or the feed gas density should be kept constant during the simulation. In a private conversation with the authors of Liu et al work [2], it was confirmed that the particle balance for H 2 O and He are not solved in their simulation. Therefore, in the developed model in PLASIMO, the density of these two species are kept constant. Figure 4 shows the comparison of the electron density and temperature between PLASIMO model and the reported results in Liu et al paper [2].
In PLASIMO, the power density for the electron energy balance is set to 20 W m −3 . Another important point is the high reported electron temperature in Liu et al study [2]. In a private conversation, the author confirms the factor 3/2 in the electron mean energy was missed in their calculation. In figure 4, this factor is applied for the results from Liu et al [2] and depicted again. In conclusion, due to the different algorithm for global model, and possible different modeled terms in particle balance equations, the agreement between two models is acceptable. In the next step, the chemistry set is refined, and the developed model is applied for the microwave plasma discharge.

H 2 O-He global model in microwave induced plasma
In order to have a complete and more accurate kinetics as much as possible, the chemistry set of Liu et al study [2] is refined as follows: • For most of the electron impact reactions, the suitable cross sections from literature are found and used. Therefore, instead of using rate coefficients based on Maxwellian EEDF, rate coefficients for these reactions can be calculated by BOLSIG+. In appendix E, one can find comments about the chosen cross sections for most of these reactions. are removed from the list because at high microwave temperature, the production and sustainability of these species are not very high. • Production of added species from the ground step, stepwise ionization from the added excited species, superelastic collisions for all added excites species, vibrationaltranslational (VT) reactions for vibrationally excited species, quenching reactions for electronically and rotationally added exited species, and some missing reactions in elastic collisions, step-wise ionization, dissociative ionization, dissociative recombination, and superelastic collisions for species that already were used in Liu et al study [2] are added to the chemistry set. The complete list of reactions can be found in the attached electronic files for the microwave model. In appendices E and D comments on the used cross sections and added reactions are explained in more detail.
In total, the refined model has 56 species and 630 reactions. Similar to the previous section, the EEDF is calculated in the beginning of the simulation with a fixed mixture composition (a fixed ratio of H 2 O in He, see section 5.1) and a fixed ionization degree 1 × 10 −5 .
The performance of the refined model is checked for the microwave reactor. The forward vortex microwave reactor of DIFFER is used in this study. The plasma is ignited and sustained by a 2.45 GHz microwave field. A rectangular waveguide (WR-340) intersects a quartz tube which contains the plasma. The flow is injected tangentially from the top of the tube to produce a swirl that confines the plasma to the center of the tube and thus prevents destructive contact of the hot plasma core with the tube wall. The rectangular microwave cavity with 2.45 GHz and 43.18 mm height is used in this experimental study [27]. The formed electromagnetic field in the cavity is TE 10 mode in a standing wave configuration, and the power can be adjusted to produce a stable plasma. More information about the setup and used waveguide can be found in [27][28][29][30].
The volume flow rate of He which is the carrier gas in this study is controlled and adjusted by a mass flow controller (MFC) and expressed in standard liters per minutes, the water mass flow rate (g hr −1 ) is also controlled by the MFC. The H 2 O and He are mixed in an evaporater and then injected into the reactor. In this study, the He flow rate is kept constant at 10 slm, H 2 O has a fixed mass flow rate 55 g hr −1 , power varies between 300-1000 W, and pressure is kept constant at 500 ± 4 mbar.
The main challenge is to convert the 3D experiment to a 0D model. In this context, several parameters should be set as input or modeled as frequencies. In the following sections, these parameters are discussed.

H 2 O ratio in He-H 2 O mixture
In order to define the ratio of H 2 O in the He + H 2 O mixture, conserved quantities like mass or number of particles-if there is no reactions-should be used. The measured values from the experiment are the mass flow rate of liquid water in g hr −1 , and the volume flow rate of helium gas in slm. Therefore, either the mole flow rate ratio or the mass flow rate ratio should be used. Mole flow rate in this study is a good choice due to the non-reactive nature of the mixing of water and helium before injection into the reactor and particle based nature of the PLASIMO global model. Mass flow rate ratio usually is the better option for mass density based solvers like typical Navier-Stokes solvers. The values that are used for unit conversion in this study are as follows: • He: molar weight 4.002 g mol −1 , mass density 0.164 kg m −3 or g L −1 (at 20 • C and 1 atm), • H 2 O: 18.015 molar weight g mol −1 .
Therefore, the molar flow rate ratio of H 2 O in the mixture can be calculated as whereṁ H 2 O andV He stand for the mass flow rate and the volume flow rate. M w and ρ m are the molar weight and mass density. For the case with 10 slm He and 55 g hr −1 H 2 O, the ratio of H 2 O in the injection mixture is 11.04%. In the following equations, unit conversions are done and each symbol presents a value in SI units.ṁ has unit kg s −1 ,V has a unit in m 3 s −1 , M w has a unit in kg mol −1 , and ρ m has a unit in kg m −3 .

Plasma volume
The volume of the plasma is important to determine the correct power density for the electron energy balance (3). In addition, to model transport and other possible terms like heat transfer mechanisms in the heavy particles energy balance, the geometry of the plasma plays role. In this study, the volume of the plasma is determined by the optical emission profiles measured by a CCD camera. Figure 5 shows the variation of optical emission at various powers for 11.04% of H 2 O in He. The plasma volume is approximated by a cylinder whose height and diameter is full width at half maximum (FWHM) of the fitted Gaussian profiles on the Abelinverted light emission intensity profiles in axial and radial lines that pass the center of the plasma. Figure 6 shows an example of this procedure. Below table lists the obtained height and diameter of plasma volume based on experimental measurements (table 1).

Inlet flow
Because the steady state results are important in this specific reactor, the global model should be run until it reaches equilibrium (no variation versus time). If only chemical source terms are taken into account, it is not possible to a reach steady state condition without keeping the density of feed gases (H 2 O-He) constant. In other words, because we do not want to keep the density of the feed gases constant in this study, transport terms are added to the model. The inlet flow into the volume of the plasma is a source for H 2 O and He particle balances. These source terms can be good representations of the convective inflow to the assumed plasma volume. They are obtained by converting the mass flow rate of H 2 O and the volume flow rate of He to the number of particles and dividing by the volume of the plasma. The used equation and for He inflow iṡ If the heavy particles energy balance is solved self-consistently with the rest of particle balances and the electron energy balance, the enthalpy that this inflow adds to the volume also should be taken into account (see equation (12)). For this   (38) and (39), proper unit conversion is used to produce the correct value for the inlet flow frequency in 1 m −3 s −1 .

Convection loss
The convection loss in this study is only considered from the outlet cross section of the cylinder that represents the volume of the plasma. This is the first term on the rhs of equation (2). This transport term is considered for all species involved in the model. In order to define the proper frequency, we assume that the bulk velocity perpendicular to the exit section, u b , is constant in the whole area of the outlet section of the plasma volume. In order to calculate u b , one can use the conservation of mass, but due to observed numerical instability in particle based global model, it is not used in this study. The alternative method to handle this situation is to use the total particle inlet flow rate (s −1 ) and the total density of the model. Then, by knowing the pressure and gas temperature, the volume flow rate loss for each species can be calculated as where N A is Avogadro constant. Again proper unit conversion is used in (40) to produce an appropriate volume flow rate loss for each species in m 3 s −1 . N t is calculated based on the updated gas temperature in each iteration and pressure of the vessel (that is assumed constant) Then, the convection loss for each species can be modeled aṡ In addition, this loss has an impact on the heavy particles energy balance, 5/2k B T g is the considered energy loss per particle in this simulation. Again, enthalpy loss is taken into account (the first term on rhs of (12)) and the effect of U int i for molecular species is not considered.

Diffusion loss
The diffusion loss in this study is considered only from the side (the fourth term on rhs of equation (2)). This loss is considered only for neutral species. Inward diffusion is neglected.
u d b i can be modeled with the help of a diffusion coefficient and a length scale that is suitable for the direction along which the diffusion occurs. In this simulation, where R is the radius of the cylinder that represents the plasma volume, and D i stand for the diffusion coefficient. In this study, the diffusion coefficients for neutral particles is approximated with following equation [26, pp 133-136] with m i mass of species i and ν m i the momentum collision frequency.
Here, the momentum collision frequency. The momentum collision frequency can be estimated as where v i , σ, and n He represent the mean velocity for species i, the momentum transfer cross section, and the particle density of He that is the dominant species in this study, respectively. If we assume that the heavy particles have a Maxwell-Boltzmann distribution, the mean velocity can be calculated as [25, pp 176-177] Due to the fact that He is the dominant back ground gas, in this study, the momentum transfer cross section of species i with He is considered. A fixed value 10 −17 m 2 is assumed as the cross section for all neutral species. Because heavy particles energy balance is self-consistently solved in this simulation, 5/2k B T g is the enthalpy per particle which is assumed to be lost due to the inter-diffusion heat flux, equation (12). Finally, the diffusion loss like the convection loss can be written asṅ

Heavy particles energy balance
In contrast to section 4 where the heavy particles temperature is assumed to be constant, in the microwave induced plasma, rising of the gas temperature is expected. Therefore, the heavy particles energy balance is solved self-consistently along with the electron energy balance and particle balances for each involved species in the model. Equation (12) in accompanied by equation (14) for estimation of the heavy particles energy is added to the ODE system of this global model. As is explained in sections 5.4 and 5.5, the first and fourth terms on the rhs of equation (12) are assumed to be transport related frequencies for the heavy particles energy balance in this model. The chemical source term is calculated as explained in section 2.2. The elastic source term has the same value as the elastic term in the electron energy balance but with the opposite sign (it is a source for the heavy particles energy). The electron dissociation reactions that have contributions to the gas heating in this model, see sections 2.2 and 2.3, are listed in table 2. The amount of energy for each reaction that has a contribution to the gas heating (after multiplying by the rate of the reaction) is written in front of each reaction. This extra energy is also considered as a sink term (on top of ΔE) for the electron energy balance.
A partition sum calculation for each species is done in order to have a proper heavy particles energy representation, see equation (14). For atomic species, only translational energy is considered. If for one species, different electronic states should be considered, each state is taken into account separately in the model. For example, O represents the ground state of oxygen atom with zero energy, and O(1D) and O(1S) represent two different electronic states of the oxygen atom with energies 4.5248 eV and 6.7472 eV. Diatomic molecules are defined in one of the following four cases: • One specific electronic state is considered, so the partition sum should be done over all rotational and vibrational levels for this electronic state (second and third summation of equation (61)) if the necessary frequencies are available. Species which are modeled under this case are H 2 , • One specific electronic state is considered, but for the internal partition sum the harmonic oscillator equation (66) and the rigid rotor equation (64) approximations are used for vibrational and rotational contributions, respectively. Species which are modeled under this case are OH, OH − , OH + , OH(A), He * 2 , and He + 2 . • Electronic state and vibrational level are specified, so the partition sum is done over all possible rotational levels (the third summation in equation (61)). Species which are modeled under this case are H 2 (v = 1-3) and O 2 (v = 1-4). • The state is defined completely (electronic state is defined with specific rotational and vibrational quantum numbers), so only translational energy is taken into account like H 2 (ν 0 , J 2 ).
If no frequency is available for the diatomic molecule, only translational energy contribution of the species is considered in the energy equation. Two examples for this condition are H 2 (b 3 Σ + u ) and HeH + . Note that for all explained cases above the translational energy are taken into account besides internal energy (14).
If the electronic state is defined for polyatomic molecules, all possible vibrational and rotational levels for this specific state are summed (second and third summation of (63)) to calculate the internal energy contribution besides the translational energy-species which are modeled under this condition are O 3 (12) for heat transfer frequency, S h . In the developed microwave induced plasma in this study, there is no need to take into account heat conduction loss to the surrounding of the assumed plasma volume because the cold gas enters the plasma by inflow and hot gas leaves the volume by convection and inter-diffusion losses, see sections 5.4 and 5.5. However, in the results and discussion part in section 6, the impact of heat conduction loss will be discussed. The conduction loss is due to the −∇ · q term in the conservative form of the total heavy particles energy equation [17, p 24]. Equation (12) is derived with the same procedure that has been applied to the continuity equation in section 2.1. If the integration over the volume is converted to the surface, for heat transfer gradient, we have In the above equation, only conduction heat transfer is taken into account. The above equation can be simplified by defining the S side and A plasma for the assumed plasma volume as follows, similar to sections 2.1 and 2.3, where L and R are the height and the radius of the plasma volume in this study, and k is the heat conductivity of the mixture. ΔT g is the difference between initial gas temperature-that is assumed to be equal to the surrounding temperature of the plasma volume-and the updated gas temperature in each iteration of the model. The initial gas temperature in all test cases in this study is set to 500 K. Volume division in (47) is due to the present of the plasma volume in lhs of the energy equation after integration and simplification similar to section 2.1 In the modeled heat conduction in equation (47), it is assumed that the heat transfer loss is S q = q c . Inter-diffusion heat transfer loss is explained in diffusion loss section 5.5. For heat conductivity, the tabular heat conductivity of He versus temperature is used in this study because He is the dominant species.

Electron density measurement
The averaged electron density in the plasma is obtained noninvasively from a quasi-optical free-space microwave interferometry [29]. This technique uses phase shift measurements of a 168 GHz diagnostic beam through the plasma medium, in comparison to free space propagation, to infer the dielectric properties of the medium (i.e. electron density). The beam is focused with a horn antenna assembly, with the beam axis aligned to the center of the discharge, perpendicular to the axis of the plasma filament. A focus waist of 2-3 mm is realized, which is smaller than the typical plasma diameter.
Using the Lorenz conductivity model, the phase shift and plasma width for a lossless medium can be related to the average electron density according to [31] as with ω d the diagnostic wave frequency, ε 0 the permittivity of free space, m e the electron mass, e the electron charge, Δφ the measured phase shift, and c the speed of light. d p is the width of the plasma that is equal to the diameter of the cylinder that is assumed as the plasma volume in this study. The spot radius ( 1 e of the peak power of the diagnostic beam) is approximately 4.5 mm at the location of the plasma based on a quasi-optical beam path assessment using transfer matrices in the paraxial approximation. The infinite slab approximation breaks down for the typical plasma diameters of under 3 mm measured, as more than half of the diagnostic beam may circumnavigates the plasma medium in extreme cases. A correction to account for the finite plasma width in infinite slab approximation calculation is incorporated by scaling the measured phase shift in (48) by a factor which is inversely proportional to the fraction of beam power that propagates through the plasma region. More information about the used correction for measured electron densities can be found in supplementary materials, part A, section 1.

Gas temperature measurement
The rotational temperature of the plasma is obtained from the rotational spectra of OH(A-X). An optical fiber is used to collect optical emission from plasma, which makes the spectra volume averaged. The relative intensities of the Q, P, and R-branches of the rotation spectrum of OH(A-X)(0, 0) are simulated and fitted to the obtained spectra (Ocean Optics HR + C2928) using the MassiveOES fitting tool [32]. The transition parameters are taken from [33][34][35]. Due to the high temperature of the plasma, which typically exceeds 2000 K, a Boltzmann distribution of the rotational states is assumed and a single temperature is warranted to describe the rotational structure. Deviation of the rotational population from the Boltzmann distribution associated with strong electronic quenching of OH(A-X) by H 2 O [36] are mitigated by the short timescales of thermalization associated with the high gas temperatures. Figure 7 shows the comparison of the measured electron density with the predicted one by the developed global model. The error bars of the experiment are determined from the phase shift measurement and alignment of the beam errors. The impact of possible errors due to H 2 O concentration, pressure, and input power value are not taken into account and considered negligible. In addition, it is assumed that the chosen methodology to determine the width of the plasma for the calculation of the averaged electron density (equation (48)) is the best possible option.

Results and discussion
As shown in figure 7, the model predicts the same trend as the experiment for electron density variations versus power. The simulated electron densities are in the same order of magnitude as the measured ones with the maximum difference of two times of corresponding experimental values. The observed discrepancy can be due to several factors like the amount of flow that enters the plasma volume, transport frequencies, and heat transfer mechanisms. In fact, the modeling of transport mechanisms is generally difficult to properly implement in 0D models. In the following parts, the impact of those parameters are tested by doing a sensitivity analysis study.
As is discussed in the last section, the conductive heat transfer loss is not taken into account in the heavy particles energy balance in the first attempt because it is believed its impact already is taken into account by cold inflow and loss of hot outflow due to convection and inter-diffusion losses. However, in order to check the impact of the overestimated heat loss on the electron density or other parameters of the model, this frequency is modeled by equation (47) and its effect on the electron density is shown in figure 8. As can be seen, the conduction heat loss has negligible impact on the electron density, although it has a higher impact on the gas temperature, see figure 9.
One of the important parameters on the electron density value is the TFR of feed gases that enter the volume of plasma. In this study, it is assumed the estimated plasma volume and subsequently estimated power density, as is explained in section 5.2, are determined properly. However, the TFRs of H 2 O and He that enter to the plasma volume can be a fraction of the total injected flow (because the volume of plasma is smaller than the volume of tube where microwave cavity crosses the tube). In order to see the impact of the flow on plasma properties, the sensitivity analysis is done by assuming 2/3, 1/2, and 1/3 of mass flow rate of H 2 O equation (38) and volume flow rate of He equation (39) enter the plasma volume.
As is shown in figure 8, the predicted electron densities when 2/3 of total flow enters the plasma volume shows the best match with experimental value. However, there is a huge deviation from experimental values, especially at higher powers, when only 1/3 of the total flow enters the plasma volume. This behavior is expected because by decreasing the amount of flow that enters the plasma volume at constant power density, production of electron density should be enhanced, and this enhancement should be more at higher powers (power densities). From both comparisons of electron densities (figure 8) and gas temperature (figure 9) with experimental values, it can be concluded for the studied cases in this paper by 0D models, 2/3 of the total injected flow rate of H 2 O and He into the reactor enter the plasma volume. Figures 9 and 10 show the variations of some important variables in the modeled microwave induced plasmas versus power for different flow rate with and without heat conduction loss. In addition, figures 11-13 show the variation of neutral, positive, and negative species versus power for the case that only 2/3 of total injected flow rate enters the plasma volume without heat conduction loss. H 3 O + , H 2 O + and O + 2 are recognized as the most important ions in this study. H 2 O + has the same trend as the electron density by variation of the power, flow, and consideration of the conduction heat loss. The main reason is that the electron density is calculated due to the quasi-neutrality in these models, and H 2 O + is one of the main positive species, and one of the main path to its production is from electron-direct ionization.
Electron temperature shows the same trend as H 3 O + in figure 9. The behavior can be justified due to the main production paths of H 3 O + which are from reactions of H 2 O + with H 2 O and OH. In other words, electron temperature is calculated with the help of the mean electron energy (3/2n e k B T e ) and electron density, electron density has the same trend as H 2 O + , and H 2 O + should be consumed for production of H 3 O + . Therefore, at the same power density, it is expected that the behavior of electron temperature becomes the same as H 3 O + .
Production of high gas temperature in figure 9 is the proof of the necessity of including mixture energy balance for modeling of these types of plasmas. As is shown, including the conduction heat loss overestimates heat loss from the plasma volume, especially by assuming the temperature of the surrounding gas equal to the injection temperature. This loss has bigger impact on the gas temperature at a higher power due to a larger ΔT g , equation (47).
Higher gas temperature is achieved by decreasing the TFR that enters the plasma volume in the same power density as expected. As is shown, when 1/3 total injected flow rate enters the plasma volume, predicted gas temperature by the model becomes unrealistic at high powers. As mentioned above, the predicted gas temperature by the model has less deviation from the measured values when 2/3 of total injected flow enters the plasma volume.
Pressure is another variable that is shown in figure 9. In calculation of the convection loss in this study (equations (40) and (41)), the impact of pressure is taken into account. However, there is no constraint imposed on any variable in the current version of the PLASIMO global model to keep the pressure constant. As is shown in figure 9, the maximum deviation of the pressure among all test cases in this study is 60 mbar. It can be accepted to some extent for these specific modeled conditions. However, for more sensitive cases, for example, when the plasma pressure is very close to the boundary of the discharge mode change (for instance, in microwave induced plasma for CO 2 [30], four modes in a pressure range 100-500 mbar is reported), even a 12% pressure change in the global model can have a big impact on the analysis of the results. A methodology to solve this issue can be found in [17, p 50-52].
Variation of some of the interesting species H 2 O 2 , OH and H 2 versus power at different flow with and without conductive heat loss are shown in figure 10. H 2 O 2 has a peak value and then by increasing the power, its density decreases. This behavior can be justified due to the higher gas temperature and subsequently more H 2 O 2 dissociation at higher power. OH and H 2 have a nonlinear behavior, but at high powers (600 W and higher), these three species have lower values when lower flow enters the plasma volume. Decrease in the amount of flow that goes into the plasma volume produces higher gas temperature and higher dissociation rates. In other words, at lower flow, the same power density produces higher temperature and subsequently higher dissociation rates for molecular species like OH, H 2 O 2 and H 2 .

Conversion and energy efficiency
For the application of hydrogen production, the model results are compared to experimental values associated with this work and reported in literature. The experimental conditions for this comparison are: 98.70% of H 2 O with a trace of He gas at pressure 190 mbar and power 900 W. The flow rate of H 2 O is 404 g hr −1 and volume flow rate of He at injection is 0.12 slm. For this case, the diameter and height of the cylinder that represents the plasma volume are set 6.894 mm and 28.786 mm, respectively. The conversion and energy efficiency are calculated as and respectively. ΔH is the enthalpy of the reaction H 2 O → H 2 + 1/2 O 2 , which is 2.4762 eV. SEI is the specific energy input eV molecule −1 , which is calculated from the power and H 2 O flow rate particle/s. H 2 O| in is calculated based on the pressure of the reactor (190 mbar) and injected temperature (500 K). H 2 | out is the value of H 2 at the steady state and also in afterglow simulations when the gas temperature reaches 300 K. In afterglow simulations, the density values and temperature of steady state are used as initial conditions, power is set 0 W, and electron temperature is assumed constant 500 K. In addition, no inlet and outlet frequencies due to convection and diffusion are set in these models. Figures 14 and  15 show the variation of H 2 density and gas temperature versus time in three different afterglow simulations. 'Conductive cooling' refers to the afterglow simulation that cooling with surroundings is considered by the conduction loss (47). The surrounding temperature is set 300 K for this simulation. A cooling mechanism for the other two afterglow simulations are set as where K c is the rate of cooling K s −1 that is set 1 × 10 7 and 1 × 10 8 . N t is the total particle density, see section 5.4, that is updated based on the pressure and updated temperature in each iteration of the model. The values for conversion and energy efficiency from modeling are reported in table 3. The H 2 conversion 0.1% is much below the expected chemical equilibrium values of around 10%-20% reported in [3]. An evaluation of the net H 2 production from the plasma core conditions depends critically on the extent of its preservation during cooling of the product mixture. However, the amount of H 2 in the discharge and afterglow seems to be in the same order of magnitude compared to experimental investigation of the reactor output in this study. Composition measurements of the cooled effluent gas in experiments related to this case had concentrations corresponding to approximately 0.8% conversion of H 2 O to H 2 and an energy efficiency of 1.2%. The measurements were conducted using gas chromatography, for further information please refer to supplementary materials, part A, section 2.
The water vapor was premixed with a trace of He and used as an internal standard to account for the concentrating effect of water condensation in the sampling lines on the measured H 2 and O 2 concentrations. A gas temperature of around 4000 K was obtained from the fits of the OH(A-X) spontaneous emission, which leads to a significant degree of thermal decomposition of H 2 O in the high temperature region.
Values of modeling and experiment correspond with varying degree to literature values. Plasma-driven water vapor dissociation for hydrogen production has been previously investigated in high power thermal arcs [37] and microwaves plasma [15,38,39]. Jung et al reported efficiencies of 7%-8% at reduced pressures in the range of 1-26 mbar. The dominant mode of dissociation in these experiments was attributed to direct dissociation of H 2 O by electron excitation rather than thermal dissociation, as the gas temperature did not exceed 1500 K while the electron temperature was relatively high at 5.6 eV. Asisov et al [15] reports an energy expenditure of 3.2 ×10 7 J m −3 (for H 2 ) equivalent to an energy efficiency of 35% for a microwave plasma experiment close to the one explored in this work. These optimal conditions were achieved at a pressure of 50 mbar and 0.52 g s −1 water vapor. Besides the difference in pressure to the plasma considered in this work, the other plasma conditions are comparable. The gas temperature of over 3000 K and the relatively low ionization degree of 10 −5 strongly indicates that thermal conversion is the dominant mechanism for heavy particle chemistry rather than charged particle-driven reactions. While the reported efficiency values in these experiments far exceed those obtained in our own experiments and modeling work, other thermal dissociation experiments are in closer agreement. Boudesocque et al [37] report an energy efficiency of 0.7% in an 5 kW thermal arc discharge at atmospheric pressure, which is in rough agreement to our own findings. Despite the different nature of thermal arcs, the discharge is also characterized by gas temperatures which are sufficiently high for thermal conversion. The relative change of hydrogen concentrations in the plasma and in the cooled effluent, which were measured respectively by optical emission peak ratio analysis and mass spectrometry, suggested that only 15% of H 2 in the plasma is conserved after quenching. In our own experiments a rapid recombination of plasma products is also evident by a red-hot glow of the discharge tube just downstream from the plasma (see figure 2 in supplementary materials, part A, section 2), which is caused by the strong local heat release form the strongly exothermic recombination reactions. The recombination of thermally dissociated species in the post-discharge region of thermal H 2 O plasmas therefore appears to be of great influence on the net H 2 production.
Fast quenching of the plasma products by rapid cooling can be essential to prevent recombination reactions which reduce H 2 output. In general, cooling rates depend to a large extent on the gas dynamics in and around the discharge and method of heat extraction. Cooling rates of 10 5 -10 6 K s −1 are common by gas-dynamic mixing, while higher rates are only considered achievable by more elaborate means such as supersonic expansion [40]. Although, based on the results of our afterglow modeling, even high cooling rate 10 8 K s −1 does not increase the conversion and energy efficiency of the hydrogen production in pure water microwave plasmas. Adding a source of carbon to the plasma for consumption of oxygen can be a better option to produce hydrogen form water.
Production of H 2 O 2 from H 2 O is another interesting application of H 2 O microwave induced plasmas. For the same  (53)) than the one used for H 2 is used to calculate energy efficiency. .
To date, there are barely any experimental nor numerical works available in the literature that study production of H 2 O 2 using H 2 O microwave induced plasmas, except the referred kinetic modeling in [3] for a microwave induced plasma in a supersonic flow. The low temperature for that specific configuration helped to stabilize the production of hydrogen peroxide in afterglow. Vasko et al [9] calculated an energy efficiency of 24.0 g kW h −1 for this type of microwave induced plasma. The calculated results in the present paper for three different afterglow simulations have the maximum energy efficiency of 2.56 g kW h −1 for a cooling rate of 1 ×10 8 K s −1 . Although, the obtained result is lower than supersonic configuration, it shows the impact of cooling mechanism in stabilization of hydrogen peroxide in the effluent of the microwave plasma reactor. In addition, Vasko et al reported maximum energy efficiency 0.12 g kW h −1 for production of H 2 O 2 in atmospheric pressure radio frequency glow discharge in helium-water mixture. While comparing to other plasma reactors, the calculated H 2 O 2 yields for a microwave plasma reactor reported here indicate that using a thermal plasma may be an alternative way. More experimental measurements and coupling with multi-dimensional models will however be required to understand and control the complex gases and heat flows interplay that will limit yields of both H 2 and H 2 O 2 species.
Note that the lack of spatial resolution inherent in global models remains a major limitation to its application in evaluating the intrinsically multi-dimensional aspects of the recombination kinetics and transport under influence of the gas flow dynamics. Even so, the elaborate global plasma chemical description of the H 2 O-He chemistry in this work provides valuable insight into the discharge kinetics under the 'warm plasma' conditions relevant for thermal decomposition of water, serving as a stepping stone for more elaborate reactor modeling of such systems. In addition, the practical feasibility of realizing the exceptionally rapid cooling rates needed to conserve most products in thermal H 2 O process under subsonic gas flow conditions remains to be investigated experimentally.

Conclusions
A detailed description of heavy particles energy balance implementation for a two-temperature global model is presented in this paper. A statistical thermodynamics approach is chosen to derive the balance equation for internal energy of heavy particles and calculate the gas temperature in each iteration of a model. Three analytical test cases are designed to assess the capability of this newly implemented functionality into the PLASIMO global model module. It is shown that the calculated internal energy of heavy particles and subsequently updated gas temperature in each iteration of models are very well consistent with analytical results of the designed test cases. This new functionality is useful to describe and predict behavior of plasmas with significant gas heating, like microwave induced plasmas.
A global model for H 2 O-He mixture is developed and validated systematically with application to the production of hydrogen from water and subsequently solar fuels from a mixture of CO 2 -H 2 O by microwave induced plasmas. The introduction of the chemical model is accompanied by an extended discussion of the frequencies that are used to represent transport phenomena and heat transfer mechanisms. First, the model is verified and benchmarked with reported values for electron density and temperature of a radio frequency plasma taken from literature [2]. Then, the kinetics are refined, and the model is extended by solving self-consistently the heavy particles energy balance next to particle and electron energy balances and applied to a H 2 O-He microwave induced plasma. The results of the model are validated with experimental measured values of the electron densities and gas temperatures. It is shown that the developed model is able to capture the trends and qualitative behavior of the microwave induced plasma properly.
Sensitivity analysis is done over the amount of flow that goes into the microwave plasma volume in order to find the proper amount of the injected feed gases that enter the plasma volume. It is shown that based on the results of 0D models under reported conditions of microwave induced plasmas in this study, when 2/3 of the injected feed gases enter the plasma volume, the predicted electron densities and gas temperature by the model are the best matched with corresponding experimental values. It is expected that less than the total amount of injected flow goes into the plasma due to the multi-dimensional nature of microwave forward vortex reactor used in this study. The effect of considering conductive heat transfer frequency in heavy particles energy balance is also investigated, and it is observed that this frequency is not needed for proper modeling of the microwave induced H 2 O-He discharge because of the chosen approach to include convection and diffusion frequencies and their energy losses in the presented model. Due to the nature of the vortex interactions with the plasma, it is actually expected that thermal conduction does not play any role in the heavy particles energy balance in the core of the plasma [41].
Production of H 2 from almost pure (98.70%) H 2 O in the microwave induced plasma is measured experimentally and modeled in this study. Both model and experiment have conversion below 1%. It leads to the conclusion that pure H 2 O kinetics has a major lost channel for H 2 irrespective of the cooling rates based on the results of the afterglow modeling. Adding a source of carbon to feed gases in order to absorb oxygen can provide a more efficient path to obtain H 2 from H 2 O by microwave induced plasmas.
Production of H 2 O 2 from 98.70% H 2 O microwave induced plasma is also modeled for three different thermal quenching rates. A maximum energy efficiency of 2 g kW h −1 is calculated for the afterglow simulation with a cooling rate of 1 × 10 8 K s −1 . Further experimental measurements and coupling with multi-dimensional models are however needed to investigate how to preserve the hydrogen peroxide in the effluent of the plasma. Interestingly, higher energy efficiencies are predicted by the model for H 2 O 2 production compared to H 2 .
It is shown for the microwave induced H 2 O-He mixture, the variation of pressure in the model is negligible. However, as an outlook, implementation of a proper method to keep pressure constant during simulation [17] is necessary, especially for plasmas for which the mode structure is sensitive to pressure changes. Such an example is the transition from diffuse to contracted regime in CO 2 discharges [30,42].
The input files, cross-sections data, and analytical validation test cases for the heavy particles energy balance presented in this work are made freely available as supplementary data to this paper. The analytical validation test cases for the heavy particles energy balance presented in this work can be used for the benchmark of zero-dimensional plasma models. In addition, the input files for H 2 O-He plasma chemistry will help efforts toward a larger availability and reproducible data for plasma modeling.

Acknowledgments
This research is supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) and Alliander N V in the framework of Project 13581 'PlasmaPower2Gas (PP2G): efficient electrical to chemical energy conversion via synergy effects in plasma and catalytic technology'.

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

Appendix A. Derivation of the internal energy (U m ) for heavy particles
As a first step, we should write the partition function for a mixture of atoms and molecules. Imagine we have a mixture with N s different species. Because we have interacting system, we should use canonical partition function [17, p 196] and [43, p 577-578], for each snapshot or iteration of the model. In each iteration we know the particle density of each species, let us at this stage we work with the number of particles instead of density. We assume that the number of particles for each species is the same for all members of the ensemble, so the energy of each member of ensembleẼ i can be written asẼ where N 1 ε i,1 represents the total energy ε i,1 of species 1 with population N 1 in member i of the ensemble and so on for other species. In order to write the canonical partition function Q for a mixture we have states Ns j e −βε j N Ns withÑ the total number of the ensemble, q i the molecular partition sum [43, p 564], [17, p 186] of species i, and β = 1/k B T g . If each molecule or atom from a single species is distinguishable from each other, equation (55) is acceptable for a mixture. Otherwise, for the case that each particle is indistinguishable from the rest of the particles of the same species, we have Now by knowing the canonical partition function, we can calculate the heavy particles internal energy as where V stands for volume. Two points should be mentioned about equation (57). First, formation energies are not taken into account because it is supposed that these terms will be calculated as a source term on the rhs of the heavy particles energy equation. Second, the number of particles is independent of the temperature in one snapshot of the model. In other words, we assume that the total number of particles in each snapshot of the whole system is constant and independent from the temperature (this is one of the assumptions underlying calculation of the Boltzmann distribution and the partition function in the forms that are presented in this paper). We look at a system and assume thermal equilibrium only in a moment of time. We have variation in time, but at each moment we have thermal equilibrium.
The unit of U m is J for equation (57), if we substitute number of particles N i with the particle density (n i ) we have the internal energy with unit J m −3 , which is represented by U m .
The molecular partition function q i for each species that are considered in the model can be written as a product of the different contributions of energy sources for a molecule, for example, energy of different independent modes [43, p 579] In equation (58), the energy of molecule i is written as a summation of translational (ε T ), rotational (ε R ), vibrational (ε V ), and electronic excitational (ε E ) contributions. Therefore, the molecular partition function can be written as Now by substituting equation (59) in equation (57), we can rewrite the heavy particles internal energy as and in per volume version as ε T i is equal to 3/2k B T g , and ε int should be calculated from the internal partition sum of each species that is included in the global model. In order to calculate the internal partition sum of each species, we need to provide an input temperature. In PLASIMO, the internal partition sum for diatomic molecules is calculated as [44, p 86] where σ is the symmetry factor, and β = 1/k B T g . n m , ν m , J m are maximum quantum numbers for electronic excitation, vibrational and rotational motion. g n stands for the degeneracy of nth electronic state of the molecule. For diatomic molecules, there is only one mode for vibration, and the degeneracy of vibrational motion is 1. (2J + 1) is the degeneracy of rotational motion. The energy of a molecular state can be written as contributions of electronic excitation, rotational and vibrational energy as The energy of ν-th vibrational level at nth electronic state (ε V nν ) and the rotational energy for nth electronic state and ν-th vibrational level (ε R nνJ ) are calculated base on explained approaches in [17, pp 191-194.], [44, pp 86-91], [45, pp 22-24], if all needed spectroscopic frequencies are available.
Polyatomic molecules have more than one vibrational mode, and the degeneracy of those modes can be greater than one. The internal partition sum for polyatomic molecules can be written as [44, p 95] where [ν] stands for the array of quantum numbers for all possible m modes of the vibrational motion, , and g R n[ν]J are degeneracy of electronic state, degeneracy of vibrational motion that depends on electronic state and vibrational quantum numbers for each mode, and degeneracy of rotational motion that depends on electronic state, vibrational quantum number array and rotational quantum number, respectively. The energy and degeneracy of different modes for polyatomic molecules are calculated based on explained approaches in [44, pp 95-97] and [17, pp 194-196], if again all needed spectroscopic frequencies are provided.
It is possible, more simplified approximations for rotational (rigid rotor) and vibrational partition (harmonic oscillator) sums are used if all necessary spectroscopic frequencies are not available.
The molecular partition function for a linear molecule with rigid rotor assumption can be calculated as where h is Plank's constant, c is the speed of light, and B is the rotational constant that can come from either spectroscopy or calculation [43, p 564]. For the condition that temperature T g is higher than characteristic rotational temperature θ R = hcB/k B (T g θ R ), the mean rotational energy is calculated as [43, p 600] The molecular partition function for a harmonic oscillator can be calculated with respect to the zero energy level (1/2hcν) as [43, p 596] whereν stands for the wavenumber of the vibrational mode. The mean vibrational energy for a harmonic oscillator can be calculated as [43, pp 600-601] For the case that temperature is higher than characteristic vibrational temperature θ V = hcν/k B (T g θ V ), the vibrational partition function can be approximated as

Appendix B. Derivation of the analytical solution of gas temperature for the first analytical test case of the heavy particles temperature validations
For the gas temperature, after some manipulations, we can reach this equation With further simplifications, the gas temperature equation is obtained as In order to solve equation (69), we should write it in the following form 2k f n H0 t + 1 3n H0 + 7n H2 0 + 7k f n H0 t n H0 + 2n H2 0 k B 2k f n H0 t + 1 3n H0 + 7n H2 0 + 7k f n H0 t n H0 + 2n H2 0 To solve the above equation, the following integral should be solved first Now we can write the integral (71) as a be Then, we need to calculate Also in order to find the analytical solution for energy equation we need to find an answer for We can write Therefore, we have .
Hence, the final result for the gas temperature equation (69) has the following form The value of g can be defined by the initial value of gas temperature Now by substitution of parameters, we can write The final form of the gas temperature can be obtained by substitution of (79) and (80) in (78). The result is presented in equation (21) inside the main part of the paper.

Appendix C. Derivation of the analytical solution of electron and gas temperatures for the second analytical test case of the heavy particles temperature validations
Based on equation (23), both the heavy and electron energy equations have the same source but with the opposite sign, so we can write Therefore, we have where T g 0 and T e 0 represent the gas and electron temperature at the initial step. The electron energy equation now can be written as The first order differential equation that should be solved is In order to solve the above equation, we need to find and sequentially Finally, the T e has the following form The D is defined based on the initial condition

Appendix D. List of reactions
For reactions without specified references, the cross section is estimated.  [ 47,48] [ 47,48]  • R10: the ionization of He from its electronic excited state is considered both in Liu et al [2] and Murakami et al studies [71]. Instead of a relation for the rate coefficient, a set of cross sections is used in this thesis that is obtained from a private communication in PLASIMO team [1]. • R11: production of H + from the ground state is considered in Liu et al model [2]. In the presented model in this thesis, the cross sections are used to calculate the EEDF and then a tabulated rate coefficient. • R12: ionization of the hydrogen atom from H(n = 2) is also considered in Liu et al model [2], instead of a relation for the rate coefficient, the cross sections for this process are extracted from Janev et al [54] (cross section for the process H(2s) + e → H + + 2e) in this study.  [2]; however, instead of a rate coefficient relation, the cross section from [47,56] is used in this model. • R16: this reaction is also considered in Liu et al model [2]. However, the cross section from [57] is used in this work instead of a rate coefficient relation. The cross section for this collision starts from 30 eV, but the energy differences between right and left-hand side of the reaction (8) is 18.0582 eV. This difference is due to the lack of resolution in the cross section measurement experiment. • R17: This reaction is stepwise ionization of H 2 from a H 2 (b 3 Σ + u ) which represents excited H 2 with energy 8.9 eV. The ground state ionization cross section is shifted by the threshold energy (6.5385 eV) for this reaction. This reaction is not considered in Liu et al model [2]. • R18: this reaction is stepwise ionization of H 2 from a H 2 (a 3 Σ + g ) that represents excited H 2 with energy 11.89 eV. The ground state ionization cross section is shifted by the threshold energy (3.5315 eV) for this reaction. This reaction is not considered in Liu et al model [2].
• R20: ionization of oxygen atom O + from the electronic excited state O( 1 D) is considered in Liu et al model [2]. However, the used rate coefficient is based on a Maxwellian electron energy [72]. The cross section for this Maxwellian EEDF is estimated based on the suggestion of Lee et al study [73], shifting the cross section of the ionization of oxygen atom from its ground state based on the new threshold energy (11.6506 eV). In this study, the shifted cross section is used. • R21: the ionization of O 2 from ground state is considered in Liu et al model [2]. However, the relation that is used for the rate coefficient is based on the work of Gudmundsson et al [72]. That rate coefficient is based on a Maxwellian EEDF. In the presented model, the cross section of this process is obtained from LXCat [47,56] and used instead of a relation for the rate coefficient. • R22: production of O + from O 2 is considered in Liu et al model [2]. The used relation for the rate coefficient is based on a Maxwellian EEDF [72]. In the presented model instead of the relation, the cross sections from Krishnakumar and Srivastava work [58] are extracted and used. The available cross sections for e +O 2 → O + + O + 2e [58] is a combination of this process and e +O 2 → O ++ 2 + 2e. Because m/z = 16 is for both O + and O ++ 2 , the mass spectrometer can not distinguish between these two species. Not that m and z represent the mass and the charge of each species, respectively. Therefore, the presented cross sections are a combination of two processes mentioned above. However, it is also mentioned in Krishnakumar and Srivastava work [58] that the cross sections for the production of O ++ 2 are less than 1% of the cross sections of O + 2 production [74]. In this study, 1% of the presented cross sections for production of O + 2 in [58] are subtracted from the combined cross sections for production of O + and O ++ 2 and used as the cross section for e +O 2 → O + + O + 2e. In the low energy region (<20 eV), the combined cross sections is at the same order of magnitude with 1% of O + 2 production cross sections. • R23-R29: the stepwise ionization for production O + 2 and O + from vibrational excited oxygen molecule is not considered in Liu et al model [2]. In this study, these reactions up to the forth vibrational quantum number are considered. Cross sections for these processes are estimated by shifting the ionization cross sections from the ground state with the new threshold values.
• R30: production of O + 2 from the electronic excited oxygen molecule O 2 (a) is considered in Liu et al model [2]; however, the used relation for the rate coefficient is based on a Maxwellian EEDF assumption [72]. In the presented model, the cross sections of the same process from the ground state of oxygen molecule is shifted based on the new threshold and used to calculate tabulated rate coefficients.
• R31: production of O + from electronic excited oxygen molecule O 2 (a) is not considered in Liu et al model [2]. In the presented model, this process is considered, and the cross sections of the same reaction but from ground state of oxygen molecule is shifted based on the new threshold. • R32: production of O + 2 from the electronic excited oxygen molecule O 2 (b) is considered in Liu et al work [2]. However, the used relation is based on a Maxwellian EEDF [75]. In this study, the cross sections of the ionization from the ground state of oxygen molecule is shifted by the new threshold and used for the calculation of tabulated rate coefficients.
• R33: production of O + from electronic excited oxygen molecule O 2 (b) is considered in Liu et al model [2]. However, the used relation is based on a Maxwellian EEDF assumption [75]. In the presented model, the cross sections of the same reaction but from the ground state of oxygen molecule are shifted based on the new threshold and used for the calculation of tabulated rate coefficients. • R34: ionization process of OH is considered in Liu et al model [2], but the used relation for the rate coefficient is calculated by assumption of Maxwellian energy distribution for electrons. In this model cross sections are extracted from the data of Riahi et al [59] study. Note that the cross sections are extracted from figure 3 of [59], and the points are selected in such a way that data becomes closer to the all references rather than following the curve of the calculated cross sections in Riahi et al [59] study. • R36: The cross sections of OH ionization from electron collision with H 2 O is extracted from Itikawa and Mason [52]. The threshold energy for this collision based on the cross sections data is 17.5 eV, which is in contrast with the needed threshold that should be 18.116 eV [52]. The discrepancy is probably due to the uncertainty in the energy of the electron beam that was used in the cross section experiment [52]. In this study 18.3016 eV is used for the threshold of the cross sections in order to be consistent with the used formation energy of OH, H 2 O, and H. • R37: the cross sections of ionization of O from electron collision with H 2 O starts from 25 eV [52], but the energy differences between right and left-hand side of the reaction (8) is 18.6517 eV. This difference is due to the lack of resolution in the cross section measurement experiment. In this study, 18.6517 eV is used as the threshold of the reaction. • R38: the cross sections of ionization of H from electron collision with H 2 O starts from 20 eV [52], but the energy differences between right and left-hand side of the reaction (8) is 18.6932 eV. This difference is due to the lack of resolution in the cross section measurement experiment. • R39: The cross sections of ionization of H 2 from electron collision with H 2 O starts from 30 eV [52], but the energy difference between right and left-hand side of the reaction (8) is 20.425 95 eV. This difference is due to the lack of resolution in the cross section measurement experiment. 20.425 95 eV is used as the threshold for the cross sections of this process in this study.
E.3. Electron impact: excitation, de-excitation, and dissociation considered in Liu et al model [2], but in the presented model, the cross sections for these processes are used [62] instead of rate coefficient relations. • R70: the dissociation of H + 2 to an ionized and a ground state hydrogen atom is also considered in Liu et al model [2]; however, in this work the cross sections are extracted and used from Tawara et al [62] instead of a rate coefficient relation.
• R71: the dissociation of H + 3 to H and H + is not considered in Liu et al model [2]. The cross sections of this process are extracted from [62] and used in this study. The threshold of this dissociation is 14.441 eV, but the difference between right and left-hand side species in the reaction is 8.9633 eV (8). This energy difference is added as an extra sink and source to the electron and the mixture energy balances in this study, respectively. • R72 and R73: The excitation of oxygen atom to O( 1 D) is considered in Liu et al model [2], but the used rate coefficient is based on a Maxwellian EEDF. In this study the cross sections from Laher and Gilmore work [76] are used (this set of cross sections is available in LXCat website [47,48]). In addition, de-excitation process is also taken into account. • R74 and R75: for excitation of oxygen atom to O( 1 S) instead of the used rate coefficient in Liu et al model [2], the cross sections from Laher and Gilmore [76] are used (this set is also available in LXCat website [47,48]). In addition, the de-excitation process is considered in the present model. • R76: the detachment process of O − is considered in Liu et al model [2]. However, in this study, the cross sections from Deutsch et al work [63] are extracted and used instead of a relation for the rate coefficient. • R77-R84: vibrational excitation of oxygen molecules are not considered in Liu et al work [2]. However, these excitations and de-excitations process up to forth vibrational quantum number are considered in the presented model. Cross sections are gathered from LXCat [47,48] and used for the calculation of tabulated rate coefficients. • R85 and R86: excitation of oxygen molecule to O 2 (a) and de-excitation process are considered in Liu et al model [2]. However, the used relation for rate coefficient is based on a Maxwellian EEDF [72]. In the presented model, the cross sections [47,48] are used for the calculation of tabulated rate coefficients. • R87 and R88: excitation of oxygen molecule to O 2 (b) and de-excitation process are considered in Liu et al [2]. However, the used relation for rate coefficient is based on a Maxwellian EEDF [75]. In the presented model, the cross sections [47,48] are used for the calculation of tabulated rate coefficients. • R89-R96: excitations to O 2 (a) and O 2 (b) from vibrational excited oxygen molecules are not considered in Liu et al model [2]. In this study, this type of excitation is considered and their cross sections are estimated by shifting the cross sections for the same process from the ground state of O 2 with the new threshold for each excitation.
• R97 and R98: transition of O 2 (a) to O 2 (b) is considered in Liu et al model [2]; however, the used rate coefficient is based on a Maxwellian EEDF [75]. In the presented model, the cross sections for this process [47,51] is used for the calculation of tabulated rate coefficients. In addition, the backward transition is not considered in [2], but this process is considered in this study and with the help of superelastic process of BOLSIG+, the tabulated rate coefficients are calculated. • R99: The dissociation of oxygen molecule to oxygen atoms is considered in Liu et al work [2]. However, the used relation for the rate coefficient is based on a Maxwellian EEDF [72]. In the presented model, this reaction is taken into account with the tabulated rate coefficients that are calculated based on cross sections [47,56]. The threshold of the cross sections for this process is 6.12 eV. However, the ΔE of this reaction is 5.115 eV. Therefore, this reaction is also considered as one of the electron dissociative reactions that has an extra sink and source for electron and the mixture energy balances, respectively. • R100: the dissociation of oxygen molecule to a ground state oxygen atom and an excited oxygen O(1D) is also considered in Liu et al [2] work. However, it is not clear from the introduced references [73,77,78] how the used relation for the rate coefficient is calculated. In the presented model a set of cross sections [47,51] is used to calculate the tabulated rate coefficients. This electron dissociation is also considered as one of the processes that has a contribution to gas heating. The ΔE of the process is 7.0823 eV, and the cross sections threshold is 8.4 eV. The extra energy 1.318 eV is considered as an extra sink and source for the electron and the mixture energy balances in this study, respectively. • R101: dissociation of O 2 (a) to oxygen atoms is considered in Liu et al model [2]. However, the used rate coefficient is based on a Maxwellian EEDF [72]. The cross sections for this process are used [51] in this study to calculate the tabulated rate coefficients. This reaction also has a contribution to the gas heating by 0.467 eV. The ΔE of the process is 4.133 eV, but the cross sections has a threshold 4.60 eV. • R102: the dissociation of O 2 (a) to O and O( 1 D) is not considered in Liu et al [2] model, but in the presented model, the cross sections of this process [47,51] are used for the calculation of tabulated rate coefficients. 0.2394 eV is the energy contribution of this reaction in gas heating. The ΔE of the reaction is 6.1006 eV, but the threshold of the cross sections is 6.34 eV. • R103: dissociation of O 2 (b) to oxygen atoms is considered in Liu et al model [2]. However, the used rate coefficient is based on a Maxwellian EEDF [75]. The cross sections for this process are used [47,51] in this study to calculate the tabulated rate coefficients. 0.471 eV is the contribution of this process to gas heating with a ΔE of 3.479 016 eV and a 3.95 eV threshold of the cross sections set. e are measured by Rapp and Briglia [65]. In addition, the threshold value for the polar dissociation is 17.1255 eV. It is assumed in this study that all of the reported cross sections for energy higher than 17.1255 eV in [65] belong to the polar dissociation. In the dissociative attachment the cross sections are reported maximum up to 10 eV [47,56], so this assumption is not too off for the estimation of the polar dissociation cross sections. • R119: the dissociative attachment of O 2 (a) to O − and O is considered in Liu et al model [2]. However, the used relation for the rate coefficient is based on a Maxwellian EEDF [72]. In the presented model, the cross sections for this process are extracted from [66] to calculate the tabulated rate coefficients. • R120: the dissociative attachment of O 2 (a) to O − and O( 1 D) is not considered in Liu et al model [2]. In the presented model in this thesis, the cross sections for this process are extracted from [67] [2] with relations for rate coefficients. In the presented model in this thesis, the cross sections for these processes are extracted from Nandi et al [68] study, but parts of the curve with energies above the threshold values are only taken into account. For production of OH − , there is a threshold due to the fact that the process is not exothermic.