Hybrid Phenomenological-Empirical Modeling of Catalyst Deactivation in an Industrial Steam Methane Reformer

1. Introduction

Hydrogen is an essential feedstock in the chemical industry and is extensively produced for large-scale industrial applications. Among the available hydrogen production routes, steam methane reforming (SMR) of natural gas remains the most widely applied technology. This process is well established and combines high conversion efficiency with relatively low production costs, making it the preferred option for industrial hydrogen production [1].

Currently, approximately 80% of global hydrogen production is derived from fossil resources, predominantly natural gas, which accounts for nearly 60% of the total [2,3,4]. This pathway, commonly referred to as gray hydrogen production, represents a mature, efficient, and economically competitive technology. SMR remains the dominant route, enabling large-scale hydrogen production at costs typically ranging from 1 to 2.5 USD kg⁻¹, thereby holding the current global hydrogen supply [5].

In industrial SMR plants, the desulfurized natural gas feed is preheated, mixed with steam, and often subjected to a pre-reforming step prior to entering the steam reformer. Inside the reformer, methane is converted into hydrogen and carbon monoxide, generating synthesis gas. Carbon monoxide is subsequently converted into carbon dioxide through the water–gas shift (WGS) reaction, producing additional hydrogen. The resulting gas stream is then purified, most commonly by pressure swing adsorption (PSA), which removes CO₂, CO, CH₄, N₂, and other residual components to produce high-purity hydrogen. Alternative purification routes include absorption processes and membrane-based separation. The purified hydrogen is then available for downstream applications [1,6,7].

The SMR process is based on the endothermic catalytic conversion of methane with steam over nickel-based catalysts, typically operated at high temperatures (750–900 °C) and moderate pressures (20–30 bar). Under these conditions, reformer performance arises from a strong coupling between chemical kinetics, radiative and convective heat transfer in the furnace, intra-particle diffusion, and pronounced axial and radial temperature gradients along the catalyst-filled tubes. In addition to these strongly interacting transport and reaction phenomena, the catalytic activity evolves continuously during long industrial campaigns due to progressive catalyst deactivation mechanisms such as sintering, carbon deposition, and poisoning. These effects alter the effective reaction rates and thermal profiles along the reformer tubes over time. Consequently, SMR constitutes a highly nonlinear and time-evolving reacting system, with direct implications for reactor design, process optimization, advanced control strategies, and tube mechanical integrity [8,9,10,11].

The mathematical modeling of steam reforming processes has evolved substantially over the past several decades, motivated by the need to describe the complex coupling between reaction kinetics, heat and mass transfer, and fluid-flow phenomena within industrial-scale reformers. Despite these advances, most reactor models assume constant catalyst activity and therefore do not explicitly represent the gradual loss of catalytic performance observed during extended industrial operation.

Early studies were largely based on pseudo-homogeneous models employing global reaction kinetics and simplified one-dimensional balances, which were mainly suitable for conceptual analyses and steady-state performance evaluation. The kinetic framework proposed by Xu and Froment [8] established the reference kinetic model for SMR and WGS reactions on nickel catalysts, incorporating competitive adsorption and equilibrium effects, and has since served as the foundation for most phenomenological models. Building upon this framework, subsequent studies adopted heterogeneous reactor models that explicitly account for intra-particle diffusion resistances, rigorous energy balances, and a more detailed description of heat transfer from the reformer furnace [9].

As modeling approaches matured, the literature expanded to address the coupling of the reformer with downstream process units, as well as process intensification and optimization strategies. Rajesh et al. [12] and Yu and Sosna [13] investigated multi-objective optimization of SMR and proposed heat-exchanger-type reformers aimed at improving thermal efficiency. In parallel, Schwaab et al. [14] and Oechsler et al. [15] introduced hybrid modeling frameworks combining CFD simulations at the particle or tube scale with reduced one-dimensional reactor models, enabling the representation of local transport phenomena at a computational cost compatible with dynamic simulation, state estimation, and industrial process control.

More recently, modeling efforts have shifted from the isolated reformer toward the full hydrogen production plant. Works such as those by Salem et al. [1] and Kumar et al. [16] developed integrated models including WGS and PSA sections, as well as alternative feedstocks such as biogas. These studies emphasized integration with commercial process simulators, exergy-based performance assessment, and energy and economic optimization, reinforcing the role of plant-wide modeling in supporting operational decision-making [9,11,17,18,19,20].

At a broader level, recent review studies have consolidated advances in both catalysis and process design for SMR and related hydrocarbon reforming routes. Ighalo and Amama [2] reviewed recent developments in steam methane reforming catalysis, highlighting the main catalyst deactivation mechanisms reported in the literature, including carbon deposition, sulfur poisoning, and thermally driven sintering of nickel particles. These mechanisms progressively reduce the number of active sites and alter catalyst effectiveness under industrial operating conditions. Complementarily, Mokheimer et al. [3] reviewed reforming routes for hydrogen production from methane, naphtha, diesel, and alcohols, addressing reactor design, pre-reforming strategies, catalyst selection, and operating conditions, and synthesizing the state of the art in kinetic modeling and CFD-based simulation of reformers. Together, these studies highlight the continuing need for modeling approaches capable of representing catalyst aging effects while maintaining sufficient phenomenological rigor and computational tractability for industrial applications.

Nevertheless, predictive modeling of catalyst deactivation in industrial fired steam methane reformers remains limited. Taji et al. [21] developed a one-dimensional heterogeneous reformer model incorporating a first-order catalyst deactivation law into a real-time optimization framework. Their study demonstrated the benefits of accounting for catalyst aging when determining optimal operating trajectories. However, catalyst activity was represented by a simple first-order decay model calibrated from methane conversion data, while the furnace thermal behavior and tube wall temperature were treated independently of the parameter estimation procedure. Consequently, the thermal response of the furnace and catalyst deactivation were not estimated simultaneously from industrial operating data.

Jokar et al. [22] subsequently proposed a dynamic phenomenological model employing a power-law deactivation expression calibrated with long-term plant data, providing a more flexible representation of catalyst aging than first-order models. Nevertheless, their work focused on an industrial membrane reformer operating under prescribed thermal conditions, where tube wall temperature remained a known boundary condition rather than an estimated variable. Moreover, the deactivation formulation assumed complete catalyst deactivation over sufficiently long times and therefore did not account for the residual catalytic activity commonly observed after the initial rapid sintering stage.

Other studies have increased the phenomenological rigor of reactor models through multidimensional and multiscale formulations. Cui et al. [11] developed highly detailed reactor models suitable for design and advanced control studies, although with substantial computational requirements that limit routine industrial application. Fabrik et al. [23] and Stoppacher et al. [24] investigated specific deactivation mechanisms, including sulfur poisoning and chemical looping, which differ from the long-term thermally driven catalyst aging typically encountered in conventional natural-gas-fired hydrogen plants.

At the plant-wide level, Lew et al. [25] proposed optimization frameworks for integrated hydrogen production networks considering economic, reliability, and safety objectives. However, these approaches do not explicitly model catalyst deactivation or the thermal behavior of individual reformer tubes. Likewise, several industrial patents rely on semi-empirical correlations, lookup tables, or simplified process models intended for process control and optimization rather than mechanistic prediction of catalyst aging [26,27,28,29,30,31].

Therefore, an important gap remains in the simultaneous estimation of catalyst activity and furnace thermal behavior using routinely available industrial operating data. Although phenomenological reformer models based on mass, energy, and momentum balances are well established, catalyst activity is generally assumed constant or described by simplified deactivation models, whereas tube wall temperature is typically imposed from plant measurements. As a result, the coupled evolution of catalyst aging and furnace thermal response over long industrial campaigns has received little attention.

The present work addresses this gap by proposing a hybrid phenomenological-empirical framework in which catalyst activity and tube wall temperature are estimated simultaneously from reconciled industrial operating data spanning more than one thousand days. The methodology combines a one-dimensional pseudo-homogeneous reformer model with a residual-activity deactivation model and an empirical wall-temperature correlation, preserving the thermal-kinetic coupling throughout the parameter estimation procedure. Unlike previous studies, the proposed framework does not require measured tube wall temperatures as model inputs, relying only on routinely monitored plant variables. This provides a computationally efficient approach for long-term performance assessment, catalyst condition monitoring, and industrial campaign analysis while maintaining the predictive capability required for practical engineering applications.

2. Process Description

2.1. Industrial Steam Methane Reforming Process

Hydrogen production via natural gas steam methane reforming is organized as an integrated sequence of unit operations, including feed pretreatment, conversion in an externally fired tubular reformer, adjustment of gas composition in water–gas shift reactors, and final purification by adsorption. A simplified process flow diagram of the industrial configuration considered in this study is shown in Figure 1.

In industrial practice, desulfurized natural gas is mixed with steam and fed to the reformer, where it is converted into synthesis gas rich in hydrogen and carbon monoxide. The reformate is subsequently processed in a water–gas shift (WGS) section and sent to a pressure swing adsorption (PSA) unit to produce high-purity hydrogen. Residual gas streams from the PSA are commonly recycled and used as fuel in the reformer furnace. Recent reviews confirm that steam methane reforming remains the dominant hydrogen production route worldwide, due to its technological maturity and strong integration with syngas-based chemical value chains [1,2,3,25].

Natural gas is first desulfurized to remove sulfur-containing compounds that irreversibly poison nickel reforming catalysts. Depending on feed composition, a pre-reforming stage may be employed to convert heavier hydrocarbons into methane, hydrogen, and carbon oxides, thereby reducing the tendency for carbon formation and protecting the main reformer [6,32,33].

After the pre-reforming stage, the process gas is fed to the primary reformer, where steam methane reforming reactions take place at high temperatures over nickel-based catalysts. Depending on the plant configuration, additional steam may be introduced upstream of the reformer to adjust the steam-to-carbon ratio. The industrial unit considered in this work employs a top-fired furnace configuration, in which burners are located above the radiant box and heat is transferred to vertically arranged catalyst-filled tubes. Tubular reformers may be designed with a wide range of tube and burner configurations [6,34,35]. The combined feed composition typically falls within molar ranges of 70–95% methane, 0–15% ethane, 0–10% propane and butanes, 0–5% CO₂, 0–2% N₂, and trace amounts of H₂S after desulfurization [1,3].

The main product of the reforming process is a hydrogen-rich gas, typically containing 70–75% hydrogen on a dry basis. The remaining components include methane (2–6%), carbon monoxide (7–10%), and carbon dioxide (6–14%). Key advantages of the SMR process include high hydrogen yield and the reuse of byproducts generated during reforming, contributing to overall process efficiency [36,37,38].

2.2. Fired Reformer Configuration

Large-scale industrial steam methane reformers consist of catalyst-filled tubes operating at high temperatures (800–950 °C) and pressures up to 30$\mathrm{atm}$. Under these conditions, the strongly endothermic reforming reactions require continuous external heat supplied by the furnace [6,34,39,40].

Steam methane reforming is the primary reaction occurring in the reformer (Eq. 1), while carbon monoxide produced during reforming can be converted into additional hydrogen through the water–gas shift (WGS) reaction (Eq. 2):

$$C{H}_4+H_{2}O\rightleftharpoons CO+3H_2$$

(1)

$$CO+H_{2}O\rightleftharpoons C{O}_2+H_2$$

(2)

Because steam reforming is strongly endothermic, high temperatures and low pressures favor the forward reaction. In contrast, the WGS reaction is mildly exothermic and is favored at lower temperatures, with negligible pressure dependence. Consequently, industrial reformers operate at high temperatures and moderate pressures to maximize methane conversion, while downstream WGS units are employed to further increase hydrogen production [32,34,39,41,42]. Operation at pressures below 20$\mathrm{atm}$ would require additional downstream compression due to the increase in the total number of moles during reforming [32,34,39,41].

The resulting thermal load has direct implications for both the material requirements and temperature distribution within the reformer tubes. In industrial units, catalyst-filled tubes are arranged inside a fired furnace and heated externally by combustion gases. In top-fired configurations, burners are located at the furnace roof, and heat is transferred predominantly by radiation from the flame and refractory surfaces to the tube walls. The furnace generally comprises a radiant section containing the burners and reformer tubes and a convection section used for heat recovery and process-stream preheating. Heat is transferred through the tube walls to sustain the endothermic reactions occurring within the catalyst bed [34,40,43].

Industrial reformers contain multiple tubes manufactured from high-temperature Cr–Ni alloys capable of withstanding prolonged operation under severe thermal conditions. The strong coupling between external heat transfer and reaction kinetics generates significant temperature gradients, making tube wall temperature a key operational variable for both process performance and equipment integrity [34].

2.3. Reforming Catalysts and Catalyst Deactivation

Nickel-based catalysts supported on alumina remain the industrial standard for steam methane reforming due to their favorable balance between catalytic activity, stability, and economic viability. Catalyst design must ensure high reforming activity, resistance to deactivation, and adequate mechanical properties for prolonged operation under severe reforming conditions, including high temperatures, elevated pressures, and large thermal gradients. In addition, catalyst geometry must provide sufficient external surface area while maintaining acceptable pressure drop within the packed bed [2,44].

Industrial reforming catalysts typically consist of nickel dispersed as small crystallites on porous supports, most commonly $\alpha$-Al₂O₃, which provides high thermal stability and mechanical strength under reforming conditions. Catalyst performance depends strongly on the availability of active nickel surface area, making the preservation of nickel dispersion a key factor for long-term operation. Promoters such as potassium and calcium are often incorporated to improve catalyst stability, modify nickel–support interactions, and suppress carbon accumulation through enhanced coke gasification [2,32,41,45,46,47,48,49].

Despite these stabilization strategies, catalyst activity inevitably declines during operation due to deactivation phenomena such as sulfur poisoning, carbon deposition, and thermal sintering of nickel particles. In industrial hydrogen plants, sulfur poisoning is largely mitigated through feed desulfurization, while operation at elevated steam-to-carbon ratios suppresses significant coke accumulation. Under these conditions, catalyst deactivation is primarily governed by temperature-dependent structural changes in the active phase, particularly nickel sintering [21,50,51].

Operating conditions strongly influence both the dominant deactivation mechanism and its severity. High temperatures accelerate nickel particle growth, reducing active surface area and catalyst effectiveness, whereas insufficient steam-to-carbon ratios may promote carbon formation through methane decomposition, hydrocarbon cracking, and the reverse Boudouard reaction. Although carbon deposition can block active sites, increase pressure drop, and induce mechanical stresses, industrial reformers are generally operated under conditions that minimize coke formation. Consequently, thermal and hydrothermal sintering frequently emerge as the prevailing mechanisms governing long-term catalyst deactivation in steam methane reformers [33,34,39,50,51].

The influence of temperature is particularly important because deactivation rates increase exponentially with temperature. While elevated temperatures may reduce carbon accumulation by promoting gasification reactions, these conditions simultaneously accelerate nickel sintering and the irreversible loss of active surface area. The combined effects of temperature and steam have been shown to enhance the mobility of nickel species and promote particle growth, resulting in progressive activity loss during extended operation [50,52,53]. Furthermore, interactions between sintering and carbon-related phenomena have been reported, as catalyst aging may increase the susceptibility of larger nickel particles to carbon formation [49,54,55].

From a modeling perspective, these observations highlight the importance of explicitly accounting for the dependence of catalyst activity on operating conditions, particularly temperature and reaction atmosphere. Under industrially relevant conditions with effective feed purification and steam-to-carbon ratios sufficiently high to suppress coke formation, sintering is expected to be the dominant mechanism controlling long-term catalyst deactivation. Accurate representation of temperature-dependent sintering kinetics is therefore essential for predicting catalyst aging, long-term reactor performance, and catalyst lifetime in industrial steam methane reformers [1,2,3,25]. These considerations form the basis for the catalyst activity model adopted in this work.

3. Process Modeling

3.1. Mathematical Modeling

The reactor was modeled using a one-dimensional pseudo-homogeneous formulation under quasi-steady-state conditions. In this approach, the reacting gas and catalyst phases are treated as a single effective phase, while intraparticle transport limitations are incorporated through an effectiveness factor. This formulation is widely employed in industrial steam reformer modeling because it offers a suitable balance between physical realism and computational tractability for reactor-scale simulations.

Catalyst deactivation was the only dynamic effect incorporated in the model, while gas composition, temperature, and pressure were assumed to reach steady state along the reactor length. This assumption is justified by the different time scales governing the system: thermal and compositional dynamics occur on the order of seconds to minutes, whereas catalyst deactivation due to sintering evolves over periods of months or years under industrial operating conditions.

The following assumptions were adopted:

• Plug-flow regime, characterized by a high Reynolds number;
• Only steam reforming and water-gas shift reactions are considered since methane decomposition and CO reduction are negligible due to the high steam-to-carbon ratio employed;
• Negligible radial gradients in the reactor tubes;
• A single tube is representative of all others;
• No intraparticle temperature or composition gradients;
• Uniform catalyst particle size and homogeneously distributed voids;
• Negligible axial mass and heat dispersion (Péclet number $\gg 1$);
• Ideal gas behavior (compressibility factor $\approx 1$);
• Constant effectiveness factor.

Mass balances for CH₄ and CO₂ are expressed in terms of conversion ($x_j$), as shown in Eqs.3–4. The reactor energy balance is given in Eq.5:

$${\displaystyle \frac{\partial x_{C{H}_4}}{\partial z}}={\displaystyle \frac{A_s}{F_{C{H}_4,0}}}{\rho }_{b}a\left({\eta }_1{r}_1+{\eta }_3{r}_3\right)$$

(3)

$${\displaystyle \frac{\partial x_{{CO}_2}}{\partial z}}={\displaystyle \frac{A_s}{F_{C{H}_4,0}}}{\rho }_{b}a\left({\eta }_2{r}_2+{\eta }_3{r}_3\right)$$

(4)

$$u_s{\rho }_g{c}_p{\displaystyle \frac{\partial T}{\partial z}}={\rho }_{b}a{\displaystyle \sum _{i=1}^3}{\eta }_i{r}_i\left(-\Delta H_{e,i}\right)-4{\displaystyle \frac{U}{d_{ti}}}\left(T-T_w\right)$$

(5)

where $A_s$ is the tube cross-sectional area, $F_{C{H}_4,0}$ is the inlet molar flow rate of methane, ${\rho }_b$ is the catalyst bed bulk density, and a is the catalytic activity. ${\eta }_i$, $r_i$, and $\Delta H_{e,i}$ are the effectiveness factor, reaction rate, and enthalpy of reaction i, respectively. Additional parameters include the interstitial velocity $u_s$, gas density ${\rho }_g$, heat capacity $c_p$, overall heat transfer coefficient U, inner tube diameter $d_{ti}$, and wall temperature $T_w$.

Pressure drop across the bed is calculated using the friction factor f from the modified Ergun correlation [56], as shown in Eqs. 6–7:

$$-{\displaystyle \frac{\partial p}{\partial z}}={\displaystyle \frac{f{\rho }_g{u}_s^2}{d_p}}$$

(6)

$$f={\displaystyle \frac{\left(1-{\epsilon }_b\right)}{{\epsilon }_b^3}}\left[{\displaystyle \frac{180\left(1-{\epsilon }_b\right)}{Re}}+4\right]$$

(7)

where $d_p$ is the particle diameter, ${\epsilon }_b$ is the bed void fraction, and $Re$ is the Reynolds number. Correlations and constants used are detailed in Appendix A.

3.2. Reaction Kinetics

Methane was considered the primary hydrocarbon participating in the reformer reactions. Heavier hydrocarbons present in the feed were assumed to be rapidly converted to methane and therefore were not explicitly represented in the kinetic model. Accordingly, the reaction network consists of the steam methane reforming reaction (Eq. 8), the water–gas shift reaction (Eq. 9), and the overall reforming reaction (Eq. 10) [8]:

$$C{H}_4+H_{2}O\rightleftharpoons CO+3H_2$$

(8)

$$CO+H_{2}O\rightleftharpoons C{O}_2+H_2$$

(9)

$$C{H}_4+2H_{2}O\rightleftharpoons C{O}_2+4H_2$$

(10)

The corresponding reaction rates are described by Langmuir–Hinshelwood expressions (Eqs. 11–13), which account for both surface reaction kinetics and competitive adsorption effects on nickel-based catalysts [57]:

$$r_1={\displaystyle \frac{k_1\frac{p_{C{H}_4}{p}_{H_{2}O}^{0,5}}{p_{H_2}^{1,25}}\left(1-\frac{p_{CO}{p}_{H_2}^3}{K_{e,1}{p}_{C{H}_4}{p}_{H_{2}O}}\right)}{{\left(1+K_{a,CO}{p}_{CO}+K_{a,H_2}{p}_{H_2}^{0,5}+\frac{K_{a,H_{2}O}{p}_{H_{2}O}}{p_{H_2}}\right)}^2}}$$

(11)

$$r_2={\displaystyle \frac{k_2\frac{p_{CO}{p}_{H_{2}O}^{0,5}}{p_{H_2}^{0,5}}\left(1-\frac{p_{C{O}_2}{p}_{H_2}}{K_{e,2}{p}_{CO}{p}_{H_{2}O}}\right)}{{\left(1+K_{a,CO}{p}_{CO}+K_{a,H_2}{p}_{H_2}^{0,5}+\frac{K_{a,H_{2}O}{p}_{H_{2}O}}{p_{H_2}}\right)}^2}}$$

(12)

$$r_3={\displaystyle \frac{k_3\frac{p_{C{H}_4}{p}_{H_{2}O}}{p_{H_2}^{1,75}}\left(1-\frac{p_{C{O}_2}{p}_{H_2}^4}{K_{e,3}{p}_{C{H}_4}{p}_{H_{2}O}^2}\right)}{{\left(1+K_{a,CO}{p}_{CO}+K_{a,H_2}{p}_{H_2}^{0,5}+\frac{K_{a,H_{2}O}{p}_{H_{2}O}}{p_{H_2}}\right)}^2}}$$

(13)

where $r_i$ is the reaction rate of reaction i, $k_i$ is the kinetic rate constant, $K_{e,i}$ is the equilibrium constant, and $p_j$ is the partial pressure of species j. The denominator accounts for competitive adsorption of CO, H₂, and H₂O on the catalyst surface, where $K_{a,j}$ denotes the adsorption equilibrium constant of species j.

The temperature dependence of the kinetic rate constants is described by the Arrhenius equation (Eq. 14) [8]:

$$k_i=k_{i,0}\exp \left({\displaystyle \frac{-E_{a,i}}{RT}}\right)$$

(14)

where $k_{i,0}$ is the pre-exponential factor, $E_{a,i}$ is the activation energy of reaction i, R is the universal gas constant, and T is the absolute temperature. Correlations and constants used are detailed in Appendix A.

3.3. Catalyst Deactivation Modeling

Under the operating conditions considered in this study, characterized by high reforming temperatures and steam-to-carbon ratios typically above 2.5, carbon formation is thermodynamically suppressed and sulfur poisoning is minimized by upstream feed desulfurization. Consequently, catalyst deactivation was assumed to be governed predominantly by nickel sintering, in agreement with previous studies on industrial steam reforming catalysts [6,50,53].

The dynamic behavior of catalyst activity was described using a generalized power-law deactivation model with residual activity (Eq. 15). This formulation extends the classical Levenspiel deactivation kinetic model by allowing the catalyst activity to approach a non-zero asymptotic value, thereby representing the residual activity commonly observed in industrial catalysts after prolonged operation [58,59]. The resulting activity decay is expressed as:

$${\displaystyle \frac{da}{dt}}=-k_d{\left(a-a_{\mathrm{eq}}\right)}^m$$

(15)

where a is the relative catalyst activity, $a_{\mathrm{eq}}$ is the residual activity at long operating times, $k_d$ is the deactivation rate constant, and m is the deactivation order.

Residual activity models have been proposed because the classical power-law formulation predicts complete catalyst deactivation, whereas many industrial catalytic systems retain measurable activity even after extended operating periods. Bain et al. [58] demonstrated that conventional power-law models are appropriate only when catalyst activity approaches zero and showed that the introduction of a residual activity term provides a more general description of catalyst aging. Monzón et al. [59] further demonstrated that neglecting residual activity may lead to biased estimates of intrinsic kinetic parameters, including activation energies and pre-exponential factors, particularly for catalysts that retain a stable fraction of active sites after prolonged operation.

The deactivation rate constant was assumed to follow an Arrhenius-type temperature dependence (Eq. 16). To reduce the correlation between the pre-exponential factor ($k_{d,0}$) and the activation energy ($E_{a,d}$), the Arrhenius equation was reparameterized using a reference temperature ($T_{\mathrm{ref}}$).

$$k_d=k_{d,0}exp\left[-{\displaystyle \frac{E_{a,d}}{R}}\left({\displaystyle \frac{1}{T}}-{\displaystyle \frac{1}{T_{\mathrm{ref}}}}\right)\right]$$

(16)

The reference temperature was defined as the harmonic mean of the outlet temperatures, because the Arrhenius equation depends on reciprocal temperature, leading to improved parameter scaling and numerical stability [60,61,62].

Catalyst activity was assumed to be uniform along the reactor length at each instant, varying only with time since catalyst sintering occurs over a much longer timescale than axial transport phenomena.

The model parameters ($k_{d,0}$, $E_{a,d}$, m, $a_{eq}$) were estimated using proprietary industrial operating data obtained from the reformer unit analyzed in this study.

3.4. Wall Temperature Modeling

The tube wall temperature is a key variable in the energy balance of the steam methane reforming process. In the studied unit, this temperature is measured manually using a pyrometer through inspection windows, targeting individual rows of tubes. Measurements are typically taken at the brightest regions of the tube surface, where temperatures are locally higher and the risk of material degradation is more critical. Consequently, in addition to uncertainties inherent to manual operation and instrument limitations, the recorded values tend to overestimate the actual average wall temperature.

To address this limitation, an empirical model was proposed to estimate a representative tube wall temperature from operating variables associated with furnace heat generation and process load. The formulation accounts for the thermal inertia of the reformer through accumulated variables obtained by first-order exponential smoothing of the fuel gas flow rate and the air-to-fuel ratio.

Steam methane reformers exhibit a slow thermal response because heat is stored and transported through refractory walls, catalyst-filled tubes, metallic tube walls, and combustion gases before reaching the reacting mixture. Consequently, the measured tube-wall temperature reflects not only the current firing conditions but also the accumulated thermal state of the furnace [63,64,65]. Exponential smoothing provides a computationally appropriate representation of this memory effect: each smoothed value is a weighted average of the current observation and all past values, with weights decaying geometrically with lag. This structure is mathematically equivalent to a first-order linear filter and is well suited to capture the dominant time scale of thermal diffusion in industrial reformers without requiring detailed heat transfer modeling.

The smoothing factor was determined from the characteristic dynamic response of the reformer. The autocorrelation function of the smoothed signal was then evaluated, and the characteristic time constant ($\tau$) was estimated as the lag at which the normalized autocorrelation decayed to $e^{-1}$. The smoothing factor $\beta$ was calculated as:

$$\beta =1-exp\left(-{\displaystyle \frac{1}{\tau }}\right)$$

(17)

The accumulated variables were then computed as:

$$E_{\mathrm{cum},k}=E_{\mathrm{cum},k-1}+\beta \left({\dot{m}}_{fuel,k}-E_{\mathrm{cum},k-1}\right)$$

(18)

$$R_{air/fuel,k}=R_{air/fuel,k-1}+\beta \left({\displaystyle \frac{{\dot{m}}_{air,k}}{{\dot{m}}_{fuel,k}}}-R_{air/fuel,k-1}\right)$$

(19)

where $E_{\mathrm{cum}}$ represents the accumulated fuel-related thermal contribution, $R_{air/fuel}$ is the smoothed air-to-fuel ratio, ${\dot{m}}_{fuel}$ is the furnace fuel mass flow rate, and ${\dot{m}}_{air}$ is the combustion air mass flow rate. The wall temperature was estimated by the following empirical relation:

$$T_w=b_1+b_2{E}_{\mathrm{cum}}+b_3{R}_{air/fuel}+b_4{\dot{m}}_{feed}$$

(20)

In this formulation, $E_{\mathrm{cum}}$ encodes the thermal energy supplied to the furnace over a time horizon governed by $\tau$, rather than reflecting only the instantaneous firing rate. The smoothed air-to-fuel ratio $R_{air/fuel}$ captures variations in combustion conditions and their effect on heat release characteristics. The feed flow rate ${\dot{m}}_{feed}$ accounts for the thermal load imposed by the endothermic reforming reactions. The parameters $b_i$ were estimated from industrial operating data.

Since the empirical wall temperature model is linear with respect to its parameters, a standardized effect was calculated for each regressor ($S_i$) by multiplying its estimated coefficient by the standard deviation of the corresponding input variable:

$$S_i=b_i\times {\sigma }_{x_i}$$

(21)

where $b_i$ is the estimated regression coefficient and ${\sigma }_{x_i}$ is the standard deviation of the i-th regressor computed from the industrial dataset. The resulting quantity has units of temperature and represents the expected change in wall temperature associated with a one-standard-deviation variation of the corresponding regressor, allowing the relative influence of variables with different units and scales to be compared.

3.5. Industrial Data Treatment

The identification and removal of outliers is essential to ensure the robustness of simulation analyses, as anomalous observations can distort statistical metrics and compromise model calibration, validation, and predictive capability. In industrial steam reforming units, process variability arises from both operational adjustments and transient disturbances, making systematic data filtering a necessary step to ensure representative steady-state behavior. The adopted methodology was based on the following criteria:

• A mandatory exclusion condition was established for days with reformer outlet temperature $T_{\mathrm{out}}<{780}^{\circ }\mathrm{C}$. Temperatures below this threshold are associated with atypical conversion levels and generally indicate shutdown or non-representative operating conditions.
• When the data point was not excluded by the primary criterion,, additional screening was performed based on the thresholds listed in Table 1, which reflect typical operating ranges for industrial steam reforming units.

Additionally, outliers were identified based on deviations from a 180-day moving median computed for key process variables, including reformer inlet and outlet temperatures, pressures, steam-to-carbon ratio, inlet natural gas molar flow rate, and furnace air-to-fuel ratio. The moving median provides a robust measure of central tendency, being less sensitive to transient fluctuations and extreme values than mean-based filters, which is particularly suitable for industrial time series affected by operational disturbances. The 180-day window was selected to represent the long-term operating envelope of the reformer, smoothing short-term process variability while preserving the slower trends associated with catalyst aging and progressive deactivation throughout the operating campaign.

4. Numerical Solution and Parameters Estimation

The proposed model consists of a coupled system of ordinary differential equations with axial spatial dependence and a time-dependent catalyst activity function. The governing equations (Eqs. 3–6) were solved using an adaptive explicit fourth- and fifth-order Runge–Kutta method, implemented in MATLAB^®, with relative and absolute tolerances of ${10}^{-6}$. Boundary conditions at the reactor inlet ($z=0$) were specified as $x_j=0$, $T=T_0$, and $p=p_0$. The reactor model was integrated simultaneously for all operating conditions in the dataset.

For each objective function evaluation, the catalyst activity profile was first computed analytically from the deactivation model presented in Section 3.3. The resulting activity profile, together with the wall temperature estimated by Eq. 20, were then supplied to the reactor model. Since both the catalyst deactivation model and the wall temperature model contain unknown parameters, all parameters were estimated simultaneously with the reactor simulation.

The model parameters were estimated using a weighted least-squares criterion. Residuals were normalized by the sample standard deviation of each response variable, yielding a formulation consistent with maximum-likelihood estimation under the simplifying assumption of independent, homoscedastic, and normally distributed residuals. This normalization accounts for the different numerical scales of the response variables and prevents either response from dominating the estimation procedure solely because of its magnitude. The corresponding objective function ($F_{\mathrm{obj}}$) was formulated as:

$$F_{\mathrm{obj}}={\displaystyle \frac{1}{N_x}}{\displaystyle \sum _{i=1}^{N_x}}{\left({\displaystyle \frac{x_{{\mathrm{CH}}_4,i}^{\exp }-x_{{\mathrm{CH}}_4,i}^{\mathrm{calc}}}{s_x}}\right)}^2+{\displaystyle \frac{1}{N_T}}{\displaystyle \sum _{i=1}^{N_T}}{\left({\displaystyle \frac{T_i^{\exp }-T_i^{\mathrm{calc}}}{s_T}}\right)}^2$$

(22)

where $s_x$ and $s_T$ are the sample standard deviations of the methane conversion and outlet temperature measurements, respectively, and $N_x$ and $N_T$ denote the corresponding numbers of observations. Since individual measurement uncertainties were unavailable for the industrial dataset, the adopted weighting should be interpreted as a scaling strategy rather than as an estimate of experimental error variance.

Parameter estimation was performed using a hybrid optimization strategy. A Particle Swarm Optimization (PSO) algorithm was first employed to identify promising regions of the parameter space, and the best solution was subsequently refined using the interior-point algorithm implemented in fmincon. The local optimization employed a step tolerance and an optimality tolerance of ${10}^{-8}$. Fixed random seeds were adopted to ensure reproducibility. To improve numerical conditioning, all decision variables were normalized to the interval $[0,1]$ during optimization and mapped back to their physical ranges before model evaluation.

Model performance was assessed separately for methane conversion and reformer outlet temperature using mean absolute percentage error (MAPE) and residuals in the original physical units were analyzed to verify the absence of systematic trends and bias.

Parameter uncertainty was assessed using residual bootstrap with 200 replicates. Synthetic datasets were generated by independently resampling the residuals of methane conversion and outlet temperature with replacement and adding them to the model predictions obtained using the optimal parameter set. For each bootstrap replicate, the model parameters were re-estimated using fmincon initialized at the optimal solution. Replicates that failed to converge or converged to solutions within a predefined tolerance of the parameter bounds were discarded. The remaining bootstrap estimates were used to calculate relative standard deviations (RSD) and correlation matrix, providing an empirical assessment of parameter uncertainty and dependence.

5. Results and Discussion

5.1. Validation of the Model under Fresh Catalyst Conditions

Before validating the complete dynamic model, the phenomenological model was evaluated independently, without the influence of the catalyst deactivation and wall temperature submodels. To this end, a fresh catalyst was assumed ($a=1$), and the tube wall temperature measured by optical pyrometry was used. The assumption of $a=1$ represents a reasonable approximation for the beginning of the operation lifespan, although a limited degree of initial deactivation cannot be completely ruled out. Nevertheless, its effect on this pointwise comparison is expected to be small. This analysis therefore assesses the predictive capability of the phenomenological equations, the kinetic model, and the auxiliary correlations independently of the parameters estimated in the subsequent stages of the study.

The largest deviation was observed for methane conversion (4.15%), followed by outlet pressure (0.96%) and outlet temperature (0.46%). These results show that the phenomenological model captures the overall behavior of the reformer under fresh catalyst conditions, providing a consistent foundation for incorporating the catalyst deactivation and wall temperature submodels used to describe long-term operation. The remaining discrepancies are consistent with the inherent simplifications of the model, the transport property correlations employed, and the uncertainties associated with industrial measurements.

5.2. Temporal Model Fitting

Figure 2 and Figure 3 present the temporal evolution of methane conversion and reformer outlet temperature throughout the industrial campaing after model fitting, comparing the model predictions with plant measurements. For confidentiality reasons, both variables are reported in dimensionless form by normalizing them with the mean of their respective measured values. This procedure preserves the temporal trends and relative variability of the data while allowing a meaningful assessment of model performance.

Figure 2 and Figure 3 illustrate that the model follows the temporal evolution of both variables over more than one thousand days of operation, reproducing both operational fluctuations and long-term trends. The goodness-of-fit metrics were calculated using the original dimensional data prior to normalization. For methane conversion, the model yielded a MAPE of 4.79% and a bias of $-2.00$ percentage points. In the industrial unit investigated, methane conversion is intentionally maintained within a narrow range around its average value throughout the campaign by adjusting operating conditions to compensate for the gradual loss of catalyst activity, and this behavior was successfully reproduced by the model. The reformer outlet temperature was predicted with a MAPE of 0.35% and a bias of 1.25$\mathrm{K}$, indicating that the proposed wall temperature formulation provide a consistent description of the thermal behavior of the unit.

The estimated wall temperature profile, also presented in Figure 3, highlights the coupling between this variable and the reformer outlet temperature, as well as between the thermal and catalyst deactivation models. The wall temperature directly influences the temperature inside the reformer tubes, which affects the catalyst activity decay rate, while the progressive loss of catalyst activity gradually modifies the thermal duty required to maintain reformer performance.

Figure 4 and Figure 5 present the parity plots comparing measured and predicted values over the entire industrial campaign. These plots provide an overall assessment of the model’s ability to reproduce simultaneously the two response variables employed during parameter estimation.

The predictions remain concentrated around the parity line for both response variables, indicating good agreement with the industrial measurements throughout the campaign. Considering the inherent variability of industrial operating data, this result suggests that the proposed formulation consistently captures both short-term operational fluctuations and the gradual changes associated with catalyst aging. As expected, the outlet temperature exhibits a tighter distribution around the 1:1 line than methane conversion, reflecting the higher variability typically associated with the measurement and operational control of conversion.

The combined analysis of the parity plots, MAPE, and bias indicates the absence of significant systematic deviations between measured and predicted values. For the outlet temperature, the bias of 1.25$\mathrm{K}$ corresponds to a small average overprediction, with no evidence of persistent over- or underestimation throughout the campaign. For methane conversion, the bias of $-2.00$ percentage points indicates a slight average underprediction, which is primarily attributed to experimental uncertainties, unmeasured operational disturbances, and the inherent simplifications of the phenomenological model. These results suggest that the proposed framework provides a consistent representation of both response variables employed in the simultaneous parameter estimation procedure.

Figure 6 presents the temporal evolution of the outlet molar fractions of H₂, CO, CO₂, and CH₄. Symbols represent industrial measurements, whereas solid lines correspond to the predictions of the phenomenological model for each species.

The species distribution provides a more stringent assessment of model performance than the analysis of overall methane conversion alone, since different combinations of the steam reforming and water–gas shift reactions may lead to similar conversion levels while producing distinct outlet compositions. Figure 6 demonstrates that the model reproduces the temporal evolution of all four species throughout the campaign, indicating that the adopted kinetic scheme consistently captures the coupling between the reforming and water–gas shift reactions even under progressively decreasing catalyst activity. The remaining discrepancies are consistent with measurement uncertainty, unmeasured operational disturbances, and the simplifying assumptions adopted in the phenomenological formulation.

In summary, the results presented in this subsection indicate that the simultaneous estimation of catalyst activity and wall temperature enables a consistent description of reformer performance throughout an extended industrial campaign while preserving the physical coupling between furnace thermal response and progressive catalyst deactivation.

5.3. Estimated Parameters of the Catalyst Activity Model

The parameters of the catalyst deactivation model (Eq. 15) and the wall temperature model (Eq. 20) were estimated simultaneously using the complete historical dataset from the industrial unit. This estimation strategy preserves the coupling between the thermal response of the reformer and the evolution of catalyst activity during parameter estimation, reducing inconsistencies arising from sequential error propagation.

The estimated parameters of the catalyst activity model and their corresponding relative standard deviations (RSD) are presented in Table 2.

The estimated deactivation activation energy ($E_{a,d}=123.12$ kJ mol⁻¹) lies within the range reported by Bartholomew [50,66] for Ni/$\gamma$-Al₂O₃ catalysts operating in hydrogen-rich atmospheres (31–159 kJ mol⁻¹). It is also consistent with values reported for industrial steam reformers (93 and 94 kJ mol⁻¹) [21,22], closely matching the value reported by Champon et al. [67] (126 kJ mol⁻¹) for the hydrothermal sintering of Ni/$\gamma$-Al₂O₃, as well as the activation energy associated with the surface mobility of Ni–OH species (137 kJ mol⁻¹) reported by Sehested et al. [68]. Taken together, these results support the hypothesis that thermal and hydrothermal sintering are the dominant deactivation mechanisms in the industrial unit investigated, consistent with the steam-to-carbon ratios above 2.5 maintained over the entire evaluation period, which suppress coke formation [43,50,69].

The deactivation rate constant at the reference temperature ($k_{d,0}=8.00\times {10}^{-3}$ d⁻¹) is of the same order of magnitude as values previously reported in the literature [22]. This parameter exhibited a very low relative standard deviation (RSD = 0.11%), whereas the activation energy showed the highest uncertainty among the activity model parameters (RSD = 25.41%). This difference should be interpreted in light of the characteristics of the industrial dataset. Because the reformer operates within a relatively narrow temperature range, as a consequence of the operating strategy adopted to preserve process stability and tube metallurgical integrity, the available thermal excitation is insufficient to completely decouple the two Arrhenius parameters (Eq. 16). Therefore, the larger uncertainty associated with $E_{a,d}$ primarily reflects the limited temperature variability present in the industrial data rather than an inherent limitation of the proposed model formulation.

The estimated deactivation order ($m=1.24$, RSD = 4.67%) remained close to unity, indicating an initially nearly linear activity decay followed by a gradual approach to the residual activity, consistent with the General Power Law Expression (GPLE). The inclusion of a residual activity term generally leads to moderate empirical reaction orders, typically between one and two, avoiding the high orders frequently obtained with Simple Power Law Expression (SPLE) models, where such values often arise from limitations of the mathematical formulation rather than reflecting the underlying deactivation mechanism [50,59,66,70,71].

From a mechanistic perspective, within the GPLE framework, values of m close to one are commonly associated with Ostwald ripening, whereas values approaching two are predominantly related to particle migration and coalescence, which is frequently reported for supported metal catalysts [72]. The value estimated in this work lies close to the lower bound of this mechanistic range, suggesting that atomic migration is likely to predominate during the industrial campaign, although the available data do not allow the simultaneous occurrence of both mechanisms to be statistically excluded.

Figure 7 presents the temporal evolution of the catalyst activity calculated from the parameters listed in Table 2. The predicted profile provides insight into the evolution of catalyst degradation throughout the industrial campaign.

A pronounced decrease in catalyst activity is observed during the first months of operation, from $a=1.0$ to values close to 0.5 within the first year of the campaign, followed by a progressively slower decline toward a residual activity level of approximately $a_{eq}=0.46$. This behavior is consistent with the sintering of freshly loaded nickel catalysts, in which initially small crystallites exhibit high surface mobility and undergo rapid coalescence. As the average particle size increases, surface mobility decreases, progressively reducing the sintering rate and leading the catalyst toward a more stable structural state [66,68,70].

The asymptotic convergence toward a nonzero residual activity is an intrinsic feature of the adopted formulation (Eq. 15), distinguishing it from SPLE, which assume complete catalyst deactivation over time [58,59]. This result is also consistent with the operating behavior of the industrial unit, where the progressive loss of catalyst activity is compensated through gradual adjustments of the reformer thermal conditions. Consequently, the catalyst activity and wall temperature models operate in an integrated manner to describe the evolution of the reforming process throughout the industrial campaign.

From a practical standpoint, the most relevant outcome is not the absolute value of catalyst activity at a particular instant, but its temporal evolution. Although the proposed formulation does not explicitly distinguish individual deactivation mechanisms, the aggregated activity factor captures their overall impact on reformer performance. This makes the estimated catalyst activity a useful indicator of the catalytic state of the unit, with potential applications in performance monitoring, catalyst aging assessment, maintenance planning, and catalyst replacement scheduling.

5.4. Estimated Wall Temperature Model Parameters

The estimated parameters of the empirical wall temperature model (Eq. 20) are presented in Table 3, together with their RSD and the relative contribution of each regressor ($S_i$), calculated from the observed variability of the corresponding input variables.

The characteristic time constant ($\tau$) estimated for the thermal response of the reformer was 27 days, supporting the assumption of high thermal inertia associated with the refractory lining, catalyst bed, and circulating gas within the furnace. This result justifies the use of accumulated variables in the wall temperature model, indicating that the thermal response of the reformer depends on the operating history of the furnace rather than solely on the instantaneous operating conditions, since disturbances in tube wall temperature produce slow responses in the steam reforming process [65]. Furthermore, empirical polynomial models similar to the one proposed here have previously demonstrated satisfactory performance in representing the thermal behavior of reformer furnaces, providing a computationally efficient reduced-order representation to detailed phenomenological models of radiative heat transfer [9,17,73,74].

The parameter $b_1$ (1051.80$\mathrm{K}$) represents the baseline wall temperature associated with the average thermal state of the furnace over the course of the operation. This interpretation is physically consistent, since tube wall temperature in industrial reformers results from the overall furnace thermal state, the heat distribution along the tube bundle, and the radiative interactions among neighboring tubes [75].

The coefficient $b_2$, associated with the accumulated fuel input ($E_{\mathrm{cum}}$), exhibited a sign consistent with the role of the thermal energy supplied by the furnace. An increase in fuel flow raises the heat flux transferred to the reformer tubes, directly affecting the tube wall temperature and, consequently, the catalyst deactivation kinetics described by the activity model. This behavior agrees with previous studies demonstrating the influence of fuel flow rate on wall temperature control in industrial reformers [76,77,78]. In terms of relative contribution, a one-standard-deviation variation in $E_{\mathrm{cum}}$ changes the predicted wall temperature by approximately 3.45$\mathrm{K}$, corresponding to the lowest, although still significant, sensitivity among the three regressors.

The coefficient associated with the smoothed air-to-fuel ratio ($b_3$) exhibited the highest sensitivity among all regressors. A one-standard-deviation variation in $R_{air/fuel}$ changes the predicted wall temperature by approximately 7.43$\mathrm{K}$, more than twice the sensitivity associated with $E_{\mathrm{cum}}$ and also greater than that associated with the reformer feed flow rate. It is also consistent with experimental observations showing that changes in the air-to-fuel ratio modify flame temperature and the spatial distribution of the radiative heat flux, affecting both the average tube temperature and the uniformity of tube heating [9,79]. Therefore, this variable provides an aggregated representation of combustion stoichiometry and heat transfer effects within the furnace.

The coefficient associated with the feed flow rate ($b_4$) presented the smallest estimated magnitude but an intermediate contribution when the typical variability of the input variable is considered ($S_i=$4.07$\mathrm{K}$). This behavior aligns with the endothermic nature of steam reforming, in which increasing the feed flow rate raises the reactor thermal demand, although combustion conditions remain the primary determinant of the external tube wall temperature [65,75,80]. Since the effects of heat supply are already partially represented by the other model variables, the coefficient $b_4$ predominantly captures the residual effect of the process load. The positive sign obtained, together with the absence of significant collinearity among the regressors and the good statistical precision of the estimate, suggests that this behavior reflects the operating strategy adopted in the industrial unit, where increases in process load are accompanied by coordinated adjustments in burner firing rate capable of maintaining, or even slightly increasing, the external tube wall temperature.

Figure 8 compares the wall temperature estimated by the model with the optical pyrometer measurements collected throughout the industrial campaign. The figure also presents the standard deviation band obtained from multiple measurements performed at different tubes and tube locations, allowing the variability inherent to the measurement procedure to be distinguished from the average trend represented by the model.

Figure 8 reveals two distinct operating regimes. During the first months of operation, the wall temperature increases gradually until reaching an approximately steady operating level. Thereafter, it remains nearly constant, exhibiting only low-amplitude fluctuations around its mean value.

Throughout most of the operating campaign, the model predictions remain within the standard deviation band of the pyrometer measurements, indicating that the model captures the average wall temperature trend despite the considerable scatter observed in the field measurements. As expected for a formulation based on smoothed operating variables, the model does not reproduce short-term fluctuations or gaps caused by the absence of measurements, such as the interval between approximately the 17th and 20th months of operation. These discrepancies primarily reflect the local and operator-dependent nature of pyrometer measurements, whereas the proposed model was developed to represent the global thermal response of the furnace.

The substantial variability observed in the experimental measurements highlights the practical value of estimating tube wall temperature from continuously monitored process variables. The proposed model provides a continuous and reproducible estimate of the reformer thermal state, reducing the influence of variability inherent to manual pyrometer measurements while capturing the long-term thermal evolution of the furnace. When coupled with the catalyst activity model, it enables an integrated description of the interaction between furnace thermal behavior and catalyst degradation throughout the industrial campaign, providing a practical framework for long-term performance monitoring and catalyst-aging assessment.

Monitoring the tube wall temperature distribution is widely recognized as essential for both reformer efficiency and operational integrity, serving as a basis for adjusting the furnace heat flux profile and preventing localized tube overheating [78,81].

5.5. Parameters Correlation

Figure 9 presents the correlation matrix obtained for the estimated parameters of the wall temperature and catalyst deactivation models. Its analysis provides insight into both parameter identifiability and the degree of coupling between these two submodels during the simultaneous parameter estimation procedure.

Figure 9 highlights how the Arrhenius reparameterization around the reference temperature ($T_{\mathrm{ref}}$) substantially reduced the correlation between the deactivation pre-exponential factor ($k_{d,0}$) and the apparent activation energy ($E_{a,d}$), whose statistical dependence is typically high because of the exponential form of the Arrhenius equation [60,61,62]. The residual correlation between these parameters was only $-0.25$, indicating a significant reduction in their statistical dependence despite the relatively narrow temperature range encountered during industrial operation.

Figure 9 also reveals relevant correlations among other parameter pairs. The deactivation order (m) exhibited a moderate negative correlation with $k_{d,0}$ ($-0.62$), a behavior commonly observed in nonlinear kinetic models, where different combinations of the rate constant and deactivation order can reproduce similar temporal activity trajectories within the resolution of the available experimental data.

Among the activity model parameters, the residual activity term ($a_{\mathrm{eq}}$) exhibited the strongest cross-correlations with the wall temperature model, particularly with $b_3$ ($-0.73$). This behavior is in accordance with the structure of the proposed framework, in which wall temperature strongly influences the deactivation kinetics, whereas the gradual loss of catalyst activity modifies the thermal duty required to maintain reformer performance. Consequently, part of the variability associated with the residual activity level also reflects systematic changes in the thermal operating conditions throughout the campaign, highlighting the coupling between the two submodels.

Within the wall temperature model, the strongest correlations were observed between $b_1$ and $b_3$ ($-0.71$), $b_1$ and $b_4$ ($-0.69$), and $b_2$ and $b_3$ ($0.66$). These relationships reflect the expected collinearity between the baseline thermal level and the operating variables included in the empirical formulation, which exhibit partially correlated behavior during plant operation. This result corroborates the empirical nature of the proposed model (Eq. 20), which prioritizes predictive capability based on continuously monitored process variables rather than strict parameter orthogonality.

Overall, the correlation matrix indicates levels of parameter dependence in agreement with a nonlinear estimation problem involving physically coupled variables. The simultaneous estimation strategy adopted in this work explicitly captures these dependencies, preventing the variability shared between the furnace thermal behavior and catalyst degradation from being attributed exclusively to either submodel. This integrated treatment contributes to a more physically consistent description of the long-term evolution of both the thermal state of the reformer and catalyst activity throughout the industrial campaign.

6. Conclusions

This work proposed and validated a hybrid phenomenological–empirical framework for describing the long-term operation of an industrial top-fired steam methane reformer. The methodology combines a one-dimensional pseudo-homogeneous reactor model with a temperature-dependent catalyst deactivation model and an empirical tube wall temperature correlation estimated from routinely monitored process variables. By simultaneously estimating the thermal and deactivation submodels using more than 1,000 days of industrial operating data, the proposed approach explicitly incorporates catalyst aging into reformer-scale simulations while preserving the coupling between furnace thermal behavior and reactor performance.

The results demonstrate that the framework is capable of consistently describing both methane conversion and reformer outlet temperature throughout an extended industrial campaign. The estimated deactivation parameters were physically meaningful and consistent with a sintering-dominated mechanism, while the predicted activity profile reproduced the expected progressive loss of catalyst performance and the approach toward a residual activity regime. Simultaneous estimation of catalyst activity and tube wall temperature proved particularly important for avoiding compensation effects between thermal and kinetic parameters, resulting in a more physically interpretable representation of the reforming process over long operating periods.

An important contribution of the proposed methodology is the estimation of a representative tube wall temperature from routinely monitored process variables, providing a continuous description of the reformer thermal state throughout the operating campaign. Rather than replacing pyrometric measurements, which remain essential for monitoring maximum tube temperatures and assessing tube integrity, the proposed thermal model complements available plant measurements by supplying information that is not continuously accessible in routine operation. This capability enables the direct integration of thermal behavior and catalyst aging within a single modeling framework. From an industrial perspective, the estimated catalyst activity can be interpreted as a macroscopic indicator of catalyst condition, supporting catalyst-aging assessment, campaign analysis, performance monitoring, and maintenance planning throughout the catalyst lifetime.

The proposed framework was intentionally developed as a computationally efficient alternative to more detailed multidimensional furnace and reactor models. Consequently, some limitations remain. Catalyst activity is represented by a lumped variable that does not distinguish individual deactivation mechanisms, the wall temperature correlation is empirical and therefore dependent on the operating envelope used for parameter estimation, and the activation energy of deactivation remains only moderately identifiable due to the relatively narrow temperature range imposed by industrial operating constraints. Nevertheless, the estimated parameters, uncertainty analysis, and model performance indicate that the adopted level of complexity is sufficient to capture the dominant interactions governing long-term reformer operation.

Overall, the proposed methodology bridges the gap between detailed reactor modeling and long-term industrial catalyst management. By integrating catalyst deactivation, furnace thermal response, and industrial operating data within a unified framework, it provides a practical tool for analyzing catalyst aging under realistic operating conditions and establishes a foundation for future applications involving online monitoring, catalyst replacement planning, operational optimization, and digital-twin implementations for industrial steam methane reformers.