Plasma–Catalyst Synergy in Methane Reforming: Simulation-Guided Optimization for Hydrogen Production

1. Introduction

Methane (CH₄) plays a central role in the global energy and chemical landscape, serving as both a key feedstock and an abundant hydrocarbon resource. Traditionally, methane is transformed into hydrogen and other value-added chemicals via processes such as steam methane reforming (SMR) or oxidative coupling. While these thermally driven pathways are industrially established, they suffer from significant drawbacks including high operational temperatures, poor selectivity for specific products, and considerable CO₂ emissions, thus limiting their sustainability and scalability in the context of global decarbonization targets [1,2]. There is an urgent need for more energy-efficient, flexible, and cleaner alternatives to conventional methane upgrading.

In this context, non-thermal plasma (NTP)-assisted reforming has emerged as a compelling pathway for low-temperature, electricity-driven methane activation. In NTP environments, highly energetic electrons dissociate CH₄ molecules into a cascade of reactive radicals (CH₃, CH₂, CH, H), enabling product formation even under ambient bulk gas temperatures [3,4]. This unique property of electron-induced dissociation, decoupled from thermal equilibrium constraints, allows non-oxidative methane conversion into hydrogen (H₂), ethylene (C₂H₄), acetylene (C₂H₂), and other hydrocarbons with rapid startup and shutdown capabilities suitable for distributed, on-demand applications [5,6,7]. However, despite the promise of NTP systems, current plasma-only configurations often exhibit moderate methane conversion rates and sub-optimal product selectivity, primarily due to inefficient radical utilization, side reactions, and non-selective recombination pathways [8,9]

To overcome these limitations, hybrid plasma–catalyst systems have been widely investigated as a means of steering product formation pathways and enhancing reaction rates [10,11,12]. The introduction of catalytic surfaces within plasma zones enables surface-mediated processes such as hydrogen abstraction, C–C coupling, and adsorption–desorption cycles, which are not accessible through gas-phase plasma alone. This synergy, often referred to as plasma catalysis, has been shown to boost hydrogen yields and significantly enhance ethylene selectivity [10,12]. The combination of fast plasma-initiated radical generation with tailored catalytic reaction routes contributes to higher overall efficiency, as reported in several studies involving nanosecond pulsed discharges and bimetallic or Ni-based catalyst surfaces [13,14,15]. Notably, recent works [3,16,17] demonstrated the use of structured catalysts and additive manufacturing to further intensify plasma-assisted C₂H₄ production in low-temperature reactors [11,13].

Understanding the fundamental plasma properties—such as electron temperature (Tₑ), electron density (nₑ), and their evolution under varying discharge conditions—is crucial for reactor design and performance optimization. In non-thermal regimes, electrons typically possess energies in the range of 1–5 eV, sufficient to activate CH₄ and initiate a sequence of chain reactions, while the neutral species remain near room temperature. Previous experimental and modeling investigations have shown that Tₑ and nₑ directly impact the dissociation rate, energy transfer efficiency, and lifetime of radical species in the plasma zone [14,17]. Moreover, the spatial and temporal behavior of these parameters governs the product selectivity and energy efficiency of the reactor, necessitating a detailed kinetic and energetic analysis of the system [6,18]

Given the complexity of plasma-induced chemistries and the presence of multiple reactive pathways, numerical modeling has become an indispensable tool in analyzing and predicting the behavior of plasma-assisted reforming reactors. Zero-dimensional (0D) global models, in particular, offer an efficient approach for studying the time evolution of species concentrations and energy balances without requiring detailed spatial discretization. These models rely on mass and energy conservation equations, coupled with electron impact reaction cross-sections and rate constants, to simulate the dynamic behavior of charged and neutral species over [14,19,20]. Prior global modeling studies have successfully been used to capture the fundamental kinetics of CH₄ conversion and to evaluate the impact of reactor parameters such as power input, gas flow rate, and operating pressure [13,16,21,22].

Despite significant progress, several research gaps persist. The interplay between plasma conditions and catalytic effects remains underexplored, especially in non-oxidative CH₄ conversion under nanosecond pulsed discharges. Additionally, parametric studies that map out optimal operating windows for hydrogen and olefin production are relatively scarce. A lack of experimental validation and sensitivity analysis in existing simulation work also limits the reliability and applicability of current models.

In this study, a detailed zero-dimensional, time-dependent simulation framework is developed using MATLAB’s stiff solvers to investigate methane conversion in a dielectric barrier discharge (DBD) reactor. The model includes a comprehensive set of electron impact and neutral reactions, validated against experimental trends reported in the literature [4,9]. Key variables such as power, pressure, and gas flow rate are systematically varied to ‘

examine their impact on species evolution, product selectivity, plasma properties, and energy efficiency. In addition, catalyst effects are modeled to capture the synergistic behavior of plasma–catalyst interactions. The simulation results provide insight into optimal operating conditions and confirm the reliability of global modeling as a tool for designing next-generation plasma reforming systems. By addressing both gas-phase kinetics and surface-enhanced pathways, this work contributes to the rational design of compact, electrified, and selective methane conversion technologies.

2. Methodology

2.1. System Overview and Simulation Objective

This study presents a zero-dimensional (0D), time-resolved plasma kinetic simulation framework developed in MATLAB to investigate the non-oxidative reforming of methane (CH₄) under non-thermal plasma conditions. The primary objective is to evaluate species evolution, product distribution, and the synergistic role of surface catalysis during plasma-assisted CH₄ conversion.

The model was developed to simulate the temporal dynamics of gas-phase and surface reactions, including the formation of hydrogen (H₂), ethylene (C₂H₄), and key intermediates. Simulations were conducted under base operating conditions representative of dielectric barrier discharge (DBD) systems reported in the literature:

• Discharge power: 70 W
• Pressure: 1 atm
• Flow rate: 50 mL/min

Both plasma-only and plasma–catalyst configurations were considered to quantify the catalytic enhancement mechanisms and provide parametric insights into CH₄ activation behavior, consistent with approaches used by [17,23].

2.2. Reaction Network and Mechanism Description

2.2.1. Gas-Phase Reaction Pathways

The gas-phase chemistry of non-oxidative methane reforming is governed by electron impact dissociation of CH₄, radical recombination, and molecular formation. The reactions were selected and adapted from published plasma studies [23,24,25]. Key pathways include:

• Electron-impact dissociation

  $$e^{-}+{CH}_4\to {CH}_3+H+e^{-}$$

  (1)
• Radical recombination

  $${CH}_3+{CH}_3\to C_2{H}_6$$

  (2)
• Hydrogen formation

  $$H+H+M\to H_2+M$$

  (3)

These steps are central to the formation of key products like H₂ and C₂ hydrocarbons under high-energy electron collisions, a process thoroughly described in studies by [24,25].

2.2.2. Surface Reaction Pathways (Catalyst Enabled)

When surface kinetics are activated, additional mechanisms are introduced, including adsorption of CHₓ species and surface-mediated C–C coupling. These reactions contribute significantly to olefin formation (e.g., C₂H₄) and are known to operate via recombination of surface-adsorbed radicals, as discussed by [9,26].

Surface radical coupling

$${CH}_3\left(ads\right)+{CH}_3\left(ads\right)\to C_2{H}_6\left(ads\right)$$

(4)

Desorption of products

$$C_2{H}_6\left(ads\right)\to C_2{H}_4+H_6$$

(5)

Table 1. Example Reaction Set Used in the Model.

No.	Reaction	Type	Rate Info / Source
R1	$e^{-}+{CH}_4\to {CH}_3+H+e^{-}$	Electron-impact	[23,24]
R2	$H+H+M\to H_2+M$	Radical recombination	[25]
R3	${CH}_3+{CH}_3\to C_2{H}_6$	Gas-phase coupling	[23,27]
R4	${CH}_3\left(ads\right)+{CH}_3\left(ads\right)\to C_2{H}_6\left(ads\right)$	Surface reaction	[9,26]

Table 1. Example Reaction Set Used in the Model.

No.	Reaction	Type	Rate Info / Source
R1	$e^{-}+{CH}_4\to {CH}_3+H+e^{-}$	Electron-impact	[23,24]
R2	$H+H+M\to H_2+M$	Radical recombination	[25]
R3	${CH}_3+{CH}_3\to C_2{H}_6$	Gas-phase coupling	[23,27]
R4	${CH}_3\left(ads\right)+{CH}_3\left(ads\right)\to C_2{H}_6\left(ads\right)$	Surface reaction	[9,26]

2.3. Plasma Kinetic Modeling

The model solves coupled ordinary differential equations (ODEs) describing the time evolution of species concentrations and plasma properties.

2.3.1. Species Mass Balance

Each species $C_i$evolves based on a summation over all reactions j:

$${\displaystyle \frac{d{C}_i}{dt}}={\sum }_j{v}_{ij}{R}_j$$

(6)

where:

i.

$C_i$ is the concentration of species i

ii.

$v_{ij}$is the stoichiometric coefficient

iii.

$R_j$is the rate of reaction j

This framework is suitable for 0D systems with spatial uniformity, as established in prior global plasma models [9,19].

2.3.2. Electron Energy Balance

The time evolution of average electron energy ${\epsilon }_e$ is governed by:

$${\displaystyle \frac{d{\epsilon }_e}{dt}}={\displaystyle \frac{P}{V{n}_e}}-{\sum }_k{R}_k∆{\epsilon }_k$$

(7)

where:

• P is the applied discharge power
• V is the reactor volume
• $n_e$ is the electron density
• $∆{\epsilon }_k$ is energy loss per reaction k

2.3.3. Electron Temperature

Assuming a Maxwellian energy distribution, the electron temperature $T_e$ is linked to energy by:

$${\epsilon }_e={\displaystyle \frac{3}{2}}{K}_B{T}_e$$

(8)

where $K_B$ is the Boltzmann constant. This assumption is commonly used in non-equilibrium plasma models [28,29]

2.3.4. Rate Coefficients for Electron-Impact Reactions

The rate constant k($T_e$) is calculated using:

$$k\left(T_e\right)={\int }_0^{\infty }\sigma \left(\epsilon \right)f\left(\epsilon ,T_e\right)v\left(\epsilon \right)d\epsilon$$

(9)

where:

• $\sigma \left(\epsilon \right)$ is the cross-section
• $f\left(\epsilon ,T_e\right)$is the energy distribution function (assumed Maxwellian)
• $v\left(\epsilon \right)$is electron velocity

This integral accounts for the effect of electron energy on reaction probabilities, and is consistent with formulations in plasma chemistry studies [13,25].

2.4. Simulation Framework and Numerical Method

The numerical solution of the coupled plasma-kinetic system was implemented in MATLAB using a stiff ODE solver (ode15s) to accommodate the broad range of reaction time constants and the inherent stiffness of plasma–catalytic systems. This approach is consistent with established methods for simulating transient kinetics in 0D plasma reactors [19,30].

Simulations were performed under isothermal conditions, assuming a fixed gas temperature and homogeneous mixing within the reactor volume. All species were initialized based on inlet composition—primarily methane at 50 mL/min—and small initial concentrations were assigned to electrons and radicals to initiate discharge. Electron density was initialized at$n_e=1.0\times {10}^{15}{m}^{-3}$, typical for dielectric barrier discharge start-up [31].

No boundary fluxes were applied due to the spatially lumped (0D) nature of the model. Wall losses were neglected, an assumption that simplifies plasma wall interactions without significantly affecting bulk kinetics in DBD-like systems [9,11].

The simulation time domain spanned 0 to 5 seconds, sufficient to capture both transient and quasi-steady-state behaviors. Convergence was verified by ensuring that all species’ net production rates dropped below ${10}^{-6}mol.m^{-3}.s^{-1}$ and that the total electron energy remained conserved within numerical tolerance.

This framework enables resolution of fast electron-driven events and slower surface-mediated pathways in a unified scheme. Such integration is critical in hybrid plasma–catalytic simulations where transient dynamics, radical propagation, and product desorption must be simultaneously addressed [14,26].

2.5. Parametric Study Design

To explore the operational sensitivity of the system, a structured parametric study was conducted by independently varying discharge power, reactor pressure, and methane flow rate while holding other parameters constant. This strategy aimed to quantify how plasma dynamics and surface interactions influence ${CH}_4$ conversion, $H_2$ yield, and product selectivity under different reactor conditions. The chosen ranges reflect experimental practices observed in recent plasma-catalytic reforming studies, particularly those by [9,11], and are consistent with the thermal-plasma hybrid models outlined by [14].

The discharge power was examined over a range from 50 to 90 W in 10 W increments, while pressure was varied between 0.8 and 1.2 atm. Flow rate conditions included 30, 50, and 70 mL/min, with 50 mL/min serving as the baseline, along with a standard power input of 70 W and pressure of 1 atm. These baseline values correspond to dielectric barrier discharge (DBD) operating regimes common in plasma-driven CH₄ reforming systems [3,8].

In addition to capturing steady-state species concentrations, the parametric variations were used to isolate the influence of plasma-only versus plasma–catalyst modes on product formation. Such comparative approaches were emphasized in catalyst synergy investigations by [26] and were instrumental in understanding the energy-resolved behavior reported by [31].

2.6. Performance Metrics

The assessment of simulation performance relied on three primary output indicators: methane conversion, hydrogen yield, and product selectivity. These metrics were evaluated across all time steps to monitor both transient behavior and final equilibrium states, following methodologies from [9,13]

Methane conversion was computed using the classical definition:

$$X_{{CH}_4}=\left({\displaystyle \frac{C_{{CH}_4,in}-C_{{CH}_4,out}}{C_{{CH}_4,in}}}\right)\times 100\%$$

(10)

Hydrogen yield was determined as the mole fraction of hydrogen produced relative to methane fed:

$$Y_{H_2}=\left({\displaystyle \frac{C_{H_2,out}}{C_{{CH}_4,in}}}\right)\times 100\%$$

(11)

Selectivity toward hydrogen and ethylene was used to estimate how efficiently each desired product formed relative to CH₄ consumed:

$$S_{H_2}=\left({\displaystyle \frac{C_{H_2,out}}{C_{{CH}_4,in}-C_{{CH}_4,out}}}\right)\times 100\%\mathrm{and}{S}_{{C_{2}H}_4}=\left({\displaystyle \frac{C_{{C_{2}H}_4,out}}{C_{{CH}_4,in}-C_{{CH}_4,out}}}\right)\times 100\%$$

(12)

These indicators enable not only quantitative tracking of system performance but also reveal the impact of catalyst inclusion on directing product pathways particularly toward H₂ and light olefins. This selectivity-based formulation is in line with frameworks proposed by [20] and validated through kinetic studies by [11,31]

2.7. Validation Strategy

Model reliability was established through qualitative comparison with experimentally reported trends. Although no direct data fitting or numerical replication was performed, model behavior was benchmarked against known directional responses of plasma-assisted methane conversion systems as described in several peer-reviewed works.

One key validation step involved confirming the expected increase in ${CH}_4$ conversion with increasing discharge power. This trend matches experimental findings reported by [9], where higher power enhanced electron impact frequencies and consequently accelerated ${CH}_4$dissociation. Similarly, the model accurately reproduced the inverse correlation between system pressure and both electron temperature and density. This effect, widely documented in dielectric plasma environments, was consistent with observations from [3,11,31].

Lastly, the effect of catalyst activation was validated by comparing cases with and without surface reaction pathways. The catalyst-enhanced model exhibited a notable increase in $H_2$and $C_2{H}_4$ yieldsan outcome that mirrors mechanistic interpretations proposed by [26] and confirmed experimentally by [14] in nanosecond pulsed discharge reactors.

Altogether, these validations support the physical consistency and predictive value of the simulation framework, demonstrating its suitability for investigating process dynamics in plasma-catalytic methane reforming.

3. Results and Discussions

3.1. Base Case Performance Analysis

3.1.1. Steady-State Species Densities

Figure 1 presents the steady-state species concentrations obtained from plasma-assisted non-oxidative CH₄ reforming under base case conditions (70 W, 1 atm, 50 mL/min). Among the observed products, molecular hydrogen (H₂) exhibited the highest concentration at approximately 5.1 × 10²² mol/m³, followed by ethylene (C₂H₄) at 1.5 × 10²² mol/m³. Despite the reaction activity, a substantial fraction of methane (CH₄) remained unconverted (3.8 × 10²² mol/m³), indicating a moderate conversion under these mild discharge conditions. Radical intermediates such as methyl (CH₃) and methylene (CH₂) appeared in significantly lower concentrations, around 6.0 × 10²¹ and 3.0 × 10²¹ mol/m³, respectively, which aligns with their transient nature and rapid consumption through radical recombination and surface adsorption pathways.

This distribution is characteristic of non-thermal plasma systems, where electron impact processes dominate initial CH₄ dissociation, while hydrogen abstraction and radical recombination reactions drive efficient H₂ production. These observations corroborate prior modeling and experimental studies that emphasized the role of energetic electrons in enhancing CH₄ activation and short-lived radical kinetics in DBD systems [12,13,14]. The relatively high ethylene yield reflects the contribution of surface-mediated C–C coupling pathways, which have been extensively studied in plasma–catalyst environments involving Ni- or Co-based active sites [26]. The low accumulation of radical species such as CH₃ and CH₂ further indicates that the reactor achieved a quasi-steady discharge regime, consistent with previously reported nanosecond-pulsed plasma kinetics [11,31]

3.1.2. Product Distribution and Selectivity

The temporal evolution of methane conversion and product yields is illustrated in Figure 2. CH₄ conversion increased rapidly during the first 2 seconds, followed by a plateau indicative of steady-state operation. The H₂ yield exhibited a synchronous trend, ultimately reaching 15.5%, suggesting that hydrogen production is tightly coupled to the rate of CH₄ dissociation. In contrast, C₂H₄ formation followed a slower rise, stabilizing at 4.4%, which reflects its dependence on surface-assisted secondary reactions—particularly C–C coupling of CHₓ intermediates [9,17].

Performance metrics across different CH₄ conversion levels are reported in Table 2. With increasing conversion (from 29.5% to 51.1%), H₂ yield improved from 11.8% to 15.5%, while its selectivity declined from 80.2% to 60.5%. This trend may be attributed to the formation of secondary products or energy losses to non-productive excitation and radical quenching, particularly at higher discharge energies [3]. By contrast, C₂H₄ selectivity remained relatively stable, within the range of 17.2% to 18.1%, indicating that its formation is governed more by catalytic surface phenomena than bulk plasma excitation—an interpretation consistent with reactor-scale selectivity models developed in hybrid plasma–catalyst systems [11,26].

3.2. Catalyst Effect Study

To investigate the influence of catalytic surface reactions on product distribution and conversion efficiency, simulations were conducted under identical plasma discharge conditions (70 W, 1 atm, 50 mL/min), both with and without enabling surface reaction pathways. The comparative results are shown in Figure 3, depicting the percentage change in species concentrations between plasma-only and plasma–catalyst configurations.

The inclusion of surface reaction mechanisms resulted in a substantial increase in product yields. Hydrogen (H₂) concentration increased by approximately 38%, while ethylene (C₂H₄) rose by over 150% when the catalyst was active. These changes are attributed to the activation of surface pathways: H₂ production is enhanced via recombination of surface-adsorbed hydrogen atoms, whereas the significant rise in C₂H₄ reflects efficient surface-mediated C–C coupling reactions—mechanisms that are largely inactive in gas-phase-only environments [9,11].

Simultaneously, the concentration of unconverted CH₄ decreased by approximately 10%, confirming improved methane utilization due to the extended reaction network provided by the catalytic surface. Radical intermediates such as CH₃ and CH₂ were found in lower concentrations in the catalyst-enabled model, consistent with their enhanced surface adsorption and subsequent transformation into more stable products—a trend also noted by [13,14] in hybrid nanosecond-pulsed plasma systems.

Further quantification of performance metrics is provided in Figure 4. The inclusion of surface reactions improved CH₄ conversion from 37% to 47%, while H₂ selectivity increased from 55.0% to 63.8%. Selectivity toward C₂H₄ nearly doubled, from 10.0% to 17.4%, confirming that catalytic surfaces not only extend the reaction pathway but also steer the plasma chemistry toward more desirable hydrocarbon products. This dual effect aligns with previous findings on plasma-catalytic synergy, where selective promotion of target molecules is achieved through tailored surface reaction engineering [8,26,31].

These results reinforce the role of surface-active sites in shaping the overall reaction kinetics and selectivity. By coupling non-thermal plasma activation with catalyst-enabled adsorption and recombination processes, the hybrid system harnesses both electronic excitation and surface thermochemistry—an approach that has been shown to improve single-pass conversion and reduce energy intensity in recent low-carbon reforming designs [10].

3.3. Plasma Properties

3.3.1. Electron Energy and Electron Density

The time-resolved evolution of electron energy (Tₑ) and electron density (nₑ) under baseline operating conditions (70 W, 1 atm, 50 mL/min) was analyzed to characterize the discharge regime and assess plasma stability. As shown in Figure 5, the simulation begins with a high initial electron energy of approximately 3.0 eV, which rapidly declines and stabilizes around 2.90 eV within the first second of operation. This stabilization of Tₑ suggests that the plasma system quickly transitions into a quasi-steady-state regime, where energy input, excitation, and loss mechanisms reach dynamic balance an essential condition for sustaining uniform, energy-efficient plasma chemistry throughout the discharge zone.

The behavior of nₑ, on the other hand, exhibits a continuous increase during the initial seconds, rising from 1.0 × 10¹⁶ m⁻³ at ignition to approximately 2.80 × 10¹⁶ m⁻³ at 2 seconds. This gradual build-up reflects the progressive ionization of CH₄ and subsequent accumulation of free electrons, which is typical in cold plasma environments where electron-impact ionization outpaces recombination during startup. The steady-state profile of nₑ beyond the two-second mark confirms that the reactor has transitioned into a fully developed plasma regime, consistent with simulation results observed in previous dielectric barrier discharge (DBD) studies [9,19,31].

This simultaneous stabilization of Tₑ and saturation of nₑ strongly supports the presence of a non-equilibrium discharge, characterized by high-energy electrons in an otherwise low-temperature bulk gas environment. Such behavior is particularly advantageous in plasma-assisted methane reforming, as it favors selective CH₄ activation through vibrational excitation and dissociative electron impact, while minimizing energy dissipation as heat [14,20] The observed electron dynamics thus align with established performance trends in DBD-driven CH₄ reforming systems, where precise control over plasma conditions plays a central role in optimizing hydrogen yield and hydrocarbon selectivity [8,11]

3.4. Parametric Study

3.4.1. Effect of Discharge Power

To evaluate the influence of discharge power on methane conversion efficiency, simulations were conducted over a range of 50 to 90 W while holding pressure and flow rate constant. As depicted in Figure 6, increasing power resulted in a proportional enhancement in CH₄ conversion, from 29.5% at 50 W to 47.0% at 90 W. This behavior is attributed to the higher electron densities and more frequent ionization events that accompany elevated power levels, facilitating the dissociation of CH₄ and subsequent radical generation. These findings are in line with experimental and modeling results from [9,11], where power modulation was shown to significantly alter both conversion rates and product profiles in non-thermal plasma reactors.

Importantly, while conversion increased with power, the rate of improvement diminished near the upper end of the range, indicating a plateau effect. This nonlinear behavior suggests the onset of saturation in electron-driven dissociation efficiency, consistent with the energy utilization limits observed in dielectric barrier discharge systems [8]. Beyond a critical input level, excess power may lead to inelastic collisions or gas heating, both of which contribute less effectively to CH₄ activation. Thus, the results emphasize the importance of identifying an optimal power range for maximizing yield without incurring energy penalties.

3.4.2. Effect of Pressure

The impact of reactor pressure was studied by varying the operating pressure between 0.8 and 1.2 atm. The results, shown in Figure 7, reveal a clear inverse relationship between pressure and CH₄ conversion, with peak performance observed at 0.8 atm (47.0% conversion). At higher pressures, both CH₄ conversion and electron density declined, with the latter falling to approximately 1.4 × 10¹⁶ m⁻³ at 1.2 atm.

This trend is consistent with prior kinetic analyses in similar non-equilibrium plasma systems, where increasing pressure was shown to decrease the mean free path of electrons, thereby limiting their ability to achieve high-impact energies required for CH₄ dissociation [19,20]. As electron impact events become less effective at higher collision frequencies, the plasma system shifts away from its optimal dissociative regime. These findings reinforce the necessity of tuning operational pressure to balance residence time and energy transfer efficiency in plasma-catalytic reforming.

3.4.3. Effect of Flow Rate

Flow rate is another critical variable influencing both chemical conversion and plasma behavior. Simulations performed at 30, 50, and 70 mL/min demonstrated that increased flow rate resulted in a notable decline in both CH₄ conversion and electron density (Figure 8). Specifically, conversion dropped from 47.0% at 30 mL/min to 31.0% at 70 mL/min, accompanied by a decrease in nₑ from 2.8 × 10¹⁶ to 1.6 × 10¹⁶ m⁻³.

This reduction is primarily due to the shortened residence time at higher flow rates, which limits the exposure of CH₄ molecules to high-energy electrons and active plasma zones. Additionally, faster flow leads to dilution of charged species and less efficient plasma-sustaining collisions. These trends agree with recent process modeling efforts that describe the inverse relationship between gas residence time and conversion in plasma-catalytic environments [14,16]. As such, while higher flow may increase reactor throughput, it compromises overall reforming efficiency, necessitating careful optimization.

3.5. Sensitivity and Validation

The trends obtained from the base case and parametric simulations demonstrate strong qualitative agreement with experimental data reported in recent literature on plasma-assisted methane reforming. In particular, the increase in CH₄ conversion and H₂ yield with higher discharge power matches the findings of [9],who showed that increased energy input enhances electron impact dissociation and radical generation in hybrid plasma-catalytic systems. Similar power-dependent performance trends have also been confirmed in the modeling work of [11] and the experimental studies of [16].

The model also accurately reproduces the inverse relationship between pressure and both CH₄ conversion and electron density. At higher pressures, the mean free path of electrons is reduced, leading to lower electron energy and limited plasma reactivity—an effect systematically observed in studies of packed-bed DBD systems and reported by [14,19,31]. Flow rate sensitivity is likewise consistent: higher gas flow rates reduced CH₄ residence time and limited conversion, reflecting well-known plasma-flow interaction mechanisms in confined reactors [20].

Additionally, the inclusion of surface reactions in the model yielded higher H₂ and C₂H₄ production, as well as reduced intermediate radical concentrations—an outcome that agrees with the synergistic effect of plasma-catalyst interfaces as described by [8,9,26]. In those studies, surface-mediated radical recombination and chain coupling were identified as critical for steering product selectivity.

Overall, the simulation framework presented here captures the essential plasma–catalyst dynamics that govern reaction performance. By reproducing multiple experimentally validated trends—across power, pressure, and catalyst interaction—it offers a reliable basis for performance prediction and reactor optimization.

3.6. Process Optimization

To explore the optimum operational window for hydrogen production, a two-dimensional parameter sweep was carried out by varying discharge power from 50 to 90 W and pressure from 0.8 to 1.2 atm. The resulting H₂ yield map, shown in Figure 9, illustrates the nonlinear interplay between power input and system pressure.

The simulations suggest that H₂ generation is maximized at power levels between 87 and 90 W when pressure is maintained near 1.0 atm. This outcome is driven by enhanced electron activation at high power and favorable energy transfer dynamics near atmospheric pressure. These trends are consistent with previous studies where electron-driven dissociation was shown to increase sharply with power [9,10,14], while pressure effects followed a non-monotonic pattern, as also discussed by [11,19].

At sub-atmospheric pressure, longer electron mean free paths improve energy transfer per electron but reduce gas-phase collision frequency, limiting reactive species generation. Conversely, at elevated pressures, frequent inelastic collisions reduce electron energy and suppress dissociation. These opposing trends define a narrow optimal window—centered around 1.0 atm—where plasma reactivity, energy delivery, and residence time are effectively balanced [20,31].

The model identifies a maximum H₂ yield of approximately 15%, achievable within this optimal regime. Outside of it, either limited excitation (at low power) or plasma inefficiency (at off-optimal pressures) constrains performance. These insights are aligned with current design recommendations for energy-efficient methane upgrading and provide a computational pathway for minimizing experimental burden during scale-up.

4. Conclusion

This study demonstrates the effectiveness of a zero-dimensional plasma kinetic model in simulating non-oxidative methane reforming within a dielectric barrier discharge (DBD) reactor. By systematically exploring the effects of power, pressure, flow rate, and catalytic surface presence, the simulation captures key phenomena governing product formation, including electron-driven radical initiation, surface-mediated selectivity, and the non-equilibrium behavior of plasma species. The findings highlight the dominant role of energetic electrons and transient radicals in initiating CH₄ dissociation, while underscoring the catalytic surface’s critical influence in stabilizing desirable products such as H₂ and C₂H₄. The synergy between plasma activation and catalytic enhancement yields a marked improvement in conversion and product selectivity, confirming that hybrid plasma–catalyst systems are superior to plasma-alone configurations in efficiency and controllability. The simulation also identifies a narrow operational window—characterized by moderate pressure, high but not excessive power, and reduced gas velocity—where energy utilization is optimized, and product yields peak. Furthermore, the model’s alignment with experimentally observed trends supports its reliability as a predictive tool. Its implementation in MATLAB using stiff solvers demonstrates not only computational efficiency but also scalability for broader parametric optimization and reactor design tasks. These insights contribute directly to the rational development of compact, electrified, and tunable methane reforming systems for sustainable hydrogen and light hydrocarbon production. Future work may extend this framework toward spatially resolved (1D or 2D) models and integrate advanced surface reaction mechanisms for catalyst-specific performance prediction, further bridging the gap between simulation and experimental reactor design.