Ab initio insights into support-induced sulfur resistance of Ni-based reforming catalysts

Abstract

Ni-based catalysts are well established for industrial H₂ production via methane steam reforming; however, their susceptibility to sulfur poisoning necessitates expensive desulfurisation and limits the development of low temperature processes using renewable feedstocks. Designing next-generation catalysts requires an atomic-level understanding of the factors that affect the catalyst sulfur tolerance, but this is difficult to obtain due to complex interactions between the Ni catalyst and non-inert metal oxide supports. In this work, we investigate the atomic-level mechanisms driving the support-induced sulfur resistance of Ni catalysts, emphasising the role of disorder in Ni-bound sulfur–oxygen adsorption complexes and support defect chemistry in promoting catalyst regeneration. The thermodynamic driving force for oxygen-mediated sulfur removal from a Ni(111) surface, which is indicative of the regenerative effects of support oxygen buffering, is investigated using grand canonical Monte Carlo (GCMC) sampling of a lattice model that is parameterised using density functional theory (DFT). The outcome is predictions of the equilibrium surface coverage and composition of co-adsorbed S and O atoms on Ni(111) at length scales that are inaccessible to DFT simulations. The GCMC predictions are validated using a fine-tuned machine learned interatomic potential to reveal entropic contributions for catalyst regeneration at experimentally relevant surface coverages, demonstrating an integrated approach for efficiently exploring the complex combinatorial space of adsorption complexes with near ab initio accuracy. Simulations of the surface chemistry of Ni(111) are complemented by predictions of the energetics of bulk defect formation in prototypical metal oxide support materials to provide insights into the proclivity for oxygen release and phase transformation during catalytic reactions. The computational modelling is correlated with experimental characterisation and methane steam reforming activity tests for H₂S-poisoned Ni nanoparticle catalysts, allowing us to rationalise the experimentally observed differences in the catalyst sulfur tolerance and establish strategies for future catalyst optimisation. The work demonstrates the integration of ab initio computational modelling, statistical sampling and machine learning, in a combined framework that complements experimental characterisation, to inform the rational design of catalyst support materials for sustainable H₂ production.

---

1 Introduction

Methane steam reforming (MSR) is an established industrial process that produces ∼95% of the global H₂ supply¹ via the conversion of natural gas (primarily CH₄, with smaller amounts of higher hydrocarbons) to syngas (mixtures of CO, CO₂ and H₂), at high temperature and pressure, in the presence of a catalyst. The commercial Ni-based catalysts are highly susceptible to sulfur poisoning by impurities in the feedstock, e.g., H₂S, SO₂, H₂SO₄ and/or COS, and therefore an expensive feed desulfurisation process is necessary to achieve sub-ppm sulfur concentrations.² The additional cost and complexity of feed desulfurisation also limits the development of biogas reforming processes for scalable H₂ production from renewable feedstocks, e.g., using solid oxide fuel cells³ or via combined steam and dry reforming.⁴ Understanding the factors that affect the catalyst sulfur tolerance is essential to enable the direct use of sulfur-containing feedstocks; a challenge that is particularly important for Ni-based catalysts as they are more economically viable than those based on platinum group metals (PGMs).

A number of strategies have been considered to enhance the sulfur tolerance of Ni-based catalysts, such as alloying with PGMs, including Au, Cu, Mn, Pd, Pt and Rh.⁵ Alloys are widely reported in the literature and are proposed to enhance the catalyst sulfur tolerance via different mechanisms, e.g., promoting sulfur scavenging by secondary metallic active phases,⁶ promoting sulfur oxidation and desorption at high temperatures^7,8 and suppressing the dissociative adsorption of feedstock poisons like H₂S.⁹ The optimisation of metal oxide supports is another effective strategy to enhance the sulfur tolerance of supported Ni nanoparticles during catalytic reforming reactions, with the mechanism widely hypothesised to involve oxygen buffering from reducible supports like CeO₂ and Y₂O₃.^10,11 In these materials, lattice oxygen is proposed to migrate from the support to the Ni active phase under reducing conditions at high temperatures, resulting in the oxidation and desorption of catalyst poisons e.g., C → CO₂ (ref. 12–16) and S → SO₂.^17–20 Similarly, a number of established chemical and electrochemical regeneration methods have been shown to restore the activity of poisoned Ni catalysts by modulating the transfer of oxygen to the poisoned Ni active sites. Chemical regeneration of sulfur-poisoned Ni catalysts can be achieved using exposure in steam, H₂ and/or O₂ depending on the degree of sulfur poisoning.^21,22 Electrochemical regeneration can also be used to control the O²⁻ spillover from both aqueous environments and Y₂O₃-stabilised ZrO₂ (YSZ) supports, towards sulfur-poisoned Pt and Ni species, enabling catalyst oxidative regeneration using a negative electrode potential.^23–25

Ab initio computational modelling methods, such as density functional theory (DFT), provide an atomic-level insight into the surface chemistry of sulfur-poisoned Ni nanoparticles. Atomic sulfur is often used to represent H₂S poisoning at low/medium surface coverage (θ_S) due to the predicted dissociative adsorption of H₂S → S on Ni(111), which does not cause surface reconstruction or sulfur penetration into the Ni bulk as observed at high θ_S.^26–29 DFT studies of oxygen-mediated sulfur removal from Ni(111) show that both atomic O and molecular O₂ (which adsorbs dissociatively) can lead to the sequential oxidation of S → SO → SO₂, which then desorbs at high temperatures.^30,31 These studies were limited to idealised adlayer representations of S, with θ_S = 0.25 monolayer (ML) and 0.5 ML, and do not account for variations in configurational entropy at intermediate coverages; therefore, whether the formation of SO₂ is thermodynamically or kinetically driven at experimentally relevant surface coverages remains unresolved. Constructing more experimentally relevant predictive models for S and O adsorption on Ni(111) requires extensive sampling of the large configurational space of adsorption complexes, which is computationally infeasible with DFT alone. Statistical sampling algorithms, such as grand canonical Monte Carlo (GCMC), must therefore be considered as they are well suited for exploring the configurational space of adsorption complexes on a lattice model of the surface, where adsorbates occupy predefined adsorption sites.^32,33 In GCMC, the ground state of the system is estimated by stochastically sampling a DFT-parameterised Hamiltonian through adsorbate perturbations such as adsorption, desorption or diffusion.³⁴ The GCMC approach allows the system to explore a wide range of chemically relevant surface configurations, producing extended models that are beyond the atomistic length scales afforded by DFT, whilst ensuring all accessible states contribute to the statistical ensemble when determining surface properties at thermodynamic equilibrium.

Lattice models simplify the sampling of the configurational space of adsorption complexes but neglect off-lattice effects, such as many-body lateral interactions and surface reconstruction, which can be non-negligible under experimental reaction conditions. To account for off-lattice effects, extended GCMC-predicted adlayers can be refined using classical interatomic potentials (IPs) to perform geometry optimisation and/or molecular dynamics simulations.^35–37 Classical simulations are a computationally efficient approach for modelling materials at the length scales unaffordable using DFT, but the accuracy of these simulations is dependent on that of the underlying IP. Modern machine learned interatomic potentials (MLIPs) offer a promising approach for balancing accuracy and computational efficiency by avoiding the predefined functional forms used in traditional IPs, enabling MLIPs to capture complex potential energy surfaces with greater flexibility. Recent advancements in neural network (e.g., SchNet,³⁸ PaiNN,³⁹ M3GNet,⁴⁰ CHGNet,⁴¹ and MACE⁴²) and Gaussian process-based (e.g., GAP⁴³) MLIPs have enabled more accurate modelling of chemical reactivity on transition metal surfaces.^44–46 Among these methods, the MACE⁴² architecture, based on message-passing neural networks (MPNNs) and the Atomic Cluster Expansion (ACE),⁴⁷ is popular as it requires less training data compared to other architectures; thus a MACE model provides a computationally tractable means for simulating off-lattice effects in extended surfaces with near ab initio accuracy.⁴⁸

Accurate simulations of poisoning and reactivity of Ni-based MSR catalysts are also very challenging to realise due to the interplay between oxygen buffering (causing catalyst regeneration) and phase transformations of the metal oxide support (causing catalyst deactivation). For example, Ni/γ-Al₂O₃ catalysts can undergo progressive Ni substitution for Al, resulting in the in situ transformation of Ni/γ-Al₂O₃ to spinel-type NiAl₂O₄.⁴⁹ Conflicting reports exist for the utility of Ni-based spinel-type oxides and whether they deactivate Ni-based catalysts⁵⁰ or enhance catalytic activity^51–56 and tolerance to S and C poisons⁵⁷ due to the facile formation of oxygen vacancies. Accurate predictions of the energetics of oxygen vacancy formation and substitutional doping for these support materials are non-trivial using DFT, particularly for reducible transition metal oxides (TMOs) e.g., TiO₂, and rare-earth metal oxides (REOs) e.g., CeO₂, which are experimentally reported to exhibit favourable oxygen buffering capacities.^58,59 The Coulomb self-interaction error (SIE) of local and semi-local DFT, when simulating materials with partially filled d or f orbitals, results in erroneous defect formation energies in TMOs and REOs;^60–62 therefore, it is necessary to use beyond-DFT methods with corrective schemes to combat the SIE. Hubbard corrected density functional theory (DFT+U) is a popular approach as it is computationally tractable for large systems (e.g., defects in large supercells) and involves an ad hoc energy correction applied selectively to localised orbitals, e.g., Ti 3d orbitals in TiO₂ and Ce 4f orbitals in CeO₂.⁶³ Despite the benefits of DFT+U in computational efficiency, the determination of appropriate simulation parameters, including the Hubbard U value and projector, is non-trivial for simulating defects in TMOs and REOs with accuracy that matches experimental observations, and care is therefore necessary in application.^64,65

In this work, a combined computational and experimental approach is adopted to investigate the enhanced sulfur tolerance of Ni nanoparticles on reducible metal oxide supports, with the aim of establishing strategies for future catalyst optimisation. We investigate the thermodynamic driving force for oxygen-mediated sulfur removal from Ni(111), indicative of the regenerative effects of support oxygen buffering, using GCMC sampling of a DFT-parameterised lattice model. The GCMC-predicted adlayers enable the prediction of the surface coverage and composition of competitively adsorbed S and O atoms as a function of temperature and the chemical potentials of S and O across an extended Ni(111) surface. The GCMC-predicted adlayers are validated using geometry optimisation simulations with a fine-tuned MACE MLIP to reveal entropic contributions and limitations to catalyst regeneration at experimentally relevant surface coverages. Simulations of the surface chemistry of Ni(111) are complemented by DFT+U predictions of the energetics of bulk defect formation (oxygen vacancies and Ni substitution) in prototypical metal oxide support materials, providing insights into the proclivity for oxygen release and phase transformation during catalytic reactions. The computational modelling is correlated with experimental characterisation (TPD-MS, XPS, ICP) and MSR activity testing of H₂S-poisoned Ni nanoparticle catalysts to rationalise the experimentally observed differences in the catalyst sulfur tolerance. The work demonstrates the integration of ab initio computational modelling, statistical sampling and machine learning to construct more realistic models of complex catalytic materials, which further complement experimental characterisation to inform future strategies for catalyst rational design.

2 Methodology

2.1 Electronic structure calculations

2.1.1 DFT. All electronic structure calculations were performed using the Fritz-Haber Institute ab initio materials simulation (FHI-aims) software package,⁶⁶ which uses an all electron numerical atom-centred orbital (NAO) basis set, interfaced with the Python-based Atomic Simulation Environment (ASE).⁶⁷ Periodic boundary conditions were applied using converged k-point sampling with the standard light basis set (2020), with equivalent accuracy to the TZVP Gaussian-type orbital basis set,⁶⁸ as decided after benchmarking of the bulk Ni vacancy formation energy (see the SI, section S1.1.1). Relativistic effects were accounted for using the zeroth order regular approximation (ZORA)⁶⁶ as a scalar correction. The system charge and spin were set to zero, given the reported quenching of Ni(111) surface magnetic moments following oxygen adsorption⁶⁹ and the temperatures of MSR far exceeding the Curie temperature of Ni (631 K), only below which long-range magnetic order is observed.⁷⁰ The mBEEF meta-GGA exchange correlation density functional was used,^71,72 as defined in Libxc,⁷³ providing the best accuracy compared to other local and semi-local functionals (see SI section S1.1.2). Dispersion corrections were not explicitly included as sulfur and oxygen bind strongly to Ni(111) through short-range chemisorption, which are well described by the mBEEF density functional.⁷¹ For such systems, long-range van der Waals interactions provide only minor contributions to adsorption energies, whilst any van der Waals correction may also be detrimental to the representation of the support material; therefore, no further dispersion corrections are included. Self-consistent field (SCF) optimisation of the electronic structure was achieved using a convergence criteria of 1 × 10⁻⁶ eV for the change in total energy, 1 × 10⁻⁴ eV for the change in the sum of eigenvalues and 1 × 10⁻⁶ e a⁻³₀ for the change in charge density. Unit cell equilibrium volumes (V₀) were calculated by fitting to the Birch–Murnaghan equation of state using ASE.⁷⁴ Geometry optimisation was performed using the quasi-Newton BFGS algorithm^75–78 with a force convergence criteria of 0.01 eV Å⁻¹. The pristine Ni(111) surface was modelled using a six layer symmetric periodic slab, of which the bottom three layers were frozen to mimic the system bulk, resulting in a converged surface energy in line with computational literature and experimental references (see SI section S1.1.3). A 20 Å vacuum gap was used in the direction perpendicular to the surface to eliminate artificial interactions between periodic images. A dipole correction was applied to compensate for the inhomogeneous electric field arising from surface adsorption. Adsorption energies were calculated as:

ΔE_Ads = E_{[Ni(111)+Ads]} − E_Ni(111) + μ_Ads (1)

where the chemical potential of the adsorbed species (μ_Ads) was calculated using the energies of isolated atomic S, atomic O, molecular SO and molecular SO₂.

2.1.2 DFT+U and defect calculations. All DFT+U calculations were performed with FHI-aims, using the on-site definition of the occupation matrix and the fully localised limit (FLL) double counting correction.⁶³ A Hubbard correction was applied to treat the Coulomb self-interaction of Ti 3d orbital electrons in tetragonal rutile TiO₂ and Ce 4f orbital electrons in cubic CeO₂. No Hubbard correction was applied for the Ni dopants or for γ-Al₂O₃. Hubbard U values for Ti 3d and Ce 4f orbital electrons were chosen as U^{Ti 3d} = 2.575 eV and U^{Ce 4f} = 2.653 eV, which are both valid with a refined atomic-like Hubbard projector, as defined in the SI section S1.2. Hubbard U values and projectors were simultaneously determined using a machine learning-based workflow, with the target of reproducing the bulk material covalency as calculated using hybrid-DFT, which results in numerically stable self-consistent simulations of point defects.⁶⁵ Defect calculations in γ-Al₂O₃, TiO₂ and CeO₂ were performed using the supercell sizes listed in the SI section S1.2, with suitable sizes to ensure a consistent defect concentration across the three systems whilst also accurately representing the dilute limit. Defect energies (ΔE_Defect) following substitution of a host metal atom (Al in γ-Al₂O₃, Ti in TiO₂ and Ce in CeO₂) with a Ni atom were calculated as:

ΔE_Defect = E_{Defective Bulk} + μ_Host − E_{Stoichiometric Bulk} − μ_Dopant (2)

where the chemical potentials μ_Host and μ_Dopant were calculated using the energy of bulk Ti (hexagonal close packed) as well as Al, Ce and Ni (all face-centred cubic). Oxygen vacancy formation energies (ΔE_OV) were calculated as:

ΔE_OV = E_{Defective Bulk} + μ_O − E_{Stoichiometric Bulk} (3)

where the chemical potential μ_O was calculated using half the energy of an isolated O₂ molecule. Defect calculations in TiO₂ and CeO₂ were performed using the occupation matrix release (OMR) method to initialise Ti³⁺ and Ce³⁺ polarons at nearest neighbour atoms to the defect. The DFT+U-predicted total energy (E^DFT+U) is pre-converged using fixed orbital occupancies until ΔE^DFT+U ≤ 0.001 eV, below which all orbital occupancies are calculated self-consistently.⁶³

2.2 Monte Carlo sampling

All lattice modelling and Monte Carlo sampling was performed using the Surface Science Modeling and Simulation Toolkit (SuSMoST) software package,³⁴ considering adsorption complexes of S, O, SO and their pairs, and the occupation of hollow HCP and hollow FCC active sites on Ni(111) motivated by our results in section 3.1 and 3.2. Full DFT geometry optimisation was performed for 70 symmetrically inequivalent pairs of adsorption complexes on either a 10 × 10 or 7 × 7 Ni(111) surface supercell within a 10 Å or 5 Å radial cutoff, respectively, as explained further in section 3.2, before calculating the energy of lateral interactions, ΔE_Lateral, using:where E_Ni(111) is the energy of the pristine surface, is the energy of a pair of adsorbates x at sites s₁ and s₂ for x ∈ {S, O} and s₁, s₂ ∈ {Hollow HCP, Hollow FCC}, is the energy of a single adsorbate x occupying site s₁ and is the energy of a single adsorbate x occupying site s₂. 35 adsorption complexes consisting of pairs of S–S, O–O and S–O atoms, with |ΔE_Lateral| ≥ 0.04 eV, were chosen for parameterising a pairwise Hamiltonian ([script letter H]) for subsequent GCMC sampling, based on the generalised lattice-gas model of adsorption monolayers by Akimenko et al.:⁷⁹where L is a set of lattice sites, σ_i is an adsorption complex at site i, ΔE_Ads(σ_i) is the adsorption energy of the adsorption complex at site i in the zero coverage limit and ΔE_lateral(σ_i, σ_j, r_ij) is the energy of lateral interactions between adsorption complexes at sites i and j, given the distance r_ij between the two sites. Geometry optimisation of S–O pairs with a short interatomic separation of 1.45 Å, corresponding to adsorption at neighbouring hollow HCP and hollow FCC active sites, resulted in atomic diffusion to other active sites and therefore these adsorption complexes were disregarded for subsequent GCMC sampling. Similarly, molecularly adsorbed SO was predicted to be less stable than individually adsorbed S and O atoms at low surface coverage, and therefore was not included in the GCMC sampling (see section 3.2).

GCMC sampling was performed on a hexagonal lattice of 30 × 30 centers with periodic boundary conditions, which was large enough to avoid finite size effects. Each Monte Carlo step involved 30 × 30 attempted moves, i.e., one attempt for each active site per step to change the state of the adsorbed layer through adsorption, desorption and surface diffusion of atomic S and O. The acceptance or rejection of a new configuration of the model adsorbed layer of S and O was determined using the Metropolis algorithm,⁸⁰ where a new configuration is accepted if the total energy ([script letter H]) is less than that of the previous configuration (i.e., Δ[script letter H] ≤ 0 eV) or, if Δ[script letter H] > 0 eV, the new configuration is accepted with the probability . One million Monte Carlo steps were used to reach thermodynamic equilibrium and then the same number of steps were used to calculate ensemble averages. The parallel tempering algorithm was used to improve convergence to equilibrium and calculate the temperature dependence of the predicted adlayer coverage and composition, while also accounting for variations in configurational entropy.⁸¹ The following temperatures were used for parallel tempering replicas: 300, 400, 600, 800, 1000, 1200, 1500 and 1700 K. Each simulation was performed with varying relative chemical potentials (μ^R) of sulfur (μ^R_S) and oxygen (μ^R_O) between −1 and 1 eV, which correspond to the adsorption energies of a single S or O atom on Ni(111) in the zero coverage limit, before geometry relaxation. Negative values of μ^R correspond to surfaces that are less likely to adsorb atoms in the zero coverage limit, whilst positive values of μ^R correspond to surfaces that are more likely to adsorb atoms in the zero coverage limit. We note that non-zero coverages are still possible for both positive and negative values of μ^R after geometry relaxation, due to entropic effects or attractive lateral interactions. To enable direct comparison with experiment, the relative chemical potentials used for GCMC sampling were mapped to gas phase partial pressures, corresponding to reservoirs of O₂ and H₂S, using ideal gas thermodynamics at the same temperature and a standard-state pressure of 1 bar:

μ^R_S(T, p) = ΔE^S_Ads + [G_H2S(T, p) − E_H2S] − [G_H2(T, p) − E_H2] (6)

where ΔE^S_Ads (ΔE^O_Ads) are the DFT-computed adsorption energies for a S (O) atom on Ni(111) in the zero-coverage limit; G_H2S, G_H2 and G_O2 are the Gibbs free energies of the isolated H₂S, H₂ and O₂ molecules, respectively, obtained from ideal gas thermochemistry using ASE; and E_H2S, E_H2 and E_O2 are the DFT-computed energies of the isolated H₂S, H₂ and O₂ molecules, respectively.

2.3 Many-body tensor representations

To quantify the differences in the GCMC-predicted spatial distribution of adsorbed S and O on Ni(111), the GCMC-predicted adlayers were encoded into structural fingerprints using many-body tensor representations (MBTRs),⁸² with the DScribe Python library.^83,84 Two-body MBTRs were used to encode pairwise interatomic distances between adsorbed S and O atoms as a smooth density distribution over a continuous grid, which was then discretised into five MBTR descriptors and reduced to a one-dimensional descriptor using principal component analysis (PCA) with the Scikit-learn Python library.⁸⁵ The principal component output from PCA (PC^MBTR) captures the most significant trends in the spatial disorder of co-adsorbed S and O in the GCMC-predicted adlayers. All hyperparameters for evaluating the MBTRs and PC^MBTR are listed in the SI section S2.

2.4 Interatomic potential training and inferencing

The GCMC predictions were validated using geometry optimisation calculations with a MACE (version 0.3.10) MLIP,⁴² providing a computationally efficient means to relax the high-coverage GCMC-predicted adlayers on the 30 × 30 Ni(111) surface (∼5800 atoms, surface area ∼50 nm²). The MACE MLIP was trained using the diverse dataset of 5921 DFT-optimised structures collected in the work, including isolated atoms and molecules (S, O, SO, SO₂ and SO₃), Ni(111) periodic slab models of different thicknesses and adsorption complexes involving S, O, SO and SO₂ at both low and high surface coverage on Ni(111). Training was performed using multihead replay fine-tuning of the off-the-shelf MACE-MPA-0 (medium) foundation model,⁴⁶ trained on approximately 146 000 unique materials in the Material Project Trajectory (MPTrj) dataset^86,87 and 3.2 million unique materials in a subset of the Alexandria dataset.⁸⁸ No dispersion correction was used and the model precision was set to float64. A randomly selected 4737 structures (80%) were used for model training, with the remaining 1184 structures (20%) used for validation. The Adam optimiser⁸⁹ was used to minimise a cost function comprised of an equally weighted average of energy and force errors, with the learning rate set to 0.01. The MACE model consists of two message-passing layers and employs a radial cutoff for learning interatomic interactions of 6 Å, resulting in a total receptive field of 12 Å, which is greater than the distance when lateral interactions between surface adsorbed pairs of S–S, O–O and S–O atoms decay to zero at low surface coverage, as computed using DFT. Fine-tuning was performed for 24 epochs, to balance cost and accuracy due to plateauing of the energy and force errors (Fig. 1(b) and (c), respectively). The fine-tuned model gave a training (validation) root mean squared error (RMSE) of 14.4 (14.2) meV per atom in total energies and 16.3 (17.2) meV Å⁻¹ in atomic forces. When inferenced on the full dataset, the pre-trained foundation model gave a RMSE of 1.43 × 10¹⁰ meV in total energies and 10.7 eV Å⁻¹ in maximum atomic forces, which were reduced by >99% upon fine-tuning the model as shown in the parity plots in Fig. 1(d) and (e).

Fig. 1 (a) Overview of the use of grand canonical Monte Carlo (GCMC) sampling and a fine-tuned MACE machine learned interatomic potential for studying the co-adsorption of S and O atoms on Ni(111) at thermodynamic equilibrium. The MACE model is fine-tuned from the MACE-MPA-0 pre-trained foundation model for 24 epochs, which results in a reduction in the (b) energy and (c) force errors until both start to plateau. When inferenced on the full dataset of DFT-optimised structures, the fine-tuned model yields a reduction in the RMSE in total energies and maximum atomic forces of >99% vs. the pre-trained foundation model, as shown in the parity plots for (d) total energies and (e) maximum atomic forces.

The fine-tuned MACE model was then used as the ASE calculator to run geometry optimisation calculations using the BFGS algorithm^75–78 with a force convergence criteria of 0.01 eV Å⁻¹. Six GCMC-predicted adlayers of differing coverages and intermixing of adsorbed S and O were validated using MACE: for μ^R_S = −1 eV, μ^R_O = −1 eV, −0.7 eV and −0.5 eV, and T = 600 K and 1200 K. The accuracy of the GCMC-predicted adlayers were validated by computing the root mean squared deviation (RMSD) of the S and O atomic positions (x and y co-ordinates) between the initial GCMC-predicted adlayers and the final MACE-optimised adlayers:

where x^GCMC_i and y^GCMC_i are the x and y coordinates of atom i (either S or O) in the initial GCMC-predicted adlayer and x^MACE_i and y^MACE_i are the corresponding coordinates in the final MACE-optimised adlayer.

2.5 Experimental characterisation

To investigate how support oxygen buffering affects the sulfur tolerance of the Ni catalyst, we selected three model supports spanning a range of reducibilities. γ-Al₂O₃ is chosen as a high surface area, structurally robust support material with negligible oxygen buffering behaviour.⁹⁰ Rutile TiO₂ is chosen as a moderately reducible support material, which can form oxygen vacancies and facilitate mild oxygen buffering at high temperatures.⁵⁹ CeO₂ is chosen as the prototypical support material for strong oxygen buffering under catalytic reaction conditions due to the ease of switching between the Ce³⁺ and Ce⁴⁺ oxidation states, and low oxygen vacancy formation energy.^58,90

The three supported catalysts of 10 wt% NiO on γ-Al₂O₃ (commercial, surface area = 140 m² g⁻¹), rutile TiO₂ (commercial, surface area = 20 m² g⁻¹) and CeO₂ (commercial, surface area = 20 m² g⁻¹) were synthesised using the standard incipient wetness impregnation method, where the support materials were first impregnated with a Ni nitrate precursor solution, then dried and calcined at 773 K for 2 hours to obtain the final catalyst samples.⁹¹ The catalysts were pelletised to a size of 250–355 μm and activated in a tube furnace, in a mixture of 10% H₂ in N₂ at 923 K for 10 hours. Scanning electron microscopy (SEM) was used to visualise the morphology of the prepared catalysts using a Zeiss Ultra 55 field emission electron microscope equipped with in-lens secondary electron and backscattered detectors. X-ray diffraction (XRD) was performed using a Bruker D8 Advance Davinci design unit to measure the NiO crystallite size in the prepared catalysts.

A 1 g portion of each catalyst was then saturated with H₂S at room temperature for 18 hours in a fixed bed reactor, using a feed gas of 100 ppm of H₂S in a mixture of 2.5% H₂ in N₂, with a relative humidity of 50% and a flowrate of 500 ml min⁻¹. The total sulfur content following room temperature saturation was quantified using inductively coupled plasma (ICP) analysis. As the focus of this work is to investigate the thermodynamic driving force for sulfur removal and catalyst regeneration, rather than the kinetics of sulfur adsorption under operating reaction conditions, the room temperature sulfur loading protocol provides a consistent baseline from which we assess the temperature-dependent catalyst regeneration behaviour. We note that the measured sulfur content for each catalyst is expected to be a high (upper bound) estimate, with reduced adsorption at higher temperatures. The surface speciation of the H₂S-poisoned catalysts, with a measurement depth of 5–10 nm, was analysed using X-ray photoelectron spectroscopy (XPS). Temperature programmed desorption-mass spectrometry (TPD-MS), using a Micromeritics Autochem II Chemisorption analyser linked with a MKS Cirrus 2 mass spectrometer, was used to track the desorption of H₂O, SO and SO₂ from the H₂S-poisoned catalysts under a fixed temperature ramp of 10 K min⁻¹, from room temperature to 1223 K, in N₂.

MSR activity testing was carried out in a low-pressure rig designed to flow dry gas mixtures of N₂, CH₄ (and higher hydrocarbons) and H₂ for catalyst pre-reduction. The dry gas composition used was 68.4% CH₄ and 3.6% C₂H₆, with a balance of N₂. The dry gas mixture is then combined with steam (following prior heating and evaporation in an oven) forming a reaction gas mixture that is flowed through a packed catalyst bed, contained in a quartz tube, within a furnace that is electrically heated up to 1223 K. The MSR activity for each H₂S-poisoned catalyst was evaluated at steady state, at temperatures of 873, 973 and 1073 K, under regulated outlet backpressures of 100, 120 and 150 mbar, respectively. During the reaction, the dry gas is combined with steam resulting in a steam to carbon ratio of 3 : 1, with a total gas flowrate of 200 ml min⁻¹. The quartz tube (diameter 0.4 cm) was loaded to a 1.5 cm bed length, equating to 0.097 g (0.094 cm³) of catalyst and 0.155 g (0.094 cm³) of SiC inert dilutant. We note that the studied support materials are chosen as model systems to investigate the key principles driving the catalyst sulfur tolerance, but are not immediately compatible with existing industrial MSR processes due to differences in the catalyst form (i.e., pellets vs. powders) and thermal instability at very high temperatures over long timescales.

3 Results and discussion

3.1 Atomic and molecular adsorption on Ni(111)

To ascertain the number of non-equivalent adsorption sites on Ni(111), atomic S and O were adsorbed at the four initial positions illustrated in Fig. 2(a): hollow HCP, hollow FCC, atop and bridge sites. Geometry optimisation of atomic S adsorbed at both atop and bridge sites resulted in S diffusion to the hollow HCP site, whilst atomic O adsorbed at atop and bridge sites diffused to hollow HCP and hollow FCC sites, respectively. The hollow HCP sites in Fig. 2(b) and (d) and the hollow FCC sites in Fig. 2(c) and (e) were therefore determined to be the relevant non-equivalent sites for adsorption.

Fig. 2 (a) The four studied adsorption sites on the Ni(111) surface: (1) hollow HCP, (2) hollow FCC, (3) atop and (4) bridge. The unit cell boundaries are denoted with black dashed lines. (b)–(i) The most stable single atom (S and O) and molecular (SO and SO₂) adsorption complexes on a 1 × 1 Ni(111) surface, calculated using DFT with the mBEEF exchange correlation density functional, where (b) and (c) correspond to S adsorption, (d) and (e) correspond to O adsorption, (f) and (g) correspond to SO adsorption and (h) and (i) correspond to SO₂ adsorption. (a)–(i) are top down views of the Ni(111) surface and the bottom row is a side view for adsorption complexes (f)–(i). The corresponding adsorption energies for the adsorption complexes (b)–(i) are listed in the SI section S1.1.4.

Both atomic S and O strongly chemisorb on the Ni(111) surface and display an energetic preference for adsorption at hollow FCC sites, by 0.05 eV for S and 0.23 eV for O. The trends in adsorption energies and site preferences are in agreement with computational literature detailed in SI section S1.1.4, although the absolute values of adsorption energies are found to vary slightly with the choice of exchange correlation density functional, as GGAs from the literature tend to underbind,⁹² and the choice of Ni(111) surface model parameters.^29,93–95 The adsorption of molecular SO was also considered, with both S and O directly bonded to the surface. At both hollow HCP and FCC sites, S-bound SO was calculated to be more energetically stable by 2.35 eV and 2.10 eV, respectively. Finally, we tested SO₂ adsorption at the four initial positions in Fig. 2(a), from which the non-equivalent adsorption sites were atop and bridge sites in Fig. 2(h) and (i), respectively. SO₂ is calculated to be most stable when S occupies the bridge site of Ni(111), as is reported experimentally,⁹⁶ with the same preferential stability as reported in the DFT study of Liu et al.⁹⁵ All calculated adsorption energies are reported in SI section S1.1.4.

3.2 Pairwise and many-body lateral interactions on Ni(111)

The four non-equivalent adsorption complexes of atomic S and O in Fig. 2(b)–(e), were used to construct new adsorption complexes of S–S, O–O and S–O pairs at low surface coverage on a 10 × 10 Ni(111) surface (for S–S and O–O pairs) and a 7 × 7 Ni(111) surface for S–O pairs (to reduce computational cost at no detriment to accuracy). Following geometry optimisation, the energies of adsorbed single atoms and pairs were then used to compute lateral energies (E_lateral, defined in section 2.2, eqn (4)), which are plotted in Fig. 3(a)–(c) for pairs of S–S, O–O and S–O, respectively. Lateral interactions are repulsive for all pairs in Fig. 3(a)–(c), indicating that the O-mediated removal of adsorbed S occurs at high surface coverage and would require a large supply of O atoms to the surface to overcome the repulsive lateral interactions between adsorbed S and O, e.g., from a reducible metal oxide support with a large oxygen buffering capacity or using a high partial pressure of O₂ gas during experimental catalyst regeneration. All adsorption complexes corresponding to |E_Lateral| ≥ 0.04 eV, i.e., green markers in Fig. 3(a)–(c), were used to parameterise the pairwise Hamiltonian ([script letter H], defined in section 2.2, eqn (5)) for GCMC sampling. Geometry optimisation of S–O pairs at low surface coverage reveals the instability of short-range interactions of ≤1.45 Å between adjacent hollow HCP and hollow FCC sites, which results in atomic diffusion to neighbouring sites in Fig. 3(d) and (e). We therefore do not include short-range S–O interactions in the GCMC sampling by assigning E_lateral = ∞ eV within the lattice model for both initial configurations in Fig. 3(d) and (e).

Fig. 3 Lateral energies between adsorbed (a) S–S, (b) O–O and (c) S–O atomic pairs, at low surface coverage on Ni(111), calculated using DFT with the mBEEF exchange correlation density functional. Green (red) markers correspond to adsorption complexes that are included (not included) in the pairwise GCMC Hamiltonian. The marker shape corresponds to the type of active site occupied by each atom in the pairs. The initial (top row) and final optimised geometries (bottom row) for DFT relaxations of short-range S–O interactions, where S occupies a hollow-HCP site and O occupies a hollow-FCC site in (d) and (f), whilst S occupies a hollow-FCC site and O occupies a hollow-HCP site in (e) and (g). Adsorption complexes (d) and (e) correspond to low surface coverage on a 7 × 7 Ni(111) surface, whilst complexes (f) and (g) correspond to high surface coverage on a 1 × 1 Ni(111) surface. The relative energy for each adsorption complex (d)–(g), calculated using eqn (9), is listed underneath each subfigure.

We investigate the validity of excluding short-range S–O interactions from the GCMC sampling, which would create the conditions necessary for the oxidation of S → SO, by considering how the S and O surface coverages affect the energetics of S oxidation. The geometry optimisation simulations in Fig. 3(d) and (e) were repeated on a smaller 1 × 1 Ni(111) surface in Fig. 3(f) and (g), respectively, corresponding to a higher surface coverage, before evaluating the relative stability (ΔE_Relative) of an adsorbed SO molecule at the most stable hollow-FCC site vs. atomic S and O, using:

ΔE_Relative = E^n×n_SO/Ni(111) − E^n×n_S,O/Ni(111) (9)

where E^n×n_SO/Ni(111) is the energy of a geometry optimised SO molecule adsorbed at a hollow-FCC site on an n × n Ni(111) surface and E^n×n_S,O/Ni(111) is the energy of a geometry optimised pair of S and O atoms adsorbed at an initial interatomic separation of 1.45 Å on an n × n Ni(111) surface.

Comparing the relative energies in Fig. 3(d)–(g), there is a significant site-dependence in the energetic feasibility of S oxidation to SO, where relaxation of S adsorbed at hollow-FCC sites and O adsorbed at hollow-HCP sites dramatically reduces ΔE_Relative compared to relaxation of S adsorbed at hollow-HCP sites and O adsorbed at hollow-FCC sites. This observation is consistent with the spin-polarised DFT study of Das and Saida, who calculated ΔE_Relative = 0.41 eV for S adsorbed at a hollow-FCC site and O adsorbed at a hollow-HCP site and ΔE_Relative = 2.98 eV for both atoms adsorbed at hollow-FCC sites, on a 2 × 2 Ni(111) surface.⁹⁷ Our results further show a strong coverage-dependence for the feasibility of S oxidation, as shown by the reduction in ΔE_Relative from 0.57 eV to 0.01 eV by increasing the surface coverage from Fig. 3(e)–(g). The pairwise GCMC Hamiltonian, which excludes short-range S–O interactions that are energetically unfavourable at low surface coverage, is concluded to be valid for simulated adlayers with low θ_S and θ_O only, shown as the lighter regions in the GCMC-predicted isotherms in Fig. 4(a) and (b), as well as regions of low intermixing between S and O shown as the lighter regions in Fig. 4(c). In these regions, strong adsorbate interactions with the Ni(111) surface exceed any attractive lateral interactions between adsorbed S and O as may be required for the formation of oxidised sulfur species.

Fig. 4 GCMC-predicted surface coverages of (a) S and (b) O at 600 K for relative chemical potentials of S (μ^R_S) and O (μ^R_O) ranging between −1 eV and 0.2 eV, as defined in section 2.2. (c) The principal component derived from two-body many-body tensor representations (PC^MBTR, discussed in the SI section S2), which encodes the pairwise interatomic distances between adsorbed S and O atoms across 10 GCMC-predicted adlayers for 441 combinations of μ^R_S and μ^R_O at 600 K. The secondary axes in (a), (b) and (c) show the equivalent gas phase thermodynamic control variables corresponding to the relative chemical potentials, including the ratio of partial pressures (p) of H₂S to H₂ (for a fixed p_H2 = 1 bar) and the partial pressure of O₂, which were obtained from ideal gas thermodynamics at the same temperature and a standard-state pressure of 1 bar. (d) The root-mean-square deviation (RMSD) in S and O x and y atomic co-ordinates, between GCMC-predicted and MACE-reoptimised adlayers. Bars represent the mean RMSD for each μ^R_O value at T = 600 K and 1200 K. Error bars represent the standard deviation of the RMSD. All bars correspond to μ^R_S = −1 eV, thereby testing the validity of adlayers with varied intermixing of adsorbed S and O atoms on Ni(111), which increases for larger values of μ^R_O.

Under sulfur-rich conditions (μ^R_S → −1 eV), the GCMC-predicted isotherm in Fig. 4(a) predicts a large sulfur coverage of up to 0.45 ML that is thermodynamically stable even at extremely low H₂S feed concentrations in a H₂S/H₂ mixture, on the order of parts per million. This reflects the strong chemisorption of atomic S to Ni(111) relative to the weak thermodynamic driving force for desorption into H₂S. In contrast, Fig. 4(b) shows that co-adsorbed oxygen can reduce sulfur coverages on Ni(111) via site competition under sufficiently oxygen-rich conditions (μ^R_O → −1 eV); although this does not occur under any realistic oxygen partial pressures at 600 K. These results suggest that a high temperature is essential for oxygen-assisted catalyst regeneration via site competition between co-adsorbed S and O.

To investigate the entropic contributions to catalyst regeneration via oxidation of S → SO, we validated six GCMC-predicted adlayers for μ^R_S = −1 eV, μ^R_O = −1 eV, −0.7 eV and −0.5 eV, and T = 600 K and 1200 K, using geometry optimisation simulations with the fine-tuned MACE model (trained on both low coverage and high coverage DFT relaxations). The mean and standard deviation of the RMSD of adsorbate atomic displacements is shown in Fig. 4(d), where the MACE relaxation trajectories do not lead to S oxidation. In all cases in Fig. 4(d), the differences in the GCMC-predicted and MACE-optimised adlayer structures are driven by surface diffusion of some adsorbed S atoms to nearest neighbour sites without any S oxidation to SO or SO₂, whilst the RMSD in atomic positions is consistently lower for adsorbed O than S (discussed in the SI section S3). The results suggest that combinations of μ^R_S and μ^R_O that lead to higher coverages and intermixing of S and O, illustrated by the dark blue regions in Fig. 4(c), create conditions that are necessary but not sufficient alone for SO formation and that thermal activation is essential for SO formation irrespective of the degree of S and O co-adsorption. As a result, the use of metal oxide support materials with a large oxygen buffering capacity can aid the regeneration of S-poisoned catalysts at high temperature, where the formation and desorption of SO and SO₂ is feasible. However, tuning the support oxygen buffering capacity is unlikely to improve the sulfur tolerance of low temperature catalysts, which requires modification of the Ni catalyst to reduce the high affinity of S, O, SO and SO₂. These findings are consistent with the kinetic modelling of S oxidation on Ni(111) by Galea et al., who combined DFT simulations with TPD experiments to investigate the removal of adsorbed S atoms using gas-phase O₂.³¹ Their TPD results showed no SO₂ formation at temperatures below 600 K for surfaces with low S coverage, indicating that direct oxidation of S atoms is not thermally accessible at these conditions. Instead, S removal was only observed above 600 K and at sufficiently high O₂ exposures, to facilitate O-assisted S diffusion and oxidation. Their DFT calculations similarly demonstrated a high activation barrier (>1 eV) for SO formation from isolated S and O atoms on Ni(111).

3.3 Reversible vs. irreversible catalyst deactivation

The results in section 3.2 can be used to rationalise the outcomes of experimental MSR activity testing of fresh and H₂S-poisoned Ni nanoparticle catalysts in Fig. 5, which shows methane conversion as a function of the reaction temperature. For both H₂S-poisoned Ni/TiO₂ and H₂S-poisoned Ni/CeO₂, catalyst regeneration and partial restoration of activity (to ∼80% and ∼50% of that of fresh Ni/TiO₂ and Ni/CeO₂, respectively) is achieved upon increasing the temperature beyond 973 K. Although H₂S-poisoned Ni/TiO₂ is restored to the highest absolute value of catalytic activity in Fig. 5(a), ICP analysis indicates a total uptake of H₂S during room temperature saturation of 0.11 weight percentage of sulfur (%_{S wt}), which is an order of magnitude lower than that of Ni/γ-Al₂O₃ (2.14%_{S wt}) and Ni/CeO₂ (2.53%_{S wt}). The reduced sulfur loading on Ni/TiO₂ likely stems from the reduced dispersion of Ni in the experimentally prepared catalyst, as evident by the SEM imaging in the SI section S4, which is consistent with the much larger XRD-determined NiO crystallite size of 17.9 nm on TiO₂ vs. 12.1 nm on CeO₂. As a result, Fig. 5(a) shows that H₂S-poisoned Ni/CeO₂ is restored to a substantially greater catalytic activity than H₂S-poisoned Ni/TiO₂, relative to its sulfur-content, which is in line with our DFT+U calculated oxygen vacancy formation energies of 3.44 eV for CeO₂ and 5.35 eV for TiO₂, i.e., oxygen from the CeO₂ lattice facilitates S oxidation. Both values are much lower than the DFT-calculated oxygen vacancy formation energy of 7.00 eV for γ-Al₂O₃, indicating support oxygen buffering may drive the enhanced sulfur resistance of Ni/CeO₂, although not in a manner to reduce the temperature required for catalyst regeneration, as discussed in section 3.2.

Fig. 5 (a) Temperature profile for MSR activity testing of fresh and H₂S-poisoned Ni catalysts supported on (b) γ-Al₂O₃, (c) TiO₂ and (d) CeO₂. The reduction in temperature from 1073 K to 873 K after t = 6 hours was only performed for the H₂S-poisoned catalysts. All fresh catalysts were subject to an additional pre-reduction in H₂ at 923 K, prior to t = 0 hours. The H₂S-poisoned catalysts contain 0.11%_{S wt}, 2.14%_{S wt} and 2.53%_{S wt} for Ni/TiO₂, Ni/γ-Al₂O₃ and Ni/CeO₂, respectively, as determined using ICP. As such, Ni/CeO₂ is regenerated substantially more than Ni/TiO₂ relative to its sulfur content.

The H₂S-poisoned Ni/γ-Al₂O₃ catalyst was found to deactivate irreversibly in Fig. 5(b), with no restoration of catalytic activity upon increasing temperature. Given the measured activity of the fresh Ni/γ-Al₂O₃ catalyst, which is subject to a pre-reduction in H₂ at 923 K, the irreversible deactivation of H₂S-poisoned Ni/γ-Al₂O₃ is likely due to the variation in the Ni oxidation state with respect to the reducibility of the reaction environment. The observed irreversible catalyst deactivation is consistent with the experimentally reported in situ transformation of Ni/γ-Al₂O₃ to spinel-type NiAl₂O₄, i.e., switching the Ni oxidation state from Ni⁰ in Ni²⁺ on the surface and in the bulk, which is inactive for MSR.^98–100 The suppression of Ni⁰ when Ni/γ-Al₂O₃ is exposed to oxidising atmospheres, e.g., when exposed to air in ambient conditions before characterisation, is further supported by the Ni 2p_3/2 XPS spectra in Fig. 6(a), where the Ni surface speciation on the different supports is distinctly different at ∼853 eV, which corresponds to Ni⁰, whilst being similar at ∼856 eV, which corresponds to Ni²⁺.¹⁰¹ Given that the relative intensity of the peak at ∼853 eV is lowest for H₂S-poisoned Ni/γ-Al₂O₃, this suggests that γ-Al₂O₃ suppresses the formation of Ni⁰ in oxidising conditions.

Fig. 6 Normalised XPS spectra for (a) Ni 2p_3/2 and (b) S 2p for the three H₂S-poisoned Ni catalysts following room temperature saturation with H₂S (before MSR activity testing). (c) Substitutional defect energies for Ni^×_Al in bulk γ-Al₂O₃ (DFT), Ni^×_Ti in bulk TiO₂ (DFT+U) and Ni^×_Ce in bulk CeO₂ (DFT+U), calculated using the mBEEF exchange correlation density functional and Hubbard parameters detailed in the SI section S1.2. The defect energies are plot alongside the corresponding Ni 3d e_g orbitals, including both 3d_z² and 3d_x²−y² orbitals. Large differences between 3d_z² and 3d_x²−y² orbital occupancies are reportedly characteristic of systems with stabilising Jahn–Teller distortions.^103,104

To investigate the driving force for irreversible catalyst deactivation further, we calculated the energetics of substitutional defect formation in the support materials using DFT and DFT+U, as outlined in section 2.1.2. As shown in Fig. 6(c), the substitutional defect energy for Ni^×_Al in γ-Al₂O₃ is calculated as 6.08 eV, which is lower than Ni^×_Ti in TiO₂ (6.67 eV) and Ni^×_Ce in CeO₂ (13.61 eV), supporting a hypothesis that the deactivating phase transformation is more favourable for Ni/γ-Al₂O₃, whereas Ni/TiO₂ and Ni/CeO₂ are more resistant to forming bulk solid solutions. Fig. 6(c) further shows that the increasing defect energies from Ni^×_Al to Ni^×_Ce correlate inversely with the polarisation of the Ni 3d e_g orbitals, comprised of the 3d_z² and 3d_x²−y² orbitals that align along the metal–oxygen bonds,¹⁰² which is characteristic of complex oxides containing divalent ions such as Ni²⁺ resulting in stabilisation via Jahn–Teller distortions that break the system symmetry.^103,104 These results indicate an energetic favourability for the initial stages of phase transformation in γ-Al₂O₃, in agreement with the DFT+U-parameterised Monte Carlo study of Elias et al., who concluded the NiAl₂O₄ can be more stable than NiO and γ-Al₂O₃ in Ni-rich conditions at high temperatures.⁴⁹ The predicted insolubility of Ni in CeO₂ is in contrast with literature-reported defect energies of ∼2–3 eV using DFT+U in a planewave basis.^105,106 Whilst the two sets of results are not directly comparable due to differences in the employed Hubbard projectors, our results align with previous work that shows self-consistent DFT+U in a NAO framework can successfully rationalise experimentally observed defect chemistry in TMOs, e.g., the varying oxidation states of Nb and W dopants in different TiO₂ polymorphs^64,65 and the energetics of Mg doping in LiCoO₂,⁶⁵ the results for which can vary ambiguously in the plane-wave DFT+U literature.^107–110 The large defect energy for Ni^×_Ce is confirmed as not an artifact of our chosen DFT+U parameters by repetition of the calculation using standalone DFT, which yields a defect formation energy of 13.81 eV.

3.4 Sulfur speciation and the role of water

To gain further insights into the mechanisms that drive sulfur removal from the H₂S-poisoned catalysts, TPD-MS was performed in N₂ to track the signals for H₂O, SO and SO₂, which correspond to measurements from mass spectrometry (Fig. 7). For H₂S-poisoned Ni/CeO₂, sulfur removal occurs partially in a low temperature regime (between 423–573 K) and also a high temperature regime (beyond 973 K), which can be attributed to lattice and surface oxygen, respectively, based on the thermogravimetric analysis of Zhu et al., who studied pure and Ni-doped CeO₂ nanorods showing surface oxygen release between 423–593 K and lattice oxygen release between 593–1073 K.¹¹¹ Liu et al. similarly used TPD-MS to investigate SO₂ release from H₂S-poisoned CeO₂, concluding that peaks between 473–673 K corresponded to the formation of SO₂ that could react with lattice oxygen above 673 K to form Ce(SO₄)₂, and then this decomposes back to SO₂ at 873 K.¹¹² The role of oxygen in facilitating sulfur removal was further supported by observations that SO₂ TPD-MS signals were greatest when the catalyst was pretreated in O₂, compared to inert Ar or reducing H₂.¹¹²

Fig. 7 Temperature-programmed-desorption-mass spectrometry (TPS-MS) spectra obtained using a fixed temperature ramp of 10 K min⁻¹ from room temperature to 1223 K in N₂ for (a) H₂O (mass = 18 g mol⁻¹) release from H₂S-poisoned γ-Al₂O₃, TiO₂ and CeO₂, (b) SO (mass = 48 g mol⁻¹) release from H₂S-poisoned γ-Al₂O₃ and CeO₂, and (c) SO₂ (mass = 64 g mol⁻¹) release from H₂S-poisoned γ-Al₂O₃ and CeO₂. The TPD-MS spectra for SO and SO₂ release from H₂S-poisoned Ni/TiO₂ were negligible (due to the lower H₂S loading as discussed in section 3.3) and therefore are not shown. TPD-MS signals for H₂S (mass = 34 g mol⁻¹) release from all catalysts were negligible, indicating H₂S desorption and/or dissociation before analysis. These catalysts were not subject to a pre-reduction in H₂ at 923 K, as discussed for the fresh catalysts in section 3.3.

Fig. 7(b) and (c) show a greater TPD-MS signal for SO and SO₂ release from H₂S-poisoned Ni/γ-Al₂O₃ at low temperatures than H₂S-poisoned Ni/CeO₂. We attribute this difference to the increased formation of surface Ni_xAl_1−xO₂ solid solutions, based on our calculated bulk defect formation energies in Section 3.3 and the H₂ temperature programmed reduction (TPR) study of Shan et al., which correlated the bimodal distribution at low temperatures in Fig. 7(b) and (c) to the existence of both Ni⁰ and Ni²⁺ on the catalyst surface.¹¹³ To rationalise the differences between the high temperature SO and SO₂ desorption behaviour from Ni/γ-Al₂O₃ and Ni/CeO₂ in Fig. 7(b) and (c), the S 2p XPS spectra in Fig. 6(b) is considered, where sulfates and sulfides (NiS) were identified as the peaks at ∼169 eV and ∼162 eV, respectively. Around 85% of all sulfur species in the three H₂S-poisoned catalysts were quantified to be sulfates using curve fitting of the S 2p XPS spectra in Fig. 6(b).

The temperature-dependent oxidation (reduction) of SO₂ to (from) sulfates is hypothesised to drive the differences in the TPD-MS spectra of Ni/γ-Al₂O₃ and Ni/CeO₂ in Fig. 7(b) and (c). The hypothesis is supported by the study of Hamzehlouyan et al., who combined TPD and diffuse reflectance infrared Fourier transform spectroscopy (DRIFTS) to investigate SO₂ release from SO₂-poisoned Pt/Al₂O₃ catalysts, concluding that SO₂-TPD peaks at ∼509 K and ∼947 K correspond to the desorption of molecularly adsorbed SO₂ and the dissociation of aluminium sulfate, respectively.¹¹⁴ Furthermore, Smirnov et al. used temperature-resolved XPS to show that water vapour inhibits SO₂ oxidation to sulfates on an Al₂O₃ thin film but enhances sulfate formation on a CeO₂ thin film, due to a Ce³⁺ redox-mediated mechanism of SO₂ oxidation.¹¹⁵ Together with our TPD-MS results in Fig. 7(a), which show orders of magnitude greater water adsorption on Ni/γ-Al₂O₃ than Ni/CeO₂ due to the 7× greater surface area, the findings of Hamzehlouyan et al. and Smirnov et al. support the hypothesis that SO and SO₂ desorb at lower temperatures from Ni/γ-Al₂O₃ as water vapour inhibits the formation and retention of thermally stable sulfates.

4 Conclusions

Understanding the atomic level mechanisms that govern the sulfur tolerance of Ni-based catalysts is essential for designing next-generation catalysts for industrial H₂ production via MSR and low-temperature processes from renewable feedstocks. In this study, a combined computational and experimental approach is adopted to investigate the enhanced sulfur tolerance of Ni nanoparticles on reducible metal oxide supports, with the aim of uncovering strategies for future catalyst optimisation. Combining DFT, GCMC and a fine-tuned MACE MLIP, we show that a high oxygen chemical potential provided via support oxygen buffering is not sufficient alone for the removal of adsorbed S from Ni(111), with thermal activation being essential. The results support experimental MSR activity tests showing that the catalytic activity of Ni supported on reducible CeO₂ can be readily restored from a poisoned state at high temperatures, compared to Ni supported on less reducible TiO₂ and γ-Al₂O₃. The results are further validated using DFT+U computed oxygen vacancy formation energies for the bulk support materials, which show the ease of oxygen vacancy formation in the order CeO₂ > TiO₂ > γ-Al₂O₃. The MSR activity testing also indicates the critical role of phase transformations into catalytically inactive phases, which is widely reported to occur for Ni/γ-Al₂O₃, and that agrees with our DFT+U computed defect energies for substitutional Ni doping, which indicate the initial stages of bulk phase transformation are more favourable in the order γ-Al₂O₃ > TiO₂ > CeO₂. TPD-MS and XPS highlight the critical role of water in the formation of thermally stable sulfate species that can increase the temperatures required for catalyst regeneration.

Overall, the combined computational and experimental investigation points to three critical aspects for the rational design of metal oxide support materials for sulfur tolerant catalysts: (1) the feasibility of bulk oxygen vacancy formation in the support; (2) the resistance of the bulk support to phase transformations into catalytically inactive solid solutions; and (3) the support- and temperature-dependent surface chemistry of SO₂ to sulfates. The integration of ab initio computational modelling, statistical sampling and machine learning further demonstrates the importance of advanced workflows for studying complex catalytic materials in a manner that faithfully bridges theory and experiment.

Author contributions

All authors contributed to the conceptualisation of the project as well as software choices and computational and/or experimental method development. AC performed the electronic structure calculations presented in this work and the analysis of results. PS performed the Monte Carlo simulations presented in this work. AC and AH performed the MACE MLIP fine-tuning and inferencing. CH performed the experimental characterisation presented in this work. All authors contributed to the preparation of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

All input/output files for electronic structure calculations, Monte Carlo sampling and MACE fine-tuning are available open-source in the GitHub repository https://github.com/amitmc1/GCMC-Adlayers and as a supplementary dataset on Figshare at the DOI: https://doi.org/10.6084/m9.figshare.29562377.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5cy01279a.

Acknowledgements

We thank Andrew Steele from Johnson Matthey for assisting with the experimental setup and characterisation, as well as Gregory Goodlet, Riho Green and Jason Raymond from the Advanced Characterisation Department at Johnson Matthey Technology Centre for performing the SEM, XPS and ICP analysis, respectively. We thank Harald Oberhofer, Matthias Kick and Maximilian Brand for valuable scientific discussions regarding the implementation of DFT+U in FHI-aims. We acknowledge funding by the Prosperity Partnership project Sustainable Catalysis for Clean Growth, funded by the UK Engineering and Physical Sciences Research Council (EPSRC), bp through the bp International Centre for Advanced Materials (bp-ICAM) and Johnson Matthey plc in collaboration with Cardiff University and The University of Manchester (EPSRC grant number EP/V056565/1). A. C. acknowledges funding by the Collaborative Computational Project Number 5 (CCP5) as part of a Postgraduate Industrial Secondment, which facilitated computational and experimental collaboration with partners at Johnson Matthey (EPSRC grant number EP/V028537/1). A. J. L. acknowledges funding by the UKRI Future Leaders Fellowship program (MR/T018372/1 and MR/Y034279/1). The authors acknowledge computational resources and support from the Supercomputing Wales project, which is part-funded by the European Regional Development Fund (ERDF) via the Welsh Government; and the UK National Supercomputing Services ARCHER and ARCHER2, accessed via membership of the Materials Chemistry Consortium, which is funded by Engineering and Physical Sciences Research Council (EP/L000202/1, EP/R029431/1, and EP/T022213/1).

Notes and references

1. J. Moore, J. Durham, A. Eijk, E. Karakas, R. Kurz, J. Lesak, M. McBain, P. McCalley, L. Moroz, Z. Mohamed, B. Pettinato, G. Phillippi, H. Watanabe and B. Williams, Machinery and Energy Systems for the Hydrogen Economy, Elsevier, 2022, ch. 9 – Compressors and expanders, pp. 333–424.
2. W. S. Jablonski, S. M. Villano and A. M. Dean, Appl. Catal., A, 2015, 502, 399–409.
3. B. Hua, M. Li, Y.-F. Sun, Y.-Q. Zhang, N. Yan, J. Chen, J. Li, T. Etsell, P. Sarkar and J.-L. Luo, J. Mater. Chem. A, 2016, 4, 4603–4609.
4. N. Schiaroli, M. Volanti, A. Crimaldi, F. Passarini, A. Vaccari, G. Fornasari, S. Copelli, F. Florit and C. Lucarelli, Energy Fuels, 2021, 35, 4224–4236.
5. J. A. Rodriguez and J. Hrbek, Acc. Chem. Res., 1999, 32, 719–728.
6. C. Brady, J. Pan and B. Xu, Catal. Sci. Technol., 2020, 10, 8429–8436.
7. C. Xie, Y. Chen, Y. Li, X. Wang and C. Song, Appl. Catal., A, 2010, 390, 210–218.
8. S. L. Lakhapatri and M. A. Abraham, Appl. Catal., A, 2009, 364, 113–121.
9. A. Cho, B. Hwang and J. W. Han, Catal. Sci. Technol., 2020, 10, 4544–4552.
10. J. Xu, J. Harmer, G. Li, T. Chapman, P. Collier, S. Longworth and S. C. Tsang, Chem. Commun., 2010, 46, 1887–1889.
11. G. Sun, K. Hidajat, X. Wu and S. Kawi, Appl. Catal., B, 2008, 81, 303–312.
12. U. Oemar, K. Hidajat and S. Kawi, Int. J. Hydrogen Energy, 2015, 40, 12227–12238.
13. Z. Li and K. Sibudjing, ChemCatChem, 2018, 10, 2994–3001.
14. D. Guo, Y. Lu, Y. Ruan, Y. Zhao, Y. Zhao, S. Wang and X. Ma, Appl. Catal., B, 2020, 277, 119278.
15. L. Pino, C. Italiano, A. Vita, M. Laganá and V. Recupero, Appl. Catal., B, 2017, 218, 779–792.
16. H. Wang, X. Dong, T. Zhao, H. Yu and M. Li, Appl. Catal., B, 2019, 245, 302–313.
17. G.-R. Hong, K.-J. Kim, S.-Y. Ahn, B.-J. Kim, H.-R. Park, Y.-L. Lee, S. S. Lee, Y. Jeon and H.-S. Roh, Catalysts, 2022, 12, 1670.
18. Y.-L. Lee, K.-J. Kim, G.-R. Hong, S.-Y. Ahn, B.-J. Kim, H.-R. Park, S.-J. Yun, J. W. Bae, B.-H. Jeon and H.-S. Roh, ACS Sustainable Chem. Eng., 2021, 9, 15287–15293.
19. S. d. S. Eduardo, J. P. Mendonca, P. N. Romano, J. M. A. R. de Almeida, G. Machado and M. A. S. Garcia, Hydrogen, 2023, 4, 493–522.
20. M. A. Ocsachoque, J. I. Eugenio Russman, B. Irigoyen, D. Gazzoli and M. G. González, Mater. Chem. Phys., 2016, 172, 69–76.
21. L. Li, C. Howard, D. L. King, M. Gerber, R. Dagle and D. Stevens, Ind. Eng. Chem. Res., 2010, 49, 10144–10148.
22. P. Wachter, C. Gaber, J. Raic, M. Demuth and C. Hochenauer, Int. J. Hydrogen Energy, 2021, 46, 3437–3452.
23. A. de Lucas-Consuegra, A. Caravaca, P. Martínez, J. Endrino, F. Dorado and J. Valverde, J. Catal., 2010, 274, 251–258.
24. S. Zha, Z. Cheng and M. Liu, J. Electrochem. Soc., 2006, 154, B201.
25. T. Morooka, T. Shishido, R. Devivaraprasad, G. Elumalai, M. Aoki, T. Shirasawa, T. Nakanishi, A. Ishikawa, T. Kondo and T. Masuda, J. Phys. Chem. C, 2024, 128, 16426–16436.
26. W. Chen, T. Li, L. Peng, G. Shen, Z. Jiang, B. Huang and H. Zuo, Comput. Theor. Chem., 2024, 1231, 114443.
27. M. Zhang, Z. Fu and Y. Yu, Appl. Surf. Sci., 2019, 473, 657–667.
28. J.-H. Wang and M. Liu, Electrochem. Commun., 2007, 9, 2212–2217.
29. C. R. Bernard Rodríguez and J. A. Santana, J. Chem. Phys., 2018, 149, 204701.
30. C.-H. Yeh and J.-J. Ho, ChemPhysChem, 2012, 13, 3194–3203.
31. N. M. Galea, J. M. Lo and T. Ziegler, J. Catal., 2009, 263, 380–389.
32. C. Schwennicke and H. Pfnür, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 56, 10558–10566.
33. C. Lazo and F. J. Keil, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 245418.
34. S. S. Akimenko, G. D. Anisimova, A. I. Fadeeva, V. F. Fefelov, V. A. Gorbunov, T. R. Kayumova, A. V. Myshlyavtsev, M. D. Myshlyavtseva and P. V. Stishenko, J. Comput. Chem., 2020, 41, 2084–2097.
35. L. Gai, Y. K. Shin, M. Raju, A. C. T. van Duin and S. Raman, J. Phys. Chem. C, 2016, 120, 9780–9793.
36. T. Demeyere, T. Ellaby, M. Sarwar, D. Thompsett and C.-K. Skylaris, ACS Catal., 2025, 15, 5674–5682.
37. T. Demeyere, H. U. Islam, T. Ellaby, M. Sarwar, D. Thompsett and C.-K. Skylaris, Phys. Chem. Chem. Phys., 2025, 27, 10011–10022.
38. K. T. Schütt, F. Arbabzadah, S. Chmiela, K. R. Muller and A. Tkatchenko, Nat. Commun., 2017, 8, 13890.
39. K. Schutt, O. Unke and M. Gastegger, Int. Conf. Mach. Learn., 2021, 9377–9388.
40. C. Chen and S. P. Ong, Nat. Comput. Sci., 2022, 2, 718–728.
41. B. Deng, P. Zhong, K. Jun, J. Riebesell, K. Han, C. J. Bartel and G. Ceder, Nat. Comput. Sci., 2023, 3, 192–202.
42. I. Batatia, D. P. Kovacs, G. Simm, C. Ortner and G. Csányi, Adv. Neural Inf. Process. Syst., 2022, 35, 11423–11436.
43. A. P. Bartók, M. C. Payne, R. Kondor and G. Csányi, Phys. Rev. Lett., 2010, 104, 136403.
44. W. G. Stark, J. Westermayr, O. A. Douglas-Gallardo, J. Gardner, S. Habershon and R. J. Maurer, J. Phys. Chem. C, 2023, 127, 24168–24182.
45. H. Jung, L. Sauerland and S. Stocker, et al., npj Comput. Mater., 2023, 9, 114.
46. I. Batatia, P. Benner, Y. Chiang, A. M. Elena, D. P. Kovács, J. Riebesell, X. R. Advincula, M. Asta, M. Avaylon and W. J. Baldwin, J. Chem. Phys., 2025, 163, 184110.
47. R. Drautz, Phys. Rev. B, 2019, 99, 014104.
48. D. P. Kovács, I. Batatia, E. S. Arany and G. Csányi, J. Chem. Phys., 2023, 159, 044118.
49. I. Elias, A. Soon, J. Huang, B. S. Haynes and A. Montoya, Phys. Chem. Chem. Phys., 2019, 21, 25952–25961.
50. P. Littlewood, S. Liu, E. Weitz, T. J. Marks and P. C. Stair, Catal. Today, 2020, 343, 18–25.
51. F. F. Tao, J. J. Shan, L. Nguyen, Z. Wang, S. Zhang, L. Zhang, Z. Wu, W. Huang, S. Zeng and P. Hu, Nat. Commun., 2015, 6, 7798.
52. L. Yu, M. Song, P. T. Williams and Y. Wei, Ind. Eng. Chem. Res., 2019, 58, 11770–11778.
53. F. Chen, C. Wu, L. Dong, A. Vassallo, P. T. Williams and J. Huang, Appl. Catal., B, 2016, 183, 168–175.
54. D. Li, Y. Li, X. Liu, Y. Guo, C.-W. Pao, J.-L. Chen, Y. Hu and Y. Wang, ACS Catal., 2019, 9, 9671–9682.
55. A. Jamsaz, N. Pham-Ngoc, M. Wang, D. H. Jeong and E. W. Shin, Chem. Eng. J., 2024, 500, 156932.
56. S. Li, J. Li, Z. He, Y. Sheng and W. Liu, Catal. Sci. Technol., 2024, 14, 5864–5873.
57. S. T. Misture, K. M. McDevitt, K. C. Glass, D. D. Edwards, J. Y. Howe, K. D. Rector, H. He and S. C. Vogel, Catal. Sci. Technol., 2015, 5, 4565–4574.
58. Z.-K. Han, W. Liu and Y. Gao, JACS Au, 2025, 5, 1549–1569.
59. T. Zhang, P. Zheng, J. Gao, X. Liu, Y. Ji, J. Tian, Y. Zou, Z. Sun, Q. Hu and G. Chen, et al., Nat. Commun., 2024, 15, 6827.
60. N. L. Nguyen, N. Colonna, A. Ferretti and N. Marzari, Phys. Rev. X, 2018, 8, 021051.
61. J. P. Perdew and A. Zunger, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 23, 5048–5079.
62. M. Reticcioli, U. Diebold and C. Franchini, J. Phys. Condens. Matter, 2022, 34, 204006.
63. M. Kick, K. Reuter and H. Oberhofer, J. Chem. Theory Comput., 2019, 15, 1705–1718.
64. A. Chaudhari, A. J. Logsdail and A. Folli, J. Phys. Chem. C, 2025, 129, 15453–15461.
65. A. Chaudhari, K. Agrawal and A. J. Logsdail, Digital Discovery, 2025, 4, 3701–3727.
66. V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter and M. Scheffler, Comput. Phys. Commun., 2009, 180, 2175–2196.
67. A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng and K. W. Jacobsen, J. Phys. Condens. Matter, 2017, 29, 273002.
68. O. Lamiel-Garcia, K. C. Ko, J. Y. Lee, S. T. Bromley and F. Illas, J. Chem. Theory Comput., 2017, 13, 1785–1793.
69. F. Passek and M. Donath, Phys. Rev. Lett., 1993, 71, 2122–2125.
70. B. Legendre and M. Sghaier, J. Therm. Anal. Calorim., 2011, 105, 141–143.
71. J. Wellendorff, K. T. Lundgaard, K. W. Jacobsen and T. Bligaard, J. Chem. Phys., 2014, 140, 144107.
72. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
73. S. Lehtola, C. Steigemann, M. J. Oliveira and M. A. Marques, SoftwareX, 2018, 7, 1–5.
74. F. Birch, Phys. Rev., 1947, 71, 809–824.
75. C. G. Broyden, IMA J. Appl. Math., 1970, 6, 76–90.
76. R. Fletcher, Comput. J., 1970, 13, 317–322.
77. D. F. Shanno, Math. Comput., 1970, 24, 647–656.
78. D. Goldfarb, Math. Comput., 1970, 24, 23–26.
79. S. S. Akimenko, V. A. Gorbunov, A. V. Myshlyavtsev and P. V. Stishenko, Phys. Rev. E, 2016, 93, 062804.
80. D. Landau and K. Binder, A guide to Monte Carlo simulations in statistical physics, Cambridge University Press, 2021.
81. D. J. Earl and M. W. Deem, Phys. Chem. Chem. Phys., 2005, 7, 3910–3916.
82. H. Huo and M. Rupp, Mach. Learn.: Sci. Technol., 2022, 3, 045017.
83. L. Himanen, M. O. J. Jäger, E. V. Morooka, F. Federici Canova, Y. S. Ranawat, D. Z. Gao, P. Rinke and A. S. Foster, Comput. Phys. Commun., 2020, 247, 106949.
84. J. Laakso, L. Himanen, H. Homm, E. V. Morooka, M. O. J. Jäger, M. Todorović and P. Rinke, J. Chem. Phys., 2023, 158, 234802.
85. F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot and E. Duchesnay, J. Mach. Learn. Res., 2011, 12, 2825–2830.
86. A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder and K. A. Persson, APL Mater., 2013, 1, 011002.
87. B. Deng, Materials Project Trajectory (MPtrj) Dataset, 2023, https://figshare.com/articles/dataset/Materials_Project_Trjectory_MPtrj_Dataset/23713842.
88. M. M. Ghahremanpour, P. J. van Maaren and D. van der Spoel, Sci. Data, 2018, 5, 180062.
89. K. D. B. J. Adam, et al., arXiv, 2014, preprint, arXiv:1412.6980,  DOI:10.48550/arXiv.1412.6980.
90. P. Promhuad, B. Sawatmongkhon, K. Theinnoi, T. Wongchang, N. Chollacoop, E. Sukjit, S. Tunmee and A. Tsolakis, ACS Omega, 2024, 9, 19282–19294.
91. J. R. Sietsma, A. Jos van Dillen, P. E. de Jongh and K. P. de Jong, in Scientific Bases for the Preparation of Heterogeneous Catalysts, ed. E. Gaigneaux, M. Devillers, D. De Vos, S. Hermans, P. Jacobs, J. Martens and P. Ruiz, Elsevier, 2006, vol. 162 of Studies in Surface Science and Catalysis, pp. 95–102.
92. M. A. Hefnawy, S. A. Fadlallah, R. M. El-Sherif and S. S. Medany, J. Mol. Graphics Modell., 2023, 118, 108343.
93. Y. Bai, D. Kirvassilis, L. Xu and M. Mavrikakis, Surf. Sci., 2019, 679, 240–253.
94. V. Alexandrov, M. L. Sushko, D. K. Schreiber, S. M. Bruemmer and K. M. Rosso, Corros. Sci., 2016, 113, 26–30.
95. L. Liu, C. Zhang, W. Wang, G. Li and B. Zhu, Molecules, 2023, 28, 6739.
96. T. Yokoyama, S. Terada, A. Imanishi, Y. Kitajima, N. Kosugi and T. Ohta, J. Electron Spectrosc. Relat. Phenom., 1996, 80, 161–164.
97. N. K. Das and W. A. Saidi, J. Chem. Phys., 2017, 146, 154701.
98. J. Zieliński, J. Mol. Catal., 1993, 83, 197–206.
99. S. Maluf and E. Assaf, Fuel, 2009, 88, 1547–1553.
100. Y. Guo, T. P. Tran, L. Zhou, Q. Zhang and H. Kameyama, J. Chem. Eng. Jpn., 2007, 40, 1221–1228.
101. D. Li, Q. Zhu, Z. Bao, L. Jin and H. Hu, Fuel, 2024, 363, 131045.
102. J. P. Allen and G. W. Watson, Phys. Chem. Chem. Phys., 2014, 16, 21016–21031.
103. X. Wang, D. Santos-Carballal and N. H. de Leeuw, J. Chem. Phys., 2024, 160, 154713.
104. R. Prasad, R. Benedek and M. M. Thackeray, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 134111.
105. Z. Chafi, N. Keghouche and C. Minot, Surf. Sci., 2007, 601, 2323–2329.
106. E. F. de Souza, C. A. Chagas, R. L. Manfro, M. M. V. M. Souza, R. Bicca de Alencastro and M. Schmal, RSC Adv., 2016, 6, 5057–5067.
107. J. A. Santana, J. Kim, P. R. C. Kent and F. A. Reboredo, J. Chem. Phys., 2014, 141, 164706.
108. K. K. Ghuman and C. V. Singh, J. Phys.: Condens. Matter, 2013, 25, 085501.
109. B. J. Morgan, D. O. Scanlon and G. W. Watson, J. Mater. Chem., 2009, 19, 5175–5178.
110. A. Raghav, K. Hongo, R. Maezono and E. Panda, Comput. Mater. Sci., 2022, 214, 111714.
111. Y. Zhu, W. Wang, G. Chen, H. Li, Y. Zhang, C. Liu, H. Wang, P. Cheng, C. Chen and G. Seong, Crystals, 2024, 14, 746.
112. B. Liu, H. Xu and Z. Zhang, Chin. J. Catal., 2012, 33, 1631–1635.
113. W. Shan, M. Luo, P. Ying, W. Shen and C. Li, Appl. Catal., A, 2003, 246, 1–9.
114. T. Hamzehlouyan, C. Sampara, J. Li, A. Kumar and W. Epling, Top. Catal., 2016, 59, 1028–1032.
115. M. Smirnov, A. Kalinkin, A. Pashis, A. Sorokin, A. Noskov, V. Bukhtiyarov, K. Kharas and M. Rodkin, Kinet. Catal., 2003, 44, 575–583.

---