Modeling Nanosecond-Pulsed Spark Discharge and Flame Kernel Evolution

Dilute combustion is considered one of the promising technologies to improve the brake thermal efficiency of an internal combustion engine [1]. The dilute combustion strategy achieved by either high excess-air ratio or high exhaust gas recirculation can raise the charge specific heat and lower the combustion temperature, both leading to a thermal efficiency gain compared to conventional spark-ignition operation. Nevertheless, under dilute operation, the ignition system faces significant challenges in terms of ignition and combustion stability. Stable ignition requires the delivery of the proper amount of electrical energy into the mixture to generate a flame kernel large enough to initiate combustion. Mixture ignitability is closely related to the cycle-to-cycle variability; therefore, a reliable ignition system with optimized operating strategies is crucial to achieving dilute combustion and support engine development effectively.

Computational fluid dynamics (CFD) can be an effective development tool that leverages high-performance computing resources. Within the last few decades, several studies have been devoted to developing spark-ignition models for engine simulation, and many of these models have been adopted in the engine multi-dimensional modeling community [2–10]. Duclos and Colin [3] have proposed an electrical circuit model to calculate the spark energy transfer rate from the secondary circuit to the gas. Spark channel elongation by the flow motion is mostly modeled with Lagrangian particles [3–5,8–10], in that Dahms et al. [4] have taken the effect of turbulent stretch and wrinkle on the arc motion into account. For a localized ignition under crossflow around the spark plug, some studies [4,9] have monitored ignitability along the spark channel based on mixture and flow conditions. The spherical flame kernel is usually initialized either in between the electrodes or at where local conditions permit, and its growth rate and temperature are computed by a set of one-dimensional (1D) governing equations. Tan and Reitz [2] have derived a plasma expansion velocity for the rate of change of kernel radius. In addition, Lucchini et al. [5] have applied a thermal expansion term for the kernel growth after the flame kernel temperature decreases below a certain threshold. Some researchers [8,10] have coupled the Eulerian-type energy deposition model with a Lagrangian-type evolution of the spark channel, instead of using 1D equations for describing the kernel growth. However, these models have not been consistently tested and evaluated at challenging operating conditions. In addition to dilute mixtures, the predictive performance of ignition models may further decrease due to the prevalence of severe in-cylinder conditions for advanced engine concepts, namely less reactive mixtures, high-speed cross flows, and intense turbulence levels around the spark plug. Moreover, complex electrical discharge features, including blow-out and re-strikes [11], pose numerical modeling challenges. To date, a comprehensive and fully predictive ignition model has not yet been developed that can be readily used by the automotive industry. To further improve model accuracy and capability, fundamental understandings of the entire ignition process are required that inform systematic model development.

Spark-ignition consists of a large number of physical processes that include electrical discharge, plasma evolution, Joule-heating of gas, flame kernel initiation, and kernel growth. Timescales of these events cover a wide spectrum from nanoseconds to milliseconds [12]. In a conventional spark system, the breakdown event ionizes the gas between the anode and cathode to open an electrically conductive channel and results in a plasma that reaches extremely high pressure and temperature values within tens of nanoseconds. Next, a shock wave propagates away from the channel within a few microseconds. The breakdown event is followed by the arc phase, where the thin plasma channel expands by the heat conduction and diffusion and initiates the combustion reaction. Finally, the glow discharge phase takes place, where the remainder circuit energy is deposited to the gas at a relatively lower temperature. As the timescale of glow discharge is of the order of milliseconds, continuous electrical energy deposition overlaps heat release by combustion and both contribute to intermediate species production as well. Thus, for a conventional ignition system, it is impossible to decouple the kernel growth from the electrical discharge.

In this study, an advanced approach is proposed to model spark-ignition energy deposition and flame kernel growth. A nanosecond-pulsed discharge (NPD) experiment is designed to lead to spark-ignition in an optically accessible spark calorimeter at Sandia National Laboratories. The experiment is purposely designed to shorten the duration of the discharge, thereby decoupling the flame kernel growth from the energy source. The flame kernel initialization model used in this study has an assumption that a plasma-chemical equilibrium state is reached in the kernel after the NPD without losses of kinetic and thermal energy by shock or radiation. The electrical energy input to the equilibrium computation uses time-resolved discharge voltage and current data along with a developed circuit model. The shape and size of the plasma channel are derived from an equilibrium plasma simulation. The proposed modeling approach is evaluated using a high-fidelity CFD procedure to compare the simulated flame kernel evolution against flame boundaries from the companion Schlieren images. Finally, a parametric study is performed for both input energy and plasma channel volume. In the following sections, a detailed description of the experimental and numerical setup is first provided. This is followed by a description of the CFD modeling methodology and relevant assumptions for the flame kernel initialization. Lastly, model assessment and a parametric study of the model inputs are presented.

All tests are performed in a custom-made (stainless steel 316), optically accessible, spark calorimeter, which is shown in Fig. 1. The experimental apparatus has been discussed in detail previously [13,14], with a summary presented here. The total internal volume of the calorimeter is 29 cm³. Three optical UV grade quartz (Corning 7980) windows are installed in the calorimeter body to visualize the ignition and flame propagation. A 20 mm diameter front viewing window enables direct imaging. Two orthogonally installed side windows with 16 mm diameter clear apertures of 12.7 mm allow for a pulsed light for Schlieren imaging of the streamer volume and ignition kernel to pass through the chamber. A pin-to-pin (P2P) opposed electrode configuration is used; the high-voltage anode is a modified non-resistor spark plug with the J-hook removed, and the anode machined to a conical shape with a rounded tip (∼70 μm radius of curvature) to maximize local electric field strengths while maintaining relatively repeatable discharge characteristics. The anode is centrally positioned at the top of the calorimeter. A 1/8 in. (3.18 mm) diameter stainless steel rod is used as the cathode. The cathode tip has a flat tip with a 200 μm plateau, as observed by a microscopic image. The cathode is installed from the calorimeter base and secured in place by a Swagelok fitting. A constant inter-electrode distance of 1.01 mm is maintained throughout the present study.

The calorimeter is filled to the desired pressure (5 bar absolute) using methane (CH₄) and desiccated air (O2: 20.95%, N2: 79.05% by volume). The mixture is heated to 343 K using resistive heat tape connected to temperature controllers with the temperature monitored using embedded K-type thermocouples. Although the temperature and pressure are lower than a typical in-cylinder condition at the time of ignition, resulting initial densities of 4.88 and 4.92 kg/m³ are comparable to engine-like conditions. For plasma discharges, density is the most relevant thermodynamic parameter. The equivalence ratios of the fuel-air mixture are varied from φ = 1 to φ = 0.55 at the given vessel condition. The lean flammability limit of the methane–air mixture at 5 bar pressure and 343 K temperature is found to be φ_lean = 0.58.

To generate the NPD within the calorimeter, a Transient Plasma Systems Inc. SSPG-101-HF high-voltage pulse generator with a ∼10 ns full width at half maximum (FWHM) pulse width and 5 ns rise time is used. A low-impedance inline attenuator probe is used to monitor pulse voltage and current for each discharge event at the upstream of the spark plug. Schlieren imaging is performed for each discharge using a pulsed green light-emitting diode, a Photron SA-Z high-speed camera, and a vertical knife-edge; a schematic of the setup is shown in Fig. 2. The camera is operated at 12 μs exposure time at 60,000 frames per second. Schlieren images are post-processed to remove noise and increase contrast.

Flame kernel growth simulations are performed using converge [15], a CFD tool that automates the mesh generation process and permits using simple orthogonal grids, locally forced grid embedding, and the adaptive mesh refinement algorithm (AMR). The turbulent flow is modeled with the dynamic structure model [16] with the large eddy simulation (LES) framework. A multi-zone well-stirred reactor model is used for the combustion, in that the multi-zone method [17] is utilized to reduce the computational cost. The GRI-Mech 3.0 chemical mechanism [18], consisting of 53 species and 325 reactions, is employed. The walls enclosing the spark calorimeter are set as non-slip wall boundaries with a constant temperature of 343 K. The Werner and Wengle wall model [19] is applied. A conjugate heat transfer (CHT) simulation is not considered due to the lack of significant contact between the burned gas and the walls.

Three equivalence ratios, one stoichiometric condition (φ = 1), and two lean conditions (φ = 0.7 and φ = 0.6) are numerically investigated. The three-dimensional computational domain covers the entire inner volume of the spark calorimeter, and the grid configuration near the electrode gap region is shown in Fig. 3. The base grid size is 1.6 mm, and the fixed grid embedding is introduced for two spherical volumes covering the electrode gap region. A 1D laminar premixed flame simulation reveals that the laminar flame thicknesses are ranging between 146 μm and 422 μm for three equivalence ratios. The AMR approach is defined to have at least five grid points across the flame thickness. So, a sub-grid scale-based AMR is activated for the mass fraction of HCO, temperature, and velocity. The minimum grid size is 25 μm for all cases by Y_HCO and temperature sub-grid scales-based AMR.

For the NPD-induced plasma, the characterization of the channel shape and size is carried out using vizspark, a multi-dimensional thermal plasma solver that computes magneto-hydrodynamic equations assuming single-temperature and equilibrium chemistry [20,21]. Figure 4 shows the axis-symmetry of the two-dimensional computational domain used in the vizspark simulation. To accurately capture the arc root on the electrode tips, the solid regions of electrodes are also included and solved via an immersed boundary method. The grid consists of the mixed triangle and quadrilateral cells, where the minimum grid size is 5 μm in the gap region, and the grid size is gradually increased up to 40 μm. To solve the current continuity equation, a potential boundary condition is imposed on the top of the anode, while a current density boundary condition is assigned to the bottom of the cathode.

Flame kernel initialization is an outcome of the electrical discharge and plasma dynamics. The majority of ignition models used for engine simulations simply neglect to solve the breakdown phase and assume an initial flame kernel with arbitrary shape and size. Also, ionized species are not taken into account. Comprehensive modeling that couples plasma physics with a flow solver has been conducted [22]. However, such an approach is still prohibitively expensive.

Plasma evolution by electrical discharge for the NPD used in this study is highly transient with multiple temporal and spatial scales, and complex physics. Relevant dynamics include streamers propagation, ionization and recombination, and joule-heating. Due to the great disparity in timescales, the result of streamer dynamics can be fed as an initial condition to the CFD solver. The output is expected to be extremely high pressure and temperature kernel with ionized species in equilibrium, which initiates the flame kernel and its early growth.

In this study, a plasma-chemical equilibrium approach is proposed for modeling the flame kernel initialization. Several assumptions are made. First, the thermodynamic properties and species composition are assumed to be in equilibrium because of the high-temperature plasma. From the experimental data, the high-temperature plasma is evident because of an enormously high current footprint. Second, ionized species from the plasma equilibrium are not introduced to the CFD domain. Adding ionized species to the CFD solver would require updated chemical kinetics that also capture multiple ionization and recombination pathways, which will significantly increase the cost of the simulation. Accordingly, ion-recombination is not taken into account. Third, the kinetic energy loss by shock is neglected in the equilibrium calculation. Heat losses due to radiative heat transfer are neglected as well.

The proposed approach calculates two equilibrium states to obtain initial kernel properties, in terms of temperature, pressure, and species composition. Only pure air is considered since the thermodynamic database for ionized hydrocarbon species is limited. Because the mole fraction of methane at the stoichiometric ratio is less than 10%, the pure air assumption is considered acceptable here. For plasma equilibrium, the air is equilibrated isochorically by maximizing the entropy with total energy deposited from the electrical discharge. A total of 14 species are considered: e−, N, N+, N−, NO, NO+, N2, N2+, N2−, O, O+, O−, O2, and O2+. Next, the methane–air mixture is equilibrated isobarically by minimizing the Gibbs free energy at the resulting temperature and pressure from the plasma equilibrium. Fourteen key species are included in the equilibrium: N, N2, O, O2, NO, C, CO, CO2, H, OH, HO2, H2O, H2, and CH4. The thermodynamic database of all ionized and ground species is obtained from the NASA website [23], and an open-source chemical kinetic software, cantera,¹ is used to calculate equilibrium. Thus, the initial kernel temperature and pressure are obtained from the plasma equilibrium, and its species compositions are taken from the second chemical equilibrium step.

To compute equilibrium, two key inputs are needed: (1) energy deposited from the electrical discharge and (2) plasma channel volume. The experimental data is used to obtain the former input, while the plasma solver is employed to compute the latter.

Voltage and current are coupled to the discharge event that takes place between the electrodes. It is difficult to measure these quantities without perturbing the electrode configuration. These quantities are typically measured somewhere upstream of the electrodes, leaving the voltage and current that pass through the discharge to be inferred. In an ideal simulation, the circuit and discharge phenomena would be coupled and manifested in the boundary conditions of the electric field and charge transport equations. The complexity of the routines that take transmission line effects into account, and their accompanying stiffness make such a detailed coupled approach unattractive. In this study, a simpler semi-empirical method to derive boundary voltage and current is proposed.

Retrieving the geometric shape and size of the plasma channel from the Schlieren images is impossible due to image saturation by the highly luminous discharge. In this study, vizspark is used to simulate the evolution of the plasma channel. Figure 6 shows the conduction current density, pressure, and temperature fields in the gap between the two electrodes. As shown earlier, the peak current density occurs at ∼90 ns, and the NPD ends at ∼110 ns. In such a short discharge duration, the electrical energy of 4 mJ is deposited into the gas. Accordingly, a very thin conductive channel is established along the electrode axis and allows the current flowing across the gap. The electrical power dissipation results in significant gas heating, and the gas temperature and pressure are extremely elevated. The mean temperature (∼25,000 K) is much greater than the adiabatic flame temperature, so it is fair to assume the equilibrium state for the gas inside the kernel. The high pressure induces a strong shock to propagate outward, and a clearly detached shock front is observed. The plasma channel is observed to be a cylindrical shape, despite some corrugations.

To obtain the diameter of the cylindrical plasma channel quantitatively, the temporal evolution of radial distributions of temperature and pressure is analyzed. Gas temperature and pressure are averaged in a number of bins in the radial direction, which is shown in Fig. 7. The sharp decrement of temperature can be seen as the interface between the high-temperature plasma channel and the fresh mixture. At the end of discharge, the diameter of the plasma channel is observed to be approximately 100 μm considering the axisymmetric domain.

By computing the plasma-chemical equilibrium with the inputs derived so far, Table 1 lists the thermodynamic conditions and species composition in the kernel. As the NPD is not varied with mixture equivalence ratio, the temperature and pressure that result from the plasma equilibrium are kept constant for all cases. Due to the high temperature, most species are dissociated to the atomic level, and only molecular N2 exists in small fractions.

A standard energy deposition approach that sources power into a prescribed shape is also tested in the CFD simulations. The power profile from Fig. 5(b) is used, and a cylindrical source shape with a diameter of 100 μm is specified. Simulations that use the standard energy deposition are not able to handle the high-power deposition within the thin cylindrical volume and result in temperature beyond the code maximum temperature limit (∼95,000 K), which forcedly terminates the simulation during the power sourcing. This is because all the energy is deposited into the ground species without consideration of the ionization. Hence, the ionization process should be considered to predict an accurate gas temperature for the flame kernel initialization.

The modeling approach for flame kernel initialization is evaluated by comparing the flame kernel evolution against the experimental data. Figure 8 shows the comparison of kernel growth between simulation results and acquired Schlieren images at φ = 1 condition. The first few frames of Schlieren images are saturated by bright light emitted from the arc and electrodes, so it is not possible to discern any meaningful comparison from those initial frames. The pseudo-Schlieren image is built by taking the gradient of the density field and plotted on the cut plane through the electrode axis. Note that it is not post-processed with a line integral to obtaining a line-of-sight picture. As the Schlieren images indicate the gradient of the density field, its volume is much larger than the burned region enclosed by the reaction layer. Based on the 1D premixed flame simulation result, the inner reaction layer temperature of 1800 K is used to show combustion activities inside the kernel.

The optical measurement reveals that the flame kernel vigorously expands in the horizontal direction. It can be inferred that the flame kernel forms a toroidal shape and surrounds the electrode without a large contact area. To investigate the effect of heat loss on the walls on the kernel evolution, a simulation with adiabatic wall boundary conditions is performed. The difference of kernel evolution between an adiabatic wall and an isothermal wall turns out to be insignificant, which implies that a CHT simulation is not needed for this type of problem. This is mostly due to the extremely short duration of discharge and the initial kernel being blown away from the electrodes and consequently not losing significant heat to them.

Simulation results are well-aligned with the experimental observations, in particular, the kernel size by initial expansion, and subsequent growth calculated by simulations match observations from the Schlieren images. As mentioned earlier, the reaction zone is embedded in the heated volume. It is interesting to note that the reaction zone is shredded due to vigorous expansion, as can be seen at 66.7 μs. However, the reactivity of the stoichiometric mixture is enough to compete with the blow-off mechanism and finally sustains the flame to propagate outward.

As the simulated initial kernel expansion matches Schlieren images, simulation results are further leveraged to analyze the dynamics that occur in the very initial period of the spark-ignition event. Figure 9 shows the temporal evolution of velocity and temperature fields. Two velocity components are used to interpret the shock and expansion dynamics. One microsecond after discharge, a shock departs away from the high-temperature kernel, and the kernel size vigorously expands from the initial diameter of 100 μm to ∼700 μm. Then, a pressure equilibrium process starts as the pressure behind the shock front becomes lower than its initial value. Balancing pressure in a free volume would assist the flame kernel in expanding spherically. However, the electrode geometry triggers the non-trivial expansion dynamic of the flame kernel.

There is no obstacle against the horizontal expansion, so the kernel consistently expands from the electrode axis after the shock departure. In contrast, vertical expansion is much more complex than horizontal one. The plateau on the cathode tip fully reflects the initial shock wave compared to the rounded tip of the anode, so a steeper pressure gradient is formed that is normal to the plateau surface. On the other hand, the bottom of the anode insulator surface plays a similar role as the cathode tip, which reflects the shock and generates a downward momentum. Two parcels of gas with oppositely directed longitudinal momenta produce a clockwise rotating vortex and assist the horizontal expansion. As a result, the flame kernel shrinks vertically but expands horizontally, which can be seen by comparing the temperature field between 1 μs and 10 μs. The upward momentum by the shock reflection at the cathode tip lasts tens of microseconds and sustains the vortex motion. Finally, the combination of all expansion dynamics results in a toroidal flame kernel shape, and then the flame propagation continues. These observations align well with the previous experiment reported in Refs. [24,25].

The flame kernel initialization model is also evaluated for a lean condition (φ = 0.7) for which experiments show successful ignition. Figure 10 shows the comparison of flame kernel evolution between simulation and experiments. Note that the NPD characteristics do not change according to the mixture equivalence ratio. The adiabatic flame temperature decreases from 2280 K (φ = 1) to 1869 K (φ = 0.7), while the temperature of the inner reaction layer is 1638 K. Accordingly, the smaller kernel volume shown in the first Schlieren image (t = 66.7 μs) is due to the lowered mixture reactivity. The simulation result indicated by pseudo-Schlieren and iso-temperature contour shows that the horizontal kernel expansion at 66.7 μs is comparable to that of the φ = 1 condition. However, as the laminar flame speed decreases from 24.1 cm/s (φ = 1) to 11.1 cm/s (φ = 0.7), the kernel growth becomes slower. Despite the lower combustion rate, the kernel grows properly and transitions into a propagating flame, which agrees with the experimental observation.

The model is further evaluated at an ultra-lean condition (φ = 0.6) which is close to the experimental lean limit, and the comparison against the experimental data is shown in Fig. 11. The laminar flame speed of φ = 0.6 mixture at the given condition is 5.9 cm/s, in which the temperature of the inner reaction layer decreases down to 1537 K. Similar to the previous two conditions, the simulation is able to capture the initial kernel expansion, but it fails to predict the later stage flame propagation. At 133.7 μs, the kernel volume represented by the pseudo-Schlieren has a comparable width, but a thinner height than that observed in the experiment. Consequently, the embedded reaction zone is packed into a small volume and broken into several pieces. The reaction is eventually quenched by blow-off, while the heated kernel volume starts to shrink.

To further analyze the combustion dynamics, Damköhler number (Da) is computed for the early stage of kernel growth, as shown in Fig. 12. Da is defined as the ratio of flow residence timescale to chemical timescale, in that the former one is the reciprocal of strain rate (1/(du/dx)) and the latter one is the rate of change of methane molar concentration ([CH4]/(d[CH4]/dt)). After the initial kernel expansion up to 66.7 μs, the Da ahead of the reacting front for φ = 0.6 condition turns out to be much smaller than that results in φ = 1 condition. At 133.7 μs, the Da of the stoichiometric case is greater than 10, while that of the lean case becomes comparable to or smaller than unity. Hence, the former case can have vigorous flame kernel growth under the reactive mixture, but the reacting front of the latter case starts to shrink because of perturbations of a chemical reaction by the flow. The kernel expansion momentum, as well as the vortex motion, reduces the residence time for the flame to be self-sustained. As a result, φ = 0.6 represents a condition beyond the lean limit for CFD simulations. The CFD lean limit lies between φ = 0.7 and φ = 0.6, and it is very close to the lean limit observed in the experiments. The remaining discrepancy between CFD and experiments is likely due to the assumptions made for this study. The heat production by ion-recombination processes, the sensitivity of near limit propagation to the chemical kinetic mechanism, and the effect of heat loss by radiation are all worth considering in a future study. Despite a fine grid resolution across the flame thickness used in this study, a higher-order direct numerical simulation may be required to unveil the flame kernel evolution in such challenging mixture conditions.

To gain more insight into the impact of boundary conditions on the CFD results, the parametric study on the two inputs (electrical energy and plasma channel size) for the plasma-chemical equilibrium calculation is conducted as follows. The aim of this study is to investigate how different flame kernel evolution results from the input variations (i.e., uncertainties in measurement or plasma simulation) in the flame kernel initialization model. For the φ = 0.6 conditions, the very first observable Schlieren image (t = 66.7 μs) is set as the reference, and the inputs are varied for energies of 1, 2, 4 mJ, and diameters of 100, 200, 400 μm, respectively.

Figure 13 presents the comparison between the acquired Schlieren image and the simulation results. Except for the case of the energy of 1 mJ and the diameter of 400 μm, a bigger diameter for kernel initialization results in a larger kernel expansion. The effect of energy input is more obvious than that of diameter; kernel expansion is largely under-predicted with half or quarter of the measured energy value. Thus, it can be concluded that energy played a key role in the proposed model. Also, it implies that accurate measurement of voltage and current signals is essential to the simulation. It is worth noting that, for the φ = 0.6 condition, neither the baseline case (E: 4 mJ, d: 100 μm) nor the other cases from the parametric study are able to sustain the flame after the expansion and lead to flame propagation. Further studies are required on the flame kernel evolution under the ultra-lean mixture condition close to the lean flammability limit.

In this study, a high-fidelity CFD procedure to simulate the flame kernel initiation and growth from a nanosecond-pulsed spark is developed, and a plasma-chemical equilibrium approach for initializing the plasma energy deposition is proposed and evaluated.

The model assumes the equilibrium state in the initial flame kernel after the NPD without losses of kinetic and thermal energy by shock and radiation. By computing the plasma equilibrium, the model takes ionization into account to properly calculate the initial temperature and pressure inside the flame kernel. In order to provide accurate energy input, a semi-empirical method is developed to post-process the measured voltage/current traces and extract the energy, power, and current being transmitted to the electrodes. Equilibrium plasma simulation is performed using vizspark to obtain the shape and size of the kernel. Combining all these inputs to the plasma-chemical equilibrium computation, the initial flame kernel is characterized by very high temperature and pressure, mostly filled with dissociated atoms of O, N, C, and H.

The proposed model is evaluated against Schlieren data from experiments performed in an optically accessible ignition test vessel. The computed flame kernel expansion matches well with the experimental measurements. The toroidal shape of the flame kernel observed from the experiment is properly reproduced by the CFD simulation. Simulations also provide insight into the mechanisms leading to the kernel evolution after the spark. Ultimately, the calculated flame kernel evolution is in good agreement with experiments for the stoichiometric (φ = 1) and one of lean conditions (φ = 0.7). Simulations fail to describe ignition success for the ultra-lean (φ = 0.6) case close to the lean flammability limit, which does not agree with experiments. Future investigations will focus on the role of chemical kinetics, ion-recombination, radiation, and grid size to bridge the gap between modeling and experiments in the flame kernel growth for the lean case.

The submitted paper has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (“Argonne”). Argonne, a US Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. This research is funded by DOE’s Vehicle Technologies Program, Office of Energy Efficiency and Renewable Energy. The authors would like to express their gratitude to Gurpreet Singh and Michael Weismiller, program managers at DOE, for their support. In addition, the authors gratefully acknowledge the computing resources provided on Bebop, a high-performance computing cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory. Finally, the authors would like to thank Anand Karpatne and Doug Breden at Esgee Technologies Inc. for their technical supports and valuable discussions on the use of the vizspark software.

There are no conflicts of interest.