Intrinsic microkinetic effects of spray-drying and SiC co-support on Mn–Na₂WO₄/SiO₂ catalysts used in oxidative coupling of methane

Abstract

This paper presents a microkinetic model to evaluate the effects of a silicon carbide (SiC) co-support and the shaping method on Mn–Na₂WO₄/SiO₂ catalysts used for the oxidative coupling of methane. The model considers mass transfer, catalytic, and gas-phase kinetics, and it is trained with experimental values (product composition) of three Mn–Na₂WO₄ catalysts for calculating the kinetic parameters using catalytic descriptors while maintaining thermodynamic consistency. The catalysts were an SiO₂-supported catalyst prepared through impregnation and two SiO₂–SiC-supported catalysts (with βSiC and α + βSiC) prepared via spray-drying. Our analysis shows how the type of SiC and preparation method affect the textural properties and result in distinct CH₃˙ radical oxidation, HO₂˙ quenching, C₂H₄ oxidation, and CO_X transformation pathways, eventually leading to CH₄ conversion and C₂ selectivity. Our approach facilitates the assessment of the effects of the promoter and support on individual and global reaction networks.

---

1. Introduction

Oxidative coupling of methane (OCM) is a promising route for natural gas valorization as C₂ chemicals are produced by the reaction of CH₄ and O₂ in a single pass in the presence of a catalyst. It has some inherent advantages over other ethylene production routes, including (i) lower greenhouse gas emissions than steam cracking of naphtha;¹ (ii) lower energy requirement and capital cost than routes involving feedstock gasification; and (iii) use of noble-metal-free catalysts, unlike ethane dehydrogenation.^2,3 However, unlike the methanol-to-olefins technology, which provides high ethylene yields and allows the adjustment of the ethylene/propylene ratio, OCM gives low ethylene and propylene yields.⁴ Thus, for OCM to be economically viable for industrial use, it is necessary to develop active and highly selective catalysts to improve the ethylene yield, preferably above 30%.⁵

Among OCM catalyst families, trimetallic combinations of Mn–Na–W/SiO₂ show remarkable CH₄ conversion (15–40%) and C₂ selectivity (55–80%). These three metals show a noteworthy synergy, which is evident from the inferior catalytic performance that results when one of them is absent.^6,7 The thermal stability of these catalysts is widely recognized. However, recent studies^8–13 have reported a concerning trend of performance drop over time on stream. Various research lines have emerged to enhance the activity/stability of Mn–Na₂WO₄/SiO₂ catalysts, including (i) the use of mixtures of Mn–Na₂WO₄/SiO₂ catalysts with alkali chlorides,¹⁴ (ii) the addition of new dopants;¹⁵ and (iii) support material.¹⁶ Our previous studies^17,18 proposed the introduction of silicon carbide (SiC) as a co-support material for Mn–Na₂WO₄/SiO₂ catalysts through spray drying. The uniform distribution of SiC and its high thermal resistance prolong the catalyst lifespan. We also found that the crystal structure of SiC plays a pivotal role in the catalyst's performance. Compared with catalysts with βSiC, α + βSiC offers enhanced metal exposure stemming from its enhanced resistance to being oxidized to SiO₂.

In this work, we investigated catalytic differences observed when SiC was introduced on Mn–Na₂WO₄/SiO₂ catalysts and the shaping method, from a microkinetic viewpoint. The developed microkinetic models of OCM represent a unique tool to analyze the interplay between radical chemistry in the gas phase and catalytic surface reactions.¹⁹ This study builds upon previous experimental findings by examining the role of SiC in the support and shaping method, now from a microkinetic standpoint that accounts for irreducible mass transport limitations, gas-phase reactions, and surface chemistry. Additionally, this work accentuates the interrelation between catalyst surface properties, crystal phases, and their consequential effects on reaction pathways and product distributions. To achieve this, we modeled two Mn–Na₂WO₄ catalysts prepared by spray-drying with different SiC phases in their co-support, using a Mn–Na₂WO₄/SiO₂ catalyst (where Mn–Na-W is impregnated onto SiO₂) as a reference.

2. Methodology

2.1 Catalyst preparation

A SiO₂-supported Mn–Na₂WO₄ catalyst, hereafter referred to as IMP SiO₂, was prepared via wetness impregnation and used as the benchmark. Precursor salts were incorporated into the support at 80 °C to obtain a target composition of 2 wt% Mn, 5 wt% Na, and 3.1 wt% W. Following impregnation, water was removed through drying at 100 °C for 6 h. With the same nominal metal loading, two different Mn–Na₂WO₄ catalysts with SiO₂/SiC supports (70/30 wt/wt) were prepared by spray-drying. Different types of SiC (i.e., SiC with different crystal phases) were introduced: one with a nanosized α + βSiC and the other with commercial porous βSiC. These two types are hereafter referred to as SD SiO₂–α + βSiC and SD SiO₂–βSiC, respectively. All catalysts were ultimately calcined in air at 800 °C for 6 h with a heating rate of 10 °C min⁻¹. Additional information regarding materials used and catalyst synthesis procedures can be found elsewhere.²⁰

2.2 Catalyst characterization

Textural properties were studied using liquid Ar (−186 °C) adsorption–desorption (Micromeritics ASAP 2040). Before the measurements, samples were degassed for 10 h at 250 °C. The Brunauer–Emmett–Teller (BET) method was employed to calculate the specific surface areas, and the Barrett–Joyner–Halenda model was used to measure the cumulative pore volume. Temperature-programmed experiments were conducted using an Altamira AMI-200ip instrument equipped with a mass spectrometer. A typical experiment began with sample pretreatment under 50 NmL min⁻¹ Ar flow at 200 °C for 2 h to remove impurities, which was followed by cooling to the initial temperature. Subsequently, the sample was exposed to O₂ with a flow rate of 50 NmL min⁻¹ for 4 h (10 vol% O₂ in N₂) at room temperature. The catalyst bed temperature ranged from room temperature to 850 °C, and then held for 30 min using a carrier gas. Following the previous procedure, when the carrier gas was pure Ar, the experiment was considered O₂ temperature-programmed desorption (O₂-TPD). In contrast, when the carrier contained different CH₄ partial pressures (0.05 bar balanced with Ar or pure CH₄), the experiment was considered a temperature-programmed surface reaction (TPSR).

2.3 Steady state kinetic data acquisition

To examine the kinetic behavior of the three OCM catalysts, we used a set of 16 tubular fixed-bed quartz microreactors arranged in parallel. Information regarding the reactor specifications, analytical methods, and definitions of performance metrics can be found in this reference,²⁰ where a data curation strategy is proposed for high-throughput kinetic data collected for the three catalysts. The range of operating conditions was as follows: temperature, 740–800 °C; pressure, 101 kPa; feed CH₄/O₂ ratio, 2.2−3.8; feed dilution, 10.8–71.2 kPa; space time, 0.68–6.15 g_c h mol_C⁻¹; CH₄ conversion, 1.1–14.1%; and O₂ conversion, 1.4−36.1%.

2.4 Reactor model

Under OCM conditions, there are constraints on the transport of radicals, which are highly reactive, that cannot be reduced. To address these irreducible mass transfer limitations, Couwenberg et al.²¹ proposed an isothermal model that considers two radial phases: the solid or intraparticle phase (characterized by the dimensionless particle coordinate ξ) and the fluid phase or interstitial phase (characterized by the dimensionless radial coordinate r). The continuity equations of each gas-phase species can be formulated for both the interstitial phase and catalyst particle or intraparticle phase (eqn (1) and (2)). In the interstitial phase, the model considers that molecular diffusion and gas-phase reactions occur. In the intraparticle phase, diffusion through pores occurs along with gas-phase and catalytic reactions.where F_v is the total gas volumetric flow rate and calculated from the total mass flow rate using ideal gas law as equation of state, ε_b is the average bed packing porosity, A_r is the reactor cross-section, C_g,i is the concentration of gas-phase species i in the interstitial phase, z is reactor bed length, D_m,i is the molecular diffusivity of gas-phase species i in the interstitial phase, r_v is the radius of the interstitial phase or average half distance between catalyst pellets, r is the dimensionless radial interstitial coordinate, R_g,i is the homogeneous net production rate of gas-phase species i, D_e,i is the effective diffusivity of gas-phase species i in the intraparticle phase, r_p is the radius of the intraparticle phase or the average pellet radius, ξ is the dimensionless radial intraparticle coordinate, C_s,i is the concentration of the gas-phase species i in the intraparticle phase, ρ_s is the catalyst density, S_s is the catalyst specific surface area, R_s,i is the heterogeneous net production rate of gas-phase species i and ε_s is the catalyst porosity. Note that R_g,i and R_s,i are defined per unit volume of gas and N_s,g represents the number of gas species in the reaction mechanism.

Both phases are coupled by the conservation of mass at the interphase, through boundary conditions that enforce equality of mass fluxes (eqn (3)) and concentrations (eqn (4)) at r = ξ = 1. The remaining boundary conditions pertain to the symmetry of the phases (eqn (5) and (6)), and the initial condition defines the feed concentration in the interstitial phase (eqn (7)). At z = 0, the concentration in the intraparticle phase is determined by eqn (2). Since no inert solid is introduced to dilute the particle, the areas of mass flux in the phases (i.e., interphase area and solid external area) are identical. The reactor is hence modeled by integrating the axial coordinate of a cylindrical gas volume with a contiguous catalyst sphere, with mass transport occurring across the phases, as schematized in Fig. 1. Assuming fully spherical particles is consistent with the previously reported SEM images of the spray-dried catalysts used in this work.^17,18 The subscript 0 indicates feed conditions.

∀z∧ξ = 1 C_g,i = C_s,i i = 1, …, N_s,g (4)

z = 0 ∧ 0 < r < 1 C_g,i = C_0,i i = 1, …, N_s,g (7)

As the catalytic net production rate of each gas species depends on the rate of individual elementary catalytic reactions and, therefore, on the coverages of surface intermediates, the conservation of surface species (eqn (8)) should be coupled with the continuity equations via the pseudo-steady-state approximation (PSSA), together with the site balance (eqn (9)).where σ is the active site density, θ_i is the fractional coverage of surface intermediate i, t is time, θ_* is the fractional coverage of the vacant sites or vacancies and N_s,s represents the number of surface intermediates in the reaction mechanism, excluding the vacancies. Species molecular diffusivities in the mixture (D_m,i) are determined from the Chapman–Enskog theory by using a mixture-averaged approach in which the molar diffusion velocity is expressed with respect to the molar average velocity and the velocities of all species j, with j ≠ i, are approximated to be equal.²²where D_i,j is the binary molecular diffusivity of gas-phase species i in gas-phase species j and X_i is the molar fraction of gas-phase species i in the mixture. The effective diffusivity in the catalyst particle (D_e,i) is determined by modifying the Dupuit law to incorporate pore constriction effects. This modification disregards the Bosanquet formula since an order-of-magnitude estimation of the average catalyst pore radius showed the effect of the Knudsen diffusion to be negligible.where τ_s is the term combining constriction and tortuosity of the catalyst.

Fig. 1 Schematic of Couwenberg et al.'s model,²¹ which was used in the present work for modeling OCM kinetics responsible for irreducible mass transfer limitations.

2.5 Homogeneous kinetic model

A reliable gas-phase mechanism is important as the coupling of CH₃˙ radicals is the primary pathway for the formation of C₂ products, which occur exclusively in the gas phase. Wang et al.²³ highlighted the importance of using detailed combustion models over ad-hoc gas-phase models developed for OCM. However, increasing the complexity of the model increases the number of equations to be solved for both phases, and therefore, this work used the most refined ad-hoc homogeneous model, which was proposed by Chen et al.²⁴ (with 39 reactions and 23 species) and validated under catalytic OCM conditions. Details of this model are available in ESI† (section S1). The net production rate of gas-phase species i (R_g,i) is computed from the forward and backward rates of each homogeneous elementary reaction step (r_j) and the stoichiometric number of i in each j (ν_j,i), where f and b denote forward and backward:where N_r,g denotes the total number of homogeneous elementary steps in the mechanism. Every gas-phase reaction is modeled as an elementary reaction that follows the law of mass action, such that, the reaction rate of step j (r_j) in the homogeneous mechanism is calculated from the product of its rate constant (k_j) and the concentration of the reactants (C_i) to the power of their stoichiometric number, applicable for both forward and backward steps. Note that ν_j,i is the stoichiometric number of reacting species i in gas-phase elementary step j (i.e., ν_j,i wherein all non-negative values are set to 0).The subscript j denotes the step number in the homogeneous reaction mechanism (Table S1†). In three-body elementary reactions, C_M denotes the concentration of an unspecified collision partner that carries away excess energy to stabilize the product molecule (forward direction) or that supplies energy to break the product molecule bond (reverse direction).Forward reaction rate constants follow the two-parameter Arrhenius equation, and their values are presented in Table S1.† The backward rate constant is computed from the forward and equilibrium constants (section 2.7).where A^f_j and E^f_a,j denote prefactor and activation energy of forward homogeneous reaction step j, respectively, R is the universal gas constant and T is temperature.

2.6 Heterogeneous kinetic model

The first step of the heterogeneous reaction network involves H˙ abstraction from CH₄ to form CH₃˙ Whether CH₄ is adsorbed on the catalyst surface or whether H˙ abstraction occurs through an Eley–Rideal step as well as the nature of O* have been discussed extensively, even via large-scale density functional theory calculations.²⁵ In this work, an Eley–Rideal CH₄ initiation with dissociative O₂ chemisorption was considered, following the model of Alexiadis et al.,²⁶ who applied this model to diverse catalysts, including Mn–Na₂WO₄/SiO₂, and considered a second-order chemisorption kinetics observed in pulsing experiments.²⁷ Nonetheless, for other catalyst families, CH₄ surface dissociation has also been reported.²⁸

A catalyst capable of abstracting H˙ from CH₄ can also abstracting H˙ from C₂H₆ owing to the lower C–H bond strength in C₂H₆, and even from C₂H₄ despite this compound's bond strength being greater than that of CH₄ by 5 kJ mol⁻¹.²⁹ Heterogeneous H˙ abstraction from C₂H₆ produces a C₂H₅˙ radical, which undergoes successive Eley–Rideal reactions with O* to form catalytic C₂H₄.³⁰ Homogeneous C₂H₄ formation from C₂H₅˙ is also possible through branched chain reactions. However, catalytic H˙ abstraction from C₂H₄ yields C₂H₃˙, a secondary source of CO in gas-phase reactions.^31,32 Thus, for a selective catalyst, effectively generating CH₃˙ radicals and facilitating their recombination to form C₂H₆ is crucial, for it is a termination step that prevents the consumption/oxidation of the radicals through CO_X-forming chain branching reactions with HO₂˙ radicals. The sequence that starts with CH₃˙ adsorption as CH₃O* and ends with formation primarily involves Langmuir–Hinshelwood reactions, with one H atom being removed from the reactant as OH* via the formation of CH₂O*.^30,31 Other relevant CO_X sources include gas-phase oxidation of C₂H₅˙ and C₂H₃˙ radicals produced through H˙ abstraction and surface C₂H₄ oxidation through adsorption, H˙ abstraction, and C–C bond cleavage.³² Therefore, effective catalysts not only activate CH₄ but also prevent deep oxidation paths, represented by HO₂˙ quenching in the surface mechanism. Finally, the regeneration of active sites involves water desorption, resulting from adsorbed OH˙ radicals yielding O*. The catalytic reaction network with 10 surface species (excluding vacancies) and 26 reversible elementary steps used in this work is presented in Table 1.

Table 1 Set of heterogeneous reversible elementary reactions used in this work

Step	Equation	Step	Equation
1	O2 (g) + 2* (s) ⇄ 2O* (s)	15	CH2O* (s) + O* (s) ⇄ CHO* (s) + OH* (s)
2	H2O (g) + * (s) ⇄ H2O* (s)	16	CH2O (g) + O* (s) ⇄ CHO˙ (g) + OH* (s)
3		17	CHO* (s) + O* (s) ⇄ CO* (s) + OH* (s)
4	CO (g) + * (s) ⇄ CO* (s)	18	CHO˙ (g) + O* (s) ⇄ CO (g) + OH* (s)
5	CH4 (g) + O* (s) ⇄ CH3˙ (g) + OH* (s)	19
6	C2H6 (g) + O* (s) ⇄ C2H5˙ (g) + OH* (s)	20	H2 (g) + O* (s) ⇄ H˙ (g) + OH* (s)
7	C2H5˙ (g) + O* (s) ⇄ C2H4 (g) + OH* (s)	21	H2O (g) + O* (s) ⇄ OH˙ (g) + OH* (s)
8	C2H4 (g) + O* (s) ⇄ C2H3˙ (g) + OH* (s)	22	OH˙ (g) + O* (s) ⇄ O˙ (g) + OH* (s)
9	CH3˙ (g) + O* (s) ⇄ CH3O* (s)	23	H2O2 (g) + O* (s) ⇄ HO2˙ (g) + OH* (s)
10	CH3O* (s) + O* (s) ⇄ CH2O* (s) + OH* (s)	24	HO2˙ (g) + O* (s) ⇄ O2 (g) + OH* (s)
11	CH3O˙ (g) + O* (s) ⇄ CH2O (g) + OH* (s)	25	HO2˙ (g) + * (s) ⇄ OH˙ (g) + O* (s)
12	C2H4 (g) + O* (s) ⇄ CH3CHO* (s)	26	OH* (s) + OH* (s) ⇄ H2O* (s) + O* (s)
13	CH3CHO* (s) + O* (s) ⇄ CH2CHO* (s) + OH* (s)
14	CH2CHO* (s) + O* (s) ⇄ CH2O* (s) + CHO* (s)

---

Heterogeneous net production rates are computed from the forward (r^f_j) and backward (r^b_j) rates of each heterogeneous reaction step.

N_r,s is the total number of heterogeneous elementary steps in the mechanism and the subscript j denotes the step number in the heterogeneous reaction mechanism (Table 1). For the surface intermediates, fractional coverages are calculated from eqn (8) and (9). Rate values of each individual surface step in both forward and backward directions are calculated using the mean field approximation, law of mass action, and the two-parameter Arrhenius law.

k_j = A_je^−E_a,j/RT j = 1, …, N_r,s (18)

where n_j is the number of sites involved in the heterogeneous elementary reaction step j. Adsorption steps are treated as sticking reactions (i.e., E^f_a,j = 0, j = 1–4, 9, 12), and their prefactors (A^f_j) are calculated from collision theory.where s_j,i is the sticking coefficient of adsorbing gas-phase species i in heterogeneous reaction step j and M_w,i is the molecular weight of gas-phase species i. Brønsted–Evans–Polanyi (BEP) relationships are used to relate the activation energy of a forward surface step to its reaction enthalpy, and reaction family-specific constants are used in the relationships.³³ BEP relationships are applicable only to the forward rate activation energy (eqn (20)), as the microscopic reversibility relates both forward and backward rate activation energies to the reaction enthalpy. The BEP parameters of the reaction families considered in this study are presented in Table 2.where E_0,f is the intrinsic energy barrier for any reaction in reaction family f, α_f is the transfer coefficient for any reaction in reaction family f and is the standard reaction enthalpy of heterogeneous reaction step j.

Table 2 Reaction families and their corresponding Polanyi parameters

Reaction family	Step (Table 1)	α f	Ref.	E 0,f	Ref.
(−) adsorptions	1–4, 9, 12, 25	0	—	0	—
(f1) H abstraction via Eley–Rideal	5–8, 11, 16, 18, 20−24	0.75	34	96.8	26
(f2) H abstraction via Langmuir–Hinshelwood	10, 13, 15, 17	0.50	35	141.3	26
(f3) OH˙ radical recombination via Langmuir–Hinshelwood	26	0.65	36	73.9	30
(f4) CO oxidation via Langmuir–Hinshelwood	19	0.26	37	67.6	26
(f5) C–C bond cleavage via Langmuir–Hinshelwood	14	0.97	37	186.7	26

---

In Table 2, for H abstraction via Eley–Rideal reactions (f₁), the intrinsic barrier is low and the transfer coefficient is high. In the BEP relationship, for reaction enthalpy values below −129 kJ mol⁻¹, highly exothermic reactions yield negative activation energy values. This applies to many exothermic reactions in that family. Here, a new scaling relationship is proposed for these reactions:

In this relationship, most of the exothermicity is attributed to the backward rather than the forward reaction.

2.7 Thermodynamic consistency

For any reversible reaction, the concentration-based thermodynamic equilibrium constant can be written as a function of the reaction entropy and enthalpy, and by combining it with the Arrhenius expressions, the enthalpic and entropic contributions to the forward (denoted with superscript f) and backward (denoted with superscript b) kinetic parameters can be isolated.where N_s,g,j is the total number of gaseous species in reaction step j, K_C,j is the concentration-based thermodynamic equilibrium constant of reaction step j, K_P,j is the pressure-based thermodynamic equilibrium constant of reaction step j, P is pressure, is the standard reaction Gibbs free energy of reaction step j and is the standard reaction entropy of reaction step j. Both enthalpy and entropy are temperature-dependent state functions. Hence, the activation energies and forward and backward prefactor values should also vary with temperature. For the homogeneous mechanism, thermodynamic consistency is ensured by directly computing the backward reaction rate constant from the forward reaction rate constant and the equilibrium constant:The equilibrium constant can be computed from the reaction Gibbs free energy, i.e., from the Gibbs energy of the formation of reactants and products. Temperature-dependent entropies and enthalpies of formation for each gaseous species are calculated via the NASA seven-coefficient polynomial parametrization of sensible heat at constant pressure within the temperature range of 200–1000 K and 1000–3500 K.where N_s,j is the total number of species in reaction step j and , and are the standard Gibbs energy, enthalpy and entropy of species i, respectively. Ensuring thermodynamic consistency in surface reaction mechanisms is more challenging because of the unavailability of Gibbs free energy values for surface species. One way to overcome this challenge is to construct the state functions for surface reactions as combinations of adsorption and analogous gas-phase reactions.³⁸ By identifying the minimum number of linearly independent surface reactions, we can categorize all reactions within the mechanism (N_r,s) into N_r,i linearly independent and N_r,d linearly dependent reactions. Once reactions are reordered such that the first N_r,i reactions are the linearly independent ones (i.e., ), any linearly dependent surface reaction can be expressed as a combination of linearly independent surface reactions.where ℜ_d is the d^th dependent reaction outside the basis set, ℜ_i is the i^th independent reaction inside the basis set and is c_d,i the coefficient of linearly decomposed reaction d onto reaction i belonging to the basis set. The set of independent reactions, known as basis set, remains fixed in size but does not limit the kinetic relevance of linearly dependent steps. The basis set size in a reaction network can be determined by finding the rank of the stoichiometric coefficient matrix; the rank often corresponds to the number of surface species. Typically, reversible adsorption–desorption steps are selected as the basis set, even if all of them are not directly included in the mechanism. In the present work, the basis set comprised the chemisorption steps of the 10 surface species, namely, O˙, OH˙, H₂O, CO₂, CO, CH₃O˙, CH₂O, CHO˙, CH₃CHO, and CH₂CHO˙. With the basis set, for the general reaction A* (s) + B* (s) ⇄ C* (s) + D* (s), the enthalpy and entropy of every catalytic surface reaction (denoted with subscript sur) in Table 1 are computed from the analogous gas-phase reaction (denoted with subscript gas) and the chemisorption functions (denoted with subscript ads) as shown in eqn (29) and (30).where N_s,s,j is the total number of surface intermediates in heterogeneous reaction step j and c_j,i is the coefficient for the adsorption of surface species i in reaction step j.

2.8 Temperature dependency of the basis set

Because the temperature dependency of the analogous gas-phase reactions is known from available thermochemistry, only temperature dependencies for chemisorption enthalpies and entropies (i.e., the basis set) are to be provided. For enthalpies, the statistical mechanics treatment of the chemisorption sensible heat considers the effect of changes in the degree of freedom (DOF). Every translational, rotational, and vibrational DOF contributes an amount of 0.5R, 0.5R, and R to the sensible heat, respectively, whereas a free, rigid, internal rotor contributes 0.5R.³⁹ The change in DOFs assumes that (i) all translational DOFs are converted to vibrational DOFs, (ii) weakly adsorbed molecules lose only one translational DOF, (iii) rotational DOFs are converted to vibrational DOFs, and (iv) for adsorbed species with a vertical axis through the adsorbed atom, one vibrational DOF gained is a free, rigid, internal rotor. With these assumptions, the temperature dependence of chemisorption enthalpy is deduced for the four cases in Table 3, and the chemisorption enthalpy at any given temperature isThe chemisorption enthalpies in the basis set at 300 K serve as model descriptors. To reduce model parameters, Su et al.⁴⁰ linked CH₃O˙ and OH˙ chemisorption enthalpies by determining of the bond energy difference between R–OH and R–OCH₃, and they set the average at 41 kJ mol⁻¹, in agreement with values in ref. 41 up to R=C₃H₇. Thus, the chemisorption enthalpy of CH₃O˙ is no longer a descriptor, but is related to the chemisorption enthalpy of OH˙ (eqn (32)). This study proposes a similar relationship for CH₂CHO˙ and CHO˙ (eqn (33)) on the basis of the average bond energy difference between R–CHO and R–CH₂CHO from aldehydes (i.e., difference in bond energy between carbonyl groups of different lengths), which is estimated to be 10 kJ mol⁻¹.⁴²where ΔQ_avg denotes the average bond energy difference between two types of hydrocarbons. The chemisorption entropy's temperature dependence is described by eqn (34), where β_i is a constant that is independent of adsorbate characteristics and binding strength. By propagating the entropic dependency over temperature evenly for the adsorption and desorption prefactors, we find that the adsorption prefactor is proportional to T^βi, while the desorption prefactor is proportional to T^−βi; β_i = 0.25 is chosen consistently.³⁸

Table 3 Temperature dependencies of chemisorption enthalpies for the steps in the basis set

Molecule	Species	Contribution of the degree of freedom				(dH/dT)i
		Translational	Rotational	Vibrational	Internal rigid rotor
Monoatomic	O˙	−3 × 0.5R		3 × R		1.5R
Diatomic	OH˙, CO	–3 × 0.5R	–2 × 0.5R	4 × R	1 × 0.5R	2R
Nonlinear polyatomic	CO2, CHO˙, CH2O, CH3O˙, CH2CHO˙	–3 × 0.5R	–3 × 0.5R	5 × R	1 × 0.5R	2.5R
Weakly bound nonlinear polyatomic	H2O, CH3CHO	−1 × 0.5R	−3 × 0.5R	4 × R		2R

---

2.9 Model descriptors

Model descriptors are primarily the enthalpies and entropies of the 10 chemisorption steps at 300 K. Empirical analogies can help reduce the enthalpic descriptors from 10 to 8, and the enthalpic descriptors are calculated at the reaction temperature. Forward and backward rate activation energies are calculated from BEP relationships and thermodynamic consistency, respectively. Adsorptions are assumed to be nonactivated, and hence, heats of chemisorption represent desorption activation energies. Temperature effects also apply to chemisorption entropies, which are linked to forward and backward rate prefactors by the even propagation of the entropic contribution to initial prefactor estimates (eqn (35) and (36)). The initial prefactor estimates (denoted with subscript init) are obtained from transition-state theory and the literature (ESI,† section S2). For adsorptions, the forward adsorption prefactor from eqn (35) is used to obtain the sticking coefficient from collision theory to ensure that the sticking coefficient does not exceed unity. If it exceeds unity, it is adjusted to unity while maintaining the prefactor ratio determined by the reaction entropy.In theory, thermodynamic properties of reactions in the basis set should suffice as descriptors. However, in practice, the density of active sites should also be considered. Under the assumption of a uniform distribution of surface intermediates and active sites, the model includes 19 descriptors: 8 chemisorption enthalpies at 300 K, 10 chemisorption entropies at 300 K, and the active site density. Fig. 2 shows the relationship between catalytic descriptors and kinetic parameters.

Fig. 2 Relationship between model descriptors and kinetic parameters of the microkinetic model.

2.10 Solution procedure and parameter estimation

The system of equations defined by the continuity equations are a set of partial differential-algebraic equations (PDEs), and they are converted into a set of differential-algebraic equations (DAEs) by using the orthogonal collocation method (ESI,† section S3), with seven and four collocation points for the intraparticle and interstitial phases, respectively. The resulting system of DAEs is solved using the SUNDIALS IDA package⁴³via the scikits.odes wrapper for Python.⁴⁴ Conversion of reactants and carbon molar selectivity of products are determined from the average mass fraction in the interstitial phase at the reactor outlet. All kinetic, thermodynamic, and transport calculations are performed in the Cantera framework⁴⁵ by using the GRI thermodynamic and transport database.⁴⁶

The vector of the 19 catalytic descriptor estimates (β) are determined through regression by minimizing the objective function in a two-step process: (i) genetic-algorithm-based optimization to conduct an order-of-magnitude search of parameters using a PFR model with the DEAP library⁴⁷ and (ii) gradientless optimization using a 1D heterogeneous reactor model with the SciPy library⁴⁸ for the refinement of the parameters. Furthermore, the associated confidence interval is calculated for each estimated descriptor value, at the 5% significance level.⁴⁹

where N_obs is the number of observations, N_res is the number of experimental responses per observation, f is the model multiresponse function, x_i is the variable representing the ith observation and ψ_i,j is the experimental (denoted with superscript “exp”) or calculated (denoted with superscript “calc”) performance metric of species j for the ith observation. To identify key descriptors in the model, after the Jacobian has been calculated, we calculate the first-order normalized sensitivity (ϕ_ij) of each descriptor (β_j) estimate for every response at each experimental condition (ψ_i), as shown in eqn (38). The Jacobian is evaluated with the numdifftools library,⁵⁰ with the spacing scaled to the order of magnitude of each descriptor to avoid round-off errors in the approximation. From the Jacobian matrix, the correlation matrix is calculated.Lastly, the global significance test of the model is performed to test the null hypothesis that all parameters would simultaneously be equal to zero. This null hypothesis is verified by comparing the regression sum of squares to the residual sum of squares divided by the number of DOFs with respect to the corresponding statistic.⁵¹where F_E is Fisher's E, n is the number of experiments and responses, p is the total number of catalytic descriptors of the model, F_inv is the inverse F distribution and α is the significance of the statistical test.

3. Results

3.1 Textural property estimation

To measure the surface area of nonporous surfaces via adsorption–desorption experiments, Ar physisorption experiments were performed. The multipoint BET surface areas are considered for the model, along with pore volume estimates that fall between the single point and Barrett–Joyner–Halenda methods. Estimated surface areas and pore volumes of the three catalysts are listed in Table 4, along with other textural properties. In the literature, only surface area values for SiO₂-supported catalysts are available, and they are consistent with those of IMP SiO₂ obtained in the present study.^26,52 Nevertheless, the values reported herein are small (<10 m_c² kg_c⁻¹) for the chosen analytical method, bearing a high uncertainty and thus demanding careful consideration, especially knowing that Mn–Na₂WO₄/SiO₂ undergo severe phase transformations under reaction condition. Differences in the surface area between IMP SiO₂ and the SiC-containing catalysts are attributed to the presence of SiC, which, unlike SiO₂, does not collapse under calcination; consequently, the SiC-containing catalysts have larger surface areas. The importance of a large surface area was shown by Wang et al.,⁵³ who linearly correlated the surface area and OCM productivity for Mn–Na₂WO₄/SiC catalysts calcined under different conditions. We set the tortuosity to 2.5 m_g² m_c⁻² based on previous results.⁵⁴ The bulk material density was estimated by weighing various catalyst volumes, and the average bed packing density was determined by loading various catalyst weights and measuring the height of the bed in the reactor. The low average bed packing density of IMP SiO₂ agrees with those reported for Mn–Na₂WO₄/SiO₂ (333 kg_c m_c⁻³ in ref. 52 and 400 kg_c m_c⁻³ in ref. 54). A key feature of SiC at the macroscopic level is that it endows the catalyst with higher density, resulting in the catalyst having an average packing density closer to the expected average packing density of a solid catalyst. From the estimates of both densities, the average bed packing porosity was obtained, and the bulk catalyst porosity was determined from the pore volume and bulk catalyst density. The intraparticle phase radius was set at 125 μm since the catalyst was sieved in the 150–300 μm range, and the interstitial phase radius was determined from the average bed packing porosity and the average particle radius.⁵⁵

Table 4 Textural properties of the catalysts studied

Property	IMP SiO2	SD SiO2–α + βSiC	SD SiO2–βSiC
Surface area (mc^2 kgc^−1)	2720	3580	4340
Pore volume (mg^3 kgc^−1)	3.4 × 10^−5	2.8 × 10^−5	2.2 × 10^−5
Bulk catalyst porosity (mg^3 mc^−3)	0.255	0.336	0.286
Catalyst tortuosity (mg^2 mc^−2)	2.5	2.5	2.5
Bulk catalyst density (kgc mc^−3)	750	1200	1300
Average bed packing density (kgc mb^−3)	450	700	800
Average bed packing porosity (mg^3 mb^−3)	0.40	0.42	0.38
Intraparticle phase radius (m)	125 × 10^−6	125 × 10^−6	125 × 10^−6
Interstitial phase radius (m)	56 × 10^−6	60 × 10^−6	51 × 10^−6

---

3.2 Temperature-programmed experiments and simulations

Temperature-programmed O₂ desorption (O₂-TPD) experiments shown in Fig. 3 revealed significant differences between the catalysts. SiC-containing catalysts exhibited single desorption peaks at lower temperatures compared with IMP SiO₂, which showed a main peak (at 792 °C) and a secondary low-temperature peak. The two types of reactive lattice oxygen are exclusive to trimetallic Mn–Na₂WO₄/SiO₂ catalysts; the strongly bonded oxygen and weakly bonded oxygen can be reversibly removed through reduction at temperatures above 700 °C and above 650 °C, respectively.²⁷ Other authors⁵⁶ have proposed that peaks at high temperatures may be linked to bulk lattice oxygen, adversely affecting OCM performance. Fig. 3 suggests that incorporating SiC into the support via spray-drying reduced the lattice oxygen strength in Mn or W species compared with that in Mn–Na₂WO₄/SiO₂. Furthermore, the breadth of peaks for the SiC-containing catalysts indicates the presence of multiple oxygen species that can be released, similar to the case of IMP SiO₂, where phases containing W and Mn interact differently with oxygen. Additionally, SD SiO₂–α + βSiC exhibited the highest value for the total amount of O₂ desorbed (64 μmol g_c⁻¹), and it was followed by IMP SiO₂ (47 μmol g_c⁻¹). These differences observed are associated with the metal to which oxygen is bonded. In other words, the preparation method, the presence of SiC, and the crystal phase of SiC influence the interaction between active sites and O₂ by modifying the electronic environment of Mn and W species, thereby altering their interaction with oxygen. This is evidenced by the observed changes in oxygen uptake, which reflect variations in oxygen mobility and dispersion. This interaction is denoted by σ in the model. This is in accord with the fact that O₂ uptake has been reported to be proportional to the Mn content.⁵⁷

Fig. 3 Evolution of temperature and measured (solid) and predicted (dashed) O₂ desorption normalized concerning the total amount of O₂ desorbed (area under the curve) over time for each of the fresh catalysts. Conditions: P = 1 bar, T_I = 25 °C, T_F = 850 °C, β = 7.5 min⁻¹, t_∞ = 30 min, F_T = 50 NmL min⁻¹, p_CH4,0 = 0 bar, W = 50 mg_c. O₂ chemisorption occurred for 6 h at 25 °C and 50 NmL min⁻¹ (10% O₂). Simulations were performed assuming plug-flow conditions and model descriptors in Table 5.

Table 5 Estimates of descriptors for each catalyst along with their 95% confidence intervals. Units: chemisorption enthalpies at 300 K, kJ mol⁻¹; chemisorption entropies at 300 K, J mol⁻¹ K⁻¹; active site density, kmol m_c⁻²

Catalyst descriptor		IMP SiO2	SD SiO2–α + βSiC	SD SiO2–βSiC
H 1	Chemisorption enthalpy, O˙	−319 ± 17	−327 ± 23	−327 ± 29
H 2	Chemisorption enthalpy, OH˙	−279 ± 14	−300 ± 20	−297 ± 25
H 3	Chemisorption enthalpy, H2O	−27 ± 3	−22 ± 3	−16 ± 6
H 4	Chemisorption enthalpy, CH2O	−138 ± 31	−154 ± 39	−146 ± 50
H 5	Chemisorption enthalpy, CHO˙	−205 ± 33	−140 ± 50	−256 ± 51
H 6	Chemisorption enthalpy, CO	−100 ± 12	−99 ± 4	−109 ± 21
H 7	Chemisorption enthalpy, CO2	−241 ± 23	−261 ± 19	−255 ± 42
H 8	Chemisorption enthalpy, CH3CHO	−27 ± 4	−37 ± 7	−70 ± 7
S 1	Chemisorption entropy, O˙	−101 ± 15	−110 ± 29	−103 ± 31
S 2	Chemisorption entropy, OH˙	−170 ± 15	−201 ± 30	−184 ± 34
S 3	Chemisorption entropy, H2O	−286 ± 15	−277 ± 24	−250 ± 29
S 4	Chemisorption entropy, CH3O˙	−130 ± 14	−148 ± 21	−138 ± 29
S 5	Chemisorption entropy, CH2O	−146 ± 34	−206 ± 48	−175 ± 57
S 6	Chemisorption entropy, CHO˙	−157 ± 25	−112 ± 50	−201 ± 48
S 7	Chemisorption entropy of CO	−216 ± 26	−228 ± 36	−231 ± 44
S 8	Chemisorption entropy, CO2	−153 ± 21	−164 ± 23	−154 ± 37
S 9	Chemisorption entropy, CH3CHO	−309 ± 37	−299 ± 66	−230 ± 41
S 10	Chemisorption entropy, CH2CHO˙	−286 ± 35	−246 ± 43	−289 ± 49
σ	Active site density	(5.02 ± 0.03) × 10^−9	(3.76 ± 0.04) × 10^−9	(2.59 ± 0.07) × 10^−9

---

The simulated desorption profiles represented by dashed lines in Fig. 3 were obtained from descriptors estimated through regression and initial coverages on the basis of the observed O₂ uptake. Since the variation of the O₂ chemisorption enthalpy across catalysts was minimal (Table 5), differences in chemisorption entropy and initial coverage could explain the experimentally observed disparities. For all three catalysts, the model predicted a low-temperature main peak around 600 °C, which was close to the experimental observation for SD SiO₂–α + βSiC and in contrast to the dual peaks of IMP SiO₂ at higher temperatures. The model predictions supported the concept that weakly bonded oxygen influences the catalytic activity in the steady-state regime.⁵⁸ However, the simulations highlighted a key model assumption, namely the mean field approximation, which may not accord with the observed O₂-active site interplay, at least for IMP SiO₂. Extending the model with two types of active sites—strong and weaker oxygen bonding sites—could provide a more realistic description, but the number of model parameters would then increase. Fleischer et al.⁵⁹ accounted for two types of oxygen species via two-step dissociative chemisorption.

Fig. 4 depicts the CH₄ temperature-programmed surface reaction (CH₄-TPSR), for which ion currents corresponding to m/z values of 2 (H₂⁺), 17 (OH⁺, CH₅⁺), 18 (H₂O⁺), 28 (CO⁺, C₂H₄⁺), 29 (C₂H₅⁺), 30 (C₂H₆⁺), 32 (O₂⁺), and 44 (CO₂⁺) were monitored. Fig. 4 shows marginal desorption of O₂ at temperatures below 600 °C, and as the temperature rises, products are formed. IMP SiO₂ exhibits two C₂H₆ peaks, at approximately 740 and 850 °C, representing a primary product, whereas SiC-containing catalysts show a single high-temperature peak, in agreement with observations in previous O₂-TPD experiments. This supports the notion that CH₄ activation occurs through strongly bound oxygen, while weakly bound oxygen leads to product formation at lower temperatures, albeit at significantly reduced rates.⁶⁰ The primary peaks of other products such as C₂H₄/CO and CO₂ is at 850 °C. Thus, low temperature C₂H₆ formation is also indicative of the prevalence of CH₃˙ radical onto the catalyst surface.⁵⁹ Furthermore, the H₂ signal reveals different kinetic behavior as SD SiO₂–α + βSiC shows activity even at low temperatures, consistent with signals such as C₂H₄/CO and H₂O whose presence is required for H₂ formation.

Fig. 4 Evolution of temperature and MS signals over time for each catalyst at different CH₄ partial pressures in the carrier gas. Conditions: P = 1 bar, T_I = 25 °C, T_F = 850 °C, β = 7.5 min⁻¹, t_∞ = 30 min, F_T = 100 NmL min⁻¹, W = 50 mg_c. O₂ chemisorption occurred for 6 h at 25 °C and 50 NmL min⁻¹ (10% O₂).

Simulations performed with the kinetic model corresponding to Fig. 4 (p_CH4,0 = 1 bar) are described in ESI† (section S4). Primary carbon-containing products, namely C₂H₆ and CO, dominated at lower temperatures (around ∼600 °C), which indicated their primary nature compared with C₂H₄ and CO₂, which appeared at higher temperatures. In particular, the absence of predicted C₂H₄ suggests that m/z = 28 in Fig. 4 could be attributed to CO. An analysis of the coverage distribution showed initial enrichments of O* and the absence of exposure to H-containing species, unlike steady-state simulations. This implies that the CH₄-TPSR experiments conducted in this work did not involve the operating regime found in steady-state conditions, which is consistent with the absence of predicted C₂H₄ controlled by the existing O* at the initial condition. O* coverage reduced over the simulation, declining noticeably upon CO₂ production. Notably, the maximum of OH* corresponded to the peak production of C₂H₆, underscoring the relationship between OH* and O* and their impact on catalyst activity. In fact, under these conditions, a positive correlation was observed between the desorbed O₂ amount and the catalyst activity (i.e., the reason for SD SiO₂–α + βSiC showing larger peaks).

3.3 Parameter estimation

Model descriptors and their significance are estimated through regression with steady-state kinetic data presented in ref. 20 are shown in Table 5.

The O˙ chemisorption enthalpy is related to the O₂ chemisorption enthalpy through the bond dissociation energy of O₂, for activating CH₄via the heterogeneous primary initiation of H abstraction. The estimated O˙ chemisorption enthalpies of all three catalysts, do not show statistically significant differences, thus implicitly suggesting a lack of tendency of SiC-containing catalysts to show stronger bonding with or affinity for O˙. O˙ chemisorption enthalpy values expressed as heats of O₂ chemisorption at 800 °C were 113, 130, and 130 kJ mol⁻¹ for IMP SiO₂, SD SiO₂–α + βSiC, and SD SiO₂–βSiC, respectively. As a result, it is clear that the O₂ chemisorption enthalpy alone does not describe the temperature-programmed experiment in Fig. 3 (i.e., the desorption is not entirely driven by the O₂ desorption barrier), wherein O* plays a pivotal role in the overall catalyst activity.

The estimated OH˙ chemisorption values were not statistically different for the three catalysts, despite the values for SD SiO₂–α + βSiC and SD SiO₂–βSiC being higher than that for IMP SiO₂. The absolute chemisorption value or heat of chemisorption of the hydroxyls indicates their degree of stability, which increases proportionally with the heat value.^26,32 These surface intermediates determine the activation barriers to many key heterogeneous steps, such as the Eley–Rideal steps, including heterogeneous CH₄ initiation, whose reaction enthalpy is proportional to (40, 27, and 30 kJ mol⁻¹ for IMP SiO₂, SD SiO₂–α + βSiC, and SD SiO₂–βSiC, respectively). Note that the OH˙ chemisorption enthalpy influences other critical routes related to deep catalytic oxidation of CH₃˙ radicals and C₂H₄, as well as the hydroxyl species regeneration step. The expected larger H₂O chemisorption entropies in Table 5, render the regeneration step more kinetically favorable for SD SiO₂–α + βSiC, especially compared with IMP SiO₂. Furthermore, weaker H₂O adsorption has been reported to result in higher CH₄ activity.³⁰

Differences in O˙ and OH˙ chemisorptions also influence the HO₂˙ quenching kinetics. In that regard, IMP SiO₂ is expected to have lower activation barriers for backward quenching reactions and overall enhanced rates. These differences are especially significant for the O*-mediated quenching reaction compared with the vacancy-mediated quenching reaction. HO₂˙ are active chain carriers in the gas phase. This interpretation cannot be decoupled from the active site density value, which influences reaction rates in linear proportion to the rate constant and quadratically for second-order reactions. The estimated active site density was significantly different across the catalysts, for IMP SiO₂ being twice that for SD SiO₂–βSiC, confirming that SiC reduced the HO₂˙ quenching capability of the Mn–Na₂WO₄ catalyst.

Another notable difference in Table 5 concerns deep oxidation routes. CH₃CHO chemisorption enthalpy showed stronger interaction with C₂H₄ for SD SiO₂–βSiC compared with SD SiO₂–α + βSiC and particularly with IMP SiO₂. Together with the sticking coefficient of CH₃˙,^26,32,61 CH₃CHO chemisorption enthalpy is crucial for C₂ selectivity.³⁰ A statistically corroborated less negative entropy indicates more kinetically prone C₂H₄ oxidation by SD SiO₂–βSiC, which is indicated by the catalyst's C₂H₄ sticking coefficient of 1.3 × 10⁻⁷ at 800 °C; the C₂H₄ sticking coefficient of SD SiO₂–α + βSiC and IMP SiO₂ were smaller, namely 6.1 × 10⁻⁹ for the latter. This trend was reversed for the C–C bond scission step, which was much more favorable for SD SiO₂–α + βSiC owing to its less negative CH₂CHO˙ chemisorption enthalpy and entropy. IMP SiO₂ behaved similar to SD SiO₂–α + βSiC with regard to C₂H₄ oxidation. The difference between the CH₃O˙ and CH₂O chemisorption entropies for SD SiO₂–α + βSiC suggests low oxidation tendency for single-carbon intermediates, in line with a significantly lower CHO˙ chemisorption enthalpy. No significant differences were observed between CO and CO₂ chemisorption descriptors for the catalysts.

The effect of model descriptors on the Arrhenius parameters of each elementary and reversible step is shown in Fig. 5. Note that the use of the 3-parameter Arrhenius expression (eqn (18)) can lead to a pronounced correlation between the activation energy and the prefactor, thereby influencing results in Fig. 5.^62,63 Based on the figure, the forward rate activation energy and prefactor values for H-abstraction steps (r₅–r₈) follow the sequence IMP SiO₂ > SiO₂–βSiC > SD SiO₂–α + βSiC, regardless of whether the hydrocarbon acts as an H source. The difference in activation energies of H abstraction from CH₄ leads to higher endothermicity, slowing CH₄ conversion and reducing the overall catalytic contribution.^26,32 This also applies to C₂H₄ yielding C₂H₃˙ and to the eventual CO_X product formation, as the latter oxidizes in the contiguous gas-phase. This trend is reversed in the case of the adsorption of C₂H₄ (r₁₂), which is more prominent in the case of SD SiO₂–βSiC. Despite higher C₂H₄ sticking coefficients, this catalyst showed smaller kinetic parameters related to C₂H₄ dehydrogenation and scission (r₁₃–r₁₄). Thus, differences in the catalytic contribution to C₂H₄ oxidation are expected between the catalysts. While the sticking rate of CH₃˙ (r₉) is similar for all catalysts, its reverse counterpart is slightly faster at lower temperatures for SD SiO₂–α + βSiC, which is indicated by the catalyst's marginally larger desorption activation energy. The sticking coefficients of r₉ are relatively high, on the order of 10⁻⁵–10⁻⁴, suggesting that the high activity resulting from the O₂ heat of chemisorption is not directed by strong inhibition of the oxidation rate of CH₃˙ radicals. Moreover, Fig. 5 shows that SD SiO₂–α + βSiC hinders the CH₂O* route more.

Fig. 5 Activation energy values of all three catalysts obtained with the microkinetic model descriptors in Table 5 for all heterogeneous (a) forward and (b) backward reactions, and prefactor values for all three catalysts obtained with the microkinetic model descriptors in Table 5 for all heterogeneous (c) forward and (d) backward reactions. The line width shows temperature-dependent variations in the 25–800 °C range. Reaction nomenclature is as presented in Table 1.

The regeneration step (r₂₆) parameters significantly differ for SD SiO₂–α + βSiC. Notably, forward and backward rate prefactors differ by two orders of magnitude unlike the other catalysts, indicating the large weight of the entropic descriptor contribution, as the OH˙ and H₂O chemisorption entropy differences are not greater than 30 J mol⁻¹ K⁻¹. This observation results from the ratio of prefactors being proportional to the exponent of the surface reaction entropy, and it is also evident from the prefactor values of H₂O chemisorption (r₂) on SD SiO₂–βSiC.

3.4 Regression assessment

Parity plots in Fig. 6 compare measured and predicted performance metrics with descriptors from Table 5. The calculated F_E values of 1329.6, 706.4, and 650.3 for IMP SiO₂, SiO₂–βSiC, and SD SiO₂–α + βSiC, respectively, reject the null hypothesis with a tabulated value of 1.6 and substantiate the global significance of the model validation. High F_E values, in order of hundreds or thousands, especially for IMP SiO₂, confirm the model's adaptability to the experimental data.⁵¹ Overall, the parity plots depict a satisfactory match for selectivities and CH₄ conversion, though the accurate capture of the experimental trends of O₂ conversion poses a challenge. This generalized disparity is because of unmeasured H₂ and H₂O yields impacting predictions. Specifically, heterogeneous steps 20–26, devoid of carbon products, are solely bound by O₂ conversion. The difficulty in capturing O₂ conversion trends is again observed when comparing the residuals to a standard normal distribution (ESI,† Fig. S3). Despite the overall linear trends, O₂ conversion deviates from linearity for all three catalysts, and variance differences across predicted values (i.e., heteroscedasticity) also exist, as indicated by slope changes in quantiles. Despite this, the relative error for O₂ conversion generally remains below the 25% relative deviation mark. Furthermore, some discrepancies in the predictions of C₂ product selectivities across the three catalysts are highly noticeable at lower CH₄ conversion rates, approximately below 5%. It is worth noting that some of these deviations may also be ascribed to heat effects causing temperature gradients^64–66 despite the flow ideality, intrinsic kinetic regime, and isothermicity of the reactor (2 mm i.d., 1 cm long) having been previously checked.²⁰

Fig. 6 Parity plots obtained by fitting experimental performance metrics to the microkinetic model with the descriptors in Table 5 for (a and b) IMP SiO₂, (c and d) SD SiO₂–α + βSiC, and (e and f) SD SiO₂–βSiC.

The binary correlation matrix of model descriptors (ESI,† Fig. S4) shows associations between certain parameters. However, the absence of strong correlations exceeding 0.95 suggests model descriptor redundancy.⁵¹ The largest binary correlation is between the chemisorption enthalpies of O˙ (H₁) and OH˙ (H₂), in the range of 0.68–0.77 across catalysts. Smaller binary correlations across the three catalysts include the CH₃O˙ chemisorption entropy (S₄), the chemisorption enthalpy of O˙ (H₁), and the chemisorption enthalpy and entropy of CO₂ (H₇ and S₈). Catalyst-specific correlations also exist, such as the chemisorption enthalpy and entropy of OH˙ (H₂ and S₂) for SiO₂–α + βSiC and SD SiO₂–βSiC, and the chemisorption enthalpy of CO and entropy of CO₂ (H₆ and S₈) for IMP SiO₂ and SD SiO₂–βSiC. Still, most descriptors exhibit absolute correlation values closer to 0 than 0.95.

3.5 Model descriptor benchmarking

Normalized sensitivity coefficients quantify the effect of model inputs on outputs at specific conditions. To identify the most influential model descriptors, we used box plots of the sensitivity coefficients of all descriptors over all the experimental conditions, and they are presented in ESI† (Fig. S5). O˙ and OH˙ chemisorption enthalpies (H₁ and H₂) had the most significant effect, consistent with being part of all Eley–Rideal steps. On the one hand, O₂ heat of chemisorption played a key role in O₂ activity and CO_X selectivity, in line with the results of Thybaut et al.³⁰ On the other hand, OH˙ chemisorption enthalpy positively influenced CH₄ conversion and C₂ product selectivity. Other remarkable effects include CO₂ chemisorption enthalpy (H₇) inhibiting both reactant conversions and aligning with global rate type-based kinetic models reported in the literature,⁶⁷ and active site density (σ) exhibiting different effects on model responses, depending on experimental conditions.

Fig. 7 compares descriptors from this study with descriptors in the literature,^{26,30–32,61} namely, H abstraction from CH₄ reaction enthalpy, O₂ chemisorption enthalpy, CH₃˙ sticking coefficient, and active site density. The significance of these four descriptors is shown in Fig. S5,† for the H abstraction from CH₄ reaction enthalpy is a function of the O₂ chemisorption enthalpy and the OH˙ chemisorption enthalpy. For instance, La–Sr/CaO and Sr/La₂O₃ in Fig. 7 exhibit OH˙ chemisorption heats of 257 kJ mol⁻¹ and 278 kJ mol⁻¹, which correlate with the respective enthalpies of H abstraction from CH₄ of 65 kJ mol⁻¹ and 44 kJ mol⁻¹.^26,32 Similar trends were observed with OH˙ chemisorption heats for IMP SiO₂, SiO₂–βSiC, and SD SiO₂–α + βSiC, resulting in H-abstraction from CH₄ reaction enthalpies of 57 kJ mol⁻¹, 44 kJ mol⁻¹ and 47 kJ mol⁻¹ at 800 °C. It is important to highlight that while these descriptor values are directly linked to the intrinsic properties of the catalyst, they are also influenced by operating conditions. For example, exposure of La–Sr/CaO to CO₂ cofeeding may alter these descriptors due to the formation of carbonates.⁶⁸

Fig. 7 Comparison of four main descriptors obtained in this work with values reported in the literature by using descriptor-based microkinetic models: (a) reaction enthalpy of H-abstraction from CH₄, (b) heat of chemisorption of O₂, (c) CH₃˙ sticking coefficient, and (d) active site density. Descriptor sources: A, Sun et al.;³¹ B, Thybaut et al.;³⁰ C, Alexiadis et al.;³² D, Kechagiopoulos et al.;⁵⁵ E, Alexiadis et al.;²⁶ F, Ahari et al.;⁵² G, Karakaya et al.;⁵⁴ H, this work. Descriptors include virtual optimal values (striped bars) and physical boundaries from Pirro et al.⁶⁹ (dashed lines). CH₃˙ sticking coefficients have been estimated from collision theory at 800 °C.

Careful examination of the descriptor values in Fig. 7 is important since they were derived under different conditions and with different considerations. For example, Ahari et al.⁵² used an isothermal plug-flow reactor, while Karakaya et al.⁵⁴ focused on nonisothermal effects in the microkinetics. Additionally, Alexiadis et al.⁶¹ demonstrated how the heat of O₂ chemisorption could vary by up to 30 kJ mol⁻¹ with changes in catalyst dilution in the bed. Hence, Fig. 7 does not depict clear trends for descriptors and catalyst families, not even between promoted and unpromoted catalysts (e.g., Li/MgO vs. Sn–Li/MgO) or within the same catalyst type (e.g., Mn–Na₂WO₄/SiO₂). In the latter case, factors such as metal loading, calcination conditions, SiO₂ type, and catalyst synthesis method significantly influence performance. Still, it appears from Fig. 7 that Mn–Na₂WO₄/SiO₂ catalysts have low CH₃˙ and C₂H₄ sticking coefficients while also exhibiting lower enthalpies of H-abstraction from the CH₄ reaction. Furthermore, Mn–Na₂WO₄/SiO₂ catalysts have low O₂ heat of chemisorption and low CH₃˙ sticking coefficients, which are crucial for improved OCM yields by SiC.

In Fig. 7, dashed lines delineate the lower and upper bounds defining feasible values,⁶⁹ and all descriptors fall within those limits. Nonetheless, chemisorption functions as descriptors have the advantage of facilitating the assessment of their physical viability. For enthalpies, the necessity for the heat of adsorption of species i to be positive as shown in Table 5. Furthermore, thermodynamics dictate that for any endothermic reaction, the forward activation energy should be larger than the heat of reaction; this is satisfied through the selection of the Polanyi parameters E_0,f and α_f to ensure that the relationship holds. For chemisorption entropies, upholding the “Langmuirian integrity” requires , which accounts for the loss of translational contribution upon adsorption. A less strict constraint, proposed by Vannice et al.,⁷⁰ is (in J mol⁻¹ K⁻¹) from the decrease of free volume upon adsorption at the standard state coverage of 0.5, with the upper limit dictated by the relationship described before. Values associated with both entropic constraints are presented in Table 6.

Table 6 Chemisorption entropies of the species in the basis set at 800 °C and their corresponding constraints. Values in bold denote violation of thermodynamic constraints. Units: chemisorption entropies, J mol⁻¹ K⁻¹

Species	IMP SiO2	SD SiO2–α + βSiC	SD SiO2–βSiC
O˙	0 < |−96.0| < 188.3	0 < |−104.7| < 188.3	0 < |−97.9| < 188.3
	41.8 < |−96.0| < 484.3	41.8 < |−104.7| < 495.3	41.8 < |−97.9| < 495.7
OH˙	0 < |−164.8| < 221.9	0 < |−195.9| < 221.9	0 < |−179.1| < 221.9
	41.8 < |−164.8| < 423.6	41.8 < |−195.9| < 452.3	41.8 < |−179.1| < 448.9
H2O	0 < |−280.7| < 235.7	0 < |−272.2| < 235.7	0 < |−244.3| < 235.7
	41.8 < |−280.7| < 65.8	41.8 < |−272.2| < 59.1	41.8 < |−244.3| < 51.1
CH3O˙	0 < |−124.9| <301.8	0 < |−143.0| < 301.8	0 < |−132.8| < 301.8
	41.8 < |−124.9| < 363.1	41.8 < |−143.0| < 391.8	41.8 < |−132.8| < 388.4
CH2O	0 < |−140.8| < 279.8	0 < |−200.3| < 279.8	0 < |−170.2| < 279.8
	41.8 < |−140.8| < 221.8	41.8 < |−200.3| < 243.9	41.8 < |−170.2| < 233.4
CHO˙	0 < |−151.8| < 276.5	0 < |−106.5| < 276.5	0 < |−195.5| < 276.5
	41.8 < |−151.8| < 314.9	41.8 < |−106.5| < 225.0	41.8 < |−195.5| < 386.9
CO	0 < |−210.9| < 236.9	0 < |−222.8| < 236.9	0 < |−225.4| < 236.9
	41.8 < |−210.9| < 173.0	41.8 < |−222.8| < 171.0	41.8 < |−225.4| < 185.4
CO2	0 < |−147.5| < 273.2	0 < |−158.9| < 273.2	0 < |−148.9| < 273.2
	41.8 < |−147.5| < 366.4	41.8 < |−158.9| < 394.3	41.8 < |−148.9| < 384.8
CH3CHO	0 < |−303.8| < 371.2	0 < |−293.4| < 371.2	0 < |−224.7| < 371.2
	41.8 < |−303.8| < 66.2	41.8 < |−293.4| < 80.5	41.8 < |−224.7| < 125.9
CH2CHO˙	0 < |−280.9| < 365.5	0 < |−240.7| < 365.5	0 < |−283.8| < 365.5
	41.8 < |−280.9| < 301.5	41.8 < |−240.7| < 211.0	41.8 < |−283.8| < 372.9

---

In Table 6, only the H₂O chemisorption entropy violates the strict constraint, which weakens model interpretability. This anomaly could result from the absence of H₂O as a model response, impacting O₂ conversion predictions. Addressing this issue may involve measuring the H₂O yields experimentally or recalibrating transition state initial estimates, such as adjusting prefactors for specific reactions. For instance, the initial estimate for the H₂O sticking coefficient was set at approximately 0.05, based on the reported value for Sn–Li/MgO.⁵⁵ However, this coefficient has been reported to be 0.5 for Li/MgO.³¹ These corrections could also resolve violations of the constraints of Vannice,⁷⁰ which serve as a heuristic guide rather than a thermodynamic constraint.

3.6 Steady-state simulations

To assess diffusional limitations for radicals and molecules, we can compare gas-phase species consumption and diffusion rates by using the diffusion length, which is determined from the square root of the effective diffusivity and the species lifetime.⁷¹ The latter is determined from the local consumption rate per unit volume of catalyst and the gas concentration in the intraparticle phase (eqn (40)). Note that only negative production rates are used as local consumption rates; otherwise, the denominator of eqn (40) is replaced by the net production rate.Fig. 8 shows the effective diffusivity and lifetime of each gas-phase species at the catalyst center (ξ = 0) and end of the catalyst bed (z = L_b). For the three catalysts, stable molecules, including reactants and products, result in the absence of concentration gradients along the radial axis of the particle since their diffusion lengths are approximately one order of magnitude larger than the average particle radius (125 μm).⁵⁵ However, H₂O₂ and CH₂O have smaller diffusion lengths, around 220 μm and 160 μm, respectively, larger than the average particle radius but smaller than the diameter. Hence, they are likely to develop concentration gradients. Previous studies⁷¹ have reported even smaller diffusion lengths for CH₂O. The radical intermediates have diffusion lengths below 10 μm, significantly smaller than the average particle size, thus validating the chosen heterogeneous reactor model. Fig. 8 indicates that the most notable expected intraparticle concentration gradients are linked to the HO₂˙ radical, which has a diffusion length of less than 1 μm. This is attributed to the rapid HO₂˙ catalytic quenching reactions despite the radical produced in both phases via homogeneous reactions.

Fig. 8 Effective intraparticle diffusivities and lifetimes of all species at the center of the particle (ξ = 0) and at the end of the catalytic reactor (z = L_b) for (a) IMP SiO₂, (b) SD SiO₂–α + βSiC, and (c) SD SiO₂–βSiC. Square symbols represent CO, CO₂, C₂H₂, and C₂H₄. Simulation conditions: T = 800 °C, P = 1 atm, feed molar ratio of CH₄/O₂/He = 3/1/0.6, W/F_CH4,0 = 4.4 g_c h mol_C⁻¹.

Fig. 9 shows the simulated concentration profiles along the catalyst bed for IMP SiO₂. Stable molecules such as CH₄ and C₂H₆ (Fig. 9a and c) do not develop discernible radial gradients, while CH₂O (Fig. 9e) exhibits a subtle radial deviation consistent with its estimated diffusion length. These species exhibit linear trends along with z, ascending (e.g., C₂H₆) or descending (e.g., CH₄). However, radicals such as CH₃˙ (Fig. 9b), produced and consumed at similar rates in the intraparticle phase, show a plateau in their concentration profiles, sharply decaying in concentration through the radial interphase because of diffusion into the interstitial phase. This underscores the catalytic role in activating the C–H bond of CH₄. Along the axial coordinate, there is an increase in interstitial phase concentrations that can be attributed to the CH₃˙ radical not being part of the initial conditions. Nevertheless, the concentration profile along the length of the bed in the interstitial phase shows a slight decrease, owing to lower CH₄ and O₂ concentrations and an increase in the CO₂ concentration, with the latter being known to hinder OCM yields.

Fig. 9 Axial and radial (interstitial and intraparticle phases) concentration profile predictions for IMP SiO₂ for (a) CH₄, (b) CH₃˙, (c) C₂H₆, (d) C₂H₅˙, (e) CH₂O, (f) CHO˙, (g) H˙, and (h) HO₂˙. Simulation conditions: T = 800 °C, P = 1 atm, feed molar ratio of CH₄/O₂/He = 3/1/0.6, W/F_CH4,0 = 4.4 g_c h mol_C⁻¹ obtained from the equivalent L_b value of 0.03 m. Simulation results: 9.7% CH₄ conversion, 19.5% O₂ conversion, 35.4% C₂H₄ selectivity, 37.6% C₂H₆ selectivity, 13.2% CO selectivity, 13.8% CO₂ selectivity.

As secondary product intermediates, the C₂H₅˙ and CHO˙ radicals (Fig. 9d and f) originate from the heterogeneous H-abstraction of C₂H₆ and diffuse at a low rate into the interstitial phase. The contribution of the interstitial phase to the gas-phase formation of C₂H₅˙ radicals via H-abstraction from C₂H₆ by OH˙ or CH₃˙ is significantly lower than that of the intraparticle phase via Eley–Rideal H-abstraction from C₂H₆. A similar analysis can be performed to produce C₂H₄. The concentration disparity for the H˙ radical between the intraparticle and interstitial phases (Fig. 9g) arises from the substantially higher total intraparticle radical concentration from heterogeneous initiations. The H˙ radical concentration also increases along the axial coordinate due to increased C₂H₆ concentration, which accelerates chain propagation and branching rates. Thus, higher rates of primary initiations by H˙ can be expected in the catalyst pores through homogeneous mechanisms, in contrast to the interstitial phase.^61,71

The radial distribution of the HO₂˙ radical (Fig. 9h) exhibits a sharp decline in the intraparticle phase, mainly because of significant heterogeneous termination or quenching reactions compared with the production rate in the interstitial phase. Interestingly, HO₂˙ concentrations develop gradients in the interstitial phase, which become more pronounced. Increased HO₂˙ concentrations positively influence CH₄ conversions, attributed to higher C₂H₆ concentrations accelerating propagation and branching rates in a branched-chain mechanism.⁷¹ However, the higher interstitial HO₂˙ concentration leads to lower C₂ selectivities. This emphasizes the importance of an effective OCM catalyst for efficiently activating CH₄ and effectively quenching HO₂˙ radicals and underscores the importance of the HO₂˙ concentration in the interstitial phase.

The simulated concentration profiles for SD SiO₂–α + βSiC are presented in ESI† (Fig. S6). Under similar conditions, SD SiO₂–α + βSiC exhibited reduced CH₄ conversion and C₂H₄ selectivity but increased C₂H₆ selectivity, resulting in similar overall C₂ selectivities. Notably, selectivity toward CO was significantly higher than that for CO₂ for SD SiO₂–α + βSiC, unlike IMP SiO₂, as evidenced by lower overall radical concentrations (crucial for OCM), such as CH₃˙, C₂H₅˙, H˙, and HO₂˙ concentrations. This can be attributed to significantly more negative OH˙ chemisorption entropy and the reduced active site density of SD SiO₂–α + βSiC compared with IMP SiO₂, which led to decreased formation of C₂H₆ and eventually lower production of C₂H₅˙, explaining the lower CH₄ conversion and C₂H₄ selectivity. More pronounced concentration profiles of CHO˙, H˙, and HO₂˙ between radial phases are also an outcome of this effect, with lower concentration peaks in their respective prevalent radial phase.

For SD SiO₂–βSiC (ESI,† Fig. S7), the predicted CH₄ conversion did not exceed that of IMP SiO₂ while similar C₂ selectivity was maintained, placing it between the latter catalyst and SD SiO₂–α + βSiC in terms of predicted activity. The CO/CO₂ ratio between the catalysts varied at the end of the catalyst bed, with SiC-containing catalysts showing molar ratios closer to 2 and IMP SiO₂ showing a ratio closer to 1. This can be key if the generated C₂H₄ undergoes further hydroformylation into propanal, where both CO and C₂H₄ are reactants and the suitability of a CO/H₂/C₂H₄ ratio of 1/1/1 for hydroformylation has been previously reported.⁷² For instance, the simulated CO/H₂/C₂H₄ ratios are 1/1.1/1.3 and 1/0.9/0.9 for IMP SiO₂ and SD SiO₂–βSiC, respectively.

The enhanced performance of SD SiO₂–βSiC, which had a lower active site density than IMP SiO₂, also highlights the influence of textural properties. The effect of catalyst porosity and surface area was observed in SD SiO₂–βSiC, which exhibited an approximately 60% larger surface area and higher catalyst porosity than IMP SiO₂. This is in accord with higher CH₄ conversions and loss of C₂ selectivity reported for larger surface areas under the same kinetics.⁵⁵ This is because higher surface areas promote initiation, generating radicals, yet homogeneous reactions in the interstitial phase do not match the consumption rate. Consequently, despite the amount of radicals in the catalyst pores, most produced radicals eventually undergo heterogeneous oxidation. Moreover, the increase in particle porosity balances the positive impact of gas-phase reactions in a CH₃˙-rich environment (where the CH₃˙ coupling rate is second order to the CH₃˙ concentration), owing to the absence of intraparticle mass transport limitations. This modulation significantly increases CH₃˙ coupling compared with deep oxidation routes.⁷¹

Fractional coverage profiles of surface intermediates for IMP SiO₂ are shown in Fig. 10; OH* is the dominant intermediate. The OH* concentration increased radially from the particle surface to the center (Fig. 10a), consistent with the production of numerous radicals such as CH₃˙ through heterogeneous initiation as a byproduct of heterogeneous Eley–Rideal steps. Despite the reactants in the primary heterogeneous initiation, namely CH₄ and O*, not showing mass transport limitations, the consumption of the CH₃˙ radical was noticeably influenced not only by its decreasing concentration but also by the intraparticle profile of OH*. This reasoning extends to other radicals primarily originating from the surface, such as C₂H₅˙. Another significant surface intermediate at z = 0 was O*, produced directly through O₂ chemisorption. This step, along with all other chemisorption steps in the reaction mechanism, is quasi-equilibrated based on the partial equilibrium indices (in the 0.46–0.5 range). Vacancies constitute a minor fraction at z = 0.

Fig. 10 Axial and radial (intraparticle phase) fractional coverage profile predictions for IMP SiO₂ for (a) OH*, (b) , (c) O*, (d) * (vacancies), (e) CH₃O*, (f) CH₃CHO*, (g) CO*, and (h) H₂O*. Simulation conditions and results are identical to those of Fig. 9.

A decrease in the OH* radical concentration was observed along the axial direction. Although expecting an increased OH* concentration due to heterogeneous initiation steps is reasonable, heightened recombination rates reduce these surface species. Subsequently, as reactions progressed, the catalyst surface became increasingly rich in (Fig. 10b). does not form any intraparticle concentration gradient, leading to an anticipated plateau in CH₄ conversion as the concentration becomes predominant; this is in accord with studies on CO₂ inhibition in various catalyst families, including Li/MgO,⁷³ La₂O₃,⁷⁴ and Mn–Na₂WO₄/SiO₂.⁷⁵ Unlike surface CO₂, CO* shows intraparticle profiles (Fig. 10g), possibly because of the low concentration of CO*. Notably, in the case of IMP SiO₂, no significant concentrations of oxygenates result from direct sticking reactions. The catalyst hence has high C₂ selectivity.

The primary difference between SD SiO₂–α + βSiC (ESI,† Fig. S8) and IMP SiO₂ lies in the faster heterogeneous oxidation pathway of CH₃˙ radicals, increasing CH₃O* and concentrations. Despite the lower CH₄ conversion and C₂ selectivity predicted for SD SiO₂–α + βSiC compared with IMP SiO₂, which are reflected in reduced OH* concentrations, an increase in coverages significantly affects CH₄ conversion. This is driven by increased CH₃O* concentrations, generating pronounced radial intraparticle gradients. In SD SiO₂–α + βSiC, concentrations of CH₃O*, CH₂O*, and CHO* are comparable, around 10⁻⁵, while concentrations of CO* and CH₃CHO* are around 10⁻¹⁰, indicating the susceptibility of CH₃˙ radicals rather than C₂H₄ to depletion.

The intraparticle species distribution of SD SiO₂–βSiC (ESI,† Fig. S9) exhibited behavior intermediate between the behavior of IMP SiO₂ and that of SD SiO₂–α + βSiC, with CH₃O* presence but a lower amount than that in SD SiO₂–α + βSiC, and the concentration was similar to that in IMP SiO₂. This could explain the CH₄ conversion of SD SiO₂–βSiC being intermediate between IMP SiO₂ and SD SiO₂–α + βSiC. Despite its high predicted C₂ selectivity, SD SiO₂–βSiC showed greater propensity to chemisorb C₂H₄.

3.7 Consumption analysis of steady-state simulations

The interpretation of concentration profiles can be further improved by using a quantitative approach involving reaction path analysis, such as the consumption analysis of Gupta and Vlachos,⁷⁶ which was used in this work for the interstitial and intraparticle phases separately; averaged concentrations along radial coordinates at the end of the catalyst bed were used in the analysis. Specifically, C-containing species were targeted, with a default equilibrium tolerance and a reaction rate cutoff value of 10⁻¹². Fig. 11 shows the interstitial phase species consumption analysis for the three catalysts.

Fig. 11 Species consumption analysis in the interstitial phase at the end of the catalyst bed (z = L_b) based on consumption rates of each carbon-containing species for (a) IMP SiO₂, (b) SD SiO₂–α + βSiC, and (c) SD SiO₂–βSiC. Simulation conditions: T = 800 °C, P = 1 atm, feed molar ratio of CH₄/O₂/He = 3/1/0.6, W/F_CH4,0 = 4.4 g_c h mol_C⁻¹. Numbers denote the molar consumption rate (in percentage) of the reactant for each specific product. The arrow width is proportional to the reaction rates at the simulation conditions.

In this phase, CH₄ is primarily consumed through reactions with H˙ and OH˙, which exclusively generate CH₃˙ radicals. The former reaction accounts for approximately 51%, 41%, and 50% of interstitial CH₃˙ production for IMP SiO₂, SD SiO₂–α + βSiC, and SD SiO₂–βSiC, respectively, while the latter contributes around 30%, 32%, and 31% for these catalysts, respectively. This is because, under OCM conditions without a catalyst, the reaction between CH₄ and O₂ is not the most important initiation step. CH₃˙ production in the interstitial phase is much lower than in the intraparticle phase, which is confirmed by the CH₃˙ profile in Fig. 9. The produced CH₃˙ undergoes various reactions in the interstitial phase, mainly CH₃˙ coupling to yield C₂H₆ (around 92% of CH₃˙ consumption rate), apart from reactions involving CH₃˙-mediated H abstraction from hydrocarbons such as CH₂O, C₂H₆, and C₂H₄.

Less relevant pathways depleting CH₃˙ radicals yield undesired products, including reactions with O₂ resulting in CH₂O and OH˙ production. CH₂O mainly originates from C₂H₃˙ radicals formed from C₂H₄ in IMP SiO₂ (52% of CH₂O compared with 29% and 40% for SD SiO₂–α + βSiC and SD SiO₂–βSiC, respectively). For SD SiO₂–α + βSiC, CH₂O in the interstitial phase mainly results from CH₃˙ reacting with O₂, which accounts for 59% of CH₂O. However, this observation is not evident from Fig. 11 as contributions are expressed in terms of species consumption. The differences observed arise from higher C₂H₃˙ concentrations in IMP SiO₂ and SD SiO₂–βSiC compared with SD SiO₂-α + βSiC owing to higher C₂H₄ concentrations, which lead to higher intrinsic hydrocarbon H-abstraction rates, including both CH₄ and C₂H₄. Additionally, concentration effects explain the lower CH₃˙ recombination rates in SD SiO₂–α + βSiC compared with the other two catalysts, as evidenced by scale differences in the legends across catalysts in Fig. 11. C₂H₄ primarily originates from pyrolytic dehydrogenation of C₂H₅˙ radical (accounting for >95%) in all three catalysts. Similarly, the primary pathway generating interstitial CO involves the pyrolytic route, accounting for more than 98% of CO-producing reactions. Notably, in the interstitial phase, CO can be consumed only via reaction with HO₂˙ to yield CO₂. Variations observed in homogeneous rates in the interstitial phase are primarily associated with the catalyst's capability to generate and consume radicals in the intraparticle phase.

In the intraparticle phase (Fig. 12), CH₄ conversion mainly occurs via heterogeneous Eley–Rideal H-abstraction, and this conversion accounts for 74%, 77%, and 77% of the CH₃˙ radical produced in IMP SiO₂, SiO₂–α + βSiC, and SD SiO₂–βSiC, respectively. Gas-phase contribution within catalyst pores involves CH₄ consumption via H˙ and OH˙ reactions that yield CH₃˙ radicals, and these reactions account for 25%, 22%, and 22% of the total intraparticle CH₄ consumption for IMP SiO₂, SiO₂–α + βSiC, and SD SiO₂–βSiC, respectively. The intraparticle phase accounts for 91.9%, 92.7%, and 93.5% of the total CH₄ net consumption rates at the end of the catalyst bed of 1.6 × 10⁻³ kmol m_b⁻³ s⁻¹, 1.3 × 10⁻³ kmol m_b⁻³ s⁻¹ and 2.1 × 10⁻³ kmol m_b⁻³ s⁻¹ for IMP SiO₂, SiO₂–α + βSiC, and SD SiO₂–βSiC, respectively, highlighting the crucial role of the catalyst and differences in radical generation between the catalysts. CH₃˙ radicals in the intraparticle phase mainly recombine to form C₂H₆ within catalyst pores (Fig. 12), with 65%, 64%, and 67% conversion rates for IMP SiO₂, SiO₂–α + βSiC, and SD SiO₂–βSiC, respectively. Although CH₃˙ radicals favor recombination in the interstitial phase, competition with oxidation routes intensifies in the intraparticle phase.

Fig. 12 Species consumption analysis for the intraparticle phase at the end of the catalyst bed (z = L_b) based on consumption rates of each carbon-containing species for (a) IMP SiO₂, (b) SD SiO₂–α + βSiC, and (c) SD SiO₂–βSiC. Simulation conditions: T = 800 °C, P = 1 atm, feed molar ratio of CH₄/O₂/He = 3/1/0.6, W/F_CH4,0 = 4.4 g_c h mol_C⁻¹. Numbers denote the molar consumption rate (in percentage) of the reactant for each specific product. The arrow width is proportional to the reaction rates at the simulation conditions.

In the intraparticle phase, catalytic oxidation pathways surpass homogeneous ones. Uncoupled CH₃˙ predominantly chemisorb CH₃O* on the catalyst surface, representing 31%, 32%, and 29% of the consumption, respectively. Gas-phase oxidation of CH₃˙ through CH₃O˙–CH₂O–CHO˙ or direct CH₂O/CHO˙ formation is limited, as most CH₃˙ is converted homogeneously to C₂H₆via H abstraction, yielding C₂H₅˙ and CH₄. Conversely, C₂H₅˙ is mainly converted to C₂H₄ (>97% for all catalysts) in the gas phase, indicating that catalyst properties should balance homogeneous-heterogeneous interactions effectively.

As noted, the origins of CO and CO₂ vary with the catalyst. IMP SiO₂ (Fig. 12a) mainly originates from CH₃˙ direct sticking, following the CH₃O*–CH₂O*–CHO*–CO*– sequence, with minimal involvement of C₂H₄ oxidation. The favored pathway from C₂H₄ to C₂H₃˙ involves gas-phase oxidation, primarily through O₂ (28%), catalytic H-abstraction by O* (24%), or reaction with H˙ (21%). Furthermore, gas-phase chain growth to C₃H₇˙ consumes 11% of C₂H₄. CH₃O* follows the expected catalytic path to CO, with 70% formed catalytically via CO desorption from CHO* through the Langmuir–Hinshelwood route.²⁸ Eley–Rideal H abstraction of CHO˙ is a minor contributor, and the remaining CO comes from gas-phase CHO˙ collision. Around half (50.3%) of CO* is desorbed, while 49.7% is oxidized to CO₂, almost entirely catalytically. Catalysts play a key role in radical quenching. HO₂˙ consumption occurs solely through heterogeneous collision with vacancies, and H-abstraction of HO₂˙ by O* partly explains the decreasing OH* concentration. Recombination consumes surface hydroxyl radicals at a rate of 1.4 × 10⁻³ kmol m_b⁻³ s⁻¹.

Regarding the origin of CO and CO₂, SD SiO₂–α + βSiC (Fig. 12b) resembles IMP SiO₂, where CH₃˙ sticking predominantly drives heterogeneous deep oxidation to surpass gas-phase CH₂O generation. Notably, in SD SiO₂–α + βSiC, up to 87% of CO is produced heterogeneously. Gas-phase CH₂O mostly reconverts CH₄ (51%) through CH₃˙ reaction, while the remaining CH₂O is converted via heterogeneous Eley–Rideal H abstraction by O* into CHO˙. In fact, from a mechanistic viewpoint, SD SiO₂–α + βSiC is comparable to SD SiO₂–βSiC (Fig. 12c). SD SiO₂–βSiC significantly influences catalytic oxidation, which accounts for 88% of CO and nearly 100% of CO₂. Furthermore, surface oxidation intermediates such as CH₂O* and CHO* originate from C₂H₄ catalytic oxidation, which accounts for 13% and 12% of the two intermediates, respectively. In particular, for all three catalysts, the HO₂˙ quenching rates were high (2.9 × 10⁻⁵ kmol m_b⁻³ s⁻¹, 4.73 × 10⁻⁵ kmol m_b⁻³ s⁻¹, and 1.6 × 10⁻⁴ kmol m_b⁻³ s⁻¹), which substantially reduced the HO₂˙ concentration (compared with ref. 55, 61 and 71) in the interstitial phase as well as the gas-phase oxidation contribution.

4. Conclusions

The present work delved into the microkinetics of OCM by using three Mn–Na₂WO₄ catalysts (IMP SiO₂, SD SiO₂–α + βSiC, and SD SiO₂–βSiC) to elucidate the effect of SiC in the catalyst support. The microkinetic analysis involved a reactor model that combined irreducible mass transfer limitations and the homogeneous–heterogeneous kinetics of OCM, while maintaining thermodynamic consistency across the kinetic parameters. In the model, textural properties, especially surface area and catalyst porosity, were pivotal factors influencing reaction outcomes. The surface area directly influenced CH₄ conversion rates and undesirable side reactions. In contrast, catalyst porosity facilitated enhanced diffusion of reactants and products, alleviating internal transport constraints while providing a confined environment for critical gas-phase reactions, such as CH₃˙ radicals coupling. Thus, incorporating SiC as a support component endows the Mn–Na₂WO₄ catalysts with increased surface area and porosity, facilitating a more balanced interplay between textural properties and kinetics.

On the kinetic side, the influence of SiC, particularly its crystal phase (βSiC vs. α + βSiC), hinges on the origin of CO_X products. While all three catalysts exhibited noteworthy selectivity toward C₂ products, SD SiO₂–βSiC showed a higher propensity for C₂H₄ chemisorption and its subsequent oxidation, unlike the other catalysts, which showed a greater inclination toward CH₃˙ oxidation, leading to the formation of CH₂O* and CHO* intermediates. SD SiO₂–α + βSiC promoted the oxidation of CH₃˙ radicals, resulting in an increased concentration of CH₃O* and an abundance of on the catalyst surface. Furthermore, the excessive accumulation of on the catalyst surface, particularly from CH₃O*, adversely affected the global reaction rate. Consequently, the catalytic contribution of CO was higher for the SiC-containing catalysts (87–88%) compared to the IMP SiO₂ catalyst (70%). This underscores the notion that in the OCM process, radicals should not only be generated at high rates but also consumed judiciously, with HO₂˙ rates contributing to high C₂ selectivities.

Overall, while SiC enhances leads to more favorable textural properties for OCM, these improvements come at the cost of a smaller active site density, leading to milder catalytic contributions. This study serves as a foundational step in the microkinetic modeling of SiC-containing catalysts for Mn–Na₂WO₄-based OCM, providing valuable insights into the catalytic role of SiC and highlighting potential areas for further model refinement, such as incorporating a more rigorous temperature-dependent characterization of chemisorption entropies, implementation of more comprehensive homogeneous models, conducting tailored experiments to propose a model that accounts for a more realistic interaction between surface and reactive oxygen and incorporating heat effects into the reactor model.

Nomenclature

Abbreviations

BEP	Brønsted–Evans–Polanyi
BET	Brunauer–Emmett–Teller
CH4-TPSR	Temperature-programmed surface reaction
CSTR	Continuously-stirred tank reactor
DAEs	Differential-algebraic equations
DOFs	Degrees of freedom
O2-TPD	Temperature-programmed O2 desorption
OCM	Oxidative coupling of methane
ODEs	Ordinary differential equations
PDEs	Partial differential-algebraic equations
PFR	Plug-flow reactor
PSSA	Pseudo-steady-state approximation

Symbols

A j	Prefactor of homogeneous or heterogeneous reaction step j, c.u.
A r	Reactor cross-section, mb^2
c d,i	Coefficient of linearly decomposed reaction d onto reaction i belonging to the basis set
c j,i	Coefficient for the adsorption of surface species i in heterogeneous reaction step j
C g,i	Concentration of the gas-phase species i in the interstitial phase, kmol mg^−3
C i	Concentration of gas-phase species i, kmol mg^−3
C M	Concentration of an unspecified collision partner in three-body reaction step j, kmol mg^−3
C s,i	Concentration of the gas-phase species i in the intraparticle phase, kmol mg^−3
D e,i	Effective diffusivity of gas-phase species i in the mixture, mg^3 mc^−1 s^−1
D i,j	Binary molecular diffusivity of gas-phase species i in gas-phase species j, mg^2 s^−1
D m,i	Molecular diffusivity of gas-phase species i in the mixture, mg^2 s^−1
E 0,f	Intrinsic energy barrier for any reaction in reaction family f, kJ mol^−1
E a,j	Activation energy of homogeneous or heterogeneous reaction step j, kJ mol^−1
f	Model multiresponse function
F E	Fisher's E, unitless
F i	Carbon molar flow rate of species i, molC s^−1
F inv	Inverse F distribution, unitless
F v	Total gas volumetric flow rate, mg^3 s^−1
F T	Total molar flow rate of species i, NmL min^−1
	Standard Gibbs energy of species i, J mol^−1
	Standard enthalpy of species i, J mol^−1
K C,j	Concentration-based thermodynamic equilibrium constant of homogeneous or heterogeneous reaction step j, c.u.
k j	Rate constant of homogeneous or heterogeneous reaction step j, c.u.
K P,j	Pressure-based thermodynamic equilibrium constant of homogeneous or heterogeneous reaction step j, unitless
L	Bed length, m or cm
M w,i	Molecular weight of gas-phase species i, kg kmol^−1
n	Number of experiments and responses
n j	Number of sites involved in the heterogeneous elementary reaction step j
N obs	Number of observations
N r,d	Total number of linearly dependent reaction steps in the surface mechanism
N r,g	Total number of homogeneous elementary steps in the mechanism
N r,i	Total number of linearly independent reaction steps in the surface mechanism
N r,s	Total number of heterogeneous elementary steps in the mechanism
N res	Number of experimental responses per observation
N s,g	Total number of gas-phase species
N s,j	Total number of species in homogeneous or heterogeneous reaction step j
N s,g,j	Total number of gaseous species in homogeneous or heterogeneous reaction step j
N s,s	Total number of surface intermediates
N s,s,j	Total number of surface intermediates in heterogeneous reaction step j
OF	Objective function, unitless
P	Pressure, bar
p	Total number of catalytic descriptors of the model
p i	Partial pressure of species i, atm
R	Universal gas constant, 8314 Pa mg^3 kmol^−1 K^−1 or 8.314 × 10^−3 kJ mol^−1 K^−1
r	Radial interstitial coordinate, dimensionless
R g,i	Homogeneous net production rate of gas-phase species i, kmol mg^−3 s^−1
r j	Rate of homogeneous or heterogeneous elementary reaction step j, kmol mg^−3 s^−1 or kmol mc^−2 s^−1
r p	Radius of the intraparticle phase or the average particle radius, mc
R s,i	Heterogeneous net production rate of gas-phase or surface species i,
ℜ d	dth dependent reaction outside the basis set
ℜ i	ith independent reaction inside the basis set
r v	Radius of the interstitial phase or average half distance between catalyst particles, mg
	Standard entropy of species i, J mol^−1 K^−1
s j,i	Sticking coefficient of adsorbing gas-phase species i in heterogeneous reaction step j, unitless
S s	Catalyst specific surface area, mc^2 kgc^−1
t	Time, s
t ∞	Holdup time, s
T	Temperature, K
W	Catalyst mass, mgc or gc
x i	Variable representing the ith observation, unitless
X i	Molar fraction of gas-phase species i in the mixture
z	Axial reactor bed coordinate, mb

Greek symbols

α	Significance of the statistical test, unitless
α f	Transfer coefficient for any reaction in reaction family f, unitless
β	β∈[Doublestruck R]^p vector of catalytic descriptor estimates of the model, c.u., or heating rate of the catalyst bed under temperature-programmed experiments, °C min^−1
β i	Constant accounting for the temperature dependence of the chemisorption entropy, unitless
	Standard reaction Gibbs free energy of homogeneous or heterogeneous reaction step j, kJ mol^−1
	Standard reaction enthalpy of homogeneous or heterogeneous reaction step j, kJ mol^−1
ΔQavg	Average bond energy difference between two types of hydrocarbons, kJ mol^−1
	Standard reaction entropy of homogeneous or heterogeneous reaction step j, J mol^−1 K^−1
ε b	Average bed packing porosity, mg^3 mb^−3
ε s	Catalyst porosity, mg^3 mc^−3
θ *	Fractional coverage of the vacant sites or vacancies
θ i	Fractional coverage of surface intermediate i
λ i	Diffusion length of gas-phase species i, m, or ith eigenvalue, dimensionless
ν j,i	Stoichiometric number of gas-phase species i in homogeneous or heterogeneous reaction step j
ξ	Radial intraparticle coordinate, dimensionless
ρ b	Average bed packing density, kgc mb^−3
ρ s	Catalyst density, kgc mc^−3
σ	Active site density, kmol mc^−2
τ i	Lifetime of gas-phase species i, s
τ s	Term combining constriction and tortuosity of the catalyst, mg^2 mc^−2
ν j,i	Stoichiometric number of gas-phase species i reacting in homogeneous or heterogeneous reaction step j (i.e., νj,i wherein all non-negative values 0)
ϕ ij	ϕ ij ∈ [Doublestruck R]^Nres first-order normalized sensitivity of the descriptor j at the ith experimental conditions, unitless
ψ i	ψ i ∈ [Doublestruck R]^Nres vector of all performance metrics for the ith observation, %
ψ i,j	Performance metric of species j for the ith observation, %

Subscripts and superscripts

0	Feed
ads	Chemisorption step
b	Backward (superscript) or reactor bed (subscript)
C	Carbon
c	Catalyst
calc	Model-based prediction
exp	Experimental
F	Final
f	Forward
g	Gas
gas	Gas-phase analogous reaction
I	Initial
init	Initial estimate
min	Minimum
s	Species
sur	Surface reaction step

Data availability

The data supporting this article have been included as part of the ESI.†

Author contributions

Gontzal Lezcano: investigation; data curation; software; formal analysis; visualization; writing – original draft. Shekhar R. Kulkarni: validation; formal analysis; writing – original draft. Vijay K. Velisoju: investigation; validation; formal analysis; writing – review & editing. Natalia Realpe: formal analysis; writing – review & editing. Pedro Castaño: conceptualization; methodology; funding acquisition; supervision; project administration; writing – review and editing.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors acknowledge the financial support, resources, and facilities provided by the King Abdullah University of Science and Technology (KAUST), BAS/1/1403.

References

1. Y. Gao, L. Neal, D. Ding, W. Wu, C. Baroi, A. M. Gaffney and F. Li, ACS Catal., 2019, 9, 8592–8621.
2. H. Saito and Y. Sekine, RSC Adv., 2020, 10, 21427–21453.
3. D. Fairuzov, I. Gerzeliev, A. Maximov and E. Naranov, Catalysts, 2021, 11, 833.
4. J. Q. Chen, A. Bozzano, B. Glover, T. Fuglerud and S. Kvisle, Catal. Today, 2005, 106, 103–107.
5. Y. Gambo, A. A. Jalil, S. Triwahyono and A. A. Abdulrasheed, J. Ind. Eng. Chem., 2018, 59, 218–229.
6. D. Kiani, S. Sourav, J. Baltrusaitis and I. E. Wachs, ACS Catal., 2019, 9, 5912–5928.
7. C. A. Ortiz-Bravo, C. A. Chagas and F. S. Toniolo, J. Nat. Gas Sci. Eng., 2021, 96, 104254.
8. N. S. Hayek, N. S. Lucas, C. Warwar Damouny and O. M. Gazit, ACS Appl. Mater. Interfaces, 2017, 9, 40404–40411.
9. N. S. Hayek, G. J. Khlief, F. Horani and O. M. Gazit, J. Catal., 2019, 376, 25–31.
10. A. Aseem, G. G. Jeba, M. T. Conato, J. D. Rimer and M. P. Harold, Chem. Eng. J., 2018, 331, 132–143.
11. T. Chukeaw, S. Sringam, M. Chareonpanich and A. Seubsai, Mol. Catal., 2019, 470, 40–47.
12. D. Matras, A. Vamvakeros, S. D. M. Jacques, N. Grosjean, B. Rollins, S. Poulston, G. B. G. Stenning, H. R. Godini, J. Drnec, R. J. Cernik and A. M. Beale, Faraday Discuss., 2021, 229, 176–196.
13. J. Ramos-Yataco and J. Notestein, Catal. Today, 2023, 416, 113770.
14. N. Hiyoshi and T. Ikeda, Fuel Process. Technol., 2015, 133, 29–34.
15. P. Wang, G. Zhao, Y. Wang and Y. Lu, Sci. Adv., 2017, 3, e1603180.
16. M. Yildiz, U. Simon, T. Otremba, Y. Aksu, K. Kailasam, A. Thomas, R. Schomäcker and S. Arndt, Catal. Today, 2014, 228, 5–14.
17. G. Lezcano, S. R. Kulkarni, V. K. Velisoju, V. E. Musteata, I. Hita, A. Ramirez, A. Dikhtiarenko, J. Gascon and P. Castaño, Mol. Catal., 2022, 527, 112399.
18. G. Lezcano, V. K. Velisoju, S. R. Kulkarni, A. Ramirez and P. Castaño, Ind. Eng. Chem. Res., 2021, 60, 18770–18780.
19. S. R. Kulkarni, G. Lezcano, V. K. Velisoju, N. Realpe and P. Castaño, ChemCatChem DOI:10.1002/cctc.202301720.
20. G. Lezcano, A. Gobouri, N. Realpe, S. R. Kulkarni, V. K. Velisoju and P. Castaño, Chem. Eng. Sci., 2024, 283, 119412.
21. P. M. Couwenberg, Q. Chen and G. B. Marin, Ind. Eng. Chem. Res., 1996, 35, 3999–4011.
22. R. J. Kee, M. E. Coltrin, P. Glarborg and H. Zhu, in Chemically Reacting Flow, John Wiley & Sons, Inc., Hoboken, NJ, USA, 2005, pp. 487–539.
23. H. Wang, C. Shao, J. Gascon, K. Takanabe and S. M. Sarathy, ACS Omega, 2021, 6, 33757–33768.
24. Q. Chen, P. M. Couwenberg and G. B. Marin, Catal. Today, 1994, 21, 309–319.
25. A. Ishikawa and Y. Tateyama, ACS Catal., 2021, 11, 2691–2700.
26. V. I. Alexiadis, M. Chaar, A. van Veen, M. Muhler, J. W. Thybaut and G. B. Marin, Appl. Catal., B, 2016, 199, 252–259.
27. Y. Gordienko, T. Usmanov, V. Bychkov, V. Lomonosov, Z. Fattakhova, Y. Tulenin, D. Shashkin and M. Sinev, Catal. Today, 2016, 278, 127–134.
28. Z. Xiong, J. Guo, Y. Deng, M. Zeng, Z. Wang, Z. Zhou, W. Yuan and F. Qi, Catal. Sci. Technol., 2024, 6882–6892.
29. M. Y. Sinev, Z. T. Fattakhova, V. I. Lomonosov and Y. A. Gordienko, J. Nat. Gas Chem., 2009, 18, 273–287.
30. J. W. Thybaut, J. Sun, L. Olivier, A. C. Van Veen, C. Mirodatos and G. B. Marin, Catal. Today, 2011, 159, 29–36.
31. J. Sun, J. W. Thybaut and G. B. Marin, Catal. Today, 2008, 137, 90–102.
32. V. I. Alexiadis, J. W. Thybaut, P. N. Kechagiopoulos, M. Chaar, A. C. Van Veen, M. Muhler and G. B. Marin, Appl. Catal., B, 2014, 150–151, 496–505.
33. A. H. Motagamwala and J. A. Dumesic, Chem. Rev., 2021, 121, 1049–1076.
34. O. V. Krylov, Catal. Today, 1993, 18, 209–302.
35. J. A. Dumesic, D. F. Rudd, L. M. Aparicio, J. E. Rekoske and A. A. Treviño, in The Microkinetics of Heterogeneous Catalysis, American Chemical Society (ACS), 1st edn, 1993.
36. N. Schumacher, A. Boisen, S. Dahl, A. Gokhale, S. Kandoi, L. Grabow, J. Dumesic, M. Mavrikakis and I. Chorkendorff, J. Catal., 2005, 229, 265–275.
37. A. Michaelides, Z.-P. Liu, C. J. Zhang, A. Alavi, D. A. King and P. Hu, J. Am. Chem. Soc., 2003, 125, 3704–3705.
38. A. B. Mhadeshwar, H. Wang and D. G. Vlachos, J. Phys. Chem. B, 2003, 107, 12721–12733.
39. A. B. Mhadeshwar and D. G. Vlachos, Combust. Flame, 2005, 142, 289–298.
40. Y. S. Su, J. Y. Ying and W. H. Green, J. Catal., 2003, 218, 321–333.
41. Y.-R. Luo, Handbook of Bond Dissociation Energies in Organic Compounds, CRC Press, 2002.
42. P. C. St. John, Y. Guan, Y. Kim, S. Kim and R. S. Paton, Nat. Commun., 2020, 11, 2328.
43. A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker and C. S. Woodward, ACM Trans. Math. Softw., 2005, 31, 363–396.
44. B. Malengier, P. Kišon, J. Tocknell, C. Abert, F. Bruckner and M.-A. Bisotti, J. Open Source Softw., 2018, 3, 165.
45. D. G. Goodwin, H. K. Moffat, I. Schoegl, R. L. Speth and B. W. Weber, Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes, Version 3.0.0., 2023, https://www.cantera.org,  DOI:10.5281/zenodo.8137090.
46. G. P. Smith, D. M. Golden, M. Frenklach, N. W. Moriarty, B. Eiteneer, M. Goldenberg, C. T. Bowman, R. K. Hanson, S. Song, Jr., W. C. Gardiner, V. V. Lissianski and Z. Qin, GRI-Mech 3.0, http://combustion.berkeley.edu/gri-mech/version30/text30.html, (accessed 16 December 2021).
47. F.-A. Fortin, F.-M. De Rainville, M.-A. Gardner, M. Parizeau and C. Gagné, Journal of Machine Learning Research, 2012, 13, 2171–2175.
48. P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, A. Vijaykumar, A. Pietro Bardelli, A. Rothberg, A. Hilboll, A. Kloeckner, A. Scopatz, A. Lee, A. Rokem, C. N. Woods, C. Fulton, C. Masson, C. Häggström, C. Fitzgerald, D. A. Nicholson, D. R. Hagen, D. V. Pasechnik, E. Olivetti, E. Martin, E. Wieser, F. Silva, F. Lenders, F. Wilhelm, G. Young, G. A. Price, G.-L. Ingold, G. E. Allen, G. R. Lee, H. Audren, I. Probst, J. P. Dietrich, J. Silterra, J. T. Webber, J. Slavič, J. Nothman, J. Buchner, J. Kulick, J. L. Schönberger, J. V. de Miranda Cardoso, J. Reimer, J. Harrington, J. L. C. Rodríguez, J. Nunez-Iglesias, J. Kuczynski, K. Tritz, M. Thoma, M. Newville, M. Kümmerer, M. Bolingbroke, M. Tartre, M. Pak, N. J. Smith, N. Nowaczyk, N. Shebanov, O. Pavlyk, P. A. Brodtkorb, P. Lee, R. T. McGibbon, R. Feldbauer, S. Lewis, S. Tygier, S. Sievert, S. Vigna, S. Peterson, S. More, T. Pudlik, T. Oshima, T. J. Pingel, T. P. Robitaille, T. Spura, T. R. Jones, T. Cera, T. Leslie, T. Zito, T. Krauss, U. Upadhyay, Y. O. Halchenko and Y. Vázquez-Baeza, Nat. Methods, 2020, 17, 261–272.
49. R. Tesser and V. Russo, Advanced Reactor Modeling with MATLAB, De Gruyter, 2020.
50. P. A. Brodtkorb and J. D'Errico, numdifftools 0.9.11, https://github.com/pbrod/numdifftools.
51. K. Toch, J. W. Thybaut and G. B. Marin, AIChE J., 2015, 61, 880–892.
52. J. Sadeghzadeh Ahari, S. Zarrinpashne and M. T. Sadeghi, Fuel Process. Technol., 2013, 115, 79–87.
53. H. Wang, R. Schmack, B. Paul, M. Albrecht, S. Sokolov, S. Rümmler, E. V. Kondratenko and R. Kraehnert, Appl. Catal., A, 2017, 537, 33–39.
54. C. Karakaya, H. Zhu, C. Loebick, J. G. Weissman and R. J. Kee, Catal. Today, 2018, 312, 10–22.
55. P. N. Kechagiopoulos, J. W. Thybaut and G. B. Marin, Ind. Eng. Chem. Res., 2014, 53, 1825–1840.
56. J. Liu, M. Huo and Y. Zhu, J. Nanosci. Nanotechnol., 2017, 17, 8818–8826.
57. J. Li, J. Chen, A. Zanina, Y. Li, C. Yu, M. Liu, G. Cui, Y. Wang, M. Zhou, E. V. Kondratenko and G. Jiang, J. Catal., 2023, 428, 115176.
58. V. I. Lomonosov, Y. A. Gordienko, M. Y. Sinev, V. A. Rogov and V. A. Sadykov, Russ. J. Phys. Chem. A, 2018, 92, 430–437.
59. V. Fleischer, R. Steuer, S. Parishan and R. Schomäcker, J. Catal., 2016, 341, 91–103.
60. W. Sun, Y. Gao, G. Zhao, J. Si, Y. Liu and Y. Lu, J. Catal., 2021, 400, 372–386.
61. V. I. Alexiadis, T. Serres, G. B. Marin, C. Mirodatos, J. W. Thybaut and Y. Schuurman, AIChE J., 2018, 64, 2603–2611.
62. M. Schwaab, L. P. Lemos and J. C. Pinto, Chem. Eng. Sci., 2008, 63, 2895–2906.
63. M. Schwaab and J. C. Pinto, Chem. Eng. Sci., 2007, 62, 2750–2764.
64. L. Pirro, P. S. F. Mendes, B. D. Vandegehuchte, G. B. Marin and J. W. Thybaut, React. Chem. Eng., 2020, 5, 584–596.
65. J. Chawla, S. Schardt, P. Lott, S. Angeli, S. Tischer, L. Maier and O. Deutschmann, Chem. Eng. J., 2024, 482, 148719.
66. L. Pirro, A. Obradović, B. D. Vandegehuchte, G. B. Marin and J. W. Thybaut, Ind. Eng. Chem. Res., 2018, 57, 16295–16307.
67. R. Petrov, S. Reshetnikov and Y. Ivanova, Fuel Process. Technol., 2021, 213, 106667.
68. Y. Cheng, P. S. F. Mendes, P. Yazdani and J. W. Thybaut, Int. J. Chem. Kinet., 2024, 703–717.
69. L. Pirro, P. S. F. Mendes, S. Paret, B. D. Vandegehuchte, G. B. Marin and J. W. Thybaut, Catal. Sci. Technol., 2019, 9, 3109–3125.
70. M. A. Vannice, S. H. Hyun, B. Kalpakci and W. C. Liauh, J. Catal., 1979, 56, 358–362.
71. P. M. Couwenberg, Q. Chen and G. B. Marin, Ind. Eng. Chem. Res., 1996, 35, 415–421.
72. H. Zhang, J. Wu, B. Xu and C. Hu, Catal. Lett., 2006, 106, 161–165.
73. S. Al-Zahrani, Q. Song and L. L. Lobban, Ind. Eng. Chem. Res., 1994, 33, 251–258.
74. S. Lacombe, Z. Durjanova, L. Mleczko and C. Mirodatos, Chem. Eng. Technol., 1995, 18, 216–223.
75. M. Ghiasi, A. Malekzadeh, S. Hoseini, Y. Mortazavi, A. Khodadadi and A. Talebizadeh, J. Nat. Gas Chem., 2011, 20, 428–434.
76. U. Gupta and D. G. Vlachos, SoftwareX, 2020, 11, 100442.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4re00403e

---