On the ignition kernel formation and propagation: an experimental and modeling approach

The next generation of advanced combustion devices is being developed to operate under ultra-high-pressure conditions. However, under such extreme conditions, flame tends to become unstable and measurement of fundamental properties such as the laminar flame speed becomes challenging. One potential method to resolve this issue is measuring the ignition-affected region during spherically expanding flame experiments. The flame in this region is more resistant to perturbations and remains smooth due to the high stretch rates (i.e. small radii). Stable flame propagation allows for improved flame measurement, however, the experimentally observed kernel propagation is a function of both inflammation and ignition plasma. Therefore, the goal of the present study is to better understand the plasma formation and propagation during the ignition process, which would allow for reliable laminar flame speed measurements. To accomplish this goal, thermal plasma operating at high pressures is studied with emphasis on the spark energy effects on the formation of the ignition kernel. The thermal effect of the plasma is experimentally observed using a high-speed Schlieren imaging system. The energy dissipated within the plasma is measured with the use of voltage and current probes with a measurement of plasma sheath voltage drop as an input to numerical modeling. The measured kernel propagation rate is used to assess the accuracy of the model. The experiments and modeling are conducted in dry air at 1, 3, and 5 atm as well as in CH4-N2 mixtures at 1 atm, and kernel radius, temperature, and mass are reported. The voltage-drop (as a non-thermal loss) is measured to be approximately 330 ± 5 V (dry air at 1 atm) for glow plasma with a large dependency on pressure, gas composition, electrode surface quality, electrode geometry, electrode shape, and current density. The same loss within the arc plasma is measured to be 15 ± 5 V, however the arc phase loss which agrees with arc propagation is significantly higher (∼45 V) which suggest additional unaccounted for phenomena occurring during the arc phase. With these losses, the modeling results are shown to predict the final kernel radius within 10%–20% of the observed kernel size. The difference found between the modeling and experimental results is determined to be a result of assuming that the primary loss mechanism (voltage drop across sheath formation) remains constant for the duration of glow discharge. The discrepancy for arc discharge is discussed with several potential sources, however, additional studies are required to better understand how the arc formation affects the kernel propagation.

The next generation of advanced combustion devices is being developed to operate under ultra-high-pressure conditions. However, under such extreme conditions, flame tends to become unstable and measurement of fundamental properties such as the laminar flame speed becomes challenging. One potential method to resolve this issue is measuring the ignition-affected region during spherically expanding flame experiments. The flame in this region is more resistant to perturbations and remains smooth due to the high stretch rates (i.e. small radii). Stable flame propagation allows for improved flame measurement, however, the experimentally observed kernel propagation is a function of both inflammation and ignition plasma. Therefore, the goal of the present study is to better understand the plasma formation and propagation during the ignition process, which would allow for reliable laminar flame speed measurements. To accomplish this goal, thermal plasma operating at high pressures is studied with emphasis on the spark energy effects on the formation of the ignition kernel. The thermal effect of the plasma is experimentally observed using a high-speed Schlieren imaging system. The energy dissipated within the plasma is measured with the use of voltage and current probes with a measurement of plasma sheath voltage drop as an input to numerical modeling. The measured kernel propagation rate is used to assess the accuracy of the model. The experiments and modeling are conducted in dry air at 1, 3, and 5 atm as well as in CH 4 -N 2 mixtures at 1 atm, and kernel radius, temperature, and mass are reported. The voltage-drop (as a non-thermal loss) is measured to be approximately 330 ± 5 V (dry air at 1 atm) for glow plasma with a large dependency on pressure, gas composition, electrode surface quality, electrode geometry, electrode shape, and current density. The same loss within the arc plasma is measured to be 15 ± 5 V, however the arc phase loss which agrees with arc propagation is significantly higher (∼45 V) which suggest additional unaccounted for phenomena occurring during the arc phase. With these losses, the modeling results are shown to predict the final kernel radius within 10%-20% of the observed kernel size. The difference found between the modeling and experimental results is determined to be a result of assuming that the primary loss mechanism (voltage drop across sheath * Author to whom any correspondence should be addressed. Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
The understanding of the electrical discharge process is important for many practical and fundamental applications including combustion [1], chemical conversion [2], flame stabilization [3], health care [4], manufacturing and surface treatment [5] to name a few. For many combustion systems, the flame is initiated using an electrical discharge to ignite air-fuel mixtures. At the spark location, a plasma channel of ionized gas (a mixture of ions, electrons, radicals, and excited particles) forms. The ignition kernel then propagates, initiating chemical reactions followed by flame front propagation [6][7][8]. This flame front propagates from the surface of the ignition kernel [9,10]. Moreover, for mixtures with a Lewis number (the thermal to mass diffusivities ratio) greater than one, there exists a minimum ignition radius [9,11,12] below which the flame can be quenched. This dependency on the kernel size makes the study of the ignition process crucial in quantifying plasma formation and propagation. A better understanding of the ignition process is also of great importance to improve fuel efficiency and reduce emissions under engine-relevant conditions [13]. However, this goal can be only achieved by a fundamental study on a well-optimized plasma (simple geometry, smooth structure) and under well-defined conditions (known pressure, temperature, gap size, electrode quality, flow motion, discharged energy). The effect of different plasma types and ignition sources on the early kernel geometry and structure has been discussed in various studies. In this context, we review two types of plasma (i.e. thermal and nonthermal) and a few common sources of ignition such as laser, spark discharge [14,15], nanosecond discharge [16,17], and dielectric barrier discharge (DBD).
Thermal and non-thermal plasmas have several differences, among them, the two important ones are: (i) the temperature of electrons in non-thermal plasmas is significantly higher than the thermal plasmas, which in the latter is equal to the temperature of heavy gas particles, and (ii) the main energy transfer mechanism in thermal plasmas is through thermal mechanisms (i.e. increasing the gas temperature via elastic collisions of electrons and heavy gas particles along with relevant transport mechanisms), while in non-thermal plasmas this occurs by both thermal (i.e. fast heat release via collisional quenching of electronically excited states) and kinetic mechanisms (i.e. rapid production of radicals and excited species via direct electron impact). The use of a laser-based ignition could provide a significant benefit from the removal of electrode interference on the ignition kernel. However, it produces a highly wrinkled early kernel and induces gas motion in the direction of the beam [18,19]. The induced gas motion is not ideal because of the potentially non-symmetric shapes and complex geometry, which makes fundamental flame studies challenging. Ignition using non-thermal plasmas such as (1) nanosecond discharge can result in additional wrinkling via generating streamwise vorticity [20][21][22], and (2) DBD may not achieve the temperatures required for ignition or be limited by the geometry of the devices [23]. On the other hand, conventional spark ignition (similar to what occurs in automotive engines) has the potential to produce the most optimized kernel in terms of shape (i.e. spherically expanding) and structure (i.e. smooth surface) [24,25]. However, this also depends on experimental factors such as electrode geometry and discharge waveform. Ziegler et al [26] showed that a relatively large breakdown gap of 2.6 mm produced unusual (non-symmetric) early kernel shapes. Kono et al [27] reported that a toroid shape was produced from the rapid expansion of breakdown [28]. Additionally, Bane et al [28] identified that the shape of the flat-tipped electrode resulted in vortices around the sharp edges of a flat-tipped electrode and consequently a non-smooth kernel. Therefore, it is required to optimize the spark ignition parameters to produce the desired kernel shape and structure.
In this work, the more conventional spark ignition system (i.e. capacitive discharge ignition [29]) is selected. This is because many of the discharge parameters, which affect the typical formation and propagation characteristics, are more easily controlled and understood. The spark ignition at high pressures (>1 atm) is a thermal plasma consisting of three distinct discharge regimes [23]: breakdown, glow, and arc. The breakdown discharge, which occurs on a time scale of 10 ns, is characterized by a sudden drop in high voltage (>1 kV) with a large rise in the current (>20 A). This is a Townsend discharge [30] where the high electric field across the gap causes the electrons to gain enough energy to cause nearly full ionization of the spark gap within a small channel. The power input results in a sudden expansion of the high-temperature channel to form a rectangular and eventually toroidal shape. The breakdown is considered insignificant for modeling the plasma discharge when the gap is sufficiently small (<0.4 mm). This is because the total energy deposited by the breakdown (and consequently the total effect on the kernel size) scales linearly with the gap size. The breakdown is of course necessary to precede and initiate the formation of the glow and arc plasmas. Glow discharge is initially assumed to be normal [26,31]. Normal glow is characterized by a minimum constant current density through the surface of the cathode (i.e. highvoltage electrode) where properties of the plasma such as the sheath formation are independent of the current [32]. Often glow discharge is referred to as non-thermal plasma, but this is highly dependent on the operating current [33]. In applications such as electronics manufacturing [5], glow discharge can be generated with low current resulting in non-thermal plasma. However, glow discharge at high currents under high pressures, which is relevant to the conditions of this work, becomes thermal plasma where local thermodynamic equilibrium (LTE) is a valid assumption [33,34]. The arc discharge is characterized by a sudden increase in current density. Often, a current threshold where the transition to arc occurs is discussed [12] but this can be misleading as this threshold is heavily dependent on electrode material and surface quality [26]. Transition to arc discharge primarily depends on the temperature of the electrode [32] where a polished electrode prevents localized heating on the surface, extending the glow range to higher currents [35]. As a result, cooled electrodes are shown to maintain glow discharge up to 14 A at pressures of 13 atm in Cobine [32]. If the surface temperature increases, the glow discharge will eventually transition to the thermionic arc discharge [36]. Experimentally the transition between modes can be observed as a discontinuity in the voltage plot with arc discharge at a much lower voltage. The transition in discharge regimes can create non-ideal conditions for kernel formation and propagation, illustrating the need for experimental control over the spark discharge. Understanding how these discharge regimes convert the electrical energy to thermal energy in the surrounding gas is important for the development of new and improved models to better describe the complex interaction between ignition and combustion.
There exist several studies for modeling high-pressure spark ignition. Maly [9] performed multi-zone modeling to study the minimum ignition energy (MIE) where the radial temperature profile was calculated. Since Maly [9] worked backward from the flame to determine the minimal plasma characteristics for ignition, the research lacks the information to produce quantitative plasma propagation data. Lim et al [37] modeled both the breakdown expansion and the diffusive glow-arc discharge. The mass diffusion for the glow regime was expressed based on a correlation between the spark energy and kernel velocity, requiring adjustable constants to scale the results. In the current work however, it is assumed that the thermal conduction from the hot plasma channel into the surrounding gas is the main mechanism responsible for mass entrainment through the kernel boundary, providing a more physical definition of kernel propagation based on available transport properties. Kravchik et al [10] and Sher et al [38] provided detailed early spark propagation models with a primary focus on the breakdown stage. Since the primary ignition energy was assumed to arise from breakdown, no information regarding the glow and arc regimes was available in their studies. Kim and Anderson [39] presented a model that predicts kernel velocity utilizing experimental data at enginerelevant conditions. Their model assumed that the characteristic spark velocity is the ratio of the change in spark voltage to the electric field at the spark gap. This suggests that the kernel velocity is only a function of the plasma size (as opposed to the thermal expansion-diffusion effects which are a result of plasma discharge). As a result, the accuracy of the results is only relevant to engine results that are cycle averaged and turbulent, not the desired time-resolved propagation under well-defined conditions. Eisazadeh-Far et al [40] developed a thermodynamic model of spark kernel propagation under a constant mass assumption. This assumption did not account for the effects of thermal conduction, resulting in a kernel temperature much higher than the experimental data. Moreover, the ignition model was heavily dependent on the initial conditions, producing highly uncertain modeling results. Recently, Lu et al [41] developed a computational fluid dynamics model for the determination of MIE. However, the experimental plasma propagation data as well as the effect of different discharge regimes were not integrated into their modeling. This is very important as the thermal efficiency of a spark can vary greatly, sometimes randomly, with time as the discharge regime changes. Further analysis is required beyond these models to provide more accurate time-resolved propagation as well as mass entrainment for the ignition kernel in the glow and arc regimes. This will also help, although not within the scope of this current work, to better understand the transition region between plasma (i.e. electrical discharge driven) and a sustained flame (i.e. chemistry driven).
The goal of the present study is to develop quantitative modeling tools for plasma formation and propagation under glow/arc discharge regimes and thus to explore a well-defined plasma propagation in non-flammable gas mixtures (e.g. dry air and nitrogen-diluted methane) at high pressures (⩾1 atm). This is accomplished through two modeling approaches; a zero-dimensional (0D) thermodynamic model and a onedimensional (1D) direct numerical simulation (DNS) model are used to predict the ignition kernel propagation. The modeling tools are compared to available experimental data to assess accuracy. In addition to modeling, extensive experimental study on the phenomena regarding the structure and properties of the plasma (sheath formation and losses) is presented as well as parametric study of significant electrode properties to better standardize and create repeatable experimental results.

Optical setup
Experiments are conducted using a constant volume combustion chamber (CVCC) shown in figure 1. The CVCC has an internal volume with a diameter and height of 13.3 cm and uses polished quartz disks on either end to provide optical access to the ignition kernel. The experimental diagnostic includes a linear Toepler-type Schlieren system where the layout is schematically shown in figure 1. The light source is a red 625 nm Thorlabs led (model M625L4), which is collected with an initial lens (Lens 1) and focused through a 1 mm diameter pinhole to remove excess stray light and make the light as close to a uniform point source as possible. After the pinhole, a planoconvex lens (Lens 2), with a 20 cm focal length, is used to pass collimated light through the plasma discharge, and another plano-convex lens (Lens 3), having a 40 cm focal length, for converging the light over the camera sensor. The Photron FASTCAM SA-Z high-speed camera is utilized to take images during the plasma formation (180 000 fps at 128 × 128-pixel resolution).

Electrical setup
The spark is created, at the center of CVCC, with two extended stainless steel wire electrodes [42] of 0.5 mm in diameter with tips polished to a point (conical tip) using a range of polishing methods. The spark is measured using a NorthStar PVM-4 high-voltage probe and a Pearson 6595 current coil. For kernel propagation experiments, the spark gap is held constant at 0.2 mm while this gap is changed between 0.2-4 mm for sheath measurements (see section 2.5). The spark discharge is generated by providing a high-voltage pulse to the electrodes initiated by passing a high current pulse through the primary windings of an automotive ignition coil (i.e. 5527 automotive coil). The result is a sinusoidal-like current shape driven by the second-order circuit (shown in figure 1) similar in nature to the referenced patent [43]. By changing the initially stored voltage in the capacitor, the magnitude of the glow/arc current is changed. The secondary coil sees a current pulse on the order of 1 A. The duration of the pulse is defined by the inductance and resistance of the ignition coil (∼0.7 ms for the ignition system used in the present study). The system can potentially produce multiple discharges with a single trigger event (because of underdamped oscillations); therefore, a diode system is used to suppress other discharges by passing the second discharge to the ground with an avalanche diode system.

Propagation measurement and uncertainty
Schlieren images are used to visualize kernel formation and growth. In the case of the ignition with a non-flammable mixture, the kernel propagation is driven by thermal conduction from the plasma region (for glow and arc discharge) [26,37]. Since the largest change in density occurs near the ambient temperature, nearly the entire volume of gas heated by the plasma (i.e. ignition kernel) can be visualized by the Schlieren method [44]. The measurement of the ignition kernel involves post-processing in MATLAB to detect the kernel edge and total volume. Because the shape of observed kernels is round (no complex overlapping geometries are considered) the radius of a sphere with the equivalent volume is considered for analysis [12].
The approach to estimating the error associated with the kernel radius measurement considers three uncertainty sources including the edge location, the pixel scaling factor, and the location of the kernel center. Each of these sources is assumed to be uncertain to the size of a single pixel. The maximum and minimum radius measurement based on the variation of these parameters is calculated while processing the nominal radius size. For example, if the scaling factor is 13 px mm −1 then the uncertainty would include the change in kernel size if the scale were changed to 13.5 and 12.5 px mm −1 . Each component of uncertainty is summed in a way that the uncertainty represents the maximum potential error in the radius measurement.

Discharge power and uncertainty
The voltage measurement is captured using a NorthStar PVM-4 probe that has a 0.4% uncertainty, and the current is captured using a Pearson 6595 current coil which has a 1% uncertainty. The probe measurements are then captured using a NI-9222 module, four-channel analog input, with 16 bits of vertical resolution and a sample rate of 500 kHz. The inductive measurement of current is affected by signal droop which causes the actual value of current to drop over the full duration of the spark. This is corrected by using equation (1) where i Meas is the signal received by the device and τ is the characteristic time constant of the current coil: The collected voltage and current data are used to find the discharge power as equation (2), where the V and I are the measured voltage and corrected current andĖ elec is the total electrical spark power within the electrodes gap, Resistive circuit losses are excluded from this measurement through the placement of the voltage probe as close to the spark location as possible and with the use of low (<0.5 ohm) resistive electrodes. Any measured resistance between the measurement location and the ground is removed with the use of ohm's law. The thermal powerĖ th which is used as an input to the predictive models, is described in equation (3), where V fall is the voltage drop across the sheath formation (see Typical discharge characteristics for voltage, current, power, and energy are shown in figure 2. The voltage waveform includes the breakdown event, which is visible as a vertical line at the start of the glow regime in figure 2(a). The glow plasma occurs for ∼0.2 ms (in this example) and at ∼350 V. Arc plasma is observed as the transition from the high voltage glow discharge to <100 V. The spark ends at ∼0.7 ms as this is the time when the current drops to 0 A. The current in this case reaches an amplitude of ∼0.8 A, which is a function of the initial capacitor energy. The power and energy of this spark event are shown for both the original and the thermal cases in figure 2(b). The thermal case, as mentioned earlier, is the case in which the voltage drop is removed, and its associated energy is released into the gas. The overall efficiency of arc plasma vs glow plasma can easily be observed in figure 2(b) where a significant ratio of glow power is lost as compared to the arc thermal/original power.
The uncertainty of the electronic measurements is determined using the total system uncertainty calculated as, where gain and offset uncertainties are 0.2% and 0.1% respectively. The u P is the uncertainty associated with the probes which are 1% for the voltage probe (u U ) and 2% for the current probe (u I ). The uncertainty of the power measurement is then calculated as, The uncertainty of the total spark energy can additionally be propagated over time; however, this measurement is not expressly used in this analysis. Additional detail on electrical uncertainty can be found in our previous publication [45].

Sheath voltage drop measurement methodology
Since the voltage drop (often called cathode fall) can change greatly [32] with the experimental conditions (e.g. composition, geometry, pressure), it is important to measure this value. The basic voltage drop measurement method for high-pressure glow plasma employs a zero-gap extrapolation [31,46] in which the inter-electrode distance is mechanically varied to observe the change in voltage in terms of the plasma channel length. The measurement conventionally consists of at least one flat plate electrode for the cathode such that the geometry  has little effect on voltage regardless of the discharge current.
In the case of the flat plate cathode, the plasma channel length is always equal to the gap distance. However, for the tip-to-tip electrode configuration, used in this paper, the plasma channel length changes with current (during the ignition process) as shown in figure 3(a).
Since the measurement of the plasma length in such a transient process using Schlieren imaging is almost impossible, it is required to have a reference case in which the plasma length can be simply obtained. The best reference case is where the plasma channel stays between the tips and its length is equal to the known value of gap size (d 0 in figure 3(a)). Although maintaining the plasma channel between the tips is experimentally challenging, it is assumed that such a condition can be satisfied when the current is very small (∼0 A). To accomplish this, an additional extrapolation (i.e. zero-current extrapolation) is required to account for the tip-to-tip electrode geometry. The zero-current extrapolation is only necessary for the glow regime as this follows a normal change in length. Therefore, the voltage drop measurement method for the tip-to-tip electrodes contains two steps. First, for each gap size, the measured voltage is plotted against the measured current. Then, the discharge voltage is extrapolated to find the zero-current voltage. Second, the zero-current voltages are plotted in terms of gap size and the voltage drop can be found by extrapolating to zero gap as seen in figure 3 To assess the voltage-drop for arc discharge, the minimum voltage (when the arc is directly between the tips) is taken for several gap lengths. At least five duplicate tests are taken for each condition to ensure repeatability and assess the variability and uncertainty in the results. The estimated error of gap voltage measurements incorporates the deviation in voltage of the redundant tests at each fitting point. The final voltage drop measurement uncertainty estimates the error from the linear polynomial fit.
The effect of surface roughness, pressure, and composition on the voltage drop is also considered experimentally. Two polished levels of 5000 and 46 000 grit are considered in this work. The 5000-level grit is sandpaper while the 46 000-level grit is in the form of 0.5 µm diamond paste. In addition to surface roughness, the effect of heavy use (100+ ignition events) on the electrode is considered. The composition will naturally have a significant impact on the sheath formation and is observed for the CH 4 -N 2 mixture (5% CH 4 and 95% N 2 ) and dry air (21% O 2 and 79% N 2 ). The effect of pressure on glow discharge is observed in the air up to 3 atm while the arc is observed up to 20 atm. Voltage drops for both glow and arc discharge are found utilizing different inter-electrode distances between 0.2-4 mm. Sheath measurement conditions are summarized in table 1. The measurement results and the effect of experimental conditions on voltage drop for both glow and arc will be discussed in section 4.2.
The range of data presented for the plasma sheath is selected to cover the experimental and modeling conditions. Arc sheath measurements were done up to 20 atm to support future research. For modeling purposes, various mixtures of arc and glow discharge are available from 1-5 atm making this range a suitable testbed for the proposed modeling (Note, measurement of sheath properties for glow discharge were only possible for pressures up to 3 atm as arc interference became too large to achieve reliable glow data above this pressure).

Modeling approaches
A thermodynamic model and DNS model were developed for predicting the radius and temperature. The governing equations and methods are explained in the following sections. For comparative, purposes both models assume that the primary non-thermal losses arise from sheath formation at the electrodes. The net thermal energy for both models is given as,Ė whereĖ th represents the net thermal power deposited within the gas for the formation of the ignition kernel from equation (8).Ė Cond andĖ rad are the heat conduction [47] to the electrode and thermal radiation loss [48] from the ignition kernel, respectively below, where ε is the net emission coefficient for air [48], A cond is the surface area of the electrode contacted by the kernel, γ is the specific heat ratio, T e is the electrode temperature, α is the thermal diffusivity, and δ is the thermal boundary layer thickness. For the conditions observed in this work, the contribution of both conduction and radiation losses are orders of magnitude smaller thanĖ th . Radiation losses only become significant at gas temperatures over 6000 K. Additionally, conduction loss may only be relevant with a significantly larger electrode [12,49]. On the other hand, Arkhipenko et al [50] indicated that thermal energy conducted through the electrode is almost entirely described by sheath losses, even for a large flat plate at a steady state discharge. Despite being negligible, in the conditions presented in this work, these terms are retained to ensure the model remains relevant in other conditions such as high pressures, high currents, and large electrode geometries in future studies. Details regarding the application and scope of both models are found in the following sections.

Thermodynamic model
The thermodynamic model solves the energy and mass conservation equations to find the radius, mass, and averaged temperature of the ignition kernel. A schematic of the model is shown in figure 4 where the dashed control volume is defined as the ignition kernel edge with energy entering and exiting the system uniformly over the entire volume.
Assumptions for the thermodynamic model include: (a) unburned and burned gases behave as an ideal gas; (b) it is assumed that the kernel is approximately spherical; (c) since the relaxation time scale of different energy modes compared to plasma expansion time is small (∼10 −9 s), all species are in LTE [33,34,38,51,52]; (d) the modeling begins shortly after the breakdown; (e) the only mass gain occurs from thermal conduction from the ignition kernel; (f) the radial temperature profile can be assumed linear; and (g) spark energy losses result from thermal radiation, heat conduction and sheath formation.
The model begins with energy balance, where m is the mass, h is the enthalpy,ṁ is the mass flow rate, subscript b represents the gas heated by the plasma channel, subscript u refers to the surrounding gas, and i is the entering conditions. By solving the energy balance for the change in kernel temperature under constant pressure conditions, we find, whereṁ b =ṁ i (the mass entering the ignition kernel) is a result of thermal conduction, T is the temperature and c p is the mixture specific heat at constant pressure. This is the form of the first governing equation that can be integrated to find the average temperature of the kernel over time using the deposited thermal energy. The ideal gas equation can be used to solve for the kernel volume. All thermodynamic properties are calculated for a wide range of temperatures from Askari et al [53,54]. It should be noted that the ideal gas equation can be differentiated and combined with equation (10) to directly solve for the kernel velocity shown in equation (11), where A is the kernel area and Ω is a tabulated thermodynamic property from [53,54]. However, it is more beneficial to solve directly for the temperature: Thermal conduction through the gas is assumed the primary mechanism for kernel expansion. Therefore, from the energy balance perspective, the energy associated with the mass gain by the kernel propagation should be equal to the energy flux by thermal conduction through the Fourier heat equation which results in, where k is the thermal conductivity of air [55], r is the kernel radius, and A is the kernel area. Since the thermodynamic model is 0D, a linear temperature change is assumed across the radius. The calculated kernel temperature T b is assumed to be the spatially averaged temperature of a spherical kernel. The radius that would correspond to the average of a linear profile would be 1/4th the radius from the kernel edge. This results in an approximate temperature gradient to use for equation (12).
Equations (10) and (12), are sufficient to fully define the kernel where the kernel radius can be found using the calculated mass and temperature with the ideal gas equation. The equations are first order discretized and since the temperature is inversely proportional to the mass, a logarithmic timescale is used to ensure small initial radii can be used without temperature becoming unbounded (a result of small initial mass).

DNS model
A DNS model was developed by modifying the existing Lagrangian transient 1D reacting flow code (LTORC) [56]. The use of Lagrangian coordinates eliminates the non-linear advection term from the conservation of energy and species equations which can cause numerical instabilities [57]. The plasma is assumed to be in thermal equilibrium so that the electrons and heavy molecules are at the same temperature. Thermodynamic and chemical kinetic properties of all species are computed using Cantera [58]. Transport properties of heavy molecules are also computed via Cantera, while electrons require special treatment. It is assumed that the electric field produced by the discharge primarily exists in the axial direction of the electrodes and that the strength of the electric field in the radial direction is negligible. A strong electric field is required to hold a net charge in the plasma, so the plasma should maintain a neutral charge in the radial direction. The absence of any significant electric field means that heavy molecules diffuse in a Fickian manner, radially. Any deviation from neutrality will produce an electric field causing charged particles to drift. The large mobility of the electron compared to the heavy molecules means electrons are most affected and will be transported to maintain charge neutrality. In effect, the diffusive flux of electrons is computed to cancel the charge flux of ions as, where the subscripts e and i are respectively the electron and species i, V i is the diffusive velocity, X i is the mole fraction, and c i is the charge. The governing equations of LTORC are given as: and the boundary conditions are given as: where r is the spatial coordinate, ψ is the lagrangian coordinate, ρ is density, P is pressure, c p is the specific heat at constant pressure, T is the temperature, t is the time, λ is the thermal conductivity,ω i is the mass production rate of species from reactions, ∆h f,i is the standard heat of formation of species i, Y i is the mass fraction of species i, andĖ th is the net thermal energy deposited per unit volume. The spatial distribution of the heat source term is defined in equation (20) as used by [59,60] where r dis is the radius of the discharge kernel, The thermal conductivity of the gas is approximated to be that of an air plasma for all cases. The thermal conductivity of air is computed as a function of temperature [55]. The method of lines is used to separate the spatial and temporal terms and the spatial terms are approximated using 2nd order finite differences. Equations (14)-(19) form a system of integro-differential equations, where equation (14) needs to be integrated to solve for r before the rest of the governing equations are evaluated. This system of integro-differential equations is solved using CVODE [61]. Additional details regarding LTORC can be found in [56].
A constant Lagrangian 1D spherical grid was used initially ranging from r = 0 to r = 2 cm with N = 1000 grid points. To resolve the plasma physics well, grid points are focused at r = 0. For a grid point j, the distance between j and j + 1 is denoted by ∆r j . The distance between the first two grid points is made to be 20 times smaller than the last two grid points which is equivalent to ∆r N−1 /∆r 1 = 20. The distance Table 2. Subset of reactions from [62]. between grid points is gradually increased with r such that ∆r j +1 /∆r j = (∆r N−1 /∆r 1 ) 1 N−1 . A mixture-averaged transport model is used to compute the diffusive velocity of heavy molecules. Two different chemical kinetic models were used, one for air and another for N 2 /CH 4 mixtures. For air, a subset of reactions from the model by Park [62] was used as listed in table 2. For N 2 /CH 4 mixture, the entirety of the FFCM1 model [63], a subset of the Park model [62] as listed in table 2, and a subset of the Eichwald et al model [64] as listed in table 3 were combined. All reaction  (21): where k (T) is the rate coefficient in cm 3 mol −1 , B is a constant, α is the temperature exponent, E a is the activation energy in cal 1 mol −1 , and R o is the universal gas constant. The LTORC simulations are initialized after breakdown and the subsequent emitted shockwave. The state of the plasma in the plasma kernel immediately after the breakdown is estimated using Cantera [58]. This is done by adding the breakdown energy to a small cylindrical volume of gas that is 0.2 mm in diameter and 0.2 mm tall and equilibrating it isochorically.
The significant discontinuity in pressure after breakdown forms a shockwave, which expands the plasma. The hydrodynamic module of AMRVAC [65] is used to solve the compressible Euler equations and simulate the change of pressure with time of the plasma kernel as it mechanically equilibrates. It is assumed that the primary cause of pressure change in the plasma during this expansion is due to the shockwave and that any change in pressure due to changes in the chemical state of the plasma is negligible. Starting from the state of the plasma after breakdown, a 0D reactor [58] is used to estimate how the thermochemical state of the plasma kernel changes during expansion by using the change of pressure vs time computed by AMRVAC. The conservation of energy and species equations of the 0D reactor are given as below, The radius and boundary layer thickness of the plasma kernel after the expansion is estimated from the final density profile of the AMRVAC simulation. These two values are used to interpolate between the post-shock state of the plasma and the state of the ambient gas to form the 1D initial conditions for LTORC. Figure 5 depicts the snapshots of plasma kernel growth under different experimental conditions (i.e. gap size, surface quality, and gas mixture composition). Each letter represents a single experiment for a total of five experiments shown. Note that, for the 'prepared' electrode (polished), such as in figures 5(b), (c) and (e) where pure glow is observed, the results (shape and size) are highly repeatable. Much of the random nature of spark discharge comes from arc discharge. The mixtures are selected to be outside of the flammability limit to only study the plasma propagation by isolation from combustion chemistry. The current flows from the cathode (high-voltage electrode) to the anode (grounded electrode) forming an ionized region in the gap between the electrodes. Depending on the conditions, this ionized region (referred to as 'plasma channel' in this paper) is represented by either pure glow, pure arc, or combined glow-arc regimes. The hightemperature plasma channel heats the surrounding gas via thermal conduction, generating a so-called 'ignition kernel' that extends beyond the plasma channel.

Spark morphology
It is important to distinguish between the plasma channel and ignition kernel as shown with purple dash-dotted and red dashed regions in figure 6, respectively. The plasma channel can be observed in some experimental conditions via the emission of bright light at the center of the ignition kernel. This bright emission is more pronounced in the absence of oxygen in the mixture as shown in figures 5(d)-(f) at t = 255 µs. The ignition kernel, which represents the density change between the heated region and surrounding gas, can be easily observed by Schlieren imaging. There also exists a thin layer (not visible in figure 5 because of the small thickness) at the surface of the cathode and anode with a large population of positive ions [66]. This layer, which is called the sheath, is responsible for a substantial electric field at the surface of the electrodes resulting in a large voltage drop. The energy dissipated within this sheath is not converted to thermal energy within the gas but is rather lost to the electrode [26,50].
Among the parameters affecting the spark morphology, electrode geometry, surface quality (i.e. roughness), gas pressure, ignition circuit (i.e. current discharge profile), and gap size are of great importance. As shown in figure 5, the gap size and surface quality can have a significant influence on plasma regime and spark morphology. Gap size, in general, increases the path-length of the plasma current and as a result the total thermal energy dissipated. However, breakdown (not included in the modeling/experimental results) sees a far greater thermal efficiency and expansion which results in the toroidal geometry in figure 5(a). The energy released during the breakdown is proportional to the gap size where gaps with a smaller size (∼0.2 mm) are found to minimize the breakdown effect as seen in figure 5(b). For glow and arc discharge the path of plasma is not only restricted to the gap and as a result is not linearly proportional as compared to breakdown for small gaps. For the case of unpolished electrodes in figure 5(d), the majority of the discharge occurs within the arc regime in contrast to the glow regime found in figure 5(e) with polished electrodes. The plasma is typically observed (via the voltage waveform and optically) to start with glow, transition to arc (if at all), then transition back to glow. If the initial glow phase is less than ∼20% of the discharge; the shape of the kernel is generally challenging to measure as apparent in figure 5(d).
Although the gas composition affects the kernel size (from changes in mixture properties and thermal efficiencies), no noticeable changes have been observed in the spark morphology (structure) as illustrated in figures 5(b), (c), and (e). The change in composition will change the characteristic discharge voltage as a result of these changed properties, however, the behavior of the kernel remains mostly the same. The kernel propagations shown in figures 5(a) and (d) are not preferred for the spherically expanding flame measurements because of the abnormal geometry. Therefore, the best practice is to adjust the gap size and surface quality to produce a more symmetric and uniform kernel shape as seen in figures 5(b), (c), and (e). This well-defined kernel morphology usually occurs in the volumetric glow discharge. However, maintaining the discharge in a pure glow regime is challenging particularly as pressure increases. Some arc discharge can retain a decent shape if arc occurs for a short duration.
The glow regime fully surrounds the cathode, unlike the arc that is highly localized. The constant current density, in a glow regime, results in a large discharge area on the cathode surface as shown in figure 5(a). The size of this area varies according to the supplied current. The formation of such discharge area results in a high electric field at the cathode surface causing the characteristic high voltages for glow plasma (>200 V). However, the arc regime is characterized by a sudden change in current density where the discharge area at the cathode is significantly decreased and forms a so-called discharge spot as depicted in figure 6(b). Such discharge transition, from a larger to a smaller area, results in a significant reduction in voltage (<100 V). The discharge spot on the cathode is in the same order of size as the anode however the location of the spot can vary widely. The ignition kernel is well defined through the schlieren method however the plasma region is only sometimes visible through intense direct radiation seen in figures 5(d) and (e).    is shown where the circles represent the average zero-current extrapolation for a given gap. These data points are then extrapolated to the zero-gap (depicted by stars) to find the voltage drop for each experimental condition.

Voltage drop
It is important to note that even for a small variation in experimental conditions such as the surface roughness of the electrode, there exists a significant change in the discharge characteristics. For instance, varying the electrode polish from 5000 grit (blue dash-dotted line) to 46 000 grit (black dotted line) results in a change in the overall discharge voltage by ∼5 V. The change produced by surface roughness is thought to be due to changing the work function of the material [67,68] where the voltage drop is expected to vary linearly with the work function [32]. This change in surface conditions will affect the field emission of electrons from the electrode surface and has been accounted for in another plasma modeling such as Venkattraman [36].
Heavy electrode usage, without re-polishing/cleaning, was also observed to result in a measurable change in the voltage drop (red dashed line). Heavy-use electrode shows an increase by ∼10 V compared with the similar 5000 grit electrode under normal-use (blue dash-dotted line). This is expected to be the combined effect of the change in the electrode surface roughness and a change in electrode surface composition (formation of oxide layer). The oxide layer, discussed in the following section, results in surface deformations as well as layers that can flake from the surface (for the thickest oxide formations) [69]. This will expose larger surface areas to the plasma. Both the mechanical roughness [67,68] and oxide composition [70] will change the electrode work function however the actual contribution of each of these effects toward the total change is beyond the scope of the work. Gas composition from dry air (blue dash-dotted line) to N 2 -CH 4 (pink solid line) creates a significant increase in the sheath voltage drop of about ∼25 V. Note that, in the N 2 -CH 4 case, the polish level of the electrodes is kept constant at 5000 grit.
The voltage drop reported in this work is around 330 V for stainless steel in the air. Note that heavy-use electrode in air is used for modeling as this appears to more closely represent electrode condition when the spark propagation data was recorded. Comparatively, in literature, the voltage drop reported by Cobine [32] for iron electrodes in the air (closest experimental condition) is 269 V. Swett [46] reported a 290-300 V drop for stainless steel electrodes (with cadmium) in air. It is suggested that the differences between these reported values are a result of many factors including electrode geometry, electrode material, electrode surface quality, as well as the current magnitude. The most significant effect is likely a result of the last factor where Cobine [32] reported a normal glow discharge while the values in Swett [46] may have abnormal discharge on account of the short pulse discharge observed. The pulse observed in this work reaches much higher currents and consequently finds a much larger value. Cobine [32] also reported a wide range of pure gas compositions, which provides insight into the effect gas composition has on the sheath. When comparing discharge in air vs N 2 , for copper electrodes, a significant reduction in the voltage drop, 160 V, is observed. However, discharge in a heavier molecule such as CO 2 is shown to increase the voltage drop over the air by 90 V. For the experimental condition observed in this research with O 2 removed, a reduction in the voltage drop is expected, however, just a small addition of CH 4 proves to reverse this effect.
The effect of pressure on voltage drop is shown in figure 7(b) for a heavy-use electrode. To have a sustained glow discharge, the maximum pressure is limited to 3 atm. Increasing pressure results in a decrease in the expected voltage drop since the current density at the cathode will increase as the length scale decreases (current is independent of pressure). This result shows that voltage drop is inversely proportional to the current density for the conditions observed. Additionally, increasing pressure leads to arc transition.
Sheath measurement for arc discharge over a wide range of pressures and different electrode surface roughness versus gap distance is shown in figure 8. As previously mentioned, the voltage drop corresponds to the zero-gap extrapolation (the circles in figure 8). As can be seen, no significant change in voltage drop is observed for different pressure and electrode roughness. Therefore, it is concluded that the value of voltage drop for arc discharge mode is around 15 ± 5 V and remains constant as expected since arc discharge is a result of thermionic emission [71] where the high thermal load on the material at the small discharge spot results in a more uniform condition (i.e. surface melting).

Effect of electrode surface morphology on voltage drop
The surface quality of the electrode was observed to effect the discharge characteristics of glow discharge (shown in figure 7). Qualitative understanding of this effect was explored by taking optical microscope images using a Keyence VHX-7000. The surface of the most polished electrode under normal-use (46 000 grit or 0.5 µms) is shown in figure 9(a), while the 5000 grit polish is shown in figure 9(b). The highly polished surface was observed to be near the optical limit (for observing surface variation) with only a few scratches visible. Comparatively the lower polish observes a significant increase in roughness. All electrode images (except for figure 9(b)) used a purely optical image. Figure 9(b) is observed with an imaging technique that uses multiple lighting angles to better indicate surface texture through shadows. The color of both figures 9(a) and (b) are observed to be uniform since they represent normal-use electrodes.
The heavy-use electrodes for both used polish levels are shown in figures 9(c) and (d). The full-color images are observed as this provides interesting information regarding the formation of oxide on the electrode surface as mentioned in the previous section. The color change observed after multiple discharges is assumed to be the buildup of oxidation on the surface. Friburg et al [69] provided data regarding the reflective colors of oxide on stainless steel. Based on this description the surface is fully coated with oxide at different layer thicknesses between 20 to 300 nm. This will completely change the composition observed by the plasma, which, alone, is enough to affect the work function and consequently the sheath formation. In addition, the roughness (while only observed comparatively) is seen to increase (especially the case for the 46 000 grit polish but also for the 5000 grit). This is inferred through observation of scanning electron microscope images provided by Friburg [69] which show the texture of the oxide layer is material that is flaking from the surface. This surface modification will likely produce an increased surface area over the freshly polished electrode.
The change from the normal-use 5000 grit ( figure 9(b)) to the heavy-use ( figure 9(d)) has a similar oxide formation to the 46 000 grit (figure 9(c)) where the color (and expected thickness) are comparable. The major difference is that some of the original surface variations from the 5000 grit are still visible where the oxide layers are thin (light blue to orange/red). On the other hand, the thickest oxide (∼300 nm thick -white/grey color) is observed to completely obscure the markings of the 5000 grit. This would suggest that after a significant number of tests both electrodes will eventually reach the same condition (roughness and oxide thickness) regardless of the initial polish level.   Table 4 provides the conditions of the presented experimental and modeling results. These conditions were selected to show a wide range of observed propagation characteristics that are a result of discharge mode (arc/glow), gas pressure, gas composition, and discharge current. The stored voltage (on the primary capacitor) represents the maximum current passed through the gap with 160 V resulting in a maximum of ∼0.8 A (see figure 2) and 100 V resulting in a max current of ∼0.5 A. The total possible energy stored by the ignition system is found to be on the order of several joules but there are many inefficiencies within the ignition circuit (i.e. transformer coils) that make this determination pointless for the spark analysis since the spark gap only sees a fraction of this energy. The voltage drop is the measured constant loss used for each spark condition, naturally the value changes depending on if the plasma is arc or plasma (section 4.2). The type of plasma for each case is also shown as a percentage, where the duration of the spark in arc discharge is presented against the total duration of the spark time (∼0.7 ms). The discharge radius is a parameter used in the DNS model which represents the size of the region in which the spark energy is dissipated (given the model is one dimensional). The initial conditions (r i and T i ) are not reported as the effect on the overall result is minimal (see the following section). In general, the measured initial radius within the uncertainty of the measurement can be used along with a physical temperature value (300-1000 K). In the following section, it is shown that the power input is the most important modeling factor. Other effects such as thermal conductivity terms or the boundary conditions show comparatively little change in the results. For this reason, the most important parameter to consider in the analysis is the sheath loss applied to the input power data. In addition, arc discharge represents a significant physical change to the experimental conditions (dimension and location of the power dissipation by plasma changes, affecting the dimensionality of the system), which is challenging to represent with these simplified models. The assumptions and results of the glow plasma are shown to be in greater agreement with the experimental results, once arc discharge is introduced, results are observed to significantly deviate from the experiment such that a 15 V sheath loss greatly overpredicts the experimental conditions (this is discussed further with the results in figure 15). Two model conditions are considered for the arc voltage drop: (1) 15 V, which would be the measured loss and on the same order of magnitude as other spark models previously discussed, and (2) N/A condition that represents the entire spark energy removal if the arc does not affect the kernel boundary (this is equivalent to a voltage drop greater than ∼80 V). This will illustrate the potential range of solutions available during arc discharge assuming the loss is greater than the assumed voltage-drop. Of course, the entire removal of the arc is not physical so the case for the pure arc is observed with a larger drop which suggests there are additional phenomena that must be accounted for during the arc beyond the initial assumptions.

Model sensitivity.
Logarithmic sensitivity analysis (S = d (ln (r final )) /d (ln (x i ))) is performed for both models to assess the uncertainty of model parameters and inputs. The sensitivity determines the relative contribution of the change of a given model input, x i , compared to the change in the final calculated kernel radius. The sensitivity is calculated using results from Case #1 and considers some of the major model inputs including; the initial radius and temperature of the kernel (T i , r i ), the conductivity of the gas, k, (given mass gain is driven primarily by thermal conduction) the boundary conditions of the kernel (T b or ∇ρ, discussed further in section 4.4.3), and the value of sheath loss (variation in input energy).
The initial inputs, and r i are shown to have an insignificant effect on the final radius, this suggests the modeling solution converges to a physical solution, where only the model results immediately following the initial conditions are affected by initial input variation. The sensitivity of thermal conductivity is important to determine since this transport property is utilized from literature for air and is not modified for alternate conditions. The sensitivity is larger but is not expected to be a major source of error. Future work may consider more detailed thermal conductivity should it be necessary.
The boundary criteria of the kernel constitute an interesting parameter to consider because there are to some degree associated with how each model physically interprets the observed kernel. In the 1D model, a temperature profile is calculated with respect to radius with the edge of the kernel represented by the max change in density ∇ρ max . This value will be physically representative of the Schlieren system but there is no guarantee that this definition is the location measured when observing the kernel. The dependence of the final radius on the selected boundary condition can be determined if the boundary is assumed to vary around the max density gradient. Naturally, the sensitivity to this condition will increase toward the end of the experiment as the temperature diffuses further outward and the temperature (and density) profile flattens. As mentioned previously the 0D model assumes a linear temperature profile such that to define a boundary edge an isotherm must be selected. Just as in the DNS model the calculated radius changes if the boundary temperature, T b , is varied, however, unlike the DNS model this parameter also defines the slope of the temperature gradient (and strength of mass gain through thermal conduction). It is suspected the sensitivity to the boundary is slightly lower than the DNS since this parameter partially defines the transport characteristics of the kernel (i.e. temperature and mass are more affected since the overall temperature profile is fixed unlike the fully solved temperature profile in the DNS model).
The voltage drop is markedly the most significant input parameter. This value represents the input power to the model and should be the first value to consider for potential errors in the result. An example of the model variation is shown in figure 10(b) with the final model radius presented for the conductivity, boundary condition, and voltage drop. Initial conditions are not displayed as they are insignificant.

Model validation.
To ascertain the validity of both models, the results of the kernel mass and temperature are examined and presented in figures 11 and 12. The results of both modeling methodologies should physically be similar to each other and be physically realistic with available experimental data. Calculated and measured temperatures of similar spark discharges are shown in figure 11 against the discharge current. As the current increases the overall intensity and temperature of the kernel should rise as shown in both sets of data.
Available literature temperature data [33,34,[72][73][74][75][76][77][78] for similar conditions is also plotted. The temperature should be carefully defined as this value can be interpreted in a few ways. The actual temperature of the ignition kernel has a profile radially, which, when measured using an optical emission method would not be capable of distinguishing emissions detected along the axis of the detector's line of sight. For this reason, when comparing the model results to experimental data, the DNS radially resolved temperature can be utilized to calculate a similar measurement (path averaged) through the center of the kernel (Green dashed line). The literature data reported are the rotational gas temperatures for similar spark conditions (1 atm, air, glow). The comparison shows the model prediction is slightly higher than the available data. There are a couple of explanations for this, (1) The experimental measurement path has some thickness: in the ideal case where the maximum temperature is precisely at the center of the kernel will result in a lower temperature if the measurement has a thickness (or a zero thickness line through the center will result in a larger value), (2) There is variation in kernel temperature through the axis of the electrode [34]: if a more detailed plasma structure is considered (positive column, negative glow, dark space), the exact measurement location within the plasma for these experiments may result in different values from this simple model, (3) The uncertainty of these reported temperatures are often several hundred degrees K. Given the many possible sources of variation between available experimental conditions, the difference of ∼500 K observed at the maximum current is suggested to be a reasonably good prediction of the kernel temperature. The temperature calculated is both physically possible and realistic while also observing the expected trend as current (and spark energy) increases. Additionally, close agreement is found between both model temperature   [33,34,[72][73][74][75][76][77][78] in atmospheric air (Case #1).
calculations. To accurately compare the temperature calculated by the Thermodynamic model (given the model is not radially resolved), a 'mass-averaged' temperature for both models is compared where the entire kernel is included in the measurement.
The mass gain in figure 12 also shows good agreement. The mass gain (through thermal conduction) is defined in the thermodynamic model using a linear temperature profile. In the following section, the DNS-calculated temperature profile is displayed which shows a non-trivial profile as conduction extends beyond the discharge region. It is these discharge profiles which might suggest the linear profile is indeed the be approximation without realizing a full 1D model. The rate of change of mass matches toward the end of discharge when diffusional effects (not input energy) are the most significant. The discrepancy in the mass profile of the two models is expected to be a result of the simplified thermal conduction with the best agreement occurring when the DNS model temperature is better approximated with a linear profile (near the middle of the discharge).

DNS results.
Additional results of the DNS model are shown in figure 13 where samples of the radial temperature profile over time and the conditions at the kernel boundary are shown. This provides some key insight into the kernel characteristics. Namely, a more accurate understanding of how the current (rises and falls over time see figure 13) affects the maximum and average kernel temperature (and consequently the rate of kernel propagation). In figure 13(a), the temperature is captured at 6 points where three are under rising current conditions and three are under a falling condition. Both the maximum and averaged temperatures will rise and fall with the current as this indicates a significant change in total power  deposited (for glow discharge). In addition, the conditions at the DNS kernel boundary are presented in figure 13(b) and are useful results for comparison to the thermodynamic model as well as to understand the assumptions used in the definition of the simplified thermal conduction for the thermodynamic model.
The first significant observation for figure 13(a) is the kernel temperature profile towards the end of discharge, has a reduced gradient. This is observed as less defined experimental kernel boundary when the kernel reaches the maximum size and will result in an increase in uncertainty of the edge location as thermal conduction (with zero additional energy) becomes the only mechanism increasing the kernel radius. The density gradient at the edge (maximum) is presented in figure 13(b) to illustrate this point.
The temperature at the edge selected from the DNS model is also presented in figure 13(b) to show that the isotherm chosen to define the kernel boundary and temperature gradient in the thermodynamic model is similar to the DNS model results (i.e. the edge of the thermodynamic model kernel can accurately be assumed to exist at an isotherm). The temperature profile calculated is compared to the results of Maly [9]. Maly calculates a glow discharge that has a maximum temperature of ∼3000 K after a long duration (2 ms). The exact magnitude and shape of the temperature profile is highly dependent on many spark characteristics (duration, energy, gap size…) however the general solution shown in figure 13(a) is in general similar. This profile in Maly [9] is however associated with the minimum ignition criteria defined in his work which is exceeded in the present work in a much shorter spark duration. Here, the temperature and overall radius are observed to vary greatly as the energy within the glow regime changes.

Model application.
Modeling results of the conditions outlined in table 4    kernel condition is pure glow discharge with results shown in figure 14(a). The model result in Case #1 is shown to predict the overall magnitude of the propagation but a discrepancy between both models and the experimental measurements is noted. This deviation in propagation profile from experimental results is expected to be a product of the assumptions made for the loss measurement values and is discussed in detail in the following section. The glow discharge conditions naturally cannot be maintained for all conditions, therefore, discharge containing some form of arc discharge is shown in figures 14(b)-(d).
The arc plasma is to some degree random, so several possible conditions are presented. In figure 14(b), the same atmospheric condition as figure 14(a) is shown with just a small fraction of arc discharge. At the moment of arc transition, the propagation rate of the kernel observes a sudden change (the arc region is shown between two vertical lines). The power measured through the voltage and current waveform of course changes with the arc transition however the expected nonthermal loss of 15 V greatly overpredicts the measured kernel. This measured value is in the same order as the theoretical voltage drops (where the cathode drop is related to the ionization potential of the electrode (∼7.9 V) and the anode drop is comparable to the work-function (∼4.4 V) of the electrode [32]. In the case of figure 14(b), both the originally proposed model loss (∼15 V drop) and the condition with zero arc energy input (high voltage drop) are shown, which indicates the full potential range of model results. The experimental results fall somewhere in between these values which suggest there are additional phenomena that are not accounted for by the present models. One explanation is the simplification of arc discharge geometry. For example, arc discharge can exist off axis that results in kernel growth that is not fully captured by the optical measurement, meaning the total effect that the spark has on the kernel is not observed.
On the other hand, the size of the plasma region can vary greatly. In figure 14(b) the cathode spot of the arc is visible and is almost the same length as the glow discharge (although the overall plasma size is limited to a much smaller diameter). The arc of this length appears to cause significant growth to the once glow-filled kernel. In the case of figure 14(c), the arc plasma is observed to exist solely between the tips of the electrode. This reduction in size compared to figure 14(b) both reduces the measured power but also appears to result in an arc that has almost no effect on the kernel size with complete kernel stagnation. It would be inappropriate with only this knowledge of the plasma size to assume the loss to the sheath formation is higher than the measured values. Rather, the energy heating of the gas must form a more complicated temperature profile with higher center temperatures, where the thermal conduction to the kernel boundary must take additional time and energy as compared to the glow discharge. A realistic modeling of this effect would require at least the resolution for the radial dimensions such as what is done in the DNS model where the size of the energy deposition can be specified. For the case with glow and arc discharges, the evidence is inconclusive to determine the accuracy of the arc loss value since the actual arc plasma dimensions as compared to the glow plasma. Further insight can be found by considering a condition of pure arc discharge since there would be no influencing effects of glow geometry.
Additional results are shown in figure 14(d) for kernel propagation at 3 atm. This condition has considerable arc propagation seen in the plateau of the data and a short burst of glow discharge towards the end seen with a final bump at 0.65 ms. Here the thermodynamic model assumes a larger loss than the 15 V. However, for the DNS model a smaller discharge radius is used suggesting that the dimensionality of arc discharge needs to be accounted for (especially when preceded by glow).
Further results are presented in figure 15 where, most interestingly, a kernel with nearly pure arc discharge is observed. Without a preceding glow discharge, the kernel boundary is known to be directly the result of the arc plasma. Despite the lack of the glow discharge kernel, a larger sheath loss (∼40 V) is used to result in a better fit. The author concludes that this is potentially a result of one of two effects: (1) the arc discharge may not be fully symmetric resulting in unmeasured heated volumes in front or behind the kernel, or (2) there is some other loss mechanism which is not accounted.
In defense of the second explanation, when the cathode spot is beyond the electrode tip, the arc plasma must still be adjacent to the highly electronegative wire. To maintain the plasma region some additional sheath may form outside of the cathode spot to isolate the plasma [66], this sheath would not be present in the arc sheath measurement since discharge between the tips was only used. As observed in Benelov [79], the sheath formation can be potentially very complicated unlike the simplistic approach taken in this work. This potential additional sheath would change with length and is highly dependent on the electrode geometry, complicating future analysis. On the other hand, temperatures are approaching the threshold where radiation becomes significant. In general, an average kernel temperature is used to calculate radiation loss, but a more detailed space-varying result may show an additional loss at the hottest part of the arc kernel near the center.
The other conditions shown in figure 15 include pure glow discharge at 3 atm and arc-glow discharge at 5 atm in the air. A lower sheath loss is selected as the pressure increases. Compared with the 1 atm at pure glow, similar errors are seen with some overprediction at the final radius. This is explained in the following section as the analysis of pure glow discharge can be more definitive since this propagation is more predictable. At 5 atm, there is significant arc discharge where, again, the loss is found to be larger than expected (∼45 V).

Interpretation of results.
The results between the modeling and experiment are noted not to be a perfect match. The difference between the model and experiment is plotted in figure 16 for selected cases. The discrepancy is observed to be on the order of 15% of the experimental radius. Additionally, a structure is noted for most of the experimental results where early time is underpredicted with the final radius overpredicted.
Several potential sources were explored to explain this difference including ideas such as heat conduction to the electrode, larger kernel surface areas from non-spherical shapes, and inaccurate edge detection as the edge density gradient ( figure 13(b)) at the end discharge reduces. However, the sensitivity (figure 10) to these parameters is not significant enough to cause the observed discrepancies. It is proposed that the primary source is most likely a result of the input power, specifically the sheath measurements.
For simplicity in measurement and modeling analysis, a constant sheath loss was assumed early in this work. This assumption has been made in other modeling works [12,26,37,39] and most empirical sheath drop values for glow discharge [32] that are often used are for a normal glow and/or have high uncertainties (especially given the many parameters which affect discharge). This assumption is likely inaccurate (at least for the discharge observed in this work), especially noticeable when observing the kernel propagation with high temporal and spatial resolution. For engine experiments, averaged and approximate values can be acceptable, and when the discharge is less likely glow.
Since the model is primarily sensitive to the input energy, the sheath measurement should be considered with greater scrutiny. A reasonable explanation is found when the current density at the cathode fall was characterized (figure 17) by measuring the size of the luminous region (where possible). For a condition with a constant sheath drop, the current density should remain constant. If the current density changes, the voltage drop should be expected to change as well.
To observe the trend of current density, the discharge area vs time was estimated for glow discharge of a tip-to-tip geometry from the available experimental results. To ensure the observed trend was accurate, the current density was also measured using a flat plate electrode which can more easily observe the sheath diameter. For both conditions, the current density was shown to decrease with time, which suggests that the glow discharge is abnormal [80]. The abnormal discharge is a result of the transient nature of the current rate of change. When the discharge area is plotted against the current, a phase shift of 0.1 ms is observed for the flat plate discharge. This is similar (although an order of magnitude larger) to the transient time (0.1 µs) of Nahemow and Wainfan [81], which observed the glow discharge at a much lower pressure of 10 Torr. Since the discharge never reaches a steady state current and the current density is shown to change, a more accurate voltage drop measurement should vary with time.
Cathode temperature is another factor that likely affects sheath formation. Arkhipenko et al [50] showed that the current density changes for uncooled electrodes which supports the theory that the voltage drop changes with time. The small mass of the electrode along with rapid surface heating would suggest the electrode temperature seen by the plasma rises rapidly.
Further analysis can be done by reversing the model calculation to make the radius the input and the spark power the output. This results in the optimal input 'perfectly match' the kernel propagation. Additionally, if the sheath formation is truly the primary mechanism for non-thermal loss, then this result could represent the real, time-varying sheath drop (assuming the models accurately represent the system, and the radius is absolutely accurate). Regardless, the observed trend of the calculated sheath-drop vs the current density appears to be related. The initial curve for the first 0.1 ms is expected to be from the initial formation of glow discharge where the current density should really asymptote to infinity (area starts at zero, where time and spatial resolution prevent higher accuracy of current density at the early time). This early time might represent the initial stabilization of the glow plasma where electrode surface temperature and ionization of the gas must stabilize and take place.
Once stabilized the calculated sheath drop increases linearly just as the current density decreases. The sheath change in current density is seen to be accentuated by the tip geometry where the plasma must grow in length. For the tip geometry observed the plot of sheath drop vs current density follows a linear trend of approximately, . This trend of course means little in terms of other geometries or discharge conditions, however, characterizing this curve under different ignition ranges may provide additional insights and improved model accuracy in future research.

Conclusions
An investigation of super-atmospheric pressure spark discharge has been conducted and involved a comprehensive study of the morphology of plasma and kernel starting with breakdown through the diffusive glow and arc plasmas. A framework for accurately reproducing the radial propagation of thermal spark plasma utilizing the input spark energy has been developed using two modeling approaches.
Reasonable agreements were found between predictions of the two models as well as between model predictions and experimental results. The models reproduce available literature kernel temperatures and the kernel radius is also shown to be within ∼15% of the experimental radius. The temperature and mass of the 0D thermodynamic model are also comparable to the DNS model.
The benefit of utilizing both modeling approaches is greater understanding and confidence in modeling. The DNS model provides additional detail regarding properties inside of the kernel and additional results to validate against experimental data. The thermodynamic model given the greater simplicity can allow for wider use and can be compared and validated against the DNS model and experimental conditions. The agreement observed here (between models) may not be the case for other conditions (i.e. higher pressure, other geometries, heavy arc discharge), but development of the DNS modeling can help produce develop accurate simplified models in such cases given additional experimental data in future studies.
Simplification of the plasma losses associated the sheath formation to a constant value is thought to be the most significant source of discrepancy between modeling and experimental data. And the sheath formation must be expected to be transient for the entire discharge if the rate of current change is significantly large.
Arc discharge was found to further deviate from the expected energy from potentially further reasons including changes in plasma geometry and additional losses based on the shape of the electrode in comparison to the plasma region.
Important conclusions from the work include: 1. The value of voltage-drop, measured with the experimental technique proposed, is expected to be at most a nominal value for the pulsed discharge (where effective length, electrode temperature, and current density are changing and not accounted for). As a result, more detailed modeling for arc discharge (or more rigorously controlled experimental arc discharge) is required to improve modeling accuracy in this regime. 2. For the small electrode in glow discharge, it is agreed that the cathode fall is the primary non-thermal loss (conduction to electrode/radiation are negligible). 3. Simplification to zero-dimension using a linear temperature profile appears to be acceptable. Additionally, the spherical approximation of the kernel is also acceptable for the glow-arc discharge. 4. The assumption of normal glow discharge is not necessarily acceptable for all discharge types. The current density should always be measured to ensure a constant loss can be applied for modeling results.
Future work to improve results will consider a changing voltage drop, V fall (T elec , t), as a function of electrode temperature and time (or current density). It is proposed that this could be done potentially with experimental data if and alternative ignition circuit design is used with a constant (or nearly constant) current profile. Alternatively, if a time-resolved plasma model (such as [82] or [83]) is implemented, it could provide insight into how sheath formation changes with time and improve the results presented here.
In summary, further understanding of arc discharge and its realization to kernel propagation deserve further studies. Also, incorporating time-resolved glow discharge (or creating a spark that has a constant cathode fall) should be accomplished.

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