Mechanistic Kinetic Modelling Framework for the Conversion of Waste Crude Glycerol to Value-Added Hydrogen-Rich Gas

Abstract

The kinetics for crude glycerol autothermal reforming was studied over S/C ratio of 2.6 and O₂/C ratio of 0.125 using 5% Ni/CeZrCa catalyst. Both power law and mechanistic kinetic models were studied. The overall power law model for crude glycerol autothermal reforming was investigated with a pre-exponential factor of 4.3 × 10¹⁰ mol/g_cat·min and activation energy of 8.78 × 10⁴ J/mol. The reaction orders with respect to crude glycerol, water and oxygen are 1.04, 0.54 and 1.78 respectively. The power law model presented an absolute average deviation of 5.84%, which showed a good correlation between the predicted and experimental rate. Mechanistic models were developed for crude glycerol autothermal reforming. For steam reforming, the Eley–Rideal approach best described the reaction rate with the surface reaction being the rate-determining step (AAD < 10%). The kinetics of the total oxidation reaction was best described by the power law model with an AAD of less than 1%, whereas for the TOR process, the molecular adsorption of crude glycerol with an AAD of 14.6% via Langmuir Hinshelwood Hougen-Watson approach was best. CO₂ methanation resulted in an AAD of 5.8% for the adsorption of carbon dioxide (CO₂) by the Eley–Rideal mechanism.

1. Introduction

The production of net-zero fuels has taken on intense focus in recent years as the world works to produce substitutes for hydrocarbon fuels, where their combustion emissions contribute to ongoing climate change. While there are a variety of fuels to consider, they are generally either zero-emission fuels (as in the case of hydrogen, which produces water vapour when combusted), or are net-zero emission fuels as may be the cases in biodiesel production (where the inputs may themselves be derived from plant-based wastes and feedstocks). Production of biodiesel in particular is an interesting net-zero fuel, since its application areas (e.g., heavy-duty construction and agricultural equipment, remote location baseload and on-demand power generation, etc.) are challenging to decarbonize through either electrification or hydrogen fuel vectors. Furthermore, upcoming legislative changes in Canada will see greater emphasis placed on low-emission fuels, where biodiesel may be used as a blending feedstock for petroleum products retailers to meet these requirements and continue serving their customers.

The production of biodiesel from plant and animal oils is relatively straightforward, in that triglycerides present in biological oils undergo transesterification in the presence of an alcohol like methanol and a catalyst to produce biodiesel and glycerol (C₃H₈O₃), which must be separated from one another in order to produce marketable biodiesel [1]. With biodiesel production expected to total 46 billion liters per year from 2023 to 2025 [2], the production of crude glycerol will total an estimated 4.6 billion liters through the same period (presuming a 10% co-production with biodiesel [3,4,5]). This presents a challenge: while there is a growing need for biodiesel (therein presenting a market opportunity for producers), by-product glycerol has limited marketable uses even at current production volumes.

The bottleneck presented by co-production of glycerol has prompted much work on the conversion and use of glycerol in value-added applications [6,7,8]. Investigations have explored the conversion of glycerol to different acids, alcohols and diols, ethers, acrolein, glycerol carbonate, and other commercial chemicals (see

Table 1. Glycerol products and their production routes [6,12,13,14,15,16].

Glycerol Products	Reaction Process	Glycerol Products	Reaction Process
Hydrogen	Oxyhydrocarbon Reforming	1,2-propanediol	Hydrogenolysis or Fermentation
Glycerol carbonate	Carboxylation or Glycerolysis	1,3-propanediol	Hydrogenolysis
Polyglycerol	Polymerization	Ethylene glycol	Hydrogenolysis
Acrolein	Dehydration	Epichlorohydrin	Halogenation
Acrylic Acid	Dehydration	Polyhydroxybutyrate (PHB)	Fermentation
1,3-Dihydroxyacetone (DHA)	Oxidation or Fermentation	Propanol	Fermentation
Glyceraldehyde	Oxidation	Butanol	Fermentation
Glyceric acid	Oxidation	Butanediol	Fermentation
Glycolic acid	Oxidation	Lactic Acid	Fermentation
Hydroxypyruvic acid	Oxidation	Propionic acid	Fermentation
Mesoxalic acid	Oxidation	Citric acid	Fermentation
Oxalic acid	Oxidation	Succinic acid	Fermentation
Tartronic acid	Oxidation	Glyceric Acid	Fermentation
Hydroxyethanoic acid	Oxidation	Mannitol	Fermentation
Formic acid	Oxidation	Erythritol	Fermentation
Glycerol Ethers	Etherification	Arabitol	Fermentation
Hydrogen	Oxyhydrocarbon Reforming	1,2-propanediol	Hydrogenolysis or Fermentation

, below). In addition, investigators have also explored the reforming of glycerol into synthesis gas [9,10,11]—a compelling option if this reforming can be completed autothermally and produced hydrogen can either be captured for fuel use or incorporated into existing hydrogen supply networks (e.g., in the case where a biodiesel facility is integrated into a petrochemical complex). A further advantage may be realized if glycerol can be reformed on a ‘feed-flexible’ basis—i.e., where crude glycerol (subject to minimal treatment) can be reacted directly to produce synthesis gas, and crude glycerol composition can be free to fluctuate as reaction conditions and feed sources change.

Crude glycerol reforming is not as straightforward as what might be assumed for the reforming of pure glycerol; following separation of crude biodiesel, crude glycerol may contain other impurities like residual methanol, soaps (from saponification side reactions [17]), water, and residual free fatty acids (FFAs) and fatty acid methyl esters (FAME) not recovered to the biodiesel stream [18]. In this case, the reforming of crude glycerol involves not just the reforming of glycerol itself, but potentially that of its contaminants and intermediate products as well.

The reforming of glycerol itself has been investigated in the literature. Sundari and Vaidya [19] studied the kinetics of pure glycerol steam reforming over Ru/Al₂O₃ catalyst. In their work, glycerol decomposition and water–gas shift were the primary reactions. In deriving the mechanistic equation, the process description included reversible glycerol adsorption, Eley–Rideal reaction of adsorbed glycerol with gas-phase steam, decomposition of the adsorbed compex to intermediates, and ultimately to gas-phase carbon dioxide and hydrogen. Their studies were conducted at a temperature of 350–500 °C and a steam-to-glycerol ratio of 3 to 12. Cheng et al. [20] studied the kinetics of glycol steam reforming using a Co-Ni/Al₂O₃ catalyst. The study parameters ranged between a steam-to-carbon ratio of 3:1 to 12:1 and a temperature of 773 to 823 K, all at a pressure of 1 atm. Power law, single- and dual-site Eley–Rideal (ER) and Langmuir–Hinshelwood–Hougen–Watson (LHHW) models were derived and parameterized, with the LHHW model describing molecular adsorption of glycerol and steam onto the catalyst surface with rate-limiting surface reaction prevailing as the most accurate model. Authayanun et al. [21] studied the thermodynamics of hydrogen production for fuel cells via the autothermal reforming of crude glycerol, where the composition of crude glycerol was dependent on the raw material used for upstream biodiesel production. The parameters analysed were operating temperature, oxygen-to-feed (crude glycerol) ratio, and steam-to-feed ratio. The crude glycerol was modeled as a mixture of pure glycerol and methanol. The product composition in their studies was solely based on thermodynamic analysis, leaving out the effects of a catalyst. They concluded that the most suitable oxygen-to-carbon ratio resided in the range of 0.4–0.7, depending on the purity of the crude glycerol, and also found that the oxygen requirement (to maintain autothermal conditions) increased as the feed flow of steam increased and the ratio of methanol to glycerol decreased. Liu and Lawal [22] developed a kinetic rate expression for the autothermal reforming of glycerol over a monolith catalyst, based on a LHHW mechanism. The standard operating conditions chosen were: 650 °C, 1 atm, steam-to-carbon ratio of 0.4, and oxygen-to-carbon ratio of 0.15. The atomic matrix approach was used to determine the number of independent reactions which could be used to describe the overall autothermal reforming process. For a catalyst weight above 0.012 g, it was found that all of the oxygen fed to the reactor was consumed—prompting only the steam reforming and water–gas shift reactions to be modelled. Both single-site and dual-site LHHW mechanisms were proposed, with the main difference between models being the assumption of whether adsorption was dissociative or non-dissociative. It was also noted, given the removal of water from glycerol, that the adsorption of OH on the active site was likely.

Reaction schemes proposed in literature are noted below in

Table 2. Proposed reaction schemes for glycerol/glycol in the literature.

Sundari and Vaidya [19]	Liu and Lawal [22]
${\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3\rightleftharpoons 4{\mathrm{H}}_2+3\mathrm{CO}$	${\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3+3{\mathrm{H}}_2\mathrm{O}\rightleftharpoons 7{\mathrm{H}}_2+3{\mathrm{CO}}_2$
$\mathrm{CO}+{\mathrm{H}}_2\mathrm{O}\rightleftharpoons {\mathrm{CO}}_2+{\mathrm{H}}_2$	${\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3+2{\mathrm{O}}_2\rightleftharpoons 4{\mathrm{H}}_2\mathrm{O}+3\mathrm{CO}$
${\mathrm{CH}}_4+2{\mathrm{H}}_2\mathrm{O}\rightleftharpoons 4{\mathrm{H}}_2+{\mathrm{CO}}_2$	$\mathrm{CO}+{\mathrm{H}}_2\mathrm{O}\rightleftharpoons {\mathrm{H}}_2+{\mathrm{CO}}_2$
Cheng et al. [20]	Authayanun et al. [21]
${\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3+3{\mathrm{H}}_2\mathrm{O}\rightleftharpoons 7{\mathrm{H}}_2+3{\mathrm{CO}}_2$	${\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3+3{\mathrm{H}}_2\mathrm{O}\rightleftharpoons 7{\mathrm{H}}_2+3{\mathrm{CO}}_2$
${\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3\rightleftharpoons 4{\mathrm{H}}_2+3\mathrm{CO}$	${\mathrm{CH}}_3\mathrm{OH}+{\mathrm{H}}_2\mathrm{O}\rightleftharpoons 3{\mathrm{H}}_2+{\mathrm{CO}}_2$
$\mathrm{CO}+{\mathrm{H}}_2\mathrm{O}\rightleftharpoons {\mathrm{CO}}_2+{\mathrm{H}}_2$	${\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3+\frac{3}{2}{\mathrm{O}}_2\rightleftharpoons 4{\mathrm{H}}_2+3{\mathrm{CO}}_2$
$2{ \mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3+{ \mathrm{H}}_2\rightleftharpoons 3{ \mathrm{CH}}_4+3{ \mathrm{H}}_2\mathrm{O}+3 \mathrm{CO}$	${\mathrm{CH}}_3\mathrm{OH}+\frac{1}{2}{\mathrm{O}}_2\rightleftharpoons 2{\mathrm{H}}_2+{\mathrm{CO}}_2$
${\mathrm{CH}}_4+{\mathrm{H}}_2\mathrm{O}\rightleftharpoons \mathrm{CO}+3{\mathrm{H}}_2$	$\mathrm{CO}+{\mathrm{H}}_2\mathrm{O}\rightleftharpoons {\mathrm{H}}_2+{\mathrm{CO}}_2$
${\mathrm{CO}}_2+{\mathrm{CH}}_4\rightleftharpoons 2\mathrm{CO}+2{\mathrm{H}}_2$	$\mathrm{CO}+3{\mathrm{H}}_2\rightleftharpoons {\mathrm{H}}_2\mathrm{O}+{\mathrm{CH}}_4$
$\mathrm{C}+{\mathrm{H}}_2\mathrm{O}\rightleftharpoons {\mathrm{H}}_2+\mathrm{CO}$	${\mathrm{CO}}_2+{\mathrm{CH}}_4\rightleftharpoons 2{\mathrm{H}}_2+2\mathrm{CO}$
${\mathrm{CO}}_2+4{\mathrm{H}}_2\rightleftharpoons 3{\mathrm{CH}}_4+2{\mathrm{H}}_2\mathrm{O}$
$2\mathrm{CO}\rightleftharpoons \mathrm{C}+{\mathrm{CO}}_2$
${\mathrm{CH}}_4\rightleftharpoons 2{\mathrm{H}}_2+\mathrm{C}$

.

Kinetic studies were completed by Ghani et al. [18] in our laboratories for the autothermal reforming of synthetic crude glycerol over a 5% Ni/CeZrCa catalyst with operating ranges of 773–873 K, 0.05–0.20 O/C ratio, and 1.6–3.6 S/C ratio. Power law modelling showed a first-order dependence on crude glycerol, a 0.5-order dependence on steam, and a second-order dependence on oxygen partial pressure.

To our present awareness, however, there does not seem to be available in literature a detailed mechanistic analysis of crude glycerol reforming kinetic data. As such, the objective of this present study (having access to the raw kinetic data produced by Ghani et al. [9]) is to postulate an appropriate set of mechanistic rate expressions that reflect the kinetics of crude glycerol reforming. In addition to this, use of the mechanistic rate expressions will be featured as part of reactor simulations to conduct a parametric study of reactor performance using this catalyst.

2. Preceding Experimental and Analytical Work

2.1. Catalyst Development and Kinetic Testing

The experimental data used were taken from the preceding study by Ghani et al. [18], which will be briefly reviewed here.

An initial catalyst screening study was conducted by the same investigators [23], testing the effect of gadolinium, magnesium, and calcium promoters on ceria–zirconia catalyst support. Surfactant-assisted support precipitation using cetyltrimethylammonium bromide (CTAB) was used, where nitrate salts of catalyst support precursors were mixed with CTAB in a 0.5:1 surfactant-to-metal molar ratio. This solution was then pH adjusted using ammonium hydroxide and aged for five days at 90 °C, where the precipitate was then filtered, washed, and dried at 120 °C for 12 h. The catalyst support was then calcined at 650 °C for 3 h, and next underwent incipient wetness impregnation to achieve a metal loading of 5% nickel (using nickel nitrate solution) prior to being calcined again at the same conditions and duration. For the prepared catalysts, a series of reactor tests were completed to perform a parametric analysis of various catalyst preparation factors including support promoter, reduction temperature, calcination temperature, steam-to-carbon ratio, and oxygen-to-carbon ratio. These parameters were evaluated for the autothermal reforming of crude glycerol, which was prepared based on the compositional analysis of industrially sourced crude glycerol received from a Canadian biodiesel producer (containing 45.6% free glycerol, 11.2% methanol, 8.3% water, 29.1% potassium palmitate soap, 3.8% oleic acid, and 1.0% sodium chloride). This screening study determined the most preferable catalyst for autothermal crude glycerol reforming to be calcium-promoted nickel–ceria–zirconia. Greater details are provided in Ghani et al. [23].

Catalyst activity testing and kinetic data collection were conducted within an electrically heated packed bed tubular reactor (12.7 mm internal diameter by ~45 mm reaction zone length), illustrated in Ghani et al. [18]. Ni/CeZrCa Catalyst (diluted with alumina) loaded into the reactor was first heated to 600 °C under pure nitrogen atmosphere, and then (once stable temperature was achieved) reduced under 100 mL/h of 5% hydrogen (in nitrogen) gas flow for two hours. Once complete, inlet gas flow was switched to compressed air supply, and the crude glycerol syringe pump (loaded with crude glycerol) was set at a rate of 0.15 mL/min [24]. Air flow was adjusted to meet the required oxygen-to-carbon ratio, while the crude glycerol water was mixed with water to achieve the required steam-to-carbon ratio. Temperature was monitored using a K-type thermocouple placed within the reactor. Reactor effluents were cooled using a water-cooled condenser, with remnant gases directed to gas chromatography for analysis. Lastly, integral reactor test data was taken to determine a power law model that could replicate the experimentally observed conversion of crude glycerol.

2.2. Preceding Reaction Rate Power Law Modelling

Ghani et al. [18] developed an empirical formula for crude glycerol (${\mathrm{C}}_{2.5}{\mathrm{H}}_7{\mathrm{O}}_2$), with the overall stoichiometric equation for autothermal reforming being:

$$\begin{array}{cc}{\mathrm{C}}_{2.5}{\mathrm{H}}_7{\mathrm{O}}_2+0.7{\mathrm{H}}_2\mathrm{O}+0.65{\mathrm{O}}_2\to 2{\mathrm{CO}}_2+0.5{\mathrm{CH}}_4+3.2{\mathrm{H}}_2\hfill & \hfill \Delta {\mathrm{H}}_{298\mathrm{K}}=-72.39\mathrm{kJ}/\mathrm{mol}\end{array}$$

(1)

According to Ghani et al. [18], a typical autothermal reaction is composed of steam reforming and partial oxidation as well as methanation, methane dry reforming, glycerol decomposition, coke formation, and other associated reactions. Some of the reactions involved are as follows:

$$\begin{array}{ccc}\begin{array}{c}\mathrm{Glycerol} \mathrm{Steam}\\ \mathrm{Reforming}\end{array}& {\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3+3{\mathrm{H}}_2\mathrm{O}\rightleftharpoons 7{\mathrm{H}}_2+3{\mathrm{CO}}_2& \Delta {\mathrm{H}}_{298\mathrm{K}}=+127.7\mathrm{kJ}/\mathrm{mol}\end{array}$$

(2)

$$\begin{array}{ccc}\begin{array}{c}\mathrm{Glycerol} \mathrm{Partial}\\ \mathrm{Oxidation}\end{array}& {\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3+\frac{3}{2}{\mathrm{O}}_2\rightleftharpoons 4{\mathrm{H}}_2+3{\mathrm{CO}}_2& \Delta {\mathrm{H}}_{298\mathrm{K}}=-597.7\mathrm{kJ}/\mathrm{mol}\end{array}$$

(3)

$$\begin{array}{ccc}\begin{array}{c}\mathrm{Glycerol} \mathrm{Total}\\ \mathrm{Oxidation}\end{array}& {\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3+2{\mathrm{O}}_2\rightleftharpoons 4{\mathrm{H}}_2\mathrm{O}+3\mathrm{CO}& \Delta {\mathrm{H}}_{298\mathrm{K}}=-715.9\mathrm{kJ}/\mathrm{mol}\end{array}$$

(4)

$$\begin{array}{ccc}\mathrm{CO} \mathrm{Methanation}& \mathrm{CO}+3{\mathrm{H}}_2\rightleftharpoons {\mathrm{H}}_2\mathrm{O}+{\mathrm{CH}}_4& \Delta {\mathrm{H}}_{298\mathrm{K}}=-206.2\mathrm{kJ}/\mathrm{mol}\end{array}$$

(5)

$$\begin{array}{ccc}\begin{array}{c}\mathrm{Methane}\\ \mathrm{Dry} \mathrm{Reforming}\end{array}& {\mathrm{CO}}_2+{\mathrm{CH}}_4\rightleftharpoons 2{\mathrm{H}}_2+2\mathrm{CO}& \Delta {\mathrm{H}}_{298\mathrm{K}}=+247.4\mathrm{kJ}/\mathrm{mol}\end{array}$$

(6)

$$\begin{array}{ccc}\mathrm{Water}–\mathrm{Gas} \mathrm{Shift}& \mathrm{CO}+{\mathrm{H}}_2\mathrm{O}\rightleftharpoons {\mathrm{H}}_2+{\mathrm{CO}}_2& \Delta {\mathrm{H}}_{298\mathrm{K}}=-41.8\mathrm{kJ}/\mathrm{mol}\end{array}$$

(7)

With the wide variety of possible reactions involved in the autothermal reforming of hydrocarbons and oxy-hydrocarbons, it is prudent to know which exact reactions define the system. As noted by Liu and Lawal [22], a matrix approach is used to obtain the rank of the matrix—indicating the number of independent reactions describing the system being studied. The matrix (shown in

Table 3. Atomic matrix for crude glycerol reforming.

Atom	Compound
	${\mathbf{H}}_2$	${\mathbf{CO}}_2$	${\mathbf{C}}_3{\mathbf{H}}_8{\mathbf{O}}_3$	${\mathbf{CH}}_4$	${\mathbf{H}}_2\mathbf{O}$	${\mathbf{O}}_2$
$\mathrm{C}$	-	1	3	1	-	-
H	2	-	8	4	2	-
O	-	2	3	-	1	2

) is constructed by identifying the atomic species C, H, and O present in each of the possible chemical species consumed or formed during the autothermal reforming of glycerol. For a total of six (6) possible species, and finding this constructed matrix to be rank three (3), taking the difference reveals the number of independent reactions—totalling three (3)—that may be used to represent this system.

Based on the final stoichiometric equation proposed by Ghani et al. [9], and based on GC analysis, it could be inferred that the three (3) main reactions are:

• Steam reforming.
• Total oxidation (where the formation of carbon monoxide was not observed).
• Carbon dioxide methanation (where the formation of methane was observed).

The power law rate expression is written in terms of reactant and/or production concentrations. It can consider the reversibility and irreversibility of the reaction. Assuming the overall reaction is elementary, the power law rate expression can be given in terms of the stoichiometric coefficients observed in Equation (1).

For the reversible power law expression, both the forward and reverse reactions are considered as follows:

$$r=k_{0,\mathrm{f}}\exp \left(-\frac{E_{\mathrm{f}}}{RT}\right)P_{\mathrm{CG}}^m{P}_{{\mathrm{H}}_2\mathrm{O}}^n{P}_{{\mathrm{O}}_2}^o-k_{0,\mathrm{r}}\exp \left(-\frac{E_{\mathrm{r}}}{RT}\right)P_{{\mathrm{CO}}_2}^q{P}_{{\mathrm{CH}}_4}^r{P}_{{\mathrm{H}}_2}^s$$

(8)

For the irreversible case, only the forward reaction is considered:

$$r=k_{0,\mathrm{f}}\exp \left(-\frac{E_{\mathrm{f}}}{RT}\right)P_{\mathrm{CG}}^m{P}_{{\mathrm{H}}_2\mathrm{O}}^n{P}_{{\mathrm{O}}_2}^o$$

(9)

2.3. Preceding Reactor Modelling in COMSOL

A straightforward reactor modelling study using the power law developed by Ghani et al. [18] was completed by Afabor et al. [25]. A transient two-dimensional axisymmetric pseudo-homogeneous reactor model was developed within COMSOL Multiphysics to replicate the experimental reactor. Comparison between experimental and calculated results indicated an AAE of 8.9%. The effect of including vs. excluding the axial dispersion term in the reactor model was determined to be negligible, while parametric simulations of crude glycerol conversion versus space velocity and reaction temperature.

As part of their conclusions, Afabor et al. [25] recommended that the next step in their line of investigation be to consider the development of a more generalized model. In the work described here, a mechanistic vs. power law type reaction rate expression is developed to help expand this generality. Furthermore, one-dimensional vs. two-dimensional reactor modelling frameworks, both respectively considering a pseudo-homogeneous and simplified heterogeneous modelling approach, are developed in MATLAB and COMSOL, respectively, and the simulation results of these reactor models are compared. Parametric study of the reactor using the mechanistic model and the reactor models developed is also provided.

3. Kinetic Analysis

3.1. Verification of Absence of Transport Limitations

In their preceding work, Ghani et al. [18] conducted qualitative experimental testing to confirm the absence of transport limitations that may lead to misinterpretation of kinetic data. This took the form of comparisons between pellet size and crude glycerol conversion, as well as reactant flow rate and crude glycerol conversion. Observing no correlations in either case, it was determined that the selected pellet size (0.8 mm) and tested flow rate conditions would not present any transport limitations.

However, greater confidence in this assertion was desired. Transport limitations may be tested quantitatively following defined criteria to evaluate the presence of transport limitations.

Internal heat transfer limitations may be assessed via the following expression [26]:

$$\mathsf{\Delta }{T}_{\max ,\mathrm{p}}=\frac{D_{\mathrm{eff}}\left(C_{A_{\mathrm{s}}}-C_A\right)\left(\Delta H_{\mathrm{r}}\right)}{{\lambda }_{\mathrm{eff},\mathrm{p}}}<1$$

(10)

In this expression, $\mathsf{\Delta }{T}_{\max ,\mathrm{p}}$ is the maximum observed thermal gradient within the catalyst pellets, $D_{\mathrm{eff}}$ is the species effective diffusivity (which may be calculated from the Bosanquet approximation between molecular diffusivity and Knudsen diffusivity [27]), $C_{A_{\mathrm{s}}}$ and $C_A$ are the species concentrations at the catalyst surface and core, $\Delta H_{\mathrm{r}}$ is the heat of reaction, and ${\lambda }_{\mathrm{eff},\mathrm{p}}$ is the catalyst pellet effective thermal conductivity. Evaluation of the kinetic data indicated $\mathsf{\Delta }{T}_{\max ,\mathrm{p}}$ was found to be 0.15 K, indicating no significant thermal gradients within the catalyst pellets.

External heat transfer limitations may be assessed by [26]:

$$\mathsf{\Delta }{T}_{\max ,\mathrm{f}}=\frac{L_{\mathrm{c}}\left(r_{\mathrm{obs}}\right)\left(\Delta H_{\mathrm{r}}\right)}{h_{\mathrm{p},\mathrm{f}}}<1$$

(11)

In this expression, $\mathsf{\Delta }{T}_{\max ,\mathrm{f}}$ is the maximum observed thermal gradient existing across the catalyst thin film (i.e., gas-phase boundary layer at the catalyst solid surface), $L_{\mathrm{c}}$ is the characteristic length (which, for a spherical geometry, will be one-third of the catalyst radius), $r_{\mathrm{obs}}$ is the observed reaction rate, and $h_{\mathrm{p},\mathrm{f}}$ is the convective heat transfer coefficient calculated for the catalyst thin film. Evaluation of the kinetic data indicated $\mathsf{\Delta }{T}_{\max ,\mathrm{f}}$ was found to be 0.6 K, indicating no significant thermal gradients existed between the catalyst pellets and the bulk gas.

A more rigorous criterion developed by Mears [28,29] was also used to evaluate the presence of thermal gradients in the catalyst thin film. The criterion is given as:

$$\frac{r_{\mathrm{obs}}{\rho }_{\mathrm{b}}{R}_{\mathrm{p}}{E}_{\mathrm{a}}\Delta H_{\mathrm{r}}}{h_{\mathrm{p},\mathrm{f}}{T}^{2}R}<0.15$$

(12)

In this expression, ${\rho }_{\mathrm{b}}$ is the catalyst bed bulk density, $R_{\mathrm{p}}$ is the catalyst pellet radius, $E_{\mathrm{a}}$ is the activation energy, $T$ is the reaction temperature, and $R$ is the universal gas constant. Calculation of this ratio reveals a value of 0.02, again indicating that thin film heat transport limitations did not exist.

Catalyst internal mass transfer limitations may be assessed using the Weisz–Prater criterion, expressed as [27,30]:

$$\frac{-r_{\mathrm{obs}}{\rho }_{\mathrm{p}}{R}_{\mathrm{p}}^2}{D_{\mathrm{eff}}{C}_{A_{\mathrm{s}}}}<0.3$$

(13)

Calculation of the Weisz–Prater criterion revealed a value of 0.24, satisfying the condition implying no presence of significant intraparticle mass transfer limitations.

Catalyst thin film mass transfer limitations may be assessed by the following expression [26]:

$$\mathrm{observed} \mathrm{reaction} \mathrm{rate}/\mathrm{film} \mathrm{mass} \mathrm{transfer} \mathrm{rate}=\frac{-r_{\mathrm{obs}}}{C_{A_{\mathrm{s}}}{k}_{\mathrm{p},\mathrm{f}}}\left(\frac{d_p}{6}\right)$$

(14)

Here, $k_{\mathrm{p},\mathrm{f}}$ represents the catalyst thin film mass transfer coefficient. This ratio was calculated to be 0.0069, indicating that there was no significant limitation caused by thin film mass transfer.

Again, the more rigorous criterion given by Mears [28] for mass transfer limitations in the catalyst thin film was also considered. The criterion is defined as:

$$\frac{-r_{\mathrm{obs}}{\rho }_{\mathrm{b}}{R}_{\mathrm{p}}n}{C_{A_{\mathrm{s}}}{k}_{\mathrm{p},\mathrm{f}}}<0.15$$

(15)

Here, $n$ represents the order of the reaction. This ratio was calculated to be 0.072, well below the limit of 0.1—indicating further that catalyst thin film mass transfer limitations were not present.

3.2. Mechanistic Model Development

Following the approach of Liu and Lawal [22], the autothermal reforming of crude glycerol was subdivided into three (3) independent reactions, where the kinetics were studied for each reaction. Reaction mechanisms following the Eley–Rideal (ER), and Langmuir–Hinshelwood–Hougen–Watson (LHHW) were developed for each of these reactions. For the sake of simplicity, crude glycerol was written out as pure glycerol in both the derived mechanistic rate expression as well as for the reversible and irreversible power law rate expressions. The reactions considered were as follows:

$$\begin{array}{cc}\begin{array}{c}\mathrm{Overall} \mathrm{Autothermal}\\ \mathrm{Reforming} \mathrm{Reaction}\end{array}& {\mathrm{C}}_{2.5}{\mathrm{H}}_7{\mathrm{O}}_2+0.7{\mathrm{H}}_2\mathrm{O}+0.65{\mathrm{O}}_2\to 2{\mathrm{CO}}_2+0.5{\mathrm{CH}}_4+3.2{\mathrm{H}}_2\end{array}$$

(16)

$$\begin{array}{cc}\mathrm{Glycerol} \mathrm{Steam} \mathrm{Reforming}& 0.8{\mathrm{C}}_{2.5}{\mathrm{H}}_7{\mathrm{O}}_2+2.4{\mathrm{H}}_2\mathrm{O}\to 2{\mathrm{CO}}_2+5.2{\mathrm{H}}_2\end{array}$$

(17)

$$\begin{array}{cc}\mathrm{Glycerol} \mathrm{Oxidation}& 0.2{\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3+0.65{\mathrm{O}}_2\to 0.5{\mathrm{CO}}_2+0.7{\mathrm{H}}_2\mathrm{O}\end{array}$$

(18)

$$\begin{array}{cc}{\mathrm{C}\mathrm{O}}_2 \mathrm{Methanation}& 0.5{\mathrm{CO}}_2+2{\mathrm{H}}_2\rightleftarrows 0.5{\mathrm{CH}}_4+{\mathrm{H}}_2\mathrm{O}\end{array}$$

(19)

3.2.1. Glycerol Steam Reforming (SR)

The steam reforming reaction occurs when steam reacts with crude glycerol to form carbon dioxide and hydrogen gas. The stoichiometric equation for this reaction is given in Equation (17). According to Sundari et al. [19] and Cheng et al. [20], a plausible mechanism for glycerol steam reforming often includes dehydrogenation and scission of C-C and C-O bonds as well as dehydration and hydration reaction steps.

In the Eley–Rideal approach, it was assumed that steam reforming did not produce methane as one of the products, and that the mechanism followed a reactant adsorption, surface reaction, and product desorption process. Vacant active sites were also accounted for in the site balance. The hypothesized mechanism is described as follows:

$$\begin{array}{cc}\begin{array}{c}\mathrm{Molecular} \mathrm{adsorption} \mathrm{of} \mathrm{crude} \mathrm{glycerol} \mathrm{onto} \mathrm{an} \mathrm{active} \mathrm{site}\end{array}& 0.8{\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3+\mathrm{S}\rightleftarrows 0.8{\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3.\mathrm{S}\end{array}$$

(20)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Dissociation} \mathrm{of} \mathrm{adsorbed} \mathrm{intermediate} \mathrm{onto} a \mathrm{vacant} \mathrm{site} \mathrm{to} \mathrm{form} \mathrm{two} \mathrm{(2)} \mathrm{different} \mathrm{oxygenated} \mathrm{hydrocarbon} \mathrm{compounds} \mathrm{(ethylene} \mathrm{glycol} \mathrm{and} \mathrm{formaldehyde)}\end{array}& {\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3.\mathrm{S}\rightleftarrows {\mathrm{C}}_2{\mathrm{H}}_6{\mathrm{O}}_2.\mathrm{S}+{\mathrm{CH}}_2\mathrm{O}.\mathrm{S}\end{array}$$

(21)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Surface} \mathrm{reaction} \mathrm{between} \mathrm{adsorbed} \mathrm{ethylene} \mathrm{glycol} \mathrm{and} \mathrm{non-adsorbed} \mathrm{steam} \mathrm{in} \mathrm{the} \mathrm{bulk} \mathrm{gas}\end{array}& {\mathrm{C}}_2{\mathrm{H}}_6{\mathrm{O}}_2.\mathrm{S}+{\mathrm{H}}_2\mathrm{O}\rightleftarrows 1.5{\mathrm{CO}}_2+4.2{\mathrm{H}}_2+\mathrm{S}\end{array}$$

(22)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Surface} \mathrm{reaction} \mathrm{between} \mathrm{adsorbed} \mathrm{formaldehyde} \mathrm{and} \mathrm{non-adsorbed} \mathrm{steam} \mathrm{in} \mathrm{the} \mathrm{bulk} \mathrm{gas}\end{array}& {\mathrm{CH}}_2\mathrm{O}.\mathrm{S}+0.2{\mathrm{H}}_2\mathrm{O}\rightleftarrows 0.5{\mathrm{CO}}_2+{\mathrm{H}}_2+\mathrm{S}\end{array}$$

(23)

In the LHHW approach, the assumption is made that both crude glycerol and steam adsorb to the catalyst active site. The approach of Cheng et al. [20] was modified to consider both molecular and atomic adsorption of the reactants in order to gain a better understanding of the mechanism. The mechanism is as follows:

$$\begin{array}{cc}\begin{array}{c}\mathrm{Molecular} \mathrm{adsorption} \mathrm{of} \mathrm{crude} \mathrm{glycerol} \mathrm{onto} \mathrm{an} \mathrm{active} \mathrm{site}\end{array}& 0.8{\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3+\mathrm{S}\rightleftarrows 0.8{\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3.\mathrm{S}\end{array}$$

(24)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Molecular} \mathrm{adsorption} \mathrm{of} \mathrm{steam} \mathrm{onto} \mathrm{an} \mathrm{active} \mathrm{site}\end{array}& 2.4{ \mathrm{H}}_2\mathrm{O}+\mathrm{S}\rightleftarrows 2.4{ \mathrm{H}}_2\mathrm{O}.\mathrm{S}\end{array}$$

(25)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Surface} \mathrm{reaction} \mathrm{between} \mathrm{adsorbed} \mathrm{glycerol} \mathrm{and} \mathrm{water}\end{array}& 0.8{\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3.\mathrm{S}+2.4{ \mathrm{H}}_2\mathrm{O}.\mathrm{S}\rightleftarrows {\mathrm{CH}}_3{\left({\mathrm{CH}}_2\right)}_{\mathrm{x}}\mathrm{COOH}.\mathrm{S}+\mathrm{S}\end{array}$$

(26)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Desorption} \mathrm{of} \mathrm{products} \mathrm{from}\\ \mathrm{the} \mathrm{active} \mathrm{site}\end{array}& {\mathrm{CH}}_3{\left({\mathrm{CH}}_2\right)}_{\mathrm{x}}\mathrm{COOH}.\mathrm{S}\rightleftarrows 2{\mathrm{CO}}_2+5.2{\mathrm{H}}_2+\mathrm{S}\end{array}$$

(27)

The proposed rate expressions for crude glycerol steam reforming are presented in

Table 4. Summary of rate equation assuming each step as a rate determining step assuming both reversible and irreversible reactions.

Equation (r.d.s.)	Rate Law Expression
ER_1 (Equation (20)) (Reversible)	$r_1=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[A\right]}^{0.8}-\frac{{\left[E\right]}^2{\left[F\right]}^{5.2}}{K_{p\left(eq\right)}\left({\left[B\right]}^{2.4}\right)}\right)}{\left(1+\frac{K_F{\left[E\right]}^2{\left[F\right]}^{5.2}}{\left({\left[B\right]}^{2.4}\right)}+\frac{K_B{\left[E\right]}^{1.5}{\left[F\right]}^{4.2}}{\left({\left[B\right]}^{2.2}\right)}+\frac{K_A{\left[E\right]}^{0.5}\left[F\right]}{\left({\left[B\right]}^{0.2}\right)}\right)}\right]$
ER_2 (Equation (21)) (Reversible)	$r_2=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[A\right]}^{0.8}-\frac{{\left[E\right]}^2{\left[F\right]}^{5.2}}{K_{p\left(eq\right)}\left({\left[B\right]}^{2.4}\right)}\right)}{{\left[1+K_A{\left[A\right]}^{0.8}+\frac{K_B\left[F\right]{\left[E\right]}^{0.5}}{\left({\left[B\right]}^{0.2}\right)}+\frac{K_W{\left[E\right]}^{1.5}{\left[F\right]}^{4.2}}{\left({\left[B\right]}^{2.2}\right)}\right]}^2}\right]$
ER_3 (Equation (22)) (Reversible)	$r_3=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left(\frac{{\left[A\right]}^{0.8}{\left[B\right]}^{2.4}}{\left[F\right]{\left[E\right]}^{0.5}}-\frac{{\left[E\right]}^{1.5}{\left[F\right]}^{4.2}}{K_{p\left(eq\right)}}\right)}{\left(1+K_1{\left[A\right]}^{0.8}+\frac{K_W{\left[A\right]}^{0.8}{\left[B\right]}^{0.2}}{\left[F\right]{\left[E\right]}^{0.5}}+\frac{K_B{\left[E\right]}^{0.5}\left[F\right]}{{\left[B\right]}^{0.2}}\right)}\right]$
ER_4 (Equation (23)) (Reversible)	$r_4=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left(\frac{{\left[A\right]}^{0.8}{\left[B\right]}^{2.4}}{{\left[E\right]}^{1.5}{\left[F\right]}^{4.2}}-\frac{\left[F\right]{\left[E\right]}^{0.5}}{K_{p\left(eq\right)}}\right)}{\left(1+K_1{\left[A\right]}^{0.8}+\frac{K_W{\left[A\right]}^{0.8}{\left[B\right]}^{2.2}}{{\left[E\right]}^{1.5}{\left[F\right]}^{4.2}}+\frac{K_C{\left[E\right]}^{1.5}{\left[F\right]}^{4.2}}{{\left[B\right]}^{2.2}}\right)}\right]$
ER_5 (Equation (20)) (Irreversible)	$r_5=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[A\right]}^{0.8}\right)}{\left(1+\frac{K_F{\left[E\right]}^2{\left[F\right]}^{5.2}}{\left({\left[B\right]}^{2.4}\right)}+\frac{K_B{\left[E\right]}^{1.5}{\left[F\right]}^{4.2}}{\left({\left[B\right]}^{2.2}\right)}+\frac{K_A{\left[E\right]}^{0.5}\left[F\right]}{\left({\left[B\right]}^{0.2}\right)}\right)}\right]$
ER_6 (Equation (21)) (Irreversible)	$r_6=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[A\right]}^{0.8}\right)}{{\left[1+K_A{\left[A\right]}^{0.8}+\frac{K_B{\left[E\right]}^{0.5}\left[F\right]}{\left({\left[B\right]}^{0.2}\right)}+\frac{K_W{\left[E\right]}^{1.5}{\left[F\right]}^{4.2}}{\left({\left[B\right]}^{2.2}\right)}\right]}^2}\right]$
ER_7 (Equation (22)) (Irreversible)	$r_7=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left(\frac{{\left[A\right]}^{0.8}{\left[B\right]}^{2.4}}{{\left[E\right]}^{0.5}\left[F\right]}\right)}{\left(1+K_A{\left[A\right]}^{0.8}+\frac{K_W{\left[A\right]}^{0.8}{\left[B\right]}^{0.2}}{{\left[E\right]}^{0.5}\left[F\right]}+\frac{K_B{\left[E\right]}^{0.5}\left[F\right]}{{\left[B\right]}^{0.2}}\right)}\right]$
ER_8 (Equation (23)) (Irreversible)	$r_8=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left(\frac{{\left[A\right]}^{0.8}{\left[B\right]}^{2.4}}{{\left[E\right]}^{1.5}{\left[F\right]}^{4.2}}\right)}{\left(1+K_A{\left[A\right]}^{0.8}+\frac{K_W{\left[A\right]}^{0.8}{\left[B\right]}^{2.2}}{{\left[E\right]}^{1.5}{\left[F\right]}^{4.2}}+\frac{K_B{\left[E\right]}^{1.5}{\left[F\right]}^{4.2}}{{\left[B\right]}^{2.2}}\right)}\right]$
LHHW_1 (Equation (24)) (Reversible)	$r_1=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[A\right]}^{0.8}-\frac{{\left[E\right]}^2{\left[F\right]}^{5.2}}{K_{p\left(eq\right)}\left({\left[B\right]}^{2.4}\right)}\right)}{\left(1+K_2{\left[B\right]}^{2.4}+K_B{\left[E\right]}^2{\left[F\right]}^{5.2}+\frac{K_W{\left[E\right]}^2{\left[F\right]}^{5.2}}{\left({\left[B\right]}^{2.4}\right)}\right)}\right]$
LHHW_2 (Equation (25)) (Reversible)	$r_2=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[B\right]}^{2.4}-\frac{{\left[E\right]}^2{\left[F\right]}^{5.2}}{K_{p\left(eq\right)}\left({\left[A\right]}^{0.8}\right)}\right)}{\left(1+K_1{\left[A\right]}^{0.8}+K_B{\left[E\right]}^2{\left[F\right]}^{5.2}+\frac{K_W{\left[E\right]}^2{\left[F\right]}^{5.2}}{\left({\left[A\right]}^{0.8}\right)}\right)}\right]$
LHHW_3 (Equation (26)) (Reversible)	$r_3=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[A\right]}^{0.8}{\left[B\right]}^{2.4}-\frac{{\left[E\right]}^2{\left[F\right]}^{5.2}}{K_{p\left(eq\right)}}\right)}{{\left(1+K_1{\left[A\right]}^{0.8}+K_B{\left[E\right]}^2{\left[F\right]}^{5.2}+K_2{\left[B\right]}^{2.4}\right)}^2}\right]$
LHHW_4 (Equation (27)) (Reversible)	$r_4=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[A\right]}^{0.8}{\left[B\right]}^{2.4}-\frac{{\left[E\right]}^2{\left[F\right]}^{5.2}}{K_{p\left(eq\right)}}\right)}{\left(1+K_1{\left[A\right]}^{0.8}+K_B{\left[A\right]}^{0.8}{\left[B\right]}^{2.4}+K_2{\left[B\right]}^{2.4}\right)}\right]$

Note: $A={\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3; B={\mathrm{H}}_2\mathrm{O}; C={\mathrm{O}}_2; D={\mathrm{CH}}_4; E={\mathrm{CO}}_2; H={\mathrm{C}}_2{\mathrm{H}}_6{\mathrm{O}}_2; I={\mathrm{CH}}_2\mathrm{O}; G={\mathrm{CH}}_3{\left({\mathrm{CH}}_2\right)}_{\mathrm{x}}\mathrm{COOH}$.

below.

3.2.2. Total Oxidative Reforming (TOR)

Oxidation reactions can be categorized into either total oxidation or partial oxidation. For partial oxidation, the obtained products are CO and H₂ [31,32,33,34], whereas for total oxidation the products formed are CO₂ and H₂O [35,36,37]. In this work, crude glycerol is reacted with oxygen to produce carbon dioxide and water. The reaction is exothermic and is given by Equation (18).

In the Eley–Rideal model, oxygen was assumed to be absorbed onto the catalyst active site with crude glycerol existing in the bulk gas phase. This assumed a single site mechanism as follows:

$$\begin{array}{cc}\begin{array}{c}\mathrm{Molecular} \mathrm{adsorption} \mathrm{of} \mathrm{crude} \mathrm{glycerol} \mathrm{onto} \mathrm{an} \mathrm{active} \mathrm{site}\end{array}& 0.65{\mathrm{O}}_2+\mathrm{S} \rightleftarrows 0.65{\mathrm{O}}_2.\mathrm{S}\end{array}$$

(28)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Surface} \mathrm{reaction} \mathrm{between} \mathrm{adsorbed} \mathrm{oxygen} \mathrm{and} \mathrm{non-adsorbed} \mathrm{crude} \mathrm{glycerol}\end{array}& 0.2{\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3+0.65{ \mathrm{O}}_2.\mathrm{S}\rightleftarrows {\mathrm{CH}}_2{\mathrm{O}}_2.\mathrm{S}\end{array}$$

(29)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Molecular} \mathrm{desorption} \mathrm{of} \mathrm{intermediate}\end{array}& {\mathrm{CH}}_2{\mathrm{O}}_2.\mathrm{S}\rightleftarrows 0.5{\mathrm{CO}}_2+0.7{\mathrm{H}}_2\mathrm{O}+\mathrm{S}\end{array}$$

(30)

A variation of this mechanism, considering dissociative vs. molecular adsorption of oxygen, was also investigated. The mechanism is as follows:

$$\begin{array}{cc}\begin{array}{c}\mathrm{Dissociative} \mathrm{adsorption} \mathrm{of} \mathrm{oxygen} \mathrm{onto} \mathrm{an} \mathrm{active} \mathrm{site}\end{array}& 0.65{\mathrm{O}}_2+\mathrm{S}\rightleftarrows 1.3\mathrm{O}.\mathrm{S}\end{array}$$

(31)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Surface} \mathrm{reaction} \mathrm{of} \mathrm{adsorbed} \mathrm{oxygen} \mathrm{and} \mathrm{non-adsorbed} \mathrm{crude} \mathrm{glycerol}\end{array}& 0.2{\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3+1.3\mathrm{O}.\mathrm{S}\rightleftarrows {\mathrm{CH}}_2\mathrm{O}.\mathrm{S}\end{array}$$

(32)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Surface} \mathrm{reaction} \mathrm{of} \mathrm{adsorbed} \mathrm{intermediate} \mathrm{with} \mathrm{remaining} \mathrm{oxygen}\end{array}& {\mathrm{CH}}_2\mathrm{O}.\mathrm{S}+\mathrm{O}.\mathrm{S}\rightleftarrows {\mathrm{CH}}_2{\mathrm{O}}_2.\mathrm{S}+\mathrm{S}\end{array}$$

(33)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Desorption} \mathrm{of} \mathrm{active} \mathrm{intermediate} \mathrm{from} \mathrm{active} \mathrm{site}\end{array}& {\mathrm{CH}}_2{\mathrm{O}}_2.\mathrm{S}\rightleftarrows 0.5{\mathrm{CO}}_2+0.7{\mathrm{H}}_2\mathrm{O}+\mathrm{S}\end{array}$$

(34)

In the Langmuir–Hinshelwood model, oxygen and crude glycerol are both adsorbed onto the active sites. The mechanism is described below.

$$\begin{array}{cc}\begin{array}{c}\mathrm{Dissociative} \mathrm{adsorption} \mathrm{of} \mathrm{oxygen} \mathrm{onto} \mathrm{an} \mathrm{active} \mathrm{site}\end{array}& 0.65{\mathrm{O}}_2+\mathrm{S}\rightleftarrows 1.3\mathrm{O}.\mathrm{S}\end{array}$$

(35)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Non-dissociative} \mathrm{adsorption} \mathrm{of} \mathrm{crude} \mathrm{glycerol} \mathrm{onto} \mathrm{an} \mathrm{active} \mathrm{site}\end{array}& 0.2{\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3+\mathrm{S}\rightleftarrows 0.2{\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3.\mathrm{S}\end{array}$$

(36)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Surface} \mathrm{reaction} \mathrm{of} \mathrm{adsorbed} \mathrm{crude} \mathrm{glycerol} (\mathrm{formation} \mathrm{of} \mathrm{carbon} \mathrm{dioxide} \mathrm{and} \mathrm{steam})\end{array}& 0.2{\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3.\mathrm{S}+1.3\mathrm{O}.\mathrm{S}\rightleftarrows 0.5{\mathrm{CO}}_2+0.7{\mathrm{H}}_2\mathrm{O}.\mathrm{S}\end{array}$$

(37)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Surface} \mathrm{reaction} \mathrm{of} \mathrm{adsorbed} \mathrm{crude} \mathrm{glycerol} (\mathrm{formation} \mathrm{of} \mathrm{formic} \mathrm{acid})\end{array}& 0.2{\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3.\mathrm{S}+1.3\mathrm{O}.\mathrm{S}\rightleftarrows {\mathrm{CH}}_2{\mathrm{O}}_2.\mathrm{S}+\mathrm{S}\end{array}$$

(38)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Dissociative} \mathrm{desorption} \mathrm{of} \mathrm{formic} \mathrm{acid}\end{array}& {\mathrm{CH}}_2{\mathrm{O}}_2.\mathrm{S}\rightleftarrows 0.5{\mathrm{CO}}_2.\mathrm{S}+0.7{\mathrm{H}}_2\mathrm{O}\end{array}$$

(39)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Molecular} \mathrm{desorption} \mathrm{of} \mathrm{carbon} \mathrm{dioxide}\end{array}& 0.5{\mathrm{CO}}_2.\mathrm{S}\rightleftarrows 0.5{\mathrm{CO}}_2+\mathrm{S}\end{array}$$

(40)

The proposed rate expressions for crude glycerol oxidation are presented in

Table 5. Rate law expression for Eley–Rideal and Langmuir–Hinschelwood models via oxidation reforming route.

Equation (r.d.s.)	Rate Law Expression
ER_1 (Equation (28))	$r_1=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[C\right]}^{0.65}-\frac{{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}}{K_{p\left(eq\right)}{\left[A\right]}^{0.2}}\right)}{\left(1+K_A{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}+\frac{K_B{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}}{{\left[A\right]}^{0.2}}\right)}\right]$
ER_2 (Equation (29))	$r_2=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[A\right]}^{0.2}{\left[B\right]}^{0.65}-\frac{{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}}{K_{p\left(eq\right)}}\right)}{\left(1+K_B{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}+K_{{ }_2}{\left[C\right]}^{0.65}\right)}\right]$
ER_3 (Equation (30))	$r_3=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[A\right]}^{0.2}{\left[B\right]}^{0.65}-\frac{{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}}{K_{p\left(eq\right)}}\right)}{\left(1+K_2{\left[C\right]}^{0.65}+K_A{\left[A\right]}^{0.2}{\left[C\right]}^{0.65}\right)}\right]$
ER_4 (Equation (31))	$r_4=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[C\right]}^{0.65}-\frac{{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}}{K_{p\left(eq\right)}{\left[A\right]}^{0.2}}\right)}{\left(1+K_C{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}+K_B{\left[C\right]}^{0.15}{\left[A\right]}^{0.2}+\frac{K_w{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}}{{\left[A\right]}^{0.2}}\right)}\right]$
ER_5 (Equation (32))	$r_5=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[A\right]}^{0.2}{\left[C\right]}^{0.15}-\frac{{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}}{K_{p\left(eq\right)}{\left[C\right]}^{0.5}}\right)}{\left(1+K_3{\left[C\right]}^{0.65}+K_B{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}+\frac{K_F{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}}{{\left[C\right]}^{0.5}}\right)}\right]$
ER_6 (Equation (33))	$r_6=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[A\right]}^{0.2}{\left[C\right]}^{0.65}-\frac{{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}}{K_{p\left(eq\right)}}\right)}{\left(1+K_3{\left[C\right]}^{0.65}+K_B{\left[C\right]}^{0.15}{\left[A\right]}^{0.2}+K_C{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}\right)}\right]$
ER_7 (Equation (34))	$r_7=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[A\right]}^{0.2}{\left[C\right]}^{0.65}-\frac{{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}}{K_{p\left(eq\right)}}\right)}{\left(1+K_3{\left[C\right]}^{0.65}+K_B{\left[C\right]}^{0.15}{\left[A\right]}^{0.2}+K_C{\left[A\right]}^{0.2}{\left[C\right]}^{0.65}\right)}\right]$
LHHW_1 (Equation (35))	$r_1=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[C\right]}^{0.65}-\frac{{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}}{K_{p\left(eq\right)}{\left[A\right]}^{0.2}}\right)}{\left(1+K_6{\left[E\right]}^{0.5}+K_1{\left[A\right]}^{0.2}+K_B{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}+\frac{K_w{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}}{{\left[A\right]}^{0.2}}\right)}\right]$
LHHW_2 (Equation (36))	$r_2=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[A\right]}^{0.2}-\frac{{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}}{K_{p\left(eq\right)}{\left[C\right]}^{0.65}}\right)}{{\left(1+K_3{\left[C\right]}^{0.65}+K_6{\left[E\right]}^{0.5}+K_B{\left[B\right]}^{0.7}+\frac{K_w{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}}{{\left[C\right]}^{0.65}}\right)}^3}\right]$
LHHW_2a (Equation (37))	$r_2=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[A\right]}^{0.2}-\frac{{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}}{K_{p\left(eq\right)}{\left[C\right]}^{0.65}}\right)}{\left(1+K_3{\left[C\right]}^{0.65}+K_6{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}+K_B{\left[E\right]}^{0.5}+\frac{K_w{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}}{{\left[C\right]}^{0.65}}\right)}\right]$
LHHW_3 (Equation (38))	$r_3=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[A\right]}^{0.2}{\left[C\right]}^{0.65}-\frac{{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}}{K_{p\left(eq\right)}}\right)}{\left(1+K_3{\left[C\right]}^{0.65}+K_1{\left[A\right]}^{0.2}+K_6{\left[E\right]}^{0.5}+K_B{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}\right)}\right]$
LHHW_4 (Equation (39))	$r_4=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[A\right]}^{0.2}{\left[C\right]}^{0.65}-\frac{{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}}{K_{p\left(eq\right)}}\right)}{{\left(1+K_3{\left[C\right]}^{0.65}+K_1{\left[A\right]}^{0.2}+K_6{\left[E\right]}^{0.5}+K_B{\left[E\right]}^{0.5}{\left[B\right]}^{0.7}\right)}^2}\right]$
LHHW_5 (Equation (40))	$r_5=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left(\frac{{\left[A\right]}^{0.2}{\left[C\right]}^{0.65}}{{\left[B\right]}^{0.7}}-\frac{{\left[E\right]}^{0.5}}{K_{p\left(eq\right)}}\right)}{\left(1+K_3{\left[C\right]}^{0.65}+K_1{\left[A\right]}^{0.2}+K_B{\left[A\right]}^{0.2}{\left[C\right]}^{0.65}+K_F\frac{{\left[A\right]}^{0.2}{\left[C\right]}^{0.65}}{{\left[B\right]}^{0.7}}\right)}\right]$
LHHW_6 Huang et al. [38]	$r_6=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left\{K_1{K}_2\left[A\right]\left[B\right]\right\}}{{\left\{1+\left(K_1\left[B\right]\right)+\left(K_2\left[A\right]\right)\right\}}^2}\right]$
LHHW_7 Huang et al. [38]	$r_7=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left\{\left(K_1{K}_2\left[A^{0.2}\right]\left[B^{0.65}\right]\right)-\frac{\left[C^{0.5}\right]\left[D^{0.7}\right]}{K_p}\right\}}{{\left\{1+\left(K_1\left[B^{0.65}\right]\right)+\left(K_2\left[A^{0.2}\right]\right)\right\}}^3}\right]$ Based on Huang et al. (2012) using a reversible reaction and the stoichiometry equation as an elementary reaction

Note: $A={\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3; B={\mathrm{H}}_2\mathrm{O}; C={\mathrm{O}}_2; D={\mathrm{CH}}_4; E={\mathrm{CO}}_2, W={\mathrm{CH}}_2{\mathrm{O}}_2; V={\mathrm{CH}}_2\mathrm{O}; S=\mathrm{active} \mathrm{site}$.

below.

3.2.3. Carbon Dioxide (CO₂) Methanation

Carbon dioxide methanation occurs between carbon dioxide and hydrogen gas to form steam and methane gas—an undesired product in this work. This reaction is exothermic and is given by equation (19). Both power law rate expressions and mechanistic rate expressions were considered in this work. For the power law rate expression, the assumption of an irreversible reaction resulted in a high AAD relative to experimental data. As such, a mechanistic model formulation similar to Koschany et al. [39], where kinetics via PLM and LHHW were studied over a NiAl(O)_x catalyst, was adopted to improve the AAD than just the conventional Arrhenius-style power law expression, where the reversibility of the methanation reaction was accounted for. Attempts were made to fit other models from the literature to the experimental data, but none were successful. Additional information on CO₂ methanation is available from Amaral [40].

Two main groupings of mechanistic Eley–Rideal models were hypothesized. The first is based on carbon dioxide onto the active site:

$$\begin{array}{cc}\begin{array}{c}\mathrm{Adsorption} \mathrm{of} \mathrm{carbon} \mathrm{dioxide} \mathrm{onto} \\ \mathrm{an} \mathrm{active} \mathrm{site}\end{array}& 0.5{\mathrm{CO}}_2+\mathrm{S}\rightleftarrows 0.5{\mathrm{CO}}_2.\mathrm{S}\end{array}$$

(41)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Surface} \mathrm{reaction} \mathrm{between} \mathrm{adsorbed} \mathrm{carbon} \mathrm{dioxide} \mathrm{and} \mathrm{hydrogen} \mathrm{gas}\end{array}& 0.5{\mathrm{CO}}_2.\mathrm{S}+2{\mathrm{H}}_2 \rightleftarrows { \mathrm{CH}}_3\mathrm{OH}.\mathrm{S}\end{array}$$

(42)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Dissociative} \mathrm{desorption} \mathrm{of} \mathrm{methanol} \mathrm{from} \mathrm{the} \mathrm{active} \mathrm{site}\end{array}& { \mathrm{CH}}_3\mathrm{OH}.\mathrm{S}\rightleftarrows 0.5{\mathrm{CH}}_4+{\mathrm{H}}_2\mathrm{O}+\mathrm{S}\end{array}$$

(43)

The second considers hydrogen adsorption onto the active site:

$$\begin{array}{cc}\begin{array}{c}\mathrm{Adsorption} \mathrm{of} \mathrm{hydrogen} \mathrm{onto} \\ \mathrm{an} \mathrm{active} \mathrm{site}\end{array}& {\mathrm{H}}_2+\mathrm{S}\rightleftarrows {\mathrm{H}}_2.\mathrm{S}\end{array}$$

(44)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Surface} \mathrm{reaction} \mathrm{between} \mathrm{adsorbed} \mathrm{carbon} \mathrm{dioxide} \mathrm{and} \mathrm{hydrogen} \mathrm{gas}\end{array}& 0.5{\mathrm{CO}}_2+2{\mathrm{H}}_2.\mathrm{S} \rightleftarrows {\mathrm{CH}}_3\mathrm{OH}.\mathrm{S}\end{array}$$

(45)

$$\begin{array}{cc}\begin{array}{c}\mathrm{Dissociative} \mathrm{desorption} \mathrm{of} \mathrm{methanol} \mathrm{from} \mathrm{the} \mathrm{active} \mathrm{site}\end{array}& { \mathrm{CH}}_3\mathrm{OH}.\mathrm{S} \rightleftarrows 0.5{\mathrm{CH}}_4+{\mathrm{H}}_2\mathrm{O}+\mathrm{S}\end{array}$$

(46)

The proposed rate expressions for carbon dioxide methanation are presented in

Table 6. Summary of rate law expression for CO₂ methanation.

Equation (r.d.s.)	Rate Law Expression
ER_1 (Equation (41))	$r_1=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[E\right]}^{0.5}-\frac{{\left[D\right]}^{0.5}\left[B\right]}{K_{p\left(eq\right)}\left[F\right]}\right)}{\left(1+K_{3\left(ads\right)}{\left[D\right]}^{0.5}\left[B\right]+\frac{K_{1\left(ads\right)}{\left[D\right]}^{0.5}\left[B\right]}{\left({\left[F\right]}^2\right)}\right)}\right]$
ER_2 (Equation (42))	$r_2=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[E\right]}^{0.5}{\left[F\right]}^2-\frac{{\left[D\right]}^{0.5}\left[B\right]}{K_{p2\left(eq\right)}}\right)}{\left(1+K_{1\left(ads\right)}{\left[E\right]}^{0.5}+K_{3\left(ads\right)}{\left[D\right]}^{0.5}\left[B\right]\right)}\right]$
ER_3 (Equation (43))	$r_3=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[E\right]}^{0.5}{\left[F\right]}^2-\frac{{\left[D\right]}^{0.5}\left[B\right]}{K_{p2\left(eq\right)}}\right)}{\left(1+K_{1\left(ads\right)}{\left[E\right]}^{0.5}+K_{3\left(ads\right)}{\left[E\right]}^{0.5}\left[F\right]\right)}\right]$

and

Table 7. Summary of Rate Law Expression for H2 adsorption on active site.

Equation (r.d.s.)	Rate Law Expression
ER_4 (Equation (42))	$r_4=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[F\right]}^2-\frac{{\left[D\right]}^{0.5}\left[B\right]}{K_{p\left(eq\right)} \ast {\left[E\right]}^{0.5}}\right)}{\left(1+K_{3\left(ads\right)}{\left[D\right]}^{0.5}\left[B\right]+\frac{K_{1\left(ads\right)}{\left[D\right]}^{0.5}\left[B\right]}{\left({\left[E\right]}^{0.5}\right)}\right)}\right]$
ER_5 (Equation (43))	$r_5=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[E\right]}^{0.5}{\left[F\right]}^2-\frac{{\left[D\right]}^{0.5}\left[B\right]}{K_{p\left(eq\right)}}\right)}{\left(1+K_{1\left(ads\right)}{\left[F\right]}^2+K_{3\left(ads\right)}{\left[D\right]}^{0.5}\left[B\right]\right)}\right]$
ER_6 (Equation (44))	$r_6=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[E\right]}^{0.5}{\left[F\right]}^2-\frac{{\left[D\right]}^{0.5}\left[B\right]}{K_{p\left(eq\right)}}\right)}{\left(1+K_{1\left(ads\right)}{\left[F\right]}^2+K_{3\left(ads\right)}{\left[D\right]}^{0.5}{\left[F\right]}^2\right)}\right]$

below, with tested literature expressions shown in

Table 8. Tested Rate Expressions from Literature.

Equation	Rate Law Expression
M1 [39]	$r=k_0{e}^{-\left(\frac{E_a}{RT}\right)}\times \left[\frac{\left({\left[E\right]}^{0.5}{\left[F\right]}^2\right)}{{\left(1+K_{1\left(ads\right)}{\left[F\right]}^2+K_{3\left(ads\right)}{\left[E\right]}^{0.5}\right)}^5}\right]$
M2 [41]	$r=\frac{{\left(\frac{K_1{K}_2{K}_3{k}_2{k}_1}{2}\right)}^{0.5}{L}^2{\left[E\right]}^{0.5}{\left[F\right]}^{0.5}}{{\left(1+{\left(\frac{2K_2{k}_1}{K_1{K}_3{k}_2}\right)}^{0.5}\frac{{\left[E\right]}^{0.5}}{{\left[F\right]}^{0.5}}+{\left(\frac{K_1{K}_2{K}_3{k}_2{k}_1}{2k_1}\right)}^{0.5}{\left[E\right]}^{0.5}{\left[F\right]}^{0.5}\right)}^2}$

Note: $A={\mathrm{C}}_3{\mathrm{H}}_8{\mathrm{O}}_3; B={\mathrm{H}}_2\mathrm{O}; C={\mathrm{O}}_2; D={\mathrm{CH}}_4; E={\mathrm{CO}}_2; F={\mathrm{H}}_2;Z={\mathrm{CH}}_3\mathrm{OH}; S=\mathrm{active} \mathrm{site}$.

.

4. Results and Discussion

4.1. Parameter Estimation

The values of the parameters were estimated using Non-Linear Regression Software (NLREG) and MATLAB 2017a. NLREG uses the Levenberg–Marquardt algorithm. In addition, parity plots were drawn to analyse the deviation of the predicted rates from the experimental ones. Models presenting an average absolute error (AAD) in excess of 20% were rejected. The parameter estimates for the successful rate expressions are shown in

Table 9. Parameter estimates for the overall power law model.

Estimated Parameter	Irreversible Expression	Reversible Expression
$k_{0f} \left[\frac{\mathrm{mol}}{{\mathrm{g}}_{\mathrm{cat}}·\min }\right]$	$4.3\times {10}^{10}$	$3.39\times {10}^5$
$E_f \left[\frac{\mathrm{J}}{\mathrm{mol}}\right]$	$8.78\times {10}^4$	$4.31\times {10}^4$
$m$	$1.04$	$-3.48\times {10}^{-5}$
$n$	$0.54$	$-0.017$
$o$	$1.78$	$9.3\times {10}^{-4}$
$k_{0r} \left[\frac{\mathrm{mol}}{{\mathrm{g}}_{\mathrm{cat}}·\min }\right]$	-	$1.38\times {10}^3$
$E_r \left[\frac{\mathrm{J}}{\mathrm{mol}}\right]$	-	$4.3\times {10}^4$
$q$	-	$4.03$
$r$	-	$-4$
$s$	-	$-8\times {10}^{-4}$

Note: $k_{0f}$ and $k_{0r}$ are the forward and reverse frequency factors; $E_f$ and $E_r$ are the activation energies for the forward and reverse reactions, respectively; $m$, $n$ and $o$ are the forward reaction rate exponents; $q$, $r$, and $s$ are the reverse reaction rate exponents.

,

Table 10. Parameter estimates for steam reforming models (with MATLAB 2017a).

Parameter	PLM	ER 3	ER 4	ER 5
$k_0 \left[\frac{\mathrm{mol}}{{\mathrm{g}}_{\mathrm{cat}}·\min }\right]$	$4.13\times {10}^{12}$ $\pm 3.69\times {10}^5$	$9.014\times {10}^{-1}$ $\pm 3.5\times {10}^{-2}$	$1.17\times {10}^5$ $\pm 9.9\times {10}^3$	$5.32\times {10}^7$ $\pm 4.1\times {10}^6$
$E_a \left[\frac{\mathrm{J}}{\mathrm{mol}}\right]$	$9.08\times {10}^4$ $\pm 2.2\times {10}^3$	$9.1208\times {10}^4$ $\pm 4.26\times {10}^3$	$9.32\times {10}^4$ $\pm 1.57\times {10}^2$	$9.11\times {10}^4$ $\pm 4.8\times {10}^3$
$m$	$2.29$ $\pm 0.07$	-	-	-
$n$	$12.8$ $\pm 0.97$	-	-	-
$K_p=K_3$	-	$1.17\times {10}^9$ $\pm 5.85\times {10}^7$	$1.13\times {10}^9$ $\pm 4.65\times {10}^6$	-
$K_B$	-	$3.85\times {10}^{-4}$ $\pm 1.99\times {10}^{-10}$	-	$6.85\times {10}^7$ $\pm 1.13\times {10}^6$
$K_A$	-	-	-	$4.2\times {10}^{-9}$ $\pm 1.6\times {10}^{-13}$
$K_1$	-	$1.79\times {10}^{-7}$ $\pm 1.6\times {10}^{-8}$	$1.56\times {10}^{-2}$ $\pm 1\times {10}^{-3}$	-
$K_c$	-	-	$2.74\times {10}^{11}$ $\pm 2\times {10}^6$	-
$K_W$	-	$4.67\times {10}^{-2}$ $\pm 1.11\times {10}^{-3}$	$1.35\times {10}^6$ $\pm 4.5\times {10}^4$	-
$K_F$	-	-	-	$2.73\times {10}^{-6}$ $\pm 2.37\times {10}^{-7}$

Note:$k_0$ and $E_{\mathrm{a}}$ are the frequency factor and activation energy, respectively; $m$ and $n$ are the power law model exponents; $K_p$ is the equilibrium constant; $K_3$, $K_B$, $K_A$, $K_1$, $K_c$, $K_W$, and $K_F$ are the adsorption constants.

,

Table 11. Parameter estimates for total oxidation models (with MATLAB 2017a).

Parameter	PLM	LHHW 2	LHHW 7
$k_0\left[\frac{{\mathrm{g}}_{\mathrm{cat}}}{\min }\right]$	$2.37\times {10}^{10}$ $\pm 2.9\times {10}^9$	$2.96\times {10}^5$ $\pm 3.78\times {10}^4$	$4.5\times {10}^{10}$ $\pm 1.8\times {10}^{-6}$
$E_a\left[\frac{\mathrm{J}}{\mathrm{mol}}\right]$	$8.9\times {10}^4$ $\pm 1.04\times {10}^3$	$7.26\times {10}^4$ $\pm 1.5\times {10}^3$	$5.46\times {10}^4$ $\pm 6.5\times {10}^3$
$m$	$-0.19$ $\pm 0.0227$	-	-
n	$3.19$ $\pm 0.24$	-	-
$K_p$	-	$1.08\times {10}^{-10}$ $\pm 2.25\times {10}^{-12}$	-
$K_B$	-	$4.56\times {10}^{-1}$ $\pm 2.05\times {10}^{-2}$	-
$K_W$	-	$5.383\times {10}^2$ $\pm 5.7\times {10}^1$	-
$K_1$	-	-	$8.45\times {10}^{-3}$ $\pm 7\times {10}^{-4}$
$K_2$	-	-	$5.07\times {10}^{-3}$ $\pm 4.68\times {10}^{-4}$
$K_3$	-	$1.7\times {10}^{-2}$ $\pm 1.23\times {10}^{-3}$	-
$K_6$	-	$3.62\times {10}^3$ $\pm 1.17\times {10}^2$	-

Note:$k_0$ and $E_{\mathrm{a}}$ are the frequency factor and activation energy, respectively; $m$ and $n$ are the power law model exponents; $K_p$ is the equilibrium constant; $K_B$, $K_W$, $K_1$, $K_2$, $K_3$, and $K_6$ are the adsorption constants.

and

Table 12. Parameter estimates for CO₂ methanation models (with MATLAB 2017a).

Parameter	PLM	ER 1
$k_0\left[\frac{{\mathrm{g}}_{\mathrm{cat}}}{\min }\right]$	$1.506\times {10}^3$ $\pm 6.20\times {10}^1$	$8.84\times {10}^6$ $\pm 7.2\times {10}^{-23}$
$E_a\left[\frac{\mathrm{J}}{\mathrm{mol}}\right]$	$8.9\times {10}^4$ $\pm 1.04\times {10}^3$	$8.12\times {10}^4$ $\pm 5.6\times {10}^3$
$m$	$-2.089$ $\pm 6\times {10}^{-2}$	-
$n$	$-0.2$ $\pm 4.87\times {10}^{-3}$	-
$K_p$	$2.79\times {10}^1$ $\pm 1.12$	$6.79\times {10}^{-5}$ $\pm 3.95\times {10}^{-6}$
$K_1$	-	$7.435\times {10}^2$ $\pm 2.76\times {10}^1$
$K_3$	-	$7.6\times {10}^4$ $\pm 8.3\times {10}^2$

Note:$k_0$ and $E_{\mathrm{a}}$ are the frequency factor and activation energy, respectively; $m$ and $n$ are the power law model exponents; $K_p$ is the equilibrium constant; $K_1$, and $K_3$ are the adsorption constants.

.

Parity plots of the various models showing the deviation between the experimental and predicted rates are shown in Figure 1, Figure 2 and Figure 3.

4.2. Thermodynamic Scrutiny

The successful mechanistic models were next subjected to thermodynamic scrutiny so as to evaluate their thermodynamic consistency. The most commonly used criteria are the Boudart–Mears–Vannice (BMV) rules [30,42] which are what is used here. The Van’t Hoff equation is the basis for thermodynamic scrutiny, and is given as follows:

$$\ln \left(K_j\right)=-\frac{\mathsf{\Delta }H}{RT}+\frac{\mathsf{\Delta }S}{R}$$

(47)

where $\mathsf{\Delta }H$ is enthalpy, $\mathsf{\Delta }S$ is entropy, $R$ is the universal gas constant, and $K_j$ is the adsorption constant, expressed as:

$$K_j=K_{j,0}\exp \left(-\frac{\mathsf{\Delta }{H}_{j,\mathrm{ads}}}{RT}\right)$$

(48)

where $K_{j,0}$ is the adsorption constant, and $\mathsf{\Delta }{H}_{j,\mathrm{ads}}$ is the heat of adsorption.

The BMV criteria are given below in Equations (49)–(51). The resulting thermodynamic scrutiny for each model is given in

Table 13. BMV Estimates for the thermodynamics scrutiny of steam reforming model.

Model	Parameter		Criterion 1	Criterion 2	Criterion 3
ER 4	$\Delta H_1$ [J mol−1]	$-1.87\times {10}^5$	$\Delta K_{J,1}<1$	$3.01\times {10}^{-8}>1.67\times {10}^{-11}$	$-17.32\le 33.00$
	$\Delta K_{J,1}$ [kPa−1]	$3.01\times {10}^{-8}$
	$\Delta S_1[$ J mol−1 K−1]	$-143.98$
	$\Delta H_2[$ J mol−1]	$-6.91\times {10}^4$	$\Delta K_{J,2}<1$	$2.35\times {10}^{-8}>1.37\times {10}^{-10}$	$-17.57\le 13.10$
	$\Delta K_{J,2}$ [kPa−1]	$2.35\times {10}^{-8}$
	$\Delta S_2[$ J mol−1 K−1]	$-146.024$

,

Table 14. BMV estimates for the thermodynamics scrutiny of total oxidation model.

Model	Parameter		Criterion 1	Criterion 2	Criterion 3
LHHW 2	$\Delta H_1$ [J mol−1]	$-6.49\times {10}^4$	$\Delta K_{J,1}<1$	$5.81\times {10}^{-9}>1.67\times {10}^{-11}$	$-18.96\le 12.40$
	$\Delta K_{J,1}$ [kPa−1]	$5.81\times {10}^{-9}$
	$\Delta S_1$ [J mol−1 K−1]	$-157.67$
	$\Delta H_{2 }$[J mol−1]	$-6.71\times {10}^4$	$\Delta K_{J,2}<1$	$4.43\times {10}^{-9}\cong 3.84\times {10}^{-9}$	$-19.23\le 12.77$
	$\Delta K_{J,2}$ [kPa−1]	$4.43\times {10}^{-9}$
	$\Delta S_2$ [J mol−1 K−1]	$-159.92$

and

Table 15. BMV estimates for the thermodynamics scrutiny of co₂ methanation model.

Model	Parameter		Criterion 1	Criterion 2	Criterion 3
ER 1	$\Delta H_1$ [J mol−1]	$-2.93\times {10}^4$	$\Delta K_{J,1}$ < 1	$1.811\times {10}^{-3}>6.80\times {10}^{-12}$	$-6.31\le 6.40$
	$\Delta K_{J,1}$ [kPa−1]	$1.811\times {10}^{-3}$
	$\Delta S_1$ [J mol−1 K−1]	−52.49
	$\Delta H_2$ [J mol−1]	$-4.75\times {10}^4$	$\Delta K_{J,2}$ < 1	$3.89\times {10}^{-3}>1.49\times {10}^{-7}$	$-5.55\le 9.47$
	$\Delta K_{J,2}$ [kPa−1]	$3.89\times {10}^{-3}$
	$\Delta S_2$ [J mol−1 K−1]	−46.14

.

$$\begin{array}{cc}I.& \exp \left(\frac{S_{j,0}^0}{R}\right)=K_{j,0}<1\end{array}$$

(49)

$$\begin{array}{cc}II.& K_{j,0}>\exp \left(-\frac{S_{j,\mathrm{g}}^0}{R}\right) \end{array}$$

(50)

$$\begin{array}{cc}III.& \ln \left(K_{j,0}\right)\le \frac{\left(12.2-0.0014\mathsf{\Delta }{H}_{j,\mathrm{ads}}\right)}{R} \end{array}$$

(51)

4.3. Comparison with Other Literature Studies

An activation energy comparison with available literature was done where the activation energy for the power law models was compared with available literature (a summary is given in

Table 16. Comparison with reforming activation energies from literature.

Reforming Process	Catalyst	Temperature (K)	Activation Energy (J/mol)	Ref.
Steam Reforming	Ni/CeZrCa	773–923	$8.91\times {10}^4$	This work
	Ni/CeO2	873	$10.34\times {10}^4$	[43]
	Ni/CeO2	673–973	$4.34\times {10}^4$	[44]
	Co-Ni/Al2O3	773–823	$6.33\times {10}^4$	[20]
Total Oxidation	Ni/CeZrCa	773–923	$8.88\times {10}^4$	This work
	N/A	N/A	$8.0\times {10}^4$ (pyrolysis)	[45]
	La2O3	773–873	$11.64\times {10}^4$ (Methane)	[46]
CO2 Methanation	Ni/CeZrCa	773–923	$8.86\times {10}^4$	This work
	La2O3	773–873	$8.78\times {10}^4$ (Methane)	[46]
	Ni/Al(O)x	555	$9.36\times {10}^4$	[39]

below). In the absence of heat and mass transfer resistances, the activation energies for steam reforming, total oxidation, and CO₂ methanation were found to be 8.09 × 10⁴ kJ/mol, 8.7 × 10⁴ kJ/mol, and 8.86 × 10⁴ kJ/mol, respectively. Variation between these activation energies may be attributed to variation between experimental conditions like steam–arbon ratio, oxygen–carbon ratio, or temperature ranges used. Other contributing factors might be the type of catalyst, feed flow rate, and the initial mole fractions of reacting species present in the solution.

5. Reactor Simulation

5.1. Modelling Approach and Solution Method

Transient one-dimensional pseudo-homogeneous and simplified heterogeneous reactor models were constructed following a cell-centred finite difference numerical discretization scheme programmed into MATLAB 2017a to draw a comparison between the suitability of the two approaches for simulating the reforming of crude glycol in the experimental reactor used for the kinetic studies. The expressions for these models are presented in

Table 17. One-dimensional pseudo-homogeneous model expressions [47].

${\epsilon }_{\mathrm{b}}\frac{\partial {\rho }_{\mathrm{g}}}{\partial t}+\frac{\partial }{\partial z}\left(u_z{\rho }_{\mathrm{g}}\right)=0$	(52)
$\frac{\partial P}{\partial z}=-\left(\frac{150\mu {\left(1-{\epsilon }_{\mathrm{b}}\right)}^2}{d_{\mathrm{p}}^2{\epsilon }_{\mathrm{b}}{}^3}+\frac{1.75{\rho }_{\mathrm{g}}\left(1-{\epsilon }_{\mathrm{b}}\right)}{d_{\mathrm{p}}{\epsilon }_{\mathrm{b}}{}^3}\left|u_z\right|\right)u_z$	(53)
$\left[\left(1-{\epsilon }_{\mathrm{b}}\right){\rho }_{\mathrm{p}}{C}_{\mathrm{p},\mathrm{s}}+{\epsilon }_{\mathrm{b}}{\rho }_{\mathrm{g}}{C}_{\mathrm{p},\mathrm{g}}\right]\frac{\partial T}{\partial t}+{\rho }_{\mathrm{g}}{C}_{\mathrm{p},\mathrm{g}}{u}_z\frac{\partial T}{\partial z}=\frac{\partial }{\partial z}\left({\lambda }_{\mathrm{e},z}\frac{\partial T}{\partial z}\right)+\frac{4U_{\mathrm{e}}}{d_{\mathrm{t}}}\left(T_{\mathrm{w}}-T\right)+S^{\prime }$ $S^{\prime }=\left(1-{\epsilon }_{\mathrm{b}}\right){\rho }_{\mathrm{p}}{\displaystyle \sum _k}\left(-\Delta H_{\mathrm{rxn},k}\right)\left({\eta }_k{R}_k\right)$	(54)
$\frac{\partial }{\partial t}\left({\epsilon }_{\mathrm{b}}{\rho }_{\mathrm{g}}{w}_j\right)+\frac{\partial }{\partial z}\left({\rho }_{\mathrm{g}}{u}_z{w}_j\right)-\frac{\partial }{\partial z}\left({\epsilon }_{\mathrm{b}}{D}_{L,j}\frac{\partial \left({\rho }_{\mathrm{g}}{w}_j\right)}{\partial z}\right)=S_j$$S_j=\left(1-{\epsilon }_{\mathrm{b}}\right){\rho }_{\mathrm{p}}{M}_j{\displaystyle \sum _k}\left({\nu }_{jk}\right)\left({\eta }_k{R}_k\right)$	(55)

and

Table 18. One-dimensional simplified heterogeneous model expressions [38].

${\epsilon }_{\mathrm{b}}\frac{\partial {\rho }_{\mathrm{g}}}{\partial t}+\frac{\partial }{\partial z}\left(u_z{\rho }_{\mathrm{g}}\right)={\displaystyle \sum _k}{a}_{\mathrm{v}}{k}_{\mathrm{g},j}{\rho }_{\mathrm{g}}\left(w_{\mathrm{s},j}-w_j\right)$	(56)
${\epsilon }_{\mathrm{b}}\frac{\partial }{\partial t}\left({\rho }_{\mathrm{g}}{w}_j\right)+\frac{\partial }{\partial z}\left(u_z{\rho }_{\mathrm{g}}{w}_j\right)-\frac{\partial }{\partial z}\left({\epsilon }_{\mathrm{b}}{D}_{L,j}\frac{\partial \left({\rho }_{\mathrm{g}}{w}_j\right)}{\partial z}\right)=a_{\mathrm{v}}{k}_{\mathrm{g},j}{\rho }_{\mathrm{g}}\left(w_{\mathrm{s},j}-w_j\right)$	(57)
${\epsilon }_{\mathrm{b}}{\rho }_{\mathrm{g}}C{p}_{\mathrm{g}}\frac{\partial T }{\partial t}+u_z{\rho }_{\mathrm{g}}C{p}_{\mathrm{g}}\frac{\partial T}{\partial z}-\frac{\partial }{\partial z}\left({\lambda }_{\mathrm{e},z}\frac{\partial T}{\partial z}\right)=h{a}_{\mathrm{v}}\left(T_{\mathrm{s}}-T\right)$	(58)
$\left(1-{\epsilon }_{\mathrm{b}}\right){\rho }_{\mathrm{p}}C{p}_{\mathrm{p}}\frac{\partial T}{\partial t}={\rho }_{\mathrm{p}}\left(1-{\epsilon }_{\mathrm{p}}\right){\displaystyle \sum _k}\left(-\Delta H_{\mathrm{rxn},k}\right)\left({\eta }_k{R}_k\right)-h{a}_{\mathrm{v}}\left(T_{\mathrm{s}}-T\right)$	(59)
$a_{\mathrm{v}}{k}_{\mathrm{g},j}{\rho }_{\mathrm{g}}\left(w_j-w_{\mathrm{s},j}\right)=\left(1-{\epsilon }_{\mathrm{b}}\right){\rho }_{\mathrm{p}}{M}_j{\displaystyle \sum _k}\left({\nu }_{jk}\right)\left({\eta }_k{R}_k\right)$	(60)

. The boundary conditions applied are summarized in

Table 19. Boundary conditions applied in one-dimensional reactor models.

Pseudo-Homogeneous		Simplified Heterogeneous
$\mathrm{Initial}\mathrm{Conditions} (t=0)$	$$\begin{gathered} \mathrm{Boundary} \\ \mathrm{Conditions} \\ (z=0,z=\mathrm{L}) \end{gathered}$$	$$\begin{gathered} \mathrm{Initial}\mathrm{Conditions} \\ (t=0) \end{gathered}$$	$$\begin{gathered} \mathrm{Boundary} \\ \mathrm{Conditions} \\ (z=0,z=\mathrm{L}) \end{gathered}$$
$P=P_{\mathrm{in}}$	$\mathrm{Inlet}: P=P_{\mathrm{in}}$ $\mathrm{Outlet}: P=P_{\mathrm{out}}$	$P=P_{\mathrm{in}}$	$P=P_{\mathrm{in}}$ $P=P_{\mathrm{out}}$
$T=T_{\mathrm{in}}$	$\mathrm{Inlet}: T=T_{\mathrm{in}}$ $\mathrm{Outlet}: \frac{\partial T}{\partial z}=0$	$T=T_{\mathrm{s}}=T_{\mathrm{in}}$	$T=T_{\mathrm{s}}=T_{\mathrm{in}}$ $\frac{\partial T}{\partial z}=0$
$W_j=W_{j,\mathrm{in}}$	$\mathrm{Inlet}: W_j=W_{j,\mathrm{in}}$ $\mathrm{Outlet}: \frac{\partial W_j}{\partial z}=0$	$W_j=W_{j,\mathrm{s}}=W_{j,\mathrm{in}}$	$W_j=W_{j,\mathrm{s}}=W_{j,\mathrm{in}}$ $\frac{\partial W_j}{\partial z}=0$

.

Additionally, two-dimensional steady-state pseudo-homogeneous and simplified heterogeneous reactor models were also constructed within COMSOL Multiphysics. The relevant expressions for these models are shown in

Table 20. Two-dimensional pseudo-homogeneous model expressions.

$D_{r,j}\left(\frac{{\partial }^2{C}_j}{\partial r^2}+\frac{1}{r}\frac{\partial C_j}{\partial r}\right)+D_{z,j}\frac{{\partial }^2{C}_j}{\partial z^2}-u_z\frac{\partial C_j}{\partial z}=\left(1-{\epsilon }_{\mathrm{p}}\right){\rho }_{\mathrm{p}}{M}_j{\displaystyle \sum _k}\left({\nu }_{jk}\right)\left({\eta }_k{R}_k\right)$

(61)

${\lambda }_{\mathrm{e},r}\left(\frac{{\partial }^{2}T}{\partial r^2}+\frac{1}{r}\frac{\partial T}{\partial r}\right)+{\lambda }_{\mathrm{e},z}\frac{{\partial }^{2}T}{\partial z^2}-C{p}_{\mathrm{g}} {\rho }_{\mathrm{g}} u_z\frac{\partial T}{\partial z}=\left(1-{\epsilon }_{\mathrm{b}}\right){\rho }_{\mathrm{p}}{\displaystyle \sum _k}\left(-\Delta H_{\mathrm{rxn},k}\right)\left({\eta }_k{R}_k\right)$

(62)

and

Table 21. Two-dimensional simplified heterogeneous model expressions.

$D_{r,j}\left(\frac{{\partial }^2{C}_j}{\partial r^2}+\frac{1}{r}\frac{\partial C_j}{\partial r}\right)+D_{z,j}\frac{{\partial }^2{C}_j}{\partial z^2}-u_z\frac{\partial C_j}{\partial z}=k_{\mathrm{g},j}{a}_{\mathrm{v}}\left(C_j-C_{j,\mathrm{s}}\right)$	(63)
${\lambda }_{\mathrm{e},r}\left(\frac{{\partial }^{2}T}{\partial r^2}+\frac{1}{r}\frac{\partial T}{\partial r}\right)+{\lambda }_{\mathrm{e},z}\frac{{\partial }^{2}T}{\partial z^2}-C{p}_{\mathrm{g}} {\rho }_{\mathrm{g}} u_z\frac{\partial T}{\partial z}=h{a}_{\mathrm{v}}\left(T-T_{\mathrm{s}}\right)$	(64)
$\frac{1}{r}{D}_{r,j}\frac{\partial }{\partial r}\left(r\frac{\partial C_{j,s}}{\partial r}\right)=\left(1-{\epsilon }_{\mathrm{p}}\right){\rho }_{\mathrm{p}}{M}_j{\displaystyle \sum _k}\left({\nu }_{jk}\right)\left({\eta }_k{R}_k\right)$	(65)
$\frac{1}{r}{\lambda }_{\mathrm{e},r}\frac{\partial }{\partial r}\left(r\frac{\partial T_{\mathrm{p}}}{\partial r}\right)=\left(1-{\epsilon }_{\mathrm{b}}\right){\rho }_{\mathrm{p}}{\displaystyle \sum _k}\left(-\Delta H_{\mathrm{rxn},k}\right)\left({\eta }_k{R}_k\right)$	(66)

, with boundary conditions presented in

Table 22. Boundary conditions applied in two-dimensional reactor models.

$$\begin{gathered} \mathbf{Inlet} \mathbf{Boundaries} \\ \left(\mathit{z}=0;0\le \mathit{r}\le \frac{{\mathit{d}}_{\mathit{t}}}{2}\right) \end{gathered}$$	$$\begin{gathered} \mathbf{Symmetry} \mathbf{Boundary} \\ \left(\mathit{r}=0;0\le \mathit{z}\le \mathit{L}\right) \end{gathered}$$	$$\begin{gathered} \mathbf{Tube} \mathbf{Wall} \mathbf{Boundary} \\ \left(\mathit{r}=\frac{{\mathit{d}}_{\mathit{t}}}{2};0\le \mathit{z}\le \mathit{L}\right) \end{gathered}$$
${P|}_{r,0}=P_{\mathrm{in}}$ ${C_j|}_{r,0}=C_{j,\mathrm{in}}$ ${T|}_{r,0}=T_{\mathrm{in}}$	${\frac{\partial C_j}{\partial r}|}_{0,z}=0$ ${\frac{\partial T}{\partial r}|}_{0,z}=0$	${\frac{\partial C_j}{\partial r}|}_{\frac{d_t}{2},z}=0$ ${\lambda }_{\mathrm{e},z}{\frac{\partial T}{\partial r}|}_{\frac{d_t}{2},z}=-U_{\mathrm{e}}\left(T_{\mathrm{w}}-T\right)$
$\mathrm{For} \mathrm{the} \mathrm{heterogeneous} \mathrm{model} \mathrm{only}, \mathrm{where} r_{\mathrm{p}}=r_{\mathrm{p}}^{\mathrm{s}}$, $D_{\mathrm{p},j}\frac{d{C}_{j,s}}{d{r}_{\mathrm{p}}}=k_{\mathrm{g},j}{a}_{\mathrm{v}}\left(C_j-C_{j,\mathrm{s}}\right), {\lambda }_{\mathrm{e},\mathrm{p}}\frac{d{T}_{\mathrm{p}}}{d{r}_{\mathrm{p}}}=h{a}_{\mathrm{v}}\left(T-T_{\mathrm{s}}\right)$

.

5.2. Review and Discussion of Simulation Results

5.2.1. Model Validation

Crude glycerol conversion predicted by the one-dimensional pseudo-homogeneous and heterogeneous reactor models developed in MATLAB were validated against the experimental work conducted by Ghani et al. [18], finding the absolute average deviation (AAD) to be 10.0% and 12.7% for each respective model. Furthermore, two-dimensional pseudo-homogeneous and heterogeneous reactor models developed in COMSOL were also validated against experimental results. The maximum absolute average deviation (AAD) between the conversions was 12.08 % and 13.12 % for the 2D PH and 2D HET models, respectively. The variations in the validation results may be as a result of the assumptions made like incompressible fluid, constant velocity and the variabilities in estimating thermo-physical properties of the reactants and products [48].

Additionally, the mole fraction was also validated with the experimental results.

Table 23. Mole fraction validation with experimental work.

Species	Absolute Average Deviation (AAD)
	1D PH	1D HET	2D PH	2D HET
O2	3.7%	16.2%	7.7%	16.2%
N2	9.1%	12.8%	3.9%	3.3%
CO2	8.6%	6.0%	18.2%	32.6%
H2	10.5%	11.6%	11.6%	20.1%
CH4	6.1%	8.4%	6.6%	0.6%

illustrates the average absolute error between the predicted and experimental mole fractions for each of the species noted. In general terms, the one-dimensional pseudo-homogeneous modelling approach produced simulated results closer to the experimental observations than other modelling approaches did. It is interesting to note varied comparisons in the errors observed for each species. In the case of oxygen, greatest accuracy was found for the 1D pseudo-homogeneous model, whereas methane outlet concentration was predicted best by the 2D heterogeneous model (despite there being significant errors using this model for oxygen, carbon dioxide, and hydrogen). Carbon dioxide appeared to be best represented by the one-dimensional models vs. two dimensional models.

Parity plots for both the one-dimensional pseudo-homogeneous and heterogeneous models are given in Figure 4. The two-dimensional pseudo-homogeneous model led to closer predictions of experimental results than the two-dimensional heterogeneous model. It was found that the numerical solution overpredicted CO₂ and N₂ and underpredicted H₂. It could be that the assumptions made for the numerical model at the early stages and exit deviated slightly from the experimental ones. In addition, the conversion from the experimental work was higher at the inlet conditions for the experimental work than it was for the numerical. Thus, the numerical model is better at predicting the conversion of crude glycerol found from the experiments between a range of 0.58–0.8.

5.2.2. Comparison between Modelling Approaches

A comparison was made between the one-dimensional models (developed in MATLAB) and the two-dimensional models (developed in COMSOL). Reactor profiles generated for temperature, species mole fraction, and crude glycerol conversion were compared between each model. In comparing the results of the one-dimensional pseudo-homogeneous and heterogeneous models, negligible difference was observed in terms of predicted species mole fractions (as seen in Figure 5). The only substantial difference noted was that of predicted temperature profile, where the heterogeneous model predicted a reactor outlet temperature of about 20 K higher than what was predicted by the pseudo-homogeneous model. Similar results were observed by Halabi et al. [49] and Lordanidis [50], and are likely the result of the incorporation of thin-film heat transfer resistance in the heterogeneous model (which is otherwise ignored in the pseudo-homogeneous model). The increased resistance to removing heat from the catalyst surface raises the catalyst surface temperature, raising the reaction rate and further increasing the rate of heat generation.

Comparison of the two-dimensional reactor models developed in COMSOL showed a higher conversion of crude glycerol in the first half of the reactor versus the one-dimensional models developed in MATLAB. This difference could be attributed to the inclusion of axial dispersion in the one-dimensional models (whereas this was omitted from the two-dimensional models). Temperature profiles generated by the two-dimensional models presented the same general trends as the one-dimensional cases. The reaction rate is high close to the entrance of the reactor due to the higher reactant concentrations. Exothermic reactions taking place at the entrance of the reactor may have contributed to the higher temperatures observed relative to the latter half of the reactor which was almost flat due to the dominance of endothermic reactions. The heat produced by the exothermic reactions was consumed by the endothermic reaction, with the overall autothermal reforming process being sustained throughout the reactor. As such, there was no substantial decrease in temperature at the outlet.

Crude glycerol conversion for each reactor modelling approach were compared in Figure 6. When the conversions were compared, it was seen that the predicted conversion was highest using the two-dimensional heterogeneous model, followed by the one-dimensional pseudo-homogeneous model, the two-dimensional pseudo-homogeneous model, and lowest using the one-dimensional heterogeneous model. This could be because of the additional terms added to the two-dimensional models, predicting stronger mixing of reactants and therein higher reaction rates and conversion. Another possible difference may be in the assumptions involved in each respective model. Noting this; however, it is important to bear in mind that each modelling approach presents similar predictions of crude glycerol conversion at the reactor outlet, with the most significant variation being seen within the first 20 mm of the reactor. Lastly, it was found that there was no significant variation in concentration in the radial direction for the two-dimensional models. In this case, combined with the better reflection of experimental data (as noted in Table 23), one-dimensional reactor modelling approaches are sufficient for modelling the autothermal reforming of crude glycerol.

5.2.3. Effects of Temperature

Temperature is an important factor that favours the kinetics of the system. From Figure 7, it can be seen that (at a particular space time) increasing the temperature increases the crude glycerol conversion of the system. According to Jimmy et al. [51], to achieve thermoneutral conditions for autothermal reforming of crude glycerol, the optimum conditions of S/C = 2.6 with an adiabatic temperature of 923 K yields maximum hydrogen. Furthermore, compared to the simulation results of Afabor et al. [25] (calculated for a space velocity of 6.0 gcat·min/mol vs. 12.71 gcat·min/mol considered here), much lower reactor outlet crude glycerol conversions are noted (e.g., for 673 K, Afabor et al. indicate ~40% conversion versus 10.9% indicated here; similarly, ~70% conversion is indicated in the previous study versus 46.3% indicated here).

5.2.4. Effects of Axial Dispersion

Axial dispersion was found to have negligible effect on the predicted exit mass fractions of each species, as shown in Figure 8. This can be attributed to the plug flow criterion ($L/d_p>50$). Additionally, at steady state, the complete mixing and the plug flow condition also ensures no back mixing. Thus, little to no variation is seen in the profile. Figure 9 also shows the effect of axial dispersion on the reactor temperature profile, revealing that the inclusion and exclusion of this term has virtually no impact on the predicted temperature profile. Overall, it is observed that the axial dispersion term has a negligible effect on model predictions.

5.2.5. Effects of Space Velocity on Reactor Performance

The overall aim of reactor design is to enhance its performance. In this work, the performance is analysed in terms of crude glycerol conversion and hydrogen yield. From Figure 10, it is seen that increasing the space time results in an increase in both crude glycerol conversion and hydrogen yield. The weight-time (W/F_A0) used were 12.71, 50.94, 76.47 and 101.94 gcat·min/mol with associated conversions of 38%, 59%, 73% and 88%, respectively. In fact, at the entrance of the reactor, a sharp increase in crude glycerol conversion and hydrogen yield was observed within the first half of the reactor length, followed by a slight continued increase in the second half. At a reactor length of 0.025 m, about 93 % of the crude glycerol conversion and hydrogen yield had been obtained and thus only 7 % more conversion and yield could be achieved by the remaining 0.02 m. Again drawing comparison with the preceding study of Afabor et al. [25], lower relative crude glycerol conversions were predicted (e.g., where Afabor et al. predicted ~74% conversion for a space velocity of 10.8 gcat·min/mol, simulation results shown here predicted 37.9%).

5.2.6. Effects of Reactor Catalyst Bed Length on Conversion

The degree of crude glycerol conversion observed in fixed-bed catalytic reactors is also dependent on the residence time, generally defined by the catalyst bed length. Logically, longer catalyst bed lengths result in greater residence time—at the cost of increased tube and catalyst materials. It is, therefore, important to understand the trade-off between catalyst bed length and observed crude glycerol conversion. Simulations were conducted at varying catalyst bed lengths of 45 mm, 75 mm and 100 mm (see Figure 11), resulting in crude glycerol outlet conversions of 84%, 88% and 90%, respectively. It can also be observed that a substantial amount of crude glycerol conversion occurs at the inlet of the reactor, within only a marginal incremental conversion observed beyond a catalyst bed length of 10 mm. Depending on the difficulty of downstream separations, it may even be reasonable to consider an 80% conversion cut-off, such that only a 10–15 mm catalytic bed length would be necessary.

6. Conclusions

The intrinsic kinetics of crude glycerol autothermal reforming (ATR) was studied over Ni/CeZrCa catalyst at conditions of 773–923 K at atmospheric pressure, steam-to-carbon ratio of 2.6, and oxygen-to-carbon ratio of 0.125. The intrinsic kinetic data provided by Ghani et al. [9] were not limited by heat or mass transfer, and conditions were assumed to be far from equilibrium. The kinetics were studied following a power law approach as well as by proposing several kinetic mechanisms following the LHHW and Eley–Rideal frameworks.

For steam reforming, it was found that the hypothesized Eley–Rideal models best described the reaction rate, with ER 4 ultimately being chosen as the preferred model following thermodynamic scrutiny, possessing an AAD of less than 10%. This model reflected a rate-limiting reaction between adsorbed intermediate and steam. The most suitable model for total oxidation was found to be the power law model. Generally, LHHW-based mechanistic rate expressions performed better than those based on an Eley–Rideal-style mechanism. The mechanistic model considering the molecular adsorption of crude glycerol was found to be the best mechanistic description, with an AAD of 14.6%. Lastly, the CO₂ methanation reaction was best described by an Eley–Rideal-type mechanism considering the adsorption of carbon dioxide, presenting an AAD of 5.8%.