Prediction of morphological changes of catalyst materials under reaction conditions by combined ab initio thermodynamics and microkinetic modelling

Abstract

In this article, we couple microkinetic modelling, ab initio thermodynamics and Wulff–Kaishew construction to describe the structural variation of catalyst materials as a function of the chemical potential in the reactor. We focus specifically on experiments of catalytic partial oxidation (CPO) of methane on Rh/α-Al₂O₃. We employ a detailed structureless microkinetic model to calculate the profiles of the gaseous species molar fractions along the reactor coordinate and to select the most abundant reaction intermediates (MARIs) populating the catalyst surfaces in different zones of the reactor. Then, we calculate the most stable bulk and surface structures of the catalyst under different conditions of the reaction environment with density functional theory (DFT) calculations and ab initio thermodynamics, considering the presence of the MARIs on the catalyst surface in thermodynamic equilibrium with the partial pressures of their reservoirs in the gas phase surrounding the catalyst. Finally, we exploit the Wulff–Kaishew construction method to estimate the three-dimensional shape of the catalyst nanoparticles and the distribution of the active sites along the reactor coordinate. We find that the catalyst drastically modifies its morphology during CPO reaction by undergoing phase transition, in agreement with spectroscopy studies reported in the literature. The framework is also successfully applied for the analysis and interpretation of chemisorption experiments for catalyst characterization. These results demonstrate the crucial importance of rigorously accounting for the structural effect in microkinetic modeling simulations and pave the way towards the development of structure-dependent microkinetic analysis of catalytic processes.

---

Introduction

Heterogeneous catalysts are dynamic systems in the sense that they can be seen as “smart materials”, which adapt their structure under the specific reaction conditions locally experienced in the reactor.¹ Such “living character” of the catalyst within the reaction environment strongly affects the surface structure, which is a crucial component in determining the observed activity and selectivity in very important metal-supported catalytic processes of industrial relevance (e.g., NH₃ synthesis, Fisher–Tropsch, hydrocarbon reforming).² Considering the complex nature of these phenomena, multiscale analysis based on microkinetic modelling is acknowledged to be the key tool for the interpretation of the experimental evidence for gaining fundamental information on the relation between the structure and the activity of heterogeneous catalysts.^3,4 However, despite the fact that the surface structure has been considered an important factor in catalysis science since the discovery of structure sensitive reactions,^5,6 the effect of the structure of the catalyst on reactivity and selectivity is at present neglected in state-of-the-art microkinetic modelling. As such, these models strongly rely on an abstract and static concept of the “catalyst material”, which is often modelled as a generic free site “*”, without explicitly accounting for the effect of the actual orientation and positioning of the atoms at the surface and their dependence on the operating conditions. All the structural effects are typically accounted for by means of fine-tuning and refinement of the kinetic parameters based on selected experimental data. This approach can lead to a satisfactory description of conversion and selectivity trends and, thus, it is a very valuable tool for the description of the macroscopic kinetic behaviour of a catalytic process. However, the material gap in the modelling of the active site intrinsically precludes the detailed analysis and understanding of the underlying mechanisms at the atomic-scale level.

A catalyst under reaction conditions generally shows numerous types of active sites, each with different activity. For instance, catalyst nanoparticles are characterized by crystal facets, edges, corners and defects, which interact differently with the reaction intermediates, yielding different reactivities. In particular, low-coordinated catalyst atoms are generally more active than high-coordinated ones, and in some conditions they give a predominant contribution to the catalytic process.^2,7 Moreover, adsorption, phase transitions and sintering processes can modify the morphology (e.g., shape and size) of the catalyst nanoparticles and consequently change the nature and the relative amount of their active sites.^8,9

An example of a system in which the catalyst undergoes drastic morphological modifications under reaction conditions is the catalytic partial oxidation (CPO) of CH₄ on a Rh/Al₂O₃ catalyst.^10,11 Grunwaldt et al. performed a 2D-mapping of the Rh oxidation state in a catalyst bed by means of X-ray absorption near edge structure (XANES) spectroscopy and a charged coupled device (CCD) camera,¹² and observed that the catalyst changes its oxidation state as a function of the O₂ chemical potential in the reactor. At the steady state, the reactor turns out to be stratified in two distinct zones: in the first zone, Rh is highly oxidized, whereas, in the second zone, Rh nanoparticles are reduced. The variation in structure observed between the two zones of the catalytic bed is sharp and it is accompanied by a drastic change in the selectivity of the process: when the catalyst is oxidized, total oxidation is observed, whereas, when Rh is in its metallic form, syngas is produced.

The modelling of the crystal habit of the catalyst material under reaction conditions and its effect on the activity easily result in a high level of complexity.¹³ This is related to the necessity of both (i) building microkinetic models not only for one facet, but for several facets (potentially including corners, edges, and defects), and (ii) coupling the prediction of the structure and shape of the nanoparticle to the local variation of the chemical potential in the reactor. In this situation, however, the first-principles description of the whole set of possible phenomena is beyond the limit of complexity accessible to even the most efficient implementation of first-principles calculations.¹³

A first step towards filling the material gap in microkinetic modelling can be achieved by decoupling the calculations of the distribution of the chemical potential in the reactor and its effect on the structure of the catalyst. We first employ a structureless microkinetic model to calculate the distribution of the chemical potential in the reactor. Structureless models allow for a reliable and detailed analysis of the macroscopic kinetic behaviour of the catalytic process. However, they do not provide any information about the catalyst structure and how it is modified by the gas phase. These phenomena are lumped together in effective kinetic parameters of the elementary steps that reflect the fine-tuning of the kinetic parameters with various experimental data. Then, at the specific conditions of the chemical potential calculated with the microkinetic model, we study how the reaction environment induces morphological changes in the catalyst structure. The description of the catalyst structural changes is achieved by using ab initio thermodynamics and Wulff–Kaishew construction.

In particular, here we focus our analysis on the CPO of CH₄ on Rh/α-Al₂O₃, for which well-established microkinetic models are reported in the literature. The structureless microkinetic model by Maestri et al.^14,15 is used for calculating the profiles of the gaseous species along the reactor coordinate and to identify the most abundant reaction intermediates (MARIs) at the catalyst surface. Then, with ab initio thermodynamics and Wulff–Kaishew construction, we calculate how the morphology of the catalyst changes when exposed to the calculated chemical potential of the gaseous species, in agreement with the experimental evidence.¹⁶ In particular, density functional theory (DFT) calculations and ab initio thermodynamics are employed to calculate the most stable bulk and surface structures of catalyst materials under different conditions of the reaction environment. For each crystal facet of the catalyst, the thermodynamically most stable surface structure (e.g. ordered or disordered adsorbate configurations, surface reconstructions) in equilibrium with the surrounding gas is selected. Then, the adhesion energy between the catalyst and the support is calculated and the 3D shape of the nanoparticles along the reactor coordinate is estimated by the Wulff–Kaishew construction method.

Methods

In this section we provide the description of the different methods employed in the analysis of the CH₄ CPO experiments on Rh catalysts.

Microkinetic modelling

The reaction environment inside the chemical reactor for CH₄ CPO experiments on Rh/α-Al₂O₃ is characterized in terms of profiles of gaseous species partial pressures and adsorbate coverages, by numerically solving mass balance equations under isothermal conditions as reported in ref. 17. Details are provided in the ESI† (section 1). The reaction rates are calculated with the detailed, thermodynamically consistent microkinetic model developed by Maestri et al.^14,15 The model is structureless, in the sense that it does not contain any explicit information about the catalyst structure and how it evolves during reaction. It consists of 82 surface reactions between 13 single-site adsorbate species and includes the effect of adsorbate–adsorbate interactions on activation energies of surface kinetics. It was derived using a hierarchical multiscale methodology based on the UBI-QEP framework and it was refined with first principles calculations to correctly account for the relevant pathways involved in water-gas shift (WGS) and reverse water-gas shift reactions.¹⁵ The resulting scheme can describe several reaction systems on a Rh/α-Al₂O₃ catalyst: CH₄ pyrolysis and oxidation, steam reforming, H₂ and CO-rich oxidation, WGS and reverse WGS.^18,19

Electronic-structure calculations

Electronic-structure calculations are performed using density-functional theory (DFT), as implemented in the Quantum Espresso²⁰ suite of codes, with PBE (Perdew, Burke, and Ernzerhof²¹) GGA ultrasoft pseudopotentials and a plane wave basis set. Rh oxide systems are described using a DFT+U approach, to correct electron self-interaction, which can be relevant for these systems.²² The value of the U parameter is selected by testing different values and comparing the band gap of the relaxed bulk Rh₂O₃ with its experimental value. A value of U of 3.5 eV is applied to the Rh atoms in order to match the experimental value of the band gap of Rh₂O₃ (1.20 eV (ref. 23)). Metallic Rh surface structures are instead described taking into account the contributions of van der Waals interactions (with the vdW-DF2 functional^24–27), which can play an important role in the description of CO* adsorption on Rh surfaces.²⁸ Further considerations on the selection of the U parameter and the employed functionals are reported in the ESI† (section 2).

Periodic surface slabs with inversion symmetry are used to characterize the catalyst surfaces in terms of specific surface free energy (per unit of area) and binding energies of adsorbates, which are positioned on both sides of the slabs. The top three and bottom three slab layers are fully relaxed until all force components acting on the atoms are below 2.6 × 10⁻² eV Å⁻¹ and the difference in energy between two calculations is lower than 1.36 × 10⁻³ eV. Convergence tests have been performed with respect to the specific surface energy, with a threshold value of 10⁻³ eV Å⁻². The resulting converged parameters are: a kinetic energy cut-off for wave functions of 35 Ry and for charge density and potential of 280 Ry, a grid of 12 × 12 × 12 k-points for bulk Rh and of 6 × 6 × 2 for bulk Rh₂O₃, 6 × 6 × 6 for bulk RhO₂, a slab height of about 12 Å and 10 Å of vacuum. For surface supercells, correspondingly smaller k-points grids are used to ensure an equivalent sampling of reciprocal space. All the parameters have been selected to achieve convergence with the different functionals or calculation set-ups herein employed. Gas phase calculations have been performed in a cubic supercell. The size of the supercell (14 Å) was selected to avoid interactions between the atoms of two neighbouring periodic cells.

The lattice constant of Rh was 3.85 Å (3.96 Å with vdW-DF2), in agreement with the experimental value (3.80 Å).²⁹ The lattice parameters of the hexagonal cell of corundum Rh₂O₃, i.e. the length of the basis vectors and the height of the cell, were 5.17 Å and 13.96 Å, respectively, consistent with the experimental values (5.21 Å and 14.15 Å).³⁰ The lattice constants of rutile RhO₂ were 4.51 Å and 3.06 Å, in good agreement with the experimental values (4.59 Å and 3.15 Å).³⁰ The geometric parameters of Rh oxides did not change after the application of Hubbard U corrections. The lattice parameters of α-Al₂O₃ were 4.80 Å and 13.09 Å, in good agreement with the experimental values (4.76 Å and 12.99 Å)³¹ and with previous theoretical studies.³²

Ab initio thermodynamics

The changes in the structure of the catalyst are predicted by means of ab initio thermodynamics³³ calculations. This method consists of selecting the most stable configuration of a system of atoms in equilibrium with a reservoir by comparison of specific thermodynamic quantities – e.g. the surface free energy – calculated by means of first-principle techniques. The main assumption is that the system is in thermodynamic equilibrium with the surrounding environment, which can be characterized by both experimental observation and microkinetic modelling.

Bulk structure. Considering that the catalyst can exist in its metallic and its oxidized form, the most stable catalyst bulk structure in specific conditions is obtained by selecting the structure that shows the lowest grand potential per Rh atom, Ω_RhxOy(T,P,P_O2):where and G^bulk_Rh(T,P) are the Gibbs free energy of the bulk structure Rh_xO_y and the bulk Rh, respectively. μ_O is the chemical potential of an oxygen atom in the catalyst structure, assumed to be equal to the chemical potential of an oxygen atom in the reservoir of the system, typically oxygen molecules in the gas phase.

The chemical potential of gaseous species is calculated in the ideal gas approximation as:

where E^ZPE_A is the zero-point energy of the molecule A, k_B is the Boltzmann constant and Δμ⁰_A(T,P⁰) is its reference chemical potential, calculated using the NIST-JANAF thermochemistry tables.³⁴

The Gibbs free energy of crystalline solid bulk systems is calculated as:

G_solid(T,P) = E^DFT_solid + F^vib_solid(T,P) (3)

where the PV term and the configurational entropy contribution are neglected.³³F^vib_solid is the vibrational free energy of the system, calculated from phonon density of states by density functional perturbation theory.²⁰

Surface structure. Concerning metallic Rh systems, for each crystalline facet that the catalyst can expose to the environment, the most stable surface structure is calculated by selecting the structure that shows the lowest specific surface energy, :where is the Gibbs free energy of the surface structure supercell, G^bulk_Rh is the energy of an atom in Rh bulk, N_Rh and N_A* are the number of Rh atoms and of the adsorbate molecules in the supercell, respectively. μ_A* is the chemical potential of the adsorbate present in the surface structure, assumed to be equal to the chemical potential of its reservoir (in the gas phase) and A_S is the area of the supercell base.

The specific surface energies of the different terminations that Rh₂O₃ can expose in equilibrium with the oxygen in the surrounding environment are calculated as:

where is the Gibbs free energy of the slab that represents the (hkil) surface termination, made of N_Rh Rh atoms and N_O oxygen atoms, is the energy of 2 metal atoms and 3 oxygen atoms in the bulk phase and μ_O is half of the chemical potential of an oxygen molecule in the gas phase.

Neglecting the contribution of the PV term,³³ we calculate the Gibbs free energy of Rh systems in the presence of adsorbates on the surfaces as:

G^surf_Rh/A = E^DFT,surf_Rh/A + F^vib,surf_Rh/A − TS^conf,surf_Rh/A (6)

The vibration contributions to the systems containing adsorbates on the Rh surfaces are calculated in the harmonic oscillator approximation, from a finite difference approximation of the Hessian matrix. The contribution of the vibrational modes of the Rh atoms below the top layer of the slabs is neglected. The configurational entropy of disordered surface structures is calculated in the lattice 2D gas approximation,³⁵ considering the adsorbates as indistinguishable, as a function of the coverages of adsorbates (ϑ_A*):

Support adhesion energy. To account for the influence of the alumina support, we investigate the morphology of α-Al₂O₃(hkil) surfaces as a function of the chemical potential of oxygen and water, which are present in the reaction environment. The specific surface energy of different terminations of Al₂O₃ is calculated as:where is the Gibbs free energy of the termination, made of N_Al aluminium atoms, N_O oxygen atoms and N_H hydrogen atoms, and is the energy of 2 metal atoms and 3 oxygen atoms in the bulk phase.

Sandwich-like symmetric slabs (shown in section 3 of the ESI†) characterized by inversion symmetry are employed to calculate the adhesion energy between the catalyst and the support surfaces. In the case of Rh and Al₂O₃ facets:

where is the Gibbs free energy of the sandwich-like slab, and 2A_S is the total contact area between the two phases, equal to twice the area of the supercell base. The supercell is constructed by ensuring that the periodicities of the two different structures (support and catalyst surfaces) are preserved (details in section 3 of the ESI†).

Wulff–Kaishew construction

The size-independent 3D shape of catalyst nanoparticles under reaction is drawn with the Wulff–Kaishew construction method. The minimum surface energy for a given volume of a crystal polyhedron is achieved if the distances of its faces from one given point are proportional to their specific surface energies:

h_hkl = λγ_hkl (10)

where h_hkl is the distance of the surface (hkl) from the centre of the nanoparticle, γ_hkl is the specific surface free energy of the surface and λ is an arbitrary constant. With ab initio thermodynamics, we calculate the specific surface energies of a set of crystal facets that the solid can expose in equilibrium with the surrounding environment, neglecting the interactions with other facets and with the support. Then, for each considered (hkl) surface, a plane normal to the vector (hkl) is drawn at a distance from the origin proportional to the specific surface energy of the plane. Planes representative of the symmetry of the bulk are also drawn. The presence of the support results in a truncation of the nanoparticle 3D shape in the direction characterized by the lowest catalyst–support adhesion energy. The structure of the catalyst–support interface is obtained by selecting the catalyst facet that shows the lowest adhesion energy, γ_adh (calculated with eqn (9)), when it is in contact with the support. Then, a plane is drawn at a distance from the centre of the Wulff–Kaishew construction proportional to the catalyst–support specific adhesion energy, h_support:³⁶

h_support = λγ_adh (11)

The space that lies inside all these planes defines the equilibrium shape for the nanoparticle in the considered conditions. The thermodynamic shape is size-independent, except in cases of exceptionally large strain effects or counting effects related to edge and corner atoms.³⁷

Results and discussion

In this section, the methods described above are employed to study the morphological changes of a 4% Rh/α-Al₂O₃ catalyst induced by the variation of chemical potential in the reactor. We focus our analysis on the experiments of Donazzi et al. in an annular reactor.³⁸

Microkinetic modelling

Microkinetic modelling is employed to calculate the gas phase composition inside the reactor and to identify the MARIs present at the catalyst surfaces. We employ a steady-state 1D heterogeneous model and the kinetics of the microkinetic model by Maestri et al.¹⁴ The geometric parameters of the reactor and the operating conditions are reported in Table 1.

Table 1 Reactor geometric parameters and operative conditions

Temperature	773 K	GHSV	2 × 10^6 Nl kgcat^−1 h^−1
Pressure	1 atm	Inlet molar fractions
Reactor length	2.2 cm	CH4	0.01
Inner diameter	0.4 cm	O2	0.01
Outer diameter	0.5 cm	N2	0.98

---

Fig. 1a shows the calculated axial profiles of partial pressures of the gaseous species that come into contact with the Rh catalyst in the reactor. The concomitant surface coverages are reported in Fig. 1b.

Fig. 1 Major species partial pressures axial profiles at the catalyst interface (a) and surface coverages (b). Reactor geometric parameters: annular geometry, length: 2.2 cm, inner diameter: 0.4 cm, outer diameter: 0.5 cm. Inlet molar fractions: CH₄ = 0.01, O₂ = 0.01, N₂ = 0.98, GHSV = 2 × 10⁶ Nl kg_cat⁻¹ h⁻¹, T = 773 K and P = 1 atm.

In agreement with previous experimental findings,¹² we observe a sharp change in the reactive environment that divides the reactor into two different zones: at a distance from the reactor inlet of about 1.2 cm, a sharp drop of O* coverage at the catalyst and of O₂ partial pressure in the gas phase occur. This drastic variation in the chemical potential of O₂ in the reactor is accompanied by a change in the selectivity of the process: syngas formation – absent in the presence of O₂ at the catalyst interface – starts to occur as soon as the amount of O₂ drops to very low values at the catalyst surface. Thus, the reactor becomes stratified in two distinct zones. In the first zone, no syngas production is observed and O* is the most abundant adsorbate at the catalyst surface. In the second zone of the reactor, syngas production is observed and CO* and H* are the MARIs. This dependence on the selectivity to syngas is in agreement with the experimental observation.¹⁰ The structureless microkinetic model relates this selectivity change to the oxygen coverage dependence of the H* and CO* oxidation pathways, which are faster at high O* coverages than their desorption, thus leading to CH₄ total oxidation.¹⁸ The experiments of Grunwaldt et al.,^10,16 however, suggest that this selectivity change is related to a change of the morphology of the catalyst in an O-rich environment. To account for this structural variation, here we associate the chemical potential of gaseous species along the reactor coordinate calculated by microkinetic modelling with particular morphologies of the catalyst. This allows us to study how the reaction environment induces changes in the catalyst structure during reaction.

Bulk phase stability

In the aim of modelling how the catalyst morphology changes as a function of the chemical potential in the reactor, we first investigate which is the thermodynamically stable bulk phase of the catalyst when it is in contact with the gaseous environment.

In their experimental work, Grunwaldt et al.¹² observed that temperature induced oxidations and reductions of the catalyst are very fast: during the stabilization period of their tests (250 °C < T < 350 °C, duration = 2 hour) they made the temperature rise and decrease multiple times and they observed that the catalyst got fully reduced at high temperatures and fully oxidized at low temperatures. Then, the oxidation state of the catalyst remained stable with time. Thus, it can be assumed that for Rh/α-Al₂O₃ systems, the morphological transformation at relatively high temperatures (T > 250 °C) is reversible and it is not hindered by metastable structures.

Hence, with eqn (1), we calculate the oxygen chemical potential at which the phase transition between Rh, Rh₂O₃ and RhO₂ bulk phases occurs. Explicating the dependence of μ_O from temperature and oxygen partial pressure with eqn (2), we derive the bulk phase diagram of Rh and its oxides, represented in Fig. 2. Under the inlet conditions of the CPO reactor (T = 773 K and P_O2 = 0.01 atm), the stable catalyst bulk phase is Rh₂O₃. At the same temperature, when the oxygen partial pressure drops below 1.2 × 10⁻³ atm, metallic Rh becomes the most thermodynamically stable bulk structure. This value of oxygen partial pressure is close to the critical value at which the sharp change in the reactive environment is predicted by reactor modelling. Therefore, the sharp variation of the selectivity of the process (no syngas vs. syngas production) is associated with a concomitant change of the oxidation state of the catalyst. In zone 2 of the reactor, where the oxygen content in the system is very low, metallic Rh remains the most stable bulk phase. RhO₂ does not become thermodynamically stable in the range of conditions of our simulation.

Fig. 2 Bulk phase diagram of Rh (in blue), corundum Rh₂O₃ (in red) and rutile RhO₂ (in purple). In dark red and dark blue are the ranges of conditions that characterize the first and the second zone of the CPO reactor at 773 K and 1 atm.

Our calculations agree with the experimental observations of Grunwaldt et al.¹² and show a modification of the structure of the catalyst as a function of the chemical potential, which is not accounted for by the microkinetic model. Thus, the structureless microkinetic model, although able to well reproduce the kinetic macroscopic behaviour of the process, cannot give insights into the underlying atomistic details at the active site.

Catalyst morphology in zone 1

After modelling how the bulk phase changes with the change of reactive environment, we now aim at predicting the morphology of the catalyst surface and its relationship with the reaction environment. We use two different models to describe the catalyst in the two regimes found in the annular reactor. In the first zone, the catalyst is oxidized, and the reservoir of oxygen atoms is gaseous oxygen in contact with the catalyst surfaces. The O₂ adsorption/desorption elementary step is at quasi-equilibrium.

Therefore, we calculate the chemical potential of O* as a function of μ_O2 as:

With the aim of calculating the catalyst morphology in the first zone of the CPO reactor, the stability of eight crystalline facets of Rh₂O₃ are investigated: (0001), (11̄02), (112̄3), (11̄01), (101̄1), (112̄0), (101̄0), and (101̄2), selected because they are the most stable surfaces of α-Al₂O₃,³⁹ which has the same bulk structure as Rh₂O₃. For each (hkil) surface, different terminations, obtained from cleaving the bulk with parallel planes, are considered. Their specific surface free energies are calculated with eqn (5), neglecting the vibrational energy and the entropy contributions, because we calculated that at reasonably low temperatures (T < 1000 K) the error introduced remains lower than 0.001 eV Å⁻².

Fig. 3 shows the plot of the specific surface energies of the Rh₂O₃ facets as a function of the oxygen chemical potential. With dashed lines are underlined the chemical potentials at which the transition to metallic Rh (Δμ_O* = −1.15 eV, in brown) and to RhO₂ (Δμ_O* = −0.55 eV, in red) occur. In the whole range of stability of bulk Rh₂O₃, the most stable surface is Rh₂O₃(11̄02), and the second is Rh₂O₃(0001), both exposing stoichiometric terminations, in agreement with the findings of Scherson et al.²³ Rh₂O₃(112̄3) turns out to be very stable at high values of oxygen chemical potential: it exposes an over-stoichiometric termination (with an excess of oxygen atoms with respect to the stoichiometric value of N_O/N_Rh = 3/2), which reduces its surface energy with the increase of oxygen chemical potential. The most stable terminations of Rh₂O₃(112̄0) and Rh₂O₃(101̄0) are stoichiometric, whereas the most stable termination of Rh₂O₃(11̄01) has an excess of oxygens with respect to the stoichiometric ratio. Rh₂O₃(101̄1) shows a stoichiometric termination for Δμ_O* < −0.75 eV, then an over-stoichiometric one. Rh₂O₃(101̄2) presents an over-stoichiometric termination with very high specific surface energy. Under-stoichiometric terminations are not stable at any value of Δμ_O* at which Rh₂O₃ is thermodynamically preferred. Representations of the structures are provided in section 4 of the ESI.†

Fig. 3 Specific surface energy plot of eight Rh₂O₃ surfaces, as a function of oxygen chemical potential: Rh₂O₃(0001) in blue, Rh₂O₃(11̄02) in orange, Rh₂O₃(112̄3) in yellow, Rh₂O₃(11̄01) in purple, Rh₂O₃(101̄1) in green, Rh₂O₃(112̄0) in cyan, Rh₂O₃(101̄0) in pink, Rh₂O₃(101̄2) in chartreuse.

The effect of the support on the three-dimensional shape of the catalyst is taken into account by calculating the specific adhesion energy of Rh₂O₃(0001) and Rh₂O₃(11̄02) facets in contact with Al₂O₃(0001) and Al₂O₃(11̄02), which are the two most exposed facets of α-Al₂O₃.^32,40 First, we characterize the Al₂O₃ surface structures: exploiting eqn (8), we investigate the stability of different terminations of the Al₂O₃ facets, with Al, O or OH as terminal groups. The stoichiometric terminations are the most stable over the entire considered range of oxygen chemical potentials, whereas hydroxylated terminations become stable for high water chemical potential (Δμ_H2O = −1.475 eV, i.e. P_H2O = 100 bar at 773 K), in agreement with the data of García-Mota et al.³² and Marmier et al.⁴⁰ The specific adhesion energy between Rh₂O₃ and α-Al₂O₃ is evaluated with eqn (9). Independent from the combination of facets, the resulting value of adhesion energy is about 0.02 eV Å⁻² (details are provided in section 3 of the ESI†). This value is very small compared to the surface free energies of the Rh₂O₃ facets, yielding Wulff plots with a large contact area between the catalyst and the support.

In Fig. 4 the 3D shapes of Rh₂O₃ nanoparticles at different values of Δμ_O*, obtained with the Wulff–Kaishew construction method (eqn (10) and (11)), are represented. At the inlet of the reactor (Fig. 1), when Δμ_O* = −1.01 eV, the most exposed surfaces are Rh₂O₃(0001), Rh₂O₃(11̄02) and Rh₂O₃(112̄3), whereas the facets Rh₂O₃(112̄0) and Rh₂O₃(101̄1) appear in the Wulff plot, even with very small areas (Fig. 4, panel a). Close to the phase transition to reduced Rh bulk, at Δμ_O = −1.11 eV, the surface area of Rh₂O₃(112̄3) is smaller, whereas that of Rh₂O₃(11̄02), Rh₂O₃(112̄0) and Rh₂O₃(101̄1) is slightly wider (Fig. 4, panel b). Considering now conditions characterized by higher oxygen chemical potential in the system, when Δμ_O = −0.55 eV (close to the bulk phase transformation to RhO₂), we find that the nanoparticle exposes only the Rh₂O₃(0001), Rh₂O₃(11̄02) and Rh₂O₃(112̄3) surfaces, with a high amount of Rh₂O₃(112̄3). The catalyst shape in these conditions is represented in Fig. 4, panel c.

Fig. 4 Percentage surface areas of the investigated crystal facets of Rh and Rh₂O₃ along the axial coordinate of the annular reactor for the CPO at 773 K. Wulff–Kaishew constructions of Rh₂O₃ nanoparticles at the inlet of the reactor, when Δμ_O = −1.01 eV (a), at the conditions at which transition from zone 1 to zone 2 occurs, Δμ_O = −1.11 eV (b) and when the bulk phase transition to rutile RhO₂ occurs, when Δμ_O = −0.55 eV (c). Wulff–Kaishew construction of Rh nanoparticles under three different operative conditions: at the beginning of zone 2 of the CPO reactor (d), at the outlet of the CPO reactor (e), under thermodynamic equilibrium conditions of the gas phase at 773 K (f).

Catalyst morphology in zone 2

As reported in Fig. 2, in the second zone of the reactor, the catalyst nanoparticles are reduced; CO* and H* are the MARIs, with coverages lower than 0.20 ML. The gas phase is not at the thermodynamic equilibrium composition, and it contains different gaseous species reservoirs of CO* and H*: H₂O, CO₂, CH₄, CO and H₂. However, only the reactions of adsorption/desorption of gaseous CO and H₂ to CO* and H* turn out to be quasi-equilibrated from the microkinetic analysis. Therefore, we calculate the chemical potentials of CO* and H* as functions of μ_CO and μ_H2:

μ_CO* = μ_CO (13)

The presence of other adsorbates does not affect the stability of Rh surfaces, as a consequence to the fact that their coverage is calculated to be lower than 0.01 ML.

With eqn (4) we calculate the specific surface free energy of all clean Rh facets characterized by Miller indices lower than 4. Then, with the Wulff construction method (eqn (10)) we calculate the equilibrium shape of clean Rh nanoparticles under vacuum and we find that the exposed facets are only six: (100), (110) (111), (211), (311) and (331). We select these six facets for the following study, as a good compromise between computational time and accuracy of results. Exploiting eqn (4), we calculate the thermodynamically stable surface structures of the six selected facets of Rh, in the presence of CO*, H* or both the adsorbates on the surfaces. Fig. 5 shows an example of a double chemical potential plot of specific surface free energy (that of the Rh(110) facet). The specific surface free energy of the facet decreases with the increase of both chemical potentials. Fig. 6 shows the double chemical potential phase diagrams of several ordered and disordered surface structures that we calculated to be the most stable on the six considered surfaces at different values of CO* and H* chemical potentials. On all the facets there are regions in which co-adsorption is preferred.

Fig. 5 Double chemical potential specific surface free energy plot of Rh(110) with CO* and/or H* adsorbed. Some of the most stable surface structures calculated by ab initio thermodynamics are represented in the figure. The coverage amount is represented with increasing the red colour for CO* and the blue for H*. The grey area represents the clean Rh(110) surface. Purple areas are characterized by co-adsorption of the two adsorbate species.

Fig. 6 Surface phase diagrams of the six Rh facets investigated in our study, in the presence of CO* and/or H* adsorbates: (a) Rh(100); (b) Rh(110); (c) Rh(111); (d) Rh(211); (e) Rh(311); (f) Rh(331). As shown in the legend, the increasing H* coverage is indicated by blue shades, and CO* coverage is represented by red shades. Purple shades represent structures with both CO* and H* co-adsorbed. Orange crosses indicate the values of CO* and H* chemical potentials characteristic of the conditions found at the beginning of zone 2 of the CPO reactor. Green crosses indicate the conditions of the reactor outlet. Black crosses indicate the CO* and H* chemical potentials at which the gas phase is in thermodynamic equilibrium at 773 K.

Our simulation of the CPO annular reactor at 773 K and 1 atm shows that, at the beginning of zone 2 (see Fig. 1), the partial pressures of CO and H₂ are 4 × 10⁻⁶ atm (Δμ_CO* = −2.40 eV) and 2 × 10⁻⁵ atm (Δμ_H* = −0.87 eV), respectively. Under these conditions, the coverage on all the facets is very low (ϑ_CO* < 0.01 ML and ϑ_H* < 0.01 ML). At the outlet of the reactor, the CO partial pressure increases to 7 × 10⁻⁴ atm (Δμ_CO* = −2.06 eV), and the H₂ partial pressure is 4.2 × 10⁻³ atm (Δμ_H* = −0.70 eV). Under these conditions, co-adsorption of CO* and H* is preferred on Rh(100). On Rh(110), Rh(211) and Rh(311), structures with adsorbed CO* and H* have similar energy to structures with only CO* adsorbed. On Rh(111), the coverages predicted are very low, with H* found to remain adsorbed more favourably than CO*. On Rh(331), CO* alone is the most stable. Considering now the conditions of the equilibrium composition of the gas phase at 773 K (P_CO = 6 × 10⁻³ atm and P_H2 = 1.4 × 10⁻² atm), we calculate that Δμ_CO* = −1.91 eV and Δμ_H* = −0.66 eV. Under these conditions, we find that CO* is the favoured adsorbate on Rh(100), Rh(110) and Rh(331). Meanwhile, on Rh(111), Rh(211) and Rh(311) the co-adsorption competes with CO* adsorbed alone. Details on the tested surface structures are reported in section 4 of the ESI.†

The effect of the support is taken into account by calculating the specific adhesion energy of Rh(100), Rh(110) and Rh(111) in contact with Al₂O₃(0001) and Al₂O₃(11̄02) surfaces. The adhesion energy is calculated with eqn (14), assuming non-coherent epitaxy between Rh and Al₂O₃ facets. Rh(111) is calculated to be the most favourable facet with which the catalyst is in contact with the support, and similar adhesion energies are obtained with the two investigated Al₂O₃ facets (0.025 eV A⁻² for Al₂O₃(0001) and 0.029 eV A⁻² for Al₂O₃(11̄02)). Details are provided in section 3 of the ESI.†

With the Wulff–Kaishew construction method (eqn (10) and (11)) we eventually calculate the equilibrium shape of the Rh catalyst nanoparticles as a function of the reactive environment. When syngas starts to be produced in the reactor, the catalyst nanoparticles (represented in Fig. 4, panel d) show a highly faceted shape, where the close-packed Rh(111) surface prevails. Rh(311) is the second most exposed surface structure, followed by Rh(331) and Rh(100). Rh(211) and Rh(110), on the other hand, present very small areas in the Wulff plot. At the reactor outlet, the exposed areas of Rh(111), Rh(211) and Rh(331) are reduced, while those of Rh(100), Rh(110) and Rh(331) are enlarged (Fig. 4, panel e). When the gas phase is at thermodynamic equilibrium, the catalyst nanoparticle is more spherical and the height to diameter ratio is increased (Fig. 4, panel f).

In Fig. 4 the profiles of the relative areas of the crystal facets exposed by the catalyst nanoparticles along the reactor coordinate are also reported. In the first zone of the CPO reactor the catalyst morphology changes with the variation of the oxygen partial pressure in the reactive environment. In particular, over-stoichiometric terminations are more abundant when the oxygen content is high. In the second zone of the reactor, the amount of stepped surfaces and Rh(110) increases with the increase of the chemical potential of H₂ and CO, whereas the relative area of Rh(111) is reduced. This clearly shows that the relative abundance of the possible active sites is influenced by the reaction environment and thus changes as a function of the chemical potential.

Comparison with chemisorption analysis

The previous analysis has been performed irrespective of the actual size of the nanoparticle, since the Wulff–Kaishew construction is size-independent. Here, we aim at testing the ability of the framework for given values of catalyst nanoparticle size. At this scope, we report an analysis and interpretation of catalyst characterization experiments on a 0.5% Rh/α-Al₂O₃ catalyst prepared by chemical vapor deposition (CVD) of Rh(acac)(CO)₂ by Beretta et al.⁴¹ These experimental data include HRTEM images obtained after aging under CPO conditions at 850 °C and H₂ and CO chemisorption experiments. HRTEM allowed for the calculation of the particle size distribution of the catalyst under the conditions of the HRTEM analysis (vacuum, room temperature). Chemisorption data are used for estimating the ratio between the number of adsorbed molecules (N_A*) and the catalyst total number of atoms (N_Rh,tot). This quantity can be employed to calculate the catalyst dispersion (Rh atoms which constitute the catalyst surfaces, N_Rh,surf, divided by N_Rh,tot), if the adsorption stoichiometry (the ratio between N_A* and N_Rh,surf) is known. Thus, the catalyst dispersion can be calculated as follows:

Here, using as an input the particle size distribution from HRTEM analysis, we test the ability of our thermodynamic model to estimate the catalyst dispersion. In particular, we calculate the structure of the catalyst and the adsorption stoichiometry under the chemisorption conditions and we compare the resulting catalyst dispersion with the data of chemisorption. In order to input the information regarding the size we proceed as follows. We assume that the number of atoms in the nanoparticles is the same under both HRTEM and chemisorption conditions. This assumption implies that sintering phenomena do not occur (room temperature), but the shape can undergo modifications induced by adsorption of probe molecules.

To this aim, with ab initio thermodynamics (eqn (4)) and Wulff–Kaishew construction (eqn (10) and (11)) we calculate the catalyst shape in equilibrium with the environment of the HRTEM analysis. We assume that the catalyst surfaces are adsorbate-free under vacuum conditions. Then, for each diameter of the particle size distribution we calculate the number of atoms in each nanoparticle by filling a Wulff–Kaishew plot with the considered diameter and counting the resulting amount of Rh atoms inside the Wulff–Kaishew plot. We consider the Wulff–Kaishew plot diameter as the distance between the closest two sides of the Wulff–Kaishew plot projected in the xy plane.

An example of a nanoparticle with a diameter of 2 nm is shown in Fig. 7 (panel a). The corresponding height from Wulff–Kaishew construction turns out to be 1.18 nm, thus leading to a catalyst height/diameter ratio of 0.59. This ratio is in agreement with the qualitative observation of Beretta et al.⁴¹ based on HRTEM images (height/diameter = 2/3).

Fig. 7 Wulff–Kaishew construction of Rh nanoparticles on Al₂O₃(0001) under the conditions of: (a) HRTEM, (b) H₂ chemisorption and (c) CO chemisorption. Rh(111) is the surface in contact with the support.

After calculating the distribution of number of atoms of the sample, we study the catalyst structure under chemisorption conditions. With ab initio thermodynamics we calculate the specific surface free energy of the Rh facets in equilibrium with the chemical potential of CO and H₂ in the gas phase. Then, with Wulff–Kaishew construction we calculate the equilibrium shape of the catalyst under the conditions at which the chemisorption pulses are performed. For each group of the distribution of number of atoms, we fill a Wulff–Kaishew plot with the corresponding number of atoms and we count the resulting Rh atoms that form the catalyst surfaces (the ones with a coordination number lower than 11). Weighting the number of surface atoms by the percentage amounts of atoms of the groups, we calculate the dispersion of the entire population of nanoparticles.

Under H₂ chemisorption conditions (P_H2 = 4.78 × 10⁻² atm, room temperature), the H* chemical potential is Δμ_H* = −0.39 eV (eqn (2) and (14)). The corresponding catalyst shape under H₂ chemisorption conditions is shown in Fig. 7 (panel b) and allows for the direct calculation of the catalyst dispersion (atoms at the surfaces divided by the total number of atoms: 0.45 Rh_surf/Rh_tot). We perform the same analysis for CO chemisorption conditions (P_CO = 5.04 × 10⁻² atm, room temperature). Under these conditions the CO* chemical potential is equal to Δμ_CO* = −0.60 eV (eqn (2) and (13)). As a result of the different chemical potential, the catalyst shape is different from the one predicted under H₂ chemisorption conditions (Fig. 7, panel e), and leads to a value of dispersion of 0.46 Rh_surf/Rh_tot.

Next, we compare the values of catalyst dispersion calculated with our thermodynamic model with the ones obtained experimentally through uptake measurements.⁴¹ A crucial information for the interpretation of such experiments is the stoichiometry of chemisorption, which can be deduced by the H₂ and CO coverages at equilibrium (ϑ_H*, ϑ_CO*). Thus, the adsorption stoichiometry can be expressed as:

where N_F is the number of catalyst crystal facets, ϑ_A*,hkl is the coverage of adsorbate A* on the facet hkl and S_rel,hkl is the relative surface area of the facet hkl.

Under H₂ chemisorption conditions, the H* coverage on the facets that the catalyst exposes to the environment is 1.00 ML (see Table 2). Therefore, the 0.47 H*/Rh_tot uptake measurement leads to the following value of dispersion:

Table 2 Relative area percentages (S_rel) and adsorbate coverages of the catalyst facets (hkl) in equilibrium with the gas phase under H₂ and CO chemisorption conditions

hkl	H2 chemisorption		CO chemisorption
	S rel [%]	ϑ H* [ML]	S rel [%]	ϑ CO* [ML]
100	15.68	1.00	91.35	0.83
110	9.32	1.00	3.01	1.00
111	58.60	1.00	5.65	0.75
211	0.00	0.67	0.00	0.67
311	16.41	1.00	0.00	0.75
331	0.00	0.67	0.00	0.67

---

This value is found to be in very good agreement with the one obtained using our thermodynamic model (0.45 Rh_surf/Rh_tot).

Concerning CO chemisorption, the surface coverage for each facet is given in Table 2, and results in a CO* adsorption stoichiometry of 0.83 CO*/Rh_surf. As such, by analogy to H₂ chemisorption, the resulting dispersion is:

which is in very good agreement with the one obtained with the thermodynamic model (0.46 Rh_surf/Rh_tot).

Conclusions

In this work, we investigated how the reactor environment influences the morphology of a catalyst, in the context of CH₄ CPO over Rh. In agreement with previous analysis, we find that the reactor is stratified in two distinct zones. The first zone is characterized by a high amount of oxygen in the gas phase in contact with the catalyst, whereas in the second zone the oxygen in the gas phase is very low. We predict that the change from the oxidized zone to the reduced zone of the reactor is accompanied by a change in the structure and morphology of the catalyst. These theoretical findings agree with previously reported operando spectroscopic data. In particular, we found Rh₂O₃ to be the thermodynamically stable bulk phase in the first zone of the reactor, whereas in the second zone, where the amount of O₂ drops to very low values at the catalyst interface, metallic Rh turns out to be the most stable bulk structure. We analysed the catalyst morphology along the reactor coordinate using Wulff–Kaishew constructions, and thus obtained 3D representations of the catalyst nanoparticles. This allowed us to estimate the relative amount of catalyst active sites that was found to change significantly along the reactor coordinate. The proposed framework was also validated via the interpretation of characterization data on a conditioned 0.5% Rh/α-Al₂O₃ catalyst and led to calculated values of catalyst dispersion that are in very good agreement with the experimental estimation via H₂ and CO chemisorption.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The project leading to this work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Project SHAPE – grant agreement No. 677423). Computational time at CINECA, Bologna (Italy) is gratefully acknowledged.

References

1. R. Schlögl, Angew. Chem., Int. Ed., 2015, 54, 3465–3520.
2. K. Honkala, A. Hellman, I. N. Remediakis, A. Logadottir, A. Carlsson, S. Dahl, C. H. Christensen and J. K. Nørskov, Science, 2005, 307, 555–558.
3. D. G. Vlachos, A. B. Mhadeshwar and N. S. Kaisare, Comput. Chem. Eng., 2006, 30, 1712–1724.
4. M. Salciccioli, M. Stamatakis, S. Caratzoulas and D. G. Vlachos, Chem. Eng. Sci., 2011, 66, 4319–4355.
5. G. Ertl, Angew. Chem., Int. Ed., 2008, 47, 3524–3535.
6. G. A. Somorjai and J. Y. Park, J. Chem. Phys., 2008, 128, 182504.
7. L. Foppa, C. Copéret and A. Comas-Vives, J. Am. Chem. Soc., 2016, 138, 16655–16668.
8. A. Berkó and F. Solymosi, J. Catal., 1999, 183, 91–101.
9. T. Avanesian, S. Dai, M. J. Kale, G. W. Graham, X. Pan and P. Christopher, J. Am. Chem. Soc., 2017, 139, 4551–4558.
10. S. Hannemann, J. D. Grunwaldt, N. van Vegten, A. Baiker, P. Boye and C. G. Schroer, Catal. Today, 2007, 126, 54–63.
11. J. Grunwaldt, J. Catal., 2001, 200, 321–329.
12. J. D. Grunwaldt, S. Hannemann, C. G. Schroer and A. Baiker, J. Phys. Chem. B, 2006, 110, 8674.
13. M. Maestri, Chem. Commun., 2017, 53, 10244–10254.
14. M. Maestri, D. G. Vlachos, A. Beretta, G. Groppi and E. Tronconi, AIChE J., 2009, 55, 993–1008.
15. M. Maestri, D. Livio, A. Beretta and G. Groppi, Ind. Eng. Chem. Res., 2014, 53, 10914–10928.
16. J. D. Grunwaldt, S. Hannemann, C. G. Schroer and A. Baiker, J. Phys. Chem. B, 2006, 110, 8674–8680.
17. M. Maestri, D. G. Vlachos, A. Beretta, G. Groppi and E. Tronconi, J. Catal., 2008, 259, 211–222.
18. M. Maestri, D. G. Vlachos, A. Beretta, P. Forzatti, G. Groppi and E. Tronconi, Top. Catal., 2009, 52, 1983–1988.
19. A. Donazzi, M. Maestri, B. C. Michael, A. Beretta, P. Forzatti, G. Groppi, E. Tronconi, L. D. Schmidt and D. G. Vlachos, J. Catal., 2010, 275, 270–279.
20. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502.
21. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
22. V. I. Anisimov, J. Zaanen and O. K. Andersen, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 44, 943–954.
23. Y. D. Scherson, S. J. Aboud, J. Wilcox and B. J. Cantwell, J. Phys. Chem. C, 2011, 115, 11036–11044.
24. T. Thonhauser, S. Zuluaga, C. A. Arter, K. Berland, E. Schröder and P. Hyldgaard, Phys. Rev. Lett., 2015, 115, 1–6.
25. T. Thonhauser, V. R. Cooper, S. Li, A. Puzder, P. Hyldgaard and D. C. Langreth, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 1–11.
26. K. Berland, V. R. Cooper, K. Lee, E. Schröder, T. Thonhauser, P. Hyldgaard and B. I. Lundqvist, Rep. Prog. Phys., 2015, 78, 66501.
27. D. C. Langreth, B. I. Lundqvist, S. D. Chakarova-Käck, V. R. Cooper, M. Dion, P. Hyldgaard, A. Kelkkanen, J. Kleis, L. Kong, S. Li, P. G. Moses, E. Murray, A. Puzder, H. Rydberg, E. Schröder and T. Thonhauser, J. Phys.: Condens. Matter, 2009, 21, 84203.
28. P. Lazić, M. Alaei, N. Atodiresei, V. Caciuc, R. Brako and S. Blügel, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 1–6.
29. P. Villars and L. D. Calvert, Pearson's Handbook of Crystallographic Data for Intermediate Phases, American Society of Metals, 1986, p. 3258.
30. M. E. Grillo, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 70, 1–6.
31. R. G. Munro, J. Am. Ceram. Soc., 1997, 80, 1919–1928.
32. M. García-Mota, M. Rieger and K. Reuter, J. Catal., 2015, 321, 1–6.
33. J. Rogal and K. Reuter, Exp. Model. Simul. Gas-Surface Interact. React. Flows Hypersonic Flights, 2007, pp. 1–18.
34. J. Malcolm and W. Chase, NIST-JANAF thermochemical tables, American Chemical Society, American Institute of Physics for the National Institute of Standards and Technology, Washington, DC, New York, 4th edn, 1998.
35. C. T. Campbell, L. H. Sprowl and L. Árnadóttir, J. Phys. Chem. C, 2016, 120, 10283–10297.
36. P. Muller and R. Kern, Surf. Sci., 2000, 457, 229–253.
37. E. Ringe, R. P. Van Duyne and L. D. Marks, Nano Lett., 2011, 11, 3399–3403.
38. A. Donazzi, A. Beretta, G. Groppi and P. Forzatti, J. Catal., 2008, 255, 259–268.
39. M. Kitayama and A. M. Glaeser, J. Am. Ceram. Soc., 2002, 85, 611–622.
40. A. Marmier and S. C. Parker, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 115409.
41. A. Beretta, A. Donazzi, G. Groppi, P. Forzatti, V. Dal Santo, L. Sordelli, V. De Grandi and R. Psaro, Appl. Catal., B, 2008, 83, 96–109.

---

Footnote

† Electronic supplementary information (ESI) available: Modelling of the annular reactor; details of the calculations of the morphology of Rh and Rh₂O₃ nanoparticles; details of the calculation of the interfacial energy between the catalyst and the support; surface structures of Rh₂O₃ and Rh. See DOI: 10.1039/c8cy00583d

---