1 Introduction

The link between Brownian motion and diffusion [1], in conjunction with the development of computers in the middle of the twentieth century, created the possibility to study diffusion processes on the computer at the molecular level. Alder and Wainwright [2] were the first to perform a molecular dynamics simulation where the trajectories of individual particles are tracked over time by integrating the Newtonian equations of motion from the interatomic forces of individual molecules over time. With increasing computational resources, tremendous methodological advances, and more accurate descriptions of molecular interactions, simulations have moved from simple model systems to increasingly accurate models of real physical systems. In this article, we provide a perspective on the connections between theoretical fundamentals, molecular simulations and experimental measurements to characterize diffusion of guest molecules in microporous and mesoporous solids with a focus on the current challenges and future opportunities for synergistic modeling and experimental approaches in designing new porous adsorbents and catalysts. We also provide examples of fundamental insights into diffusion processes obtained from the combination of simulation and physical experiments. Apart from a brief overview in Sect. 2, this article does not include a detailed discussion of molecular simulation techniques and best practices, as there are numerous excellent resources on these topics for the curious reader [3,4,5,6,7].

The importance of molecular diffusion in porous solids is well established in applications such as heterogeneous catalysis and adsorption separations, where molecules must be able to travel rapidly in and out of the porous material. One example from heterogeneous catalysis, where a detailed understanding of diffusion mechanisms is imperative, is the methanol to olefins (MTO) process. The MTO process works particularly well in microporous catalysts containing large cavities, where catalytic aromatic molecules can accumulate, connected by narrow channels or windows in the hydrocarbon pool (HP) mechanism [8, 9]. When the HP mechanism is at play, small pore (8-ring) zeolites are preferred, as ethylene and propylene products can diffuse out of the cavities, while larger HP molecules are trapped, allowing for continued catalytic turnovers. The cavities help control the size of the HP molecules and product olefins to prevent coking due to oligomerization, but the catalyst still quickly deactivates (ca. 3 h), requiring frequent regeneration cycles [10, 11]. The size of the HP molecules, and, consequently, the stability of the catalyst are intrinsically tied to the diffusion of guest molecules, which cannot be ameliorated by simply increasing the pore size without disturbing the size and composition of the HP molecules. For this, and other processes where constrained diffusion is observed, like catalytic cracking, designing nanoporous solids becomes paramount to achieve the desired performance.

While we primarily focus on applications in heterogeneous catalysis and adsorption based separations, the fundamental concepts of molecular diffusion have broader application to membrane science [12], chromatography [13], gas exploration [14], environmental remediation [15], etc. Greater cross-community collaborations can help provide new insights into diffusion processes that are valuable across disciplines. This article is organized by material type, from the most amenable to simulation techniques to the most complex. First, crystalline microporous solids are examined, including zeolites and metal–organic frameworks (MOFs). Then, we examine mesoporous amorphous solids, including carbons, silicas, and other porous oxides. Finally, we consider hierarchical porous materials with well-defined micro/meso/macropore combinations, including hierarchical zeolites with crystalline microporous regions, zeolite pellets with meso/macropores, and crystalline hierarchical MOFs.

2 Simulation methods overview

Techniques for simulating diffusion using atomistic models can be generally placed into two families: those that follow molecular trajectories over time and those that sample the free energy surface and use transition state theory (TST) to calculate the diffusion rate. Generating molecular trajectories is the more intuitive and easily understood approach. In this technique, known as molecular dynamics (MD) simulation [4, 5], the forces acting on each atom l in the system are calculated at each time step and the atomic positions are propagated forward in time according to Newton’s equations of motion: Fl = ml al in its simplest form. Accurate expressions for the potential energy are needed for calculating the forces. In classical MD simulations, the potential energy is typically described as a sum of terms describing bond stretching, bond angle bending, torsional potentials, and non-bonded interactions, such as dispersion, repulsion, and Coulombic forces. For adsorbates in a porous solid, the non-bonded host–guest and guest–guest dispersion and repulsion interactions are typically described by a Lennard–Jones potential between all pairs of atoms. Similarly, Coulomb interactions are calculated from partial charges placed on framework and adsorbate atoms. The accuracy of a simulation is typically more sensitive to the accuracy of the Lennard–Jones parameters and partial charges than to the bonded interactions, such as bond stretching.

The Newtonian equations can be integrated numerically using a variety of different integration schemes. Typical time steps are ~ 1 fs, and with modern computers MD can simulate times scales of 10s of ns routinely. In most MD simulations, molecular trajectories are generated for a system at equilibrium (i.e., without concentration gradients), and the self-diffusion coefficient for species i, \({D}_{s,i},\) is calculated from the mean-squared displacement of adsorbates over time using the Einstein relation [6]:

$$D_{{s,i}} = \mathop{\lim }\limits_{{t \to \infty }} \frac{1}{{2dt}}\left\langle {\frac{1}{{N_{i} }}\sum \limits_{{l = 1}}^{{N_{i} }} \left[ {\user2 {r}_{{\user2{l},\user2{i}}} \left( t \right) - ~{\user2{r}}_{{{\user2{l}},\user2{i}}} \left( 0 \right)} \right]^{2} } \right\rangle ,$$
(1)

where d is the dimensionality of the system, \({N}_{i}\) is the number of adsorbates of species i, and \({r}_{l,i}\)(t) is the position of adsorbate l and species i at time t. The self-diffusion coefficient \({D}_{s,i}\) can be calculated from the slope of the mean-squared displacement versus time plot. Care must be taken to simulate long enough and to consider only the linear diffusive regime, where the mean-squared displacements are not correlated with short-timescale ballistic motion. In porous media, there also exists a sub-diffusive region that is correlated with pore wall collisions. For slit pores and fractal media, diffusion may become anomalous and the diffusivity time-dependent, as discussed in Sect. 4. An alternative approach to extract the self-diffusion coefficients from MD simulations is to use the velocity autocorrelation function (VACF) and the Green–Kubo relation [6]

$$D_{{s,i}} = \frac{1}{d}\frac{1}{{N_{i} }}\mathop \sum \limits_{{l = 1}}^{{N_{i} }} \mathop \smallint \limits_{0}^{\infty } \left\langle {\user2{v}_{{\user2{l},\user2{i}}} \left( 0 \right) \cdot \user2{v}_{{\user2{l},\user2{i}}} \left( t \right)} \right\rangle dt$$
(2)

where \({{\varvec{v}}}_{{\varvec{l}},{\varvec{i}}}\) is the velocity of adsorbate l of species i. Self-diffusivities calculated from MD simulations can be compared directly to experimental measurements from pulsed field gradient nuclear magnetic resonance (PFG NMR) and quasi-elastic neutron scattering (QENS). It should be noted that QENS probes the same time and length scales as MD, while PFG NMR probes longer time and length scales. If large crystals are used, however, PFG NMR can measure the intracrystalline self-diffusivity without effects from intercrystalline diffusion. QENS can also be used to obtain transport diffusivities, which we discuss next [16,17,18].

In practical applications where concentration gradients occur, it is often the Fickian or transport diffusion coefficient that is of interest. In general, the transport and self-diffusion coefficients are different quantities, although they coincide in the limit of very low adsorbate loadings. Transport diffusivities can be computed from equilibrium properties by multiplying the so-called “corrected” diffusivity (Dc) and a thermodynamic factor (Γ) [19],

$$D_{t} = D_{c} {{\Gamma }}$$
(3)

Γ can be obtained from the derivative of the sorption isotherm \(\Gamma = \left( {\frac{\partial \ln f}{{\partial \ln c}}} \right)_{T}\), if the isotherm is represented in terms of the fugacity f and the quantity adsorbed c. The adsorption isotherm can be obtained from a grand canonical Monte Carlo (GCMC) simulation using the same atomistic model that is used in MD. The corrected diffusivity of species i (Dc,i) can be calculated from equilibrium MD as follows [7]:

$${D}_{c,i}=\frac{1}{2d{N}_{i}}\underset{t\to \infty }{\mathrm{lim}}\frac{d}{dt}\left\langle {\left(\sum_{l=1}^{{N}_{i}}\left[{{\varvec{r}}}_{{\varvec{l}},{\varvec{i}}}\left(t\right)- {{\varvec{r}}}_{{\varvec{l}},{\varvec{i}}}\left(0\right)\right]\right)}^{2}\right\rangle$$
(4)

Note that this Einstein equation does not average over all particles (as in Eq. 1), but rather tracks the collective adsorbate translational motion for a species i. Combining Eqs. 3 and 4 provides a purely equilibrium route to the transport diffusivity, Dt,i. Note that at very low adsorbate loadings, the thermodynamic factor \(\Gamma\) is close to unity and the corrected, transport, and self-diffusivities all coincide. Non-equilibrium MD techniques, which incorporate a gradient of concentration under steady-state or non-steady-state conditions, can also be used to obtain transport diffusion coefficients [20,21,22,23].

There are other situations where diffusion is anomalous and the mean-squared displacement is not proportional to time. An example is single-file diffusion, where the mean-squared displacement varies with the square root of time [24]. This may occur in situations where one-dimensional pores are so constricting that adsorbates cannot pass each other. In this situation, the single-file mobility factor F has been proposed, which has the dimensions of m2 s−1/2 [25,26,27,28].

When considering mixtures of adsorbates (or even the implicit mixture of a single adsorbate and the host), the Maxwell–Stefan diffusion formulism is useful to quantify the exchange between adsorbates in a mixture. Maxwell–Stefan diffusion coefficients (Đ) can be calculated for multicomponent mixtures using MD simulations, providing rich detail into how the presence of co-adsorbed species influences diffusion. One drawback of using Maxwell–Stefan diffusion coefficients is that they only directly correspond to the Fickian diffusion coefficients at low loadings, making direct comparison with physical experiments difficult. For zeolites and other microporous materials, where molecules cannot directly pass each other, they also fail to account for the topology of the pore network or heterogeneity in the zeolite structure [29, 30].

For molecules that diffuse slowly (~ 10–11 m2 s−1 or slower), molecular dynamics simulations become intractable, as the probability of observing diffusive events becomes small on the time scale accessible to the simulations. In this case, sampling the free energy along a diffusion pathway provides a route to obtain the diffusion coefficient using transition state theory (TST). This sampling can be accomplished using a variety of techniques, such as particle insertion (or Widom insertion) [31], the multiple histograms technique [32], or umbrella sampling [33]. For example, Dubbeldam et al. [34, 35] used configurational bias Monte Carlo to insert adsorbates at different locations in various cage-like zeolites. The free energy of the adsorbate was then averaged across a reaction coordinate, which, for this example, was the path between two zeolite cages separated by a window.

For diffusion of adsorbates in microporous materials, in which the free energy profile is characterized by a landscape of valleys (energy wells deeper than ~ 3kbT) with barriers in between, TST can be applied. The valleys are associated with adsorption sites. If the free energy for a diffusion event along a reaction coordinate is known, the rate constant for hopping from state a to state b can be calculated from TST [36]:

$${k}_{a\to b}=\sqrt{\frac{{k}_{b}T}{2\pi }}\frac{\int \delta \left[C\left(q\right)\right]\left|{\nabla }_{q}C(q)\right|{e}^{-\beta F(q)}dq}{\int {e}^{-\beta F(q)}dq}$$
(5)

Here, F(q) is the free energy along the diffusion coordinate (q), C(q) denotes the dividing surface between neighboring states, kb is the Boltzmann constant, T is temperature, and β is 1/ kbT. Rate constants such as those described in Eq. (5) overestimate diffusion, unless corrected for dynamical re-crossings [37]. TST is easiest to apply at low loading, but it can be applied at finite adsorbate loadings, allowing comparisons with experimentally relevant conditions [38].

With the rate constants in hand for all pairs of connected states (adsorption sites), the self-diffusion coefficient can be obtained by solving the so-called master equation:

$$\frac{d{p}_{A}}{dt}=-\sum_{B}{k}_{A\to B}{p}_{A}+\sum_{B}{k}_{B\to A}{p}_{B}$$
(6)

where pA(t) is the probability of finding the system in state A at time t. The master equation is typically solved using a coarse-grained stochastic simulation referred to as kinetic Monte Carlo. In the simulation, molecules move between adjacent sites with frequencies determined by the rate constants kA→B. Since each possible event A → B has a first-order rate constant associated with it, the time between events of this type is an exponentially distributed random variable. The chain of events constitutes a Poisson process, and the average time between A → B moves is 1/kA→B.

3 Microporous crystalline materials

3.1 Introduction

It may seem straightforward to predict the diffusion of molecules in microporous crystalline materials using MD simulation, as it is just a matter of calculating the forces on all atoms at each time step and integrating the Newtonian equations of motion forward in time. There are, however, complexities in the models for calculating the forces, the time scales that can be accessed with straightforward MD, and our ability to characterize materials that may not be perfectly crystalline. We illustrate these complexities in this section, starting with the simplest systems and moving toward increasing complexity. We focus on two types of crystalline microporous materials: zeolites and metal–organic frameworks (MOFs). Both classes of materials are of practical interest, and, due to their regular crystalline structures, have been used as model systems to study the fundamentals of diffusion under confinement. For both classes of materials there are many examples of studies that combine experimental measurements and computational modeling, with several comprehensive reviews on diffusion in zeolites [3, 7, 39,40,41,42,43,44,45,46,47] and MOFs [44, 48]. Here, we focus on recent comparisons of modeling and experiments, many of which are fueled by applications in separation and catalytic processes in which diffusion plays an important role, such as CO2 capture, kinetic separation processes, and catalysis in small-pore zeolites. We discuss different aspects of simulating diffusion in microporous zeolites, including intracrystalline diffusion, external surface barriers, effects of adsorbate loading and framework flexibility, and force field development. We show how results from simulation have provided molecular-level insights into experimental measurements and how experiments have been used to validate modeling. We then turn our attention to microporous MOFs, where the effects of framework flexibility, metal node identity, and adsorbate loading are highlighted. We conclude each discussion with our perspectives on how advances in computer simulations can address some of the unresolved facets of diffusion in microporous materials.

3.2 Zeolites

3.2.1 Introduction

Zeolites are microporous, crystalline aluminosilicates, which are widely used as industrial catalysts and adsorbents due to their tunable and well-defined pore sizes and shapes. The International Zeolite Association recognizes over 252 zeolite framework types, and many yet un-synthesized hypothetical structures have been proposed. The high crystallinity and, therefore, well-defined pores in zeolites have made them model systems for fundamental research on diffusion in nanopores.

In this section, we start with diffusion in medium- and large-pore zeolites, where MD simulations have shown good agreement with self-diffusivities measured by PFG NMR. We then turn to diffusion in small-pore zeolites, which have proven to be much more difficult to model. The state-of-the-art in force field development, consideration of framework flexibility, and modeling of surface barriers are discussed in the context of recent examples from the literature.

3.2.2 Intracrystalline self-diffusion in medium- and large-pore zeolites

The most straightforward way to perform an MD simulation is to simulate a system at equilibrium in a central “box” of atoms surrounded by periodic images of the central box (so-called periodic boundary conditions). This effectively models an infinite system, with no external crystal boundaries, and thus only intracrystalline diffusion is considered. By tracking the molecular positions over time—taking into account that molecules initially in the central box may end up far from their initial positions—the mean-squared displacement can be calculated, and thus the self-diffusion coefficient can be obtained from Eq. (1). Early MD simulations in zeolites took this approach, and to save computational time many studies considered the framework atoms to be fixed at their crystallographic coordinates. Force fields were often taken from previous work that had adjusted the Lennard–Jones parameters of the framework atoms to fit experimental data for thermodynamic quantities, such as the adsorption isotherm and enthalpy of adsorption. June et al. [49], for example, took this approach in modeling the self-diffusion of methane and xenon in the zeolite silicalite (all-silica MFI). As shown in Fig. 1, the predicted self-diffusion coefficients are in excellent agreement with PFG NMR data from the literature. This agreement was an early indication of the power of MD simulations to predict diffusion properties in zeolites, and it also provided support for the PFG NMR results. At the time, there was a vigorous discussion in the literature about the disagreement of several orders of magnitude when comparing diffusivities from PFG NMR and uptake measurements. The MD simulations of June et al. also provided information on the anisotropy of diffusion in the MFI channel system and indicated that methane and xenon prefer to adsorb in the straight and sinusoidal channels and avoid the intersections. Early work on MD simulations of self-diffusion studies in zeolites has been reviewed by Demontis and Suffritti [50]. In a more recent example, Guliants and Huth [51] found that their self-diffusion coefficients for nitrogen and methane in silicalite from MD were similar to those from PFG NMR at 573 K, and, interestingly, the MD results were similar when using three different force fields.

Fig. 1
figure 1

Copyright 1990 American Chemical Society

Self-diffusivities for methane in silicalite at various temperature and loadings from MD simulation. Predictions are compared against the experimental PFG NMR data of Caro et al. [52]and Jobic et al. [53] Reprinted with permission from June et al. [49].

The importance of extraframework cations on guest diffusion is illustrated in a study by Newsome and Coppens [54], who reported MD simulations for CO2 and N2 in Na-ZSM-5 (MFI). They found that the strength of the adsorbate–cation interactions define an initial effect, at low loading, where diffusivities increase with loading because the first molecules have strong interactions with the cations and then—when the cation sites are filled—additional molecules interact only with the silica portions of the zeolite surface. This happens for systems with strong adsorbate–cation interactions, such as CO2 in Na-ZSM-5. At the highest loading, the mobility decreases with loading, due to crowding in the pores as observed in all-silica MFI, for example in Fig. 1. For N2, which interacts less strongly with the sodium cations, the first effect at low loading (increasing diffusivity with increasing loading) is almost negligible, and instead a monotonic decrease in diffusivity is observed with increased loading. The effect of temperature is also important, especially with strongly interacting adsorbates, since temperature can be used to fine tune the threshold loading until which the effect of increased diffusivity at low loading is observed. The higher the temperature, the lower the threshold loading.

One advantage of PFG NMR for studying diffusion is that one can obtain separate signals for the different species in a mixture and thus obtain self-diffusivities for individual components under mixture conditions. Snurr and Kärger [55] measured self-diffusivities for CH4/CF4 mixtures in silicalite by monitoring the 1H signal from CH4 and the 19F signal from CF4. They performed MD simulations for the same mixtures and found good agreement between simulation and experiment over a range of conditions using the same force field as in the early work by June et al. [49] At a constant loading of 12 molecules per unit cell, the self-diffusivities of both species decreased as the mole fraction of the larger, slower CF4 was increased. A study by Fernandez et al. [56] using MAS and PFG NMR experiments on n-butane/iso-butane mixtures in silicalite at 363 K for loadings between 0–4 molecules/u.c. reports ratios of experimental versus MD calculated diffusion coefficients for n-butane of 6/3, 3.8/1.8, 0.4/0.01 for the respective n-C4/i-C4 loadings of 4/0, 3/1, 2/2. Hence, the values are good enough to reproduce the order of magnitude, which is the usual target, although, as the iso-butane content increases, the n-butane predictions differ increasingly from the experimental results. However, the simulations are able to qualitatively capture the fact that at a n-C4/i-C4 ratio of 2/2 iso-butane tends to locate preferentially at the channel intersections, retarding by two orders of magnitude the diffusivity of n-butane compared to the 4/0 mixture.

Faujasite is an important large-pore zeolite used in processes such as petroleum refining. It is known as zeolite X or Y, depending on the Si/Al ratio. NaY, with the same structure but fewer cations (Si136Al56Na56O384) than NaX (Si105Al87Na87O384), has been investigated for the effect of loading on the diffusion of methane, similar to the work above for MFI. A study by Déroche et al. [57] summarizes previous research and provides explanations for their results using QENS experiments and MD calculations. For temperatures above 225 K, the diffusion of methane monotonically decreases with loading. However, at lower temperatures, some methane molecules remain adsorbed to the sodium cations, reducing their self-diffusivity, but as loading increases further, higher self-diffusivities are observed, due to incoming molecules feeling a weaker interaction with the “solvated” cations. Then, after a threshold loading is passed, the expected decreasing trend of diffusivity with loading is observed. Some MD simulations were not able to reproduce these features, probably due to insufficient accuracy in describing methane–cation interactions. Granado et al. [58] made a careful comparison between experimental results and MD simulations for propane and propylene in NaX (Si104Al88Na88O384, also called 13X). For loadings up to 32 molecules/u.c., all guest molecules are adsorbed into sodium cations and hence a monotonic decrease in the diffusion coefficient is observed, with the MD results giving an excellent agreement with PFG NMR experiments.

Zeolite crystal quality and potential defects must be considered when making comparisons between experiments and simulations, since perfect crystallinity is almost always assumed in simulations and defects are difficult to characterize experimentally. For instance, disagreements between diffusivities from MD and PFG NMR for pentane isomers in NaX are claimed to be due to crystal defects of the NaX samples [59]. A similar case is reported for diffusivities of methanol in a dealuminated HY sample [60], in which disagreements between MD and QENS were attributed to defects of the sample. Defects can be expected when samples are obtained by dealumination, with faujasite being a well-documented case. Agreement between calculations and experiments in those cases is more difficult, but could be obtained, such as in the study by Wongthong et al. [61] in a sample of pure silica faujasite in which MD of small alkanes was found to give results (13 × 10–9 m2/s for propane) comparable to PFG NMR (9 × 10–9 m2/s) at low loading (8 molecules/u.c.) and with the calculated heat of adsorption being 26 kJ/mol. Another study of alkanes in pure silica faujasite gives for propane (at 8 molec./u.c.) a very similar value, 10–8 m2/s, with a calculated heat of adsorption of 20 kJ/mol [62]. Of course, when the sample is of high quality, which is typical for NaY, it is easier to find agreement of the MD simulations and the experiments, such as in the case of propane rotational constants in NaY [63].

3.2.3 Intracrystalline self-diffusion in small-pore zeolites

There has recently been a proliferation of computational studies of diffusion in small pore zeolites due to the industrial importance of the MTO reaction, which is recognized as being considerably controlled by diffusion. Gao et al. [11] performed MD simulations of methane, ethane, and propane in three cavity-type small-pore silicoaluminophosphate molecular sieves (LEV, CHA, RHO) with similar 8-ring window size. As shown in Fig. 2, the simulated self-diffusion coefficients agree reasonably well with results from PFG NMR, and the qualitative loading dependences are correctly reproduced. The parameters for the force field were taken from the work of Jee and Sholl, who used a steeper repulsion Lennard–Jones potential for the interactions of methane with zeolite oxygens in order to reproduce experimental diffusivities in silica DDR zeolites [64]. The self-diffusivity increases with loading in LEV and CHA but decreases in RHO. The diffusion was quantitatively described as an inter-cavity hopping process, and the jump frequency and jump length were extracted from the MD simulations.

Fig. 2
figure 2

Reproduced from Gao et al. [11] with permission from Elsevier

Diffusion characteristics of methane in cavity-type silicoaluminophosphate molecular sieves. (left) Experimental (PFG NMR) and simulated (MD) self-diffusion coefficients of methane in LEV, CHA and RHO at the loading of two molecules per cavity at 298 K. (right) Loading dependence of experimental and simulated self-diffusion coefficients for methane in cavity-type molecular sieves at 298 K.

A detailed MD study of the diffusion of ethene and propene in three small-pore zeolites (AEI, CHA, AFX), without acid sites (pure silica and aluminophosphate) and with acid sites (aluminosilicate and silicoaluminophosphate), allowed Ghysels et al. to establish correlations for how pore size is affected not only by topology but also by chemical composition [65]. A single newly defined parameter, named the accessible window area of the 8-ring, was demonstrated to correlate with the number of cavity-to-cavity jumps per unit time, which is directly related to the activation barrier and, hence, to the self-diffusivity.

MD simulations of ethane diffusion (using a united atom model) in zeolite CHA from two different groups show somewhat contrasting results. Dubbeldam et al., [36] who used a rigid zeolite model, reported self-diffusivities of 0.5 and 0.9 × 10–9 m2/s at low loading (1 molecule per cavity) and 300 and 600 K, respectively (Fig. 3), whilst Schüring et al., [66] who used a flexible zeolite model, found slightly higher diffusivities of 2.0 and 1.6 × 10–9 m2/s for the same loadings and temperatures. However, the former found an increase in diffusivity with temperature (as expected), whilst the latter found a decreased diffusivity with temperature, which was explained through entropic considerations: heating the system moves ethane away from eight-ring windows, on average, increasing the entropic barrier for cage-to-cage motion. This example demonstrates that for tightly fitting molecules simulations of diffusion can be sensitive to the inclusion of the framework flexibility. While these authors were able to obtain self-diffusivities from MD, sometimes the diffusive events may occur on time scales longer than those accessible to MD, which is for diffusivities typically below ~ 10–12 m2/s. As described above, in such cases transition-state theory (TST) can be used to extract diffusion coefficients from the same molecular-level model. Dubbeldam et al. used their MD results to test the TST approach, and Fig. 3 illustrates excellent agreement between the two computational methods [36].

Fig. 3
figure 3

Copyright 2006 American Chemical Society

Anisotropic diffusion of ethane in CHA-type zeolite computed by MD simulation and by dynamically corrected TST at 300 and 600 K. Reprinted with permission from Dubbeldam et al. [36].

For small-pore zeolites, even small molecules may fit very tightly in the narrowest windows that control diffusion through the material. In such cases, small changes in the force field parameters, such as the size of the framework oxygen atoms, may have a large effect on the predicted diffusion coefficients. Also, as discussed in more detail in Sect. 3.2.4 it may be essential to account for motion of the framework atoms in the MD simulation or the TST calculations. The effect of the force field on diffusion of CO2 in LTA zeolite has been studied by García-Sánchez et al. [67]. They found that small differences in the force field result in significant differences in the adsorption isotherms and very large differences in the loading dependence of the diffusivities. A similar effect was found later for methane in flexible, compared to rigid LTA zeolite [68], in which the authors found a small effect of the flexibility on adsorption and a much larger effect on diffusion. A wider study by Lim et al., [69] also with CO2 and 156 zeolites, explored how different force fields can yield very different results for adsorption enthalpies and isotherms, especially at high loading. For a good parameterization, it is suggested that diffusion coefficients be considered, as well as adsorption isotherms.

In many systems where the pore size is relatively larger than the sorbate, the self-diffusivity decreases monotonically with loading [70, 71], as shown, for example, in Fig. 1. However, for small-pore zeolites that consist of cages separated by narrow windows, the self-diffusivity often increases with loading, reaches a maximum, and then decreases (sometimes dramatically) as the saturation loading is approached, as shown, for example, in Fig. 3. An important related question is how diffusion coefficients vary with loading and composition in multicomponent systems? MD simulations can predict not only self-diffusion coefficients but also Fickian (or “transport”) diffusion coefficients and the Maxwell–Stefan diffusivities [7]. For binary mixtures, the Maxwell–Stefan equations provide a convenient framework for estimating the mixture diffusion coefficients from the single component values as well as cross terms that depend on the interaction between different components. The self-exchange Maxwell–Stefan (MS) diffusion coefficient (Ɖii) can be approximated in terms of the MS diffusion coefficient of the single component (Ɖi), the self-diffusion coefficient (Ds,i) and the fractional occupancy (θi) as θi·Ɖi/Ɖii = Ɖi/Di−1, and, similarly, for the second component in a binary mixture (i = 1, 2). Ɖi and Ds,i can both be calculated from MD simulations, through either the mean square displacements or velocity autocorrelation functions. The higher the value of θi·Ɖi/Ɖii, the stronger the correlation effect, and the weaker the confinement, the stronger the correlation. The interpretation for this is that strongly confined molecules cross the narrow windows individually, without strong intermolecular interaction, whilst, for weakly confined molecules, there are strong intermolecular interactions at all steps in the diffusion pathway. Weakly confined molecules show Ɖ increasing monotonically as the loading increases, whilst strongly confined molecules show a decreasing trend of Ɖ with loading, followed by an increasing trend at larger loadings up to saturation. Correlation effects in binary mixtures, in general, cause the more mobile species to be slowed down and the tardier species to be sped up. The model by Reed and Ehrlich [72], derived for surface diffusion (which we shall return to in Sect. 4.4.1), has been used by Krishna and coworkers in the context of zeolites, which allows an interpretation of the loading and temperature dependence of the MS diffusion coefficient, Ɖ(θ, T). The normalized diffusion coefficient, Ɖi(θi, T)/Ɖi(0, T), can be expressed as a function of loading and a parameter, δEi, which represents the reduction in the free energy barrier for inter-cage hopping of molecules across windows due to the effect of loading and temperature. For a large reduction in free energy barriers, δEi/RT > 1, there is a significant increase in Ɖi with increasing fractional loading (θi), up to a limiting value near saturation. If, on the other hand, the intermolecular interactions do not affect the free energy barrier, then δEi = 0 and Ɖi decreases linearly with (1−θi). For each temperature, with MD calculations of Ɖi at different loadings, it is possible to obtain a fit which gives the value of δEi. For different temperatures, the variation of δE i with T can also be obtained. A review by Krishna contains data for the loading dependence of the free energy of methane diffusion in LTA, CHA and FAU [73], showing different trends that can be qualitatively explained with the results outlined above. The free energy activation barrier decreases with loading, due to the higher energy of the adsorbed molecules in the cages, whilst the energy of the transition state, described as a hopping of isolated molecules through the windows, remains constant.

The degree of confinement can be quantified as the relative size of the molecule with respect to the critical dimension for diffusion. The diameter of the largest free sphere that can diffuse in a given direction or the pore-limiting diameter can be used for the critical dimension for diffusion. Krishna explains that the zero-loading MS diffusivity decreases monotonically (although not linearly) with the degree of confinement. There are considerable limitations to this concept: molecules are not spheres; molecular size is a somewhat ambiguous quantity; it may be difficult to define a unique critical pore size in a zeolite; and the critical zeolite windows may not be circular. Despite the limitations, this concept is useful and can be used to develop correlations, as shown in Fig. 4 [73]. This reasonably accurate formulation has as one of its main goals to establish equations that relate diffusivity with loading so that it becomes possible to rationalize the different behaviors in terms of the physico-chemical properties of the zeolite host and the molecular guests.

Fig. 4
figure 4

Copyright 2009 American Chemical Society

MD data for the M-S diffusivity at zero-loading, Ɖi(0), for Ne, Ar, CH4 and Kr, in a wide variety of zeolites, MOFs, COFs, and cylindrical silica pores, as a function of the degree of confinement, defined as σ/dp where σ is the Lennard–Jones size parameter of the guest and dp is the pore diameter. Reprinted with permission from Krishna [73].

Another way to rationalize specific diffusivity trends with loading is the so-called “commensurate effect,” introduced by Smit and coworkers [34], which states that molecules that fit particularly well in micropores diffuse slower. This leads to an apparent anomalous trend of the diffusivities in the series of linear n-alkanes in, for instance, chabazite, in which instead of a monotonic decrease of diffusivity with increasing chain length, a decreasing diffusivity is found up to octane (the most commensurate with CHA cavities), but with larger n-alkanes (up to n = 12) diffusing faster than n-octane. This interpretation has to be made with care that other factors are not affecting the comparison of experiments with calculations [74], and this was soon realized as an explanation of the so-called 'window effect', in which the cracking rate of n-alkanes was related to a complex bimodal function [75].

3.2.4 Framework flexibility

As noted above, for some systems the flexibility of the zeolite framework can have a significant effect on diffusion. For example, O'Malley and Catlow [76] performed a comparison of united-atom + fixed-framework with full-atom + flexible framework for diffusion of long n-alkanes in silicalite. The latter shows better agreement with experiments, which suggests that pore breathing enhances the diffusivities. Schüring et al. [77], using a rigid framework, calculated the self-diffusivity of n-hexane in silicalite as 1.8 × 10–10 m2/s, which compares with 3.8 × 10–10 m2/s obtained by Leroy et al. [78] from MD. The latter study makes a comparison of results using rigid and flexible frameworks. Self-diffusion is enhanced using a flexible silicalite for the lowest loadings and the shortest alkanes, methane and n-butane, while no effect is observed for n-hexane and n-octane diffusion. In both studies, the predicted diffusivities were in very good agreement with PFG NMR result of 1.2 × 10–10 m2/s. A small effect of framework flexibility on diffusion is due to the MFI pore diameter being considerably larger than the cross section of a hexane molecule. Bermúdez and Sastre [79] studied the pore breathing of four different structures of silicalite using four different force fields and showed differences of up to 0.4 Å in pore channel diameters, which is very significant for channels with an average size of 5.1 × 5.6 Å (MFI).

The effect of framework flexibility in zeolites is important for tight-fitting molecules, and for frameworks that exhibit considerable breathing. Otherwise, MD simulations using rigid frameworks can be used. It should be noted that in some systems, for instance smooth-walled carbon nanotubes, incorporation of framework flexibility effects is necessary to simulate diffusive behavior at low loadings (see Sect. 4.4.2) [80,81,82]. This is not the case for zeolites, because pore wall corrugations introduce sufficient randomness in molecular collisions, and diffusive behavior is observed at low densities [80].

Framework flexibility in zeolites is particularly crucial for tight-fitting molecules near the repulsive region in the host–guest intermolecular interaction. However, since force fields are often parameterized on adsorption data (isotherms, enthalpies of adsorption [83]), this is the most difficult region for a force field to perform effectively, since it is far from the region in which parameters are usually obtained. Another complication is that molecular dynamic simulations do not perform well for slow diffusion, typically below 10–12 m2/s, and so more complicated TST or other free-energy based calculations may be required. Also, simulations with a fully flexible framework are significantly more time consuming. Finally, we stress that the issue of framework flexibility cannot be separated from the effect of force field quality. If the force field is not of sufficient quality, the flexibility simulated will not correspond to the real framework dynamics. We address the quality of the forcefields and new developments in the next Sect. 3.2.5. An excellent discussion is found in the study by García-Sánchez et al. [84] and the subsequent comment by Krishna and van Baten [85]. The former argues that none of the current force fields is able to capture the average structure exactly and specifically the 8-ring windows separating the cavities of LTA zeolite, since a small deviation of the window from the crystal structure leads to very different diffusivities. The latter argues that the lattice flexibility, per se, has no influence on the self-diffusivity of methane in cavity-type zeolites with 8-ring windows. Zimmerman et al. [86] used a dynamically corrected TST method, with a detailed analysis of dynamics versus static pore size, to show that, with a particular choice of force field, the average values coincide. This allows for accurate comparisons between fixed and flexible frameworks. For zeolites with relatively wide channels, such as AFI and LTL, a factor of two difference has been found in the flexible versus rigid frameworks for diffusion coefficients of methane at large loading. At low loadings, the effect was almost negligible.

The ideas above seem to suggest that rather than the dynamic breathing, the most important aspect—and perhaps the only necessary aspect—is to accurately reproduce the dimensions of the limiting windows during the diffusion process. Krishna and van Baten suggested that the framework flexibility can be effectively taken into account by performing fixed framework simulations with a time-averaged structure obtained from the simulations of an empty framework [85]. However, for tight fitting molecules diffusion should be more sensitive to wider pore conformations, making simple averaging less reliable [87]. Awati et al. performed a systematic investigation of methane diffusion in 8-ring zeolites, and demonstrated that the simple time-averaging recipe is not generally accurate [88]. It depends on the number of distinct directions in a fluctuating 8-ring window and the degree of harmonicity of the atomic vibrations. For example, for a given forcefield, CHA exhibited almost harmonic vibrations, the time-averaged structure was similar to the rigid structure, and simulations with the rigid time-averaged structure did not capture any effects of flexibility. For other structures with anharmonic vibrations (e.g., ERI) the flexibility effects could be partially captured, but were still not in full agreement with the flexible framework simulations.

Awati et al. [88] introduced two new methods in which the flexible structure is approximated by a set of discrete rigid structures (snapshots) obtained from a MD simulation of an empty framework. In the first “changing snapshot” method, MD simulations of diffusion are performed in a rigid framework, and the rigid structure is replaced with a new random snapshot at a frequency corresponding to the breathing motion of the windows. In the second method, TST was used, averaged over a distribution of snapshots. Both methods gave excellent predictions of CH4 diffusivities in small-pore zeolites as compared to fully flexible simulations, while providing 2–3 orders of magnitude speed-up in computational cost. In cases when reliable force fields for performing flexible framework simulations are not available, the above methods can be used with the snapshots obtained from short ab initio MD simulations (see Sect. 3.2.8).

The changing snapshot methods assume that the framework dynamics are decoupled from the dynamics of diffusing molecules. Boulfelfel et al. simulated diffusion of linear alkanes in fully flexible silica LTA zeolite using Transition Path Sampling [89]. Simulations confirmed that the framework dynamics are not affected during diffusion of methane, ethane, and ethylene. However, for larger molecules (propane, propylene, n-butane) a coupling between the zeolite vibrations and molecular diffusion was observed—the molecules exhibited longer residence times at the zeolite window, causing its distortion by  ~0.1–0.2 Å. In this case the changing snapshot methods will not fully capture the widening of the pore aperture, and the influence of this coupling on the diffusion rate.

Kolokathis et al. introduced a new curvilinear umbrella sampling free-energy method to obtain diffusivities of para-xylene and benzene in sinusoidal channels of silica MFI [91]. This system is another example of tightly fitting sorbates, and small changes in the zeolite atom positions have been shown to significantly affect adsorption energetics [90]. Using a flexible force field model, the resulting diffusivities were in reasonably good agreement with the QENS measurements [91]. In agreement with experiments, the authors found that p-xylene diffuses  ~100 times faster than benzene when sorbed at low occupancy in silicalite. The simulated activation energy for p-xylene was ~25 vs ~10 kJ/mol in experiments. Benzene encounters strong entropic barriers to translational motion at channel intersections, where it can adopt a variety of orientations. The corresponding barriers for p-xylene are much lower, reflecting the inability of its major axis to reorient within channel intersections.

Most isotherms of hydrocarbons in silicalite calculated using Monte–Carlo methods employ the monoclinic structure by van Koningsveld et al. [92] In any case it is important to note that the results can be significantly affected by the silicalite structure used, with differences of up to 20% in the calculated uptake. The differences in critical window diameters have been highlighted recently [79], and confirmed, as stated above, the importance of using accurate structures for the calculations of diffusivity regardless of whether they correspond to fixed or flexible frameworks. Moreover, this work determines the average diameter that corresponds to a flexible framework calculation, allowing comparison with a fixed framework calculation using an averaged structure obtained by selecting appropriate snapshots from the molecular dynamics calculations. Since the pore diameter fluctuation is faster than the fastest diffusing molecules, molecules influence the pore breathing as they diffuse, making it possible that the diffusing molecule contributes to widening of the pore up to the maximum width in the case of tight-fitting molecules.

3.2.5 Quality of forcefields: siliceous zeolites

Detailed comparisons of simulations with experiments require well-characterized experimental crystals without defects that influence intracrystalline diffusion. Siliceous zeolites, especially zeolites synthesized using the fluoride method, which minimize the number of defects, have been very useful. In recent years, high-silica 8-ring, small pore zeolites have been shown as reliable materials for intracrystalline diffusion studies. There have been numerous studies of LTA (ITQ-29), Si-CHA, and DDR (ZSM-58) zeolites using PFG NMR [93,94,95], membrane measurements [96, 97], gravimetric uptake [98], frequency response [93, 96, 99], breakthrough [100], Zero Length Chromatography [101], and microimaging [3, 102, 103]. Despite some variability between different materials [104], an overall good agreement obtained between the microscopic and macroscopic experimental techniques for diffusion of small molecules (e.g., methane and ethane) indicates that transport is controlled predominantly by intracrystalline diffusion, with negligible contributions from the external surface resistance. The diffusivities of small molecules in these structures scale with the size of the window aperture [93, 99]. For example, the methane diffusivities at low loadings are ~10–10, ~10–11, and ~10–12 m2/s for ITQ-29, Si-CHA and siliceous DDR, respectively, and roughly an order of magnitude lower for ethane [94, 103]. Thus, these small-pore silica zeolites provide rigorous tests for molecular simulations of diffusion, as transport is less likely to be affected by surface barriers.

Molecular simulations showed that adsorption and diffusion of mixtures in window-cage zeolites like DDR is segregated: CO2 and N2 molecules occupy window sites, and hinder diffusion of methane and other hydrocarbons, which predominantly adsorb within the cages [64, 105]. To predict equilibrium loadings of mixtures the Ideal Adsorbed Solution Theory had to be modified to account for the presence of different adsorption sites [64]. TST-based kinetic Monte Carlo simulations for mixtures [106, 107] showed that CO2 diffusion rates in DDR are only weakly affected by the presence of methane at different loadings [64]. This situation is very different from other materials in which slow molecules usually retard transport of faster molecules. The segregated diffusion has been confirmed by calculating binary Onsager coefficients—off diagonal components were very small compare to the diagonal components. These calculations allowed to predict the selectivity of the DDR membrane in good agreement with the experimental measurements over a wide range of pressures [108]. Jakobtorweihen et al. used reactive Monte Carlo to investigate confinement on propene metathesis in DDR and other zeolites; this study indicated that equilibrium compositions depend strongly on the type of zeolite [109].

Diffusion of tightly fitting CH4 molecules presents a challenge to the existing force fields, due to the high sensitivity to small changes in the intermolecular interactions. To simulate diffusion of CO2 and CH4 in DDR, Jee and Sholl used a rigid framework and adjusted the potential of interactions with the zeolite to reproduce experimental diffusivities [64]. They used a steeper repulsive potential than the standard 12-6 Lennard–Jones. Huth et al. simulated diffusivities of CO2, N2 and CH4 molecules in DD3R silica zeolite using MD simulations with a rigid framework and standard LJ potentials and obtained reasonable agreement with the experimental PFG NMR measurements [110]. Inclusion of the framework flexibility did not significantly change CO2 and N2 diffusivities but produced a tenfold increase in the CH4 diffusion, in disagreement with experiments. Bonilla et al. [111] compared infrared microscopy uptake rates of methane, ethane and propene on silicoaluminophosphate SAPO-34 (CHA) with MD and TST simulations and obtained generally good agreement with experiments when using a rigid framework model and the Dubbeldam force field. Inclusion of the framework flexibility, however, produced roughly an order of magnitude faster diffusivities, indicating an inconsistency between the force fields for framework flexibility and the methane–zeolite interactions.

For larger molecules, Boulfelfel et al. simulated diffusion of linear alkanes in fully flexible silica LTA zeolite using Transition Path Sampling [89]. These calculations used force field parameters for hydrocarbon–zeolite interactions that were developed to reproduce experimental equilibrium adsorption isotherms [112] in combination with the TraPPE force field for hydrocarbons [113] and the Hill–Sauer force field for framework flexibility [114]. When compared to PFG NMR measurements [94], the calculated diffusivities of methane, ethane and ethylene were about an order of magnitude lower, while the diffusivity of propylene was two orders of magnitude higher [89]. These discrepancies demonstrate limitations of the forcefields that were developed to reproduce equilibrium properties for predicting diffusion of tightly fitting molecules.

Accurate prediction of the experimental window aperture is very important. Bermudez and Sastre [79] demonstrated significant differences in calculated zeolite channel dimensions when using different forcefields for the framework flexibility. To improve the predictions for small pore (8-ring) zeolites Boulfelfel et al. modified the popular Hill–Sauer force field for silicas and zeolites [117] using fully periodic dispersion-corrected DFT calculations for reference silica zeolites [116]. The new force field provides reliable unit cell and window dimensions and even predicts the monoclinic-to-orthorhombic phase transition in siliceous MFI zeolite [118]. For the zeolite–adsorbate interactions Fang et al. developed a first-principles-derived force field for CH4 in silica zeolites by using periodic DFT calculations and the DFT/CC (coupled cluster) method [115]. This force field is fully predictive, without using experimental data. It reproduces very accurate DFT/CC energies not only in the most energetically favorable configurations, but also in the transition states, such as, e.g., 8-ring windows. Figure 5 shows that diffusion of tightly fitting molecules is very sensitive to the intermolecular interactions in the transition state configurations. Very small changes can lead to diffusivities that differ by orders of magnitude. Simulations with the new CCFF force field and fully flexible framework predict reliable diffusivities for the reference 8-ring zeolites (Fig. 5). Recent extensions of the DFT/CC approach to larger alkanes and alkenes show good agreement with experiments for both adsorption and diffusion, and transferability across a range of siliceous zeolites [119]. Force fields based on high-level first-principles calculations hold considerable promise for making molecular simulations of adsorption and diffusion more predictive, and obtaining better agreement with experimental data.

Fig. 5
figure 5

Copyright 2018 American Chemical Society

CH4 self-diffusion coefficients at low loadings and 300 K in all-silica zeolites ITQ-29, CHA and DDR computed using a first-principles derived CCFF force field for methane–zeolite interactions [115] and modified Hill–Sauer force field for flexible zeolite frameworks [116], compared to 15 other force fields as a function of the methane–zeolite oxygen Lennard–Jones interaction parameter. Experimental self-diffusion coefficients are shown as horizontal lines [94]. Adapted with permission from Fang et al. [115].

3.2.6 External surface barriers

The external surface of zeolite crystals can also play an important role in transport into and out of a zeolite. While control over synthesis has made it possible to produce crystals that are largely free of internal defects, particularly for fluoride-assisted siliceous zeolite synthesis, this is not the case for the external surfaces [120]. The importance of external surface barriers for diffusion is well established experimentally. For example, multiscale simulations of ethylation of benzene in ZSM-5 composites would only agree with experiments, both qualitatively and quantitatively, if external surface barriers were included [121]. Nevertheless, the nature of such barriers at the molecular level remains elusive. The influence of surface barriers on the overall mass transfer in some cases is so large that it appears equivalent to almost complete blockage of most pore entrances to the diffusing molecules [122,123,124].

Surface barriers have also been invoked as another mechanism for explaining discrepancies in diffusivities from different measurements. Zeolite A is one of such systems, for which many experimental observations are still poorly understood on the molecular level. For example, corrected diffusivities for CO2 in 4A zeolite crystals of different origin differ widely (10–15–10–13 m2/s at room temperature), while the activation energy is essentially constant [125]. Possible explanations are differences in the cation distributions and formation of surface barriers, due to different dehydration conditions of the crystals [125]. MD simulations do show that cations can migrate upon diffusion of other molecules. For example, migration of sodium cations has been observed in MD simulations of water in mordenite [126] and 4A zeolites [127], methanol in NaX and NaY zeolites [128], and CO2 in 4A zeolite [129]. MD simulations using the COMPASS force field [130] predicted CO2 self-diffusivities in pure NaA zeolite of about 10–10 m2/s, and decreasing with increased K+ content for NaKA zeolites [129]. This agrees with experimental trends for the NaKA zeolites, but the simulated self-diffusivities were much higher than typically measured CO2 diffusivities in 4A crystals [125].

Water self-diffusion has been simulated in cationic LTA-type zeolites (e.g. 4A, 5A), and generally good agreement with microscopic PFG NMR and QENS measurements was obtained using different forcefields [131,132,133,134,135,136,137]. Water diffusion in 4A is fairly fast (~ 10–12–10–10 m2/s) with substantial loading dependence [138]. The agreement with various simulations was usually within an order of magnitude or better for the self-diffusivities, and much better for the activation energy. Simulations showed that interfacial charge distribution and different terminations of zeolite LTA can affect the intracrystalline diffusivities [139]. The influence of surface chemistry on zeolite membrane permeability has been explored in nonequilibrium MD simulations of water transport through an LTA membrane [135]. About three times lower fluxes were predicted for the hydrophobic membrane, while for the hydrophilic surface the predicted permeabilities were in good agreement with the experimental values [135].

At the same time a spectacular disagreement between equilibrium MD simulations (~ 10–9 m2/s) and experimental uptake measurements [140] (~ 10–15 m2/s) has been reported for water diffusion in hydrophobic silica MFI with various levels of hydrophilic defects [141]. The disagreement has been attributed to surface barriers. While the surface barriers can certainly influence transport, an underestimation of water diffusivities in gravimetric vapor uptake measurements [140] is likely a major factor as well. Earlier PFG NMR measurements and MD simulations reported similarly fast water diffusivities (~ 10–10–10–9 m2/s) in MFI [142, 143].

External surface barriers are not the only mechanism that slows transport. Another possibility is internal barriers in polycrystalline materials. Simulations of xenon diffusion at the zeolite crystal interfaces of NaY were performed to understand the differences between QENS and PFG NMR measurements [144].

The importance of surface resistance versus intracrystalline diffusion increases with decreasing crystal size. Methods for modeling transport that include interfacial resistances continue to be an active topic [145,146,147,148,149]. Simulations of hexane diffusion in hierarchical pentasil zeolite nanosheets showed decreased diffusion rates due to tortuous paths of molecules in microporous regions at low loadings [148], in agreement with ZLC measurements.

A recent study by Knio et al. has used equilibrium MD simulations on 2 nm MFI zeolite nanosheets to measure surface resistances for light gases and found that they correlate with the heats of adsorption of those gases [149]. Furthermore, they showed that adsorbates with high heats of adsorption, like CO2, have lower diffusivity in nanosheets, whereas the diffusivity of H2, which has a lower heat of adsorption, was similar to that in the bulk zeolite.

3.2.7 Large-scale screening of materials

Screening of zeolites for diffusion based separations using MD methods is computationally demanding, and standard MD is limited to fast diffusing molecules (> ~10–12 m2/s) [105, 150]. Moreover, realistic predictions of separation performance of adsorbents and membranes often require transport diffusivities of mixtures at high loadings, rather than the more easily calculated self-diffusivities. Rao and Coppens predicted mixture diffusion of CH4/C3H8 through a silicalite membrane and the permeation of N2/CH4 through a SAPO-34 membrane, using a theory based on the MS equations, but without need of intermolecular exchange coefficients or empirical relations, in contrast to the typical use of these equations [151]. Several simplified approaches have been developed that replace expensive predictions of transport (or corrected or MS) diffusivities with predictions of self-diffusivities [105, 152]. Further successful approximations, in some cases, employed lattice kinetic Monte Carlo methods to simulate diffusion [40, 64]. These approaches were successful to predict the performance of zeolite membranes for CO2/CH4 separations in very good agreement with experiments [105, 108].

For large-scale screening, TST-based methods with approximately determined energy landscapes have been applied [153]. Novel algorithms for automated determination of diffusion pathways may open possibilities of using kinetic Monte Carlo with TST-based rate constants, to explore diffusivities in a much wider range of structures [154].

3.2.8 First principles quantum chemical methods

Above, we discussed the use of quantum chemical calculations (usually DFT) to aid in force field development. Owing to ever-increasing computing power, another approach is to use ab initio (or first-principles) molecular dynamics (AIMD), in which forces are calculated by DFT at every time step. AIMD simulations are currently limited to typically 5–50 ps long trajectories, which is not enough to calculate diffusivities in zeolites from the molecular mean-square displacements, but sufficient to characterize possible effects of the flexibility of the framework, provide insights into diffusion mechanisms [155,156,157], and calculate free energy barriers [158,159,160].

AIMD showed that significant window distortion of almost ~ 1 Å is the mechanism of kinetic separation of ethylene and ethane in siliceous ITQ-55 zeolite. ITQ-55 zeolite has very narrow windows, with the smallest windows only ~ 2.1 Å wide [157]. AIMD calculations demonstrated how diffusing molecules can “brace the window open” (Fig. 6). Recent first-principles MD simulations showed that xylene molecules create significant structural distortions (up to ~ 0.5 Å) when passing through channels of silica MFI and MEL zeolites [160].

Fig. 6
figure 6

Reproduced from Bereciartua et al. [157] with permission from AAAS

Distributions of the minimal 8-ring window size for the empty structure of siliceous zeolite ITQ-55 (left) and the ITQ-55 structure with ethylene molecules constrained to the center of the 8-ring window (right), as calculated from DFT molecular dynamics simulations at 300 K.

Recently, it has become feasible to use AIMD to calculate free energy barriers to diffusion. Mace et al. used AIMD and the “blue moon” constraint method [161] to calculate free energy profiles for diffusion of CO2 and N2 in NaKA zeolite [158]. Free energy barriers calculated from constrained AIMD were used in kinetic Monte Carlo simulations to calculate uptake rates in NaKA zeolite crystals, including the effects of the Na+ /K+ cation ratio and the presence of surface skin-layer defects. [162]. The simulated uptake rates were in qualitative agreement with the role of surface defects in A zeolites [125], yet they were still much faster than those measured experimentally. Kumar et al. calculated free energy barriers for selective diffusive transport of xylenes in MFI, MEL zeolites, and their intergrowth using AIMD and umbrella sampling [160].

Recently, a molecular “trap door” mechanism has been proposed to explain CO2 diffusion in small-pore cationic zeolites with cations blocking 8-ring windows [163]. Strong temperature-dependent adsorption of small guest molecules in such zeolites has been studied by (static) DFT calculations using Nudged Elastic Band energy profiles [163]. A proposed model explains this behavior by thermal fluctuations of cations blocking the pore windows and acting as a “pore-keeping” group [164]. Coudert and Kohen used AIMD to study these systems [159]. They demonstrated that a large-amplitude thermal motion of the cations is enough to allow CO2 diffusion in Na-RHO zeolite, while blocking larger nonpolar molecules such as methane [159].

3.2.9 Outlook

Molecular modeling and simulation have been very successful in unraveling molecular-level details of diffusion in zeolites, providing insights into the loading dependence of diffusion and multicomponent behavior, and predicting diffusion for a far greater number of zeolite structures and molecules than it has been (and will be) possible to do experimentally. Detailed quantitative comparisons with experiments, however, are limited by the availability of validated experimental measurements and the complexity of real zeolite structures, their compositional variations and the level of defects, such as surface barriers. Complexities increase as the size of the diffusing molecules becomes comparable to pore sizes. For the experimentally most studied system, siliceous MFI, the agreement is quantitative for reasonably sized linear hydrocarbons, both for self as well as for the transport diffusivities. Tightly fitting molecules, such as aromatics, still present considerable challenges, although recent simulations of self-diffusion agree within a factor of two or so with the available microscopic experiments. For smaller pore size windows, such as 8-ring zeolites, the results are heavily dependent on the size of the window (silica < aluminophosphates < aluminosilicates), and the details of the force field, especially the distance between the diffusing molecule and the oxygen atoms in the zeolite window. Most of the force fields in the zeolite literature were developed to model equilibrium adsorption, which in most cases requires good sampling of the most energetically favorable configurations. Very few force fields account for higher energy configurations, characteristic of diffusing molecules in the transition state. Framework flexibility plays a role for tightly fitting molecules, and it tends to accelerate diffusion. Recent advances in modeling diffusion include the use of first principles-based methods. DFT and AIMD calculations hold considerable promise for the development of more accurate and transferable force fields and for first-principles-based, free energy methods. More reliable calculations of energy landscapes, either by improved force fields or by DFT, in combination with an automated detection of diffusion pathways and kinetic Monte Carlo methods should make it possible to predict reliable diffusivities of simple molecules in a much wider range of zeolites. At the same time, quantitative comparisons with experiments may still be hampered by the poorly understood nature of the surface barriers, internal grain boundaries and other defects in real materials, which tend to be specific to particular zeolite topologies and compositions and difficult to generalize.

3.3 Microporous metal–organic frameworks

3.3.1 Introduction

Metal–organic frameworks (MOFs) are porous crystalline solids with tunable pore environments. MOFs are constructed from metal or metal-oxide clusters, referred to as nodes, that are connected by organic species referred to as linkers. The diversity of nodes and linkers, in addition to their various connectivities and spatial arrangements can lead to micropores, mesopores, and cages of many shapes and sizes. High internal surface areas and adsorbate uptakes make MOFs promising materials for gas storage and separation applications, while metal oxide nodes and other opportunities for incorporating catalytic sites create potential applications in catalysis. The tremendous design space of MOFs and their prevalent syntheses lead to many possible materials [165]. Like zeolites, the high crystallinity of many MOFs has made them model systems for fundamental studies of molecular diffusion under confinement. Quasi-elastic neutron scattering (QENS), pulsed field gradient NMR (PFG-NMR), infrared microscopy (IRM), and other techniques have been used to measure diffusion in MOFs, while molecular simulations have been used to predict and corroborate experimental measurements, as well as provide molecular level insight into the diffusion process.

Here, we provide a perspective on connections between simulation and experiment for diffusion in microporous MOFs. There have been numerous computational contributions to this topic, and we highlight a subset of these that connect directly to experimental diffusion measurements. We first consider single component diffusion in MOFs that can be approximated as rigid, followed by a discussion of multicomponent diffusion in MOFs, with subsequent discussion on the role of framework flexibility on diffusion. We conclude with an outlook on the interaction of simulation and experiment.

3.3.2 Single component diffusion in rigid metal–organic frameworks

The first measurements of adsorbate diffusion in a MOF were published in 2006 by Stallmach et al., [166] but these measurements were predated by a set of predictive MD studies from Sarkisov et al., [167] Skoulidas [168], Skoulidas and Sholl [169], and Yang and Zhong [170]. These simulation studies drew from previous literature for modeling diffusion in other microporous materials to select reasonable force field parameters. All of these early studies considered the MOF frameworks as rigid and placed Lennard–Jones sites on the framework atoms with Lennard–Jones parameters taken from the DREIDING [171], UFF [172], or OPLS-AA [173] force fields. While the choice of force field in each of these studies was influenced by their accuracy in predicting adsorption isotherms, for example by Vishnyakov et al. [174] and Düren et al., [175] their accuracy in predicting transport properties would not be established until later comparison with experiments.

Early simulation studies mainly focused on the set of isoreticular MOFs (IRMOFs) and HKUST-1, which had been shown to have high gas storage capacity [176, 177]. IRMOF-1 is constructed from Zn4O nodes and 1,4-benzenedicarboxylate (BDC) linkers [176] with cage diameters of 10.9 and 14.3 Å [167]. By varying the choice of linkers while retaining the Zn4O nodes, a series of IRMOFs have been synthesized with tunable pore volumes [176]. HKUST-1 is comprised of Cu3(benzene-1,3,5-tricarboxylate)2(H2O)3 monomers that form a 1 nm diameter 3-D channel system [177]. Sarkisov et al. performed MD simulations for methane, n-pentane, n-hexane, n-heptane, cyclohexane, and benzene at low loadings (~ 1.25 molecule per cage) and calculated their self-diffusion coefficient. They compared the results to the diffusivities of these molecules in the zeolite silicalite as a benchmark [167]. The simulated diffusion coefficients for methane and n-hexane were similar to those in silicalite (3.08 × 10–8 m2 s−1 in IRMOF-1 compared to 0.91 × 10–8 m2 s−1 in silicalite) [167]. Stallmach et al. later compared their measured n-hexane self-diffusivities in IRMOF-1 with those of Sarkisov et al. and found good quantitative agreement (3.2–4.1 × 10–9 m2 s−1 from PFG NMR compared to 2.22 × 10–9 m2 s−1 from simulation), providing the first evidence that MD simulations in MOFs could correspond to experimental measurements. Sarkisov et al. noted that while benzene diffusion was not observed on the timescale of MD simulations, potential energy contours indicated that benzene could still diffuse in IRMOF-1 on experimental time scales [167]. This prediction was corroborated by Stallmach et al. who were able to measure a benzene diffusivity slower than n-hexane [166]. In addition to predicting diffusion coefficients, the MD simulations were able to quantify structural and energetic profiles for adsorbates during the diffusion process that were consistent with experiments.

The loading of adsorbates in zeolite micropores has been shown to affect diffusion, and the first studies of loading dependent diffusion coefficients in MOFs were performed by Skoulidas [168] for Ar in HKUST-1 and Skoulidas and Sholl [169] for Ar in HKUST-1 and IRMOFs as well as light gases in IRMOF-1. Skoulidas and Sholl simulated light gases in IRMOF-1 at various loadings, and the simulated methane self-diffusion coefficients at low loadings were consistent with self-diffusion coefficients of methane at similar conditions by Sarkisov et al. [167] Simulated methane self-diffusion coefficients in both MOFs were approximately an order of magnitude lower than measured with PFG NMR which has been attributed to methane being sensitive to defects in the micropore structure [166]. Skoulidas reported argon self-diffusivities in HKUST-1 from infinite dilution to loadings corresponding to 104 bar and found that the diffusivities decreased monotonically above 100 bar and were similar to those in the zeolite ITQ-7 [168]. Skoulidas and Sholl were able to further compare argon self-diffusivities across a series of IRMOFs and HKUST-1 which all decreased precipitously at high loadings [169].

In contrast to the IRMOF series and HKUST-1 which have distinct cages, MIL-53(Cr) [178] and MIL-47(V4) [179] form 1-D diamond shaped channels. MIL-53(Cr) and MIL-47(V) contain either Cr or V oxide nodes connected by terephthalate linkers. MIL-53 (Cr) has hydroxyl groups on each linker, while MIL-47(V) does not [178, 179]. Rosenbach et al. used a combination of molecular dynamics simulations and quasi-elastic neutron scattering (QENS) experiments to study methane diffusion in these MOFs as a function of loading. The MOFs were kept rigid in the simulations, based on XRD evidence in MIL-53(Cr) indicating minimal structural changes upon adsorption of methane [180]. Diffusivities were larger in MIL-47(V) than in MIL-53(Cr), which was attributed to the hydroxyl groups in MIL-53(Cr) adsorbing methane, which diffuses as 1-D hops between nodes [181]. Methane self-diffusion coefficients are qualitatively similar between QENS and MD, with the most significant discrepancies seen at low loading.

Dihydrogen diffusion was studied in the same MIL MOFs by Salles et al. using MD and QENS. They found a similar diffusion mechanisms as methane, where adsorption to hydroxyl groups in MIL-53(Cr) increased diffusion activation energies and introduced 1-D hops that were not observed in MIL-47(V) [182, 183].The importance of the force fields on simulating diffusion was shown by Salles et al., who used a one-site H2 model, a two-site H2 model, and a one-site H2 model with the hydrogen-MOF interactions modified [182]. Figure 7 shows that self-diffusion coefficients calculated with the modified one-site model correlate the best with experimental QENS measurements. This modified force field performed poorly in predicting H2 adsorption isotherms, cautioning the reader that a force field that works well for diffusion does not necessarily guarantee accurate adsorption predictions (and vice versa) [184]. In contrast to the monotonically decreasing dihydrogen diffusivities from Salles et al., [182, 183] dihydrogen diffusion in a suite of IRMOFs simulated by Yang and Zhong showed that hydrogen diffusivity increased slightly with loading, but then sharply decreased [170]. Simulated dihydrogen diffusivities in the [Zn(bdc)(ted)0.5]n MOF decreased by less than an order of magnitude from 0 to 100 bar, and quantum diffraction corrections systematically increased the simulated diffusion coefficients [185]. Dihydrogen diffusion coefficients were similar to those simulated by Skoulidas and Sholl in IRMOF-1, especially at high pressures [169, 185]. Differences in observed dihydrogen diffusion profiles with loading in different MOFs indicate that node identity and pore structure can influence the qualitative trends of the diffusion coefficients.

Fig. 7
figure 7

Copyright 2009 American Chemical Society

Comparison of MD (empty symbols) and QENS (filled symbols) for H2 self-diffusion coefficients in MIL-47(V) (a) and MIL-53(Cr) (b). Each colored empty symbol represents a different force field used to model the framework nonbonded interactions. Reprinted with permission from Salles et al. [182].

Alkane diffusion in MIL-47(V) was considered by Jobic et al., [186] Deroche et al., [187] and Rives et al. [188] for various adsorbate loadings and variations in chain length. C2–C5 n-alkane diffusion was measured with QENS and simulated with MD, with simulated alkane diffusion coefficients predicted to decrease with chain length, as well as adsorbate loading [186]. QENS measurements found that n-butane diffusivities were larger than propane at low loadings, which was not observed in MD, and this was attributed to the inability of n-butane to freely rotate within the MIL-47 pore [186]. The authors commented that this may be due in part to the “absence of momentum transfer from the alkanes of MIL framework… because (the framework) is considered fixed” [186]. Deroche compared the orientation of diffusing C5–C9 alkanes in MIL-47, where longer chain alkanes lie parallel along the pore [187]. Diffusivities of alkanes up to C16 measured by QENS and simulated by MD show good quantitative agreement up to C10. However, MD overpredicts diffusion coefficients of C12 and C16 by over an order of magnitude and does not show as severe of a decrease with chain length as measured with QENS [188].

HKUST-1 is a MOF with interconnected ~ 9 Å cages that are surrounded by smaller 6 Å chambers and 4.6 Å pockets [189, 190]. Alkane diffusion in HKUST-1 was simulated with MD in conjunction with infrared microscopy [191]. The loading of alkanes was shown to strongly influence diffusion, as light alkanes mostly occupy the smaller pockets at low loadings and do not adsorb preferentially in the large cages until high loadings. At low loadings, Maxwell–Stefan diffusivities increased with loading, and then sharply decreased. Adsorption isotherms from grand canonical Monte Carlo simulations were used to calculate thermodynamic factors, which were used to corroborate the observed inflections in the transport diffusivities [189]. PFG NMR experiments were also performed by Wehring et al. for small alkanes in HKUST-1 to compare with interference microscopy [192]. In addition to comparisons based on diffusion coefficients, NMR relaxation times were also simulated to compare with experiment, providing additional insight into rotational motions of adsorbates. For n-butane in HKUST-1, simulated self-diffusion coefficients correlate well with PFG NMR, but both are more than an order of magnitude larger than the Maxwell–Stefan diffusion coefficients obtained from IRM [192]. At higher loadings, results from PFG NMR and IRM agree with each other, but simulated self-diffusivities are still more than an order of magnitude higher than either experimental technique [192].

Additional comparisons between MD simulations and PFG NMR for alkanes in MOFs include work by Ford et al. [193] in which the loading of adsorbates was carefully controlled to represent experimental conditions. Methane, ethane, n-hexane, n-decane, n-dodecane, and benzene loadings were fixed between 5–9 carbon atoms per IRMOF-1 cage, which corresponds to the 20–25% loadings used in the PFG NMR experiments. Self-diffusion coefficients decreased by two orders of magnitude with increasing chain length and were similar to bulk liquid diffusion coefficients for the same adsorbate. The comparisons between MD and PFG NMR from Ford et al. [193] are summarized in Fig. 8. For alkanes larger than methane, there is very good agreement between simulation and PFG NMR. MD simulations from Sun et al. found similarly decreasing alkane diffusivities in IRMOF-1 up to n-butane [194]. The largest discrepancy in Fig. 8 is observed for methane, where simulated self-diffusion coefficients are nearly an order of magnitude lower than experimentally measured ones. This was attributed to the rapid exchange of methane between inter- and intra-crystalline spaces in the PFG NMR measurements. Standard MD with periodic boundary conditions simulates only the intracrystalline diffusion. Ford et al. also studied higher loadings of molecules and observed that self-diffusion coefficients were independent of loading at moderate conditions, but sharply decreased at higher loadings, as the cages reached their saturation loading [193].

Fig. 8
figure 8

Copyright 2012 American Chemical Society

Self-diffusion coefficients in IRMOF-1 from MD simulations and PFG-NMR at similar fractional loadings. Reprinted with permission from Ford et al. [193].

Kolokathis et al. [195] and Splith et al. [196] considered water diffusion in MIL-100, where hydrogen bonding leads to strong guest–guest and host–guest interactions. MIL-100 has trimers of octahedral metals and 1,3,4-benzenetricarboxylate linkers. The cages are approximately 25 and 29 Å, with windows of approximately 5.2 and 8.8 Å [196]. Initial MD simulations were performed in MIL-100(Fe) by Kolokathis et al. [195], where water diffusivities approached the bulk water diffusivity at higher loadings, similar to what Ford et al. [193] observed for alkanes in IRMOF-1. The presence of strongly bound (structural) water molecules that cap open metal sites was shown to have a significant stabilizing effect on host–guest interactions, but diffusivities with and without structural water were similar across all loadings. Water self-diffusivities in MIL-100(Al) from PFG NMR experiments and MD simulations showed good agreement, both in trend with loading and quantitative values. Splith et al. compared simulated diffusion coefficients in MIL-100(Al) and MIL-100(Fe), and found that the diffusivities in both MOFs approached similar values at high loading, with an order of magnitude difference at low loadings [196]. The difference in the two MOFs at low loadings could be attributed to differences in water adsorption at the metal nodes. Rudenko et al. used quantum-chemically derived force fields to study water diffusion in Mg-MOF-74 and showed that the more accurate description of water-MOF interactions was significant in improving the uptake predictions [197].

Transition state models have been used to understand concentration-dependent experimental diffusion in MOFs. Chmelik and Kärger utilized the window-hopping diffusion mechanism in the ZIF-8 MOF to predict concentration-dependent transport and corrected diffusivities of ethane and propane validated with IRM measurements [198]. They noted that because the windows between adjacent cages are small, jumps through windows are both rare events and independent of correlated molecular interactions. This implies that the hopping rate from a cage to a window is proportional to the number of molecules per cage. Chmelik and Kärger go on to show that the hopping rate is proportional to the quantity adsorbed and can be related to the external pressure from the equilibrium isotherm. The transport diffusivity can also be related to the hopping rate, allowing the authors to correlate transport diffusivity to the equilibrium thermodynamics of the system using transition state theory. This construction is not general for all systems, but for microporous solids where the diffusion mechanism is dominated by uncorrelated hops through windows, it is possible to predict the concentration dependence of the transport and correlated diffusivities. This technique has also been used by Lauerer et al. in ZSM-58 and DDR zeolites to understand the diffusion of binary mixtures [199].

3.3.3 Multicomponent diffusion in metal–organic frameworks

With their high degree of tunability, MOFs have received considerable attention as adsorbents for separation of mixtures. Separations using adsorbents can be based on differences in diffusion rates of different species or differences in their equilibrium adsorption. For the latter, fast diffusion is important in reaching equilibrium quickly. Few multicomponent experimental diffusion measurements have been performed in MOFs, yet the unique pore shapes that have been explored with molecular simulations hint at rich possibilities for kinetic separations. We focus here on alkane isomer separations, CO2/CH4 separations, and CO2/H2 separations, where most of the MOF-related multicomponent diffusion studies have been performed.

Alkane isomer separations can be difficult for distillation when the boiling points of the components are similar, and these separations are important in petroleum refining, so there is strong interest in finding adsorbents that can separate alkane isomers. Babarao et al. [200] as well as Krishna and van Baten [201] performed MD simulations to predict diffusion coefficients of alkane mixtures in MOFs and relate their relative mobilities with structural properties of the pores. Due to their ubiquity in petroleum refining, comparisons were made between zeolites and MOFs to establish benchmarks for separations performance. Babarao et al. [200] considered two PCN MOFs (PCN-6 and PCN-6′), and two IRMOFs (IRMOF-14 and IRMOF-13) for C4 and C5 separations. PCN-6′ has a similar topology to HKUST-1 with square pores that are 15.2 and 21.4 Å in diameter. PCN-6 is constructed by concatenating PCN-6′ and has smaller 9.2 Å triangular channels. Babarao showed that in PCN-6′ n-butane and i-butane mixture self-diffusivities were similar at the pressures and compositions predicted from GCMC, with a noticeable peak in diffusivities at moderate loadings in contrast to PCN-6 where n-butane diffusivities were greater than those of i-butane for all loadings considered [200]. Similar trends were observed for neo-pentane, n-pentane, and i-pentane, where mixture self-diffusivities of the isomers were similar in PCN-6′ compared to PCN-6 where n-pentane diffusion coefficients were larger than those of the other isomers. These results correspond with simulations performed by Krishna and van Baten where differences in diffusivities of the components were due to clustering of adsorbates in different pore regions. Krishna and van Baten calculated Maxwell–Stefan exchange coefficients for alkane mixtures in different zeolites and MOFs, which also exhibited significant cluster formations for binary mixtures [201].

The separation of CO2 from CH4 is important for purification of natural gas. Zeolites such as NaX have been considered for this separation, and the development of MOFs has led to renewed interest in this separation. Yang et al. used MD simulations and QENS to quantify both pure methane diffusion coefficients, and binary diffusion coefficients as a function of CO2 loadings in UiO-66 [202, 203]. UiO-66 is constructed from Zr6 nodes connected by benzene dicarboxylate linkers, and has alternating-size cages connected by windows. Both QENS and MD showed that when the CO2 loading is increased, the CH4 diffusivity increases before decreasing at higher loadings, as is summarized in Fig. 9 [203]. The agreement in predicted mixture self-diffusion coefficients is accurate within a factor of ~ 2–3 and correctly captures the QENS trend of higher CH4 diffusivity with increased CO2 loading [203]. Additional pure component CO2 and CH4 QENS measurements and MD simulations showed that CO2 preferentially occupies the smaller cages, which reduces the barrier for CH4 hopping, resulting in a higher diffusivity for CH4 in mixtures [202].

Fig. 9.
figure 9

a Self-diffusion coefficients of CH4 from QENS for pure gas (red) and CH4/CO2 mixtures at a fixed loading of 4 CH4 per unit cell (blue). b Self-diffusion coefficients of CH4 at a fixed methane loading of 4 CH4 per unit cell and variable CO2 loading from MD (empty symbols) and QENS (solid symbols). Reprinted with permission from Yang et al. [203]. Copyright 2018 American Chemical Society.

Recently, Solanki and Borah screened 4764 MOFs from the Computational Ready Experimental (CoRE) MOF database for separation of hexane isomers, and picked 11 of the most promising MOFs for more detailed pure component diffusivity simulations at infinite dilution to determine kinetic selectivity coefficients [204]. Solanki and Borah found that, among the different MOFs examined, the diffusivities of hexane isomers varied by more than an order of magnitude, indicating large kinetic selectivity to specific isomers. Interestingly, there were MOFs for which n-hexane does not have the fastest diffusivity. Krishna and van Baten have also used screening of diverse MOF topologies to compare diffusion processes [105, 205] where pure component self-diffusivities and Maxwell–Stefan diffusivities can be used to extrapolate diffusion behavior to mixtures. Krishna and van Baten have shown how different pore environments in the same MOF can produce “molecular traffic control” effects that influence kinetic separations [105].

3.3.4 Flexible microporous metal–organic frameworks

Much of the previous discussion focused on diffusion simulations where framework atoms were constrained to their equilibrium lattice positions. These rigid framework simulations have been shown to accurately predict diffusivities for a variety of adsorbates in different MOFs, as evidenced in the previous sections; however, some MOFs have flexible linkages, such that pore sizes and shapes can change dramatically due to changes in temperature or adsorbate loading. MOFs can have dramatic thermal expansion coefficients [206] that impact simulated activation energies of diffusion, and some MOFs can restructure their pores at higher adsorbate loadings [207,208,209]. In some cases, local fluctuations of the framework structure may affect window sizes and influence diffusion as in zeolites. In this section we highlight reports where MOF flexibility has been incorporated into the models and the results have been compared with experimental diffusion measurements.

Models that incorporated MOF flexibility have been implemented with the goal of better reproducing experimentally measured diffusivities, but the effect of flexibility in general depends on the specific MOF and adsorbates [210]. Amirjalayer et al. [211] used an extension of the MM3 forcefield to include the effect of framework deformation upon benzene adsorption in IRMOF-1 and compared activation energies and diffusion coefficients to rigid framework models. For the same system, Greathouse and Allendorf extended the CVFF force field to include potentials for the IRMOF-1 framework and reproduced the negative thermal expansion in IRMOF-1 as well as the power spectra, but benzene diffusion coefficients were at least an order of magnitude lower than those predicted by Amirjalayer et al. [211] despite both groups reporting nearly identical activation energies. Ford et al. [212] compared flexible and rigid forcefields for benzene in IRMOF-1 and found that diffusion coefficients for both models were similar to within a factor of 2 and showed excellent agreement with PFG NMR measurements from Stallmach et al. [166] (2.8 × 10–9 m2 s−1 MD flexible, 1.9 × 10–9 m2 s−1 MD rigid, 1.9 × 10–9 m2 s−1 PFG NMR). Comparisons between flexible and rigid IRMOF-1 frameworks for alkanes were within ~ 20–50% of each other, indicating that framework flexibility has a minimal effect on the diffusivity of small guest species in IRMOF-1, which does not exhibit conformational changes or gate-opening mechanisms at the windows [212].

One MOF with prominent flexibility is Zn(tbip) with tetrahedral Zn(II) cations linked by 5-tert-butyl isophthalate producing 4.5 Å 1-D pores [213]. Heinke et al. performed IR microscopy to track transient ethane concentration profiles including tracer exchange experiments to operate under quasi-equilibrium [214]. Seehamart et al. performed MD simulations of ethane diffusion in Zn(tbip) with a flexible forcefield and reported significant quantitative differences in the loading dependence of ethane diffusion depending on force field choice [208]. They showed that for a rigid framework, ethane diffusivity decreased with loading, but for the flexible force field, once the loading exceeded 6 ethane molecules per unit cell there was a sharp increase in diffusivity which was attributed to restructuring of the pores at higher loadings [208]. The role of framework flexibility in reshaping the pore environment was further established by Seehamart et al. where density maps of ethane, Fig. 10, show how diffusion bottlenecks at higher loadings widen, resulting in the faster observed diffusion [209]. Simulated ethane self-diffusion coefficients are roughly an order of magnitude higher than transport diffusivities measured by Heinke et al., which should coincide at infinite dilution [214]. Binary mixtures of CO2/ethane, methane/ethane, and CO2/methanol in Zn(tbip) simulated by Seehamart et al. [207] indicate that pore restructuring at different loadings may lead to improved loading-dependent kinetic separations.

Fig. 10
figure 10

Reproduced from Seehamart et al. [209] with permission from Elsevier

Ethane molecule location in Zn(tbip)for different framework models.

The MIL materials can also have significant flexibility effects depending on the node identity and adsorbate identity [215]. MIL-53(Cr) without guest molecules has square pores, but upon low water loading, a narrower diamond-shaped pore structure is formed. At higher loadings of water, the pores return to a square shape [216]. Simulated diffusion coefficients differ by an order of magnitude, depending on the loading and pore structure, with QENS experiments indicating that self-diffusion coefficients are lower than 5 × 10–10 m2 s−1, which agrees with simulation [216]. MD and QENS studies of CO2 in MIL-53(Cr) show a similar structural change at higher loadings, and simulated transport diffusivities agree with experiment at low loadings, but deviate at moderate and higher CO2 loadings. Similar to the alkane studies in MIL materials [181, 186, 187], 1-D hopping diffusion was observed for CO2 in MIL-53(Cr) [217].

Single-file diffusion was observed for neopentane in MIL-47(V) in simulations by Ghoufi and Maurin [218]. In single-file diffusion, as discussed in Sect. 2, the mean-squared displacement of molecules scales with the square root of time [24]. For MIL-47(V) and its narrow pores, neopentane mean-squared displacements at long time scales in a >400 Å supercell possessed all the characteristics of single-file diffusion. They then observed that the obtained mobilities corresponded to QENS measurements.

MOFs with imidazolate linkers and tetrahedral metal centers (referred to as zeolitic imidazolate frameworks or ZIFs) have tunable pore environments that resemble zeolite topologies [219]. Early MD studies by Liu et al. [220] and Pantatosaki et al. [221] simulated diffusivities in ZIF materials with rigid frameworks while Battisti et al. [222] considered flexible frameworks for light gas diffusion in a suite of ZIF materials. Hertag et al. used both MD simulations and IR microscopy to study methane diffusivity in ZIF-8, which has 3.4 Å windows. Measured H2/CH4 selectivities by Bux et al. showed that although the kinetic diameter of CH4 is larger than the window diameter, CH4 was still able to diffuse through the material [223]. Hertag et al. found that modeling the opening of windows in ZIF-8 using an appropriate force field allowed for impressive agreement between simulated self-diffusion coefficients and IR microscopy transport diffusivities at low loadings [224]. Haldoupis et al. also considered CH4 as well as CO2 diffusion in ZIF-8 using different force field parameters and found large differences in predicted diffusion coefficients depending on the choice of force field [225], which was also reported by Hertag et al. [224] Pantatosaki et al. used MD simulations to study CO2 and CH4 diffusion in ZIF-8 where both transport diffusivities at higher loadings and self-diffusivities at infinite dilution agreed with IR microscopy and PFG NMR measurements [226]. Chokbunpiam et al. also reported on the importance of using force fields that include framework flexibility to properly account for window opening when studying diffusion of ethane in ZIF-8 [227]. Parkes et al. and Ford et al. showed that the framework flexibility had a minimal effect on predicting diffusion for molecules much smaller than the window size [212, 228].

Diffusion of larger molecules have also been performed in ZIF-8 and related materials with flexible windows. Gee et al. looked at methanol and ethanol diffusion in ZIF-8 and ZIF-90 and made comparisons to PFG NMR [229]. Alcohols were found to adsorb at carbonyl groups in the frameworks, resulting in diffusion coefficients of ~ 10–12 m2 s−1 at loadings close to 50% of saturation. These diffusivities were at least an order of magnitude lower than those measured by PFG NMR, but MD and PFG NMR both predicted more rapid diffusion in ZIF-90 that ZIF-8. Krokidas et al. predicted single component diffusion coefficients for methane, ethane, propane, and propylene in ZIF-8 for propane/propylene separations, obtaining excellent agreement with experimental studies (measured methane self-diffusion coefficients ~ 1–1.5 × 10−10 m2 s−1 compared to 1.3 × 10–10 m2 s−1 from MD) [230]. Good agreement of the activation energies was observed as well (measured methane diffusion activation energy 26.6–38.8 kJ mol−1 compared to 30.0 kJ mol−1 from MD) [230].

To simulate diffusivities for molecules with larger kinetic diameters, where the barrier to pass through a ZIF-8 window is more pronounced, the probability that molecules diffuse through windows may be too infrequent to sample on the typical timescale of MD. This regime is typically < 10–12–10–11 m2 s−1. For slowly diffusing molecules, transition state sampling can access infrequent diffusion events, such as for larger molecules in ZIF-8. Verploegh et al. [231] used transition state theory to calculate diffusion coefficients at infinite dilution for 15 different adsorbates in ZIF-8. Figure 11 summarizes these diffusivities as a function of adsorbate size, with comparisons to experimental measurements. The agreement between simulation and experiment is good for most adsorbates, but discrepancies are seen for the three largest adsorbates. Self-diffusion coefficients of alkanes were mildly dependent on loading, and transport diffusivities increased by more than an order of magnitude with loading for C3 and C4 hydrocarbons [231]. A recent article by Verploegh et al. has extended these correlations to additional ZIFs and found a similar relationship between diffusivity and molecular diameter that also depends on the average window size [232]. These calculated diffusivities exhibited good agreement with the measurements by Chmelik for C3 alkanes [233]. Han et al. considered point defects in ZIF-8 and found that linker vacancies and dangling linkers can increase hopping rates, so that the effect of these defects must be considered [234]. A recent study by Zheng and Maurin [235] on propane/propylene kinetic separations in ZIF-8 used high temperatures to simulate diffusion coefficients and extrapolated down to standard conditions using an Arrhenius plot; these results also compared well to the measurements by Chmelik et al. [233] Zheng and Maurin showed that mechanical pressure exerted on ZIF-8 could improve propane/propylene separations by reducing the window size and significantly decreasing propylene diffusion relative to propane.

Fig. 11
figure 11

Copyright 2015 American Chemical Society

Infinite dilution self-diffusion coefficients in ZIF-8 and LTA at 35 °C. Reprinted with permission from Verploegh et al. [231].

3.3.5 Outlook

Molecular simulations of diffusion in MOFs have contributed key physical insights that would have been difficult to obtain by experimental measurements alone. Comparisons of simulation and experimental measurements have proven essential for developing accurate force fields that realistically model diffusion in MOFs. This is an ever-evolving effort that necessitates continued cooperation between modeling and experiment. We foresee that with decreasing computational cost of quantum mechanical simulations, force field development will be accelerated, especially for systems that are not well described by “generic” force fields, such as UFF. Advances in data sciences and the field of machine learning provide additional techniques for force field development [236, 237], which we expect to only accelerate in the future. Most collaborative modeling/experimental studies have used simulation to provide additional molecular-level insights about the diffusion process, but simulations can also be used in a predictive mode and there are now early studies using simulation to screen large numbers of hypothetical MOFs for their diffusion properties. This approach has been tremendously successful for gas storage [238], and we anticipate more efforts in this direction for diffusion. Recently, Bukowski and Snurr [239] performed MD simulations in over 30 Zr MOFs at varying fractional loadings for alkanes. Combined modeling and experimental efforts to tackle multicomponent diffusion in MOFs are relatively sparse, despite the significant interest in MOFs for separations. We foresee a significant increase in the number of collaborative multicomponent studies of diffusion in MOFs in the near future.

4 Amorphous Mesoporous Materials

4.1 Introduction

4.1.1 General features of diffusion in mesoporous materials

According to IUPAC, mesoporous materials have pores between 2 and 50 nm in size, although pore size has a precise meaning only when the geometrical shape is well defined, e.g., the diameter of a cylindrical pore or the pore width of a slit, defined as the distance between the two opposite walls. This nomenclature is based on the way in which nitrogen at its normal boiling point (77 K) adsorbs on, and desorbs from, pore networks. Another property is the porosity, defined as the ratio of the void volume to the total volume of voids and solids.

Such mesoporous materials are ubiquitous in industrial catalysis and separation processes, such as pressure-swing adsorption, chromatography, and membrane separation. The principal reason for their use is that they have a high specific surface area, often several 100 m2 g−1, combined with effective diffusivities that are, typically, only moderately reduced from their values in the bulk (one or two orders of magnitude, at most), compared to microporous materials (usually many orders of magnitude). Still, the confinement effects induced by the narrow pores can be significant and non-trivial, due to a complex combination of factors, described in this Section. Therefore, diffusion limitations are common in mesoporous catalysts, because of their high surface area that needs to be serviced by diffusion. Thus, the ability to model diffusion in mesoporous materials and understand confinement effects is crucial to describe diffusion-driven or limited processes and, ultimately, to design new materials with optimized pore spaces that enhance the overall catalytic or separation performance. In practice, industrial mesoporous materials are often bimodal, containing both nano- and macropores or wide mesopores, all of which need to be optimized in hierarchical porous catalysts and adsorbents. This is discussed further in Sects. 4.3 and 5.

The term “nanopores” is increasingly used for pores that are smaller than ~100 nm. As the term “nano” implies properties different from single molecules and the bulk, the term nanopore might be most useful in the context of describing phenomena in pores that even qualitatively deviate from those expected in the bulk. This is the definition we adopt here, in the context of modeling diffusion in mesopores. Otherwise, we can simply consider bulk diffusion in the individual pores, and then use techniques, described below, to evaluate the overall, effective diffusivity, \({D}_{eff}\), to account for the porosity and the pore network structure.

Here, we will concentrate on the challenge of modeling diffusion in mesopores in which there are nano-confinement effects on diffusion. The nature of these nano-confinement effects generally differs from that in micropores. As shown in Sect. 3, diffusion in micropores tends to be activated, following the Arrhenius law as a function of temperature, \(T: D\sim {e}^{-{E}_{act}/(RT)}\), with species “hopping” from one adsorption site to another with activation energy, \({E}_{act}\), akin to defect diffusion in solids (\(R\) is the ideal gas constant). In mesopores, this is not the case; a fluid state is maintained, in which molecules can pass each other continuously, and, thus, the temperature dependence is a mild power law, as in the bulk, where, generally: \(D\sim {T}^{3/2}\). While fundamentally dissimilar from diffusion in micropores, there are, nevertheless, mechanistic differences to diffusion in the bulk as well, due to the presence of the enclosing mesopore walls. This is illustrated in Fig. 12.

Fig. 12
figure 12

Reproduced from Coppens [253] with permission from Taylor & Francis

Schematic representation of the different diffusion mechanisms in a porous medium: a molecular diffusion, Knudsen diffusion, and surface diffusion, all relevant in mesoporous materials; b configurational diffusion, which is typical in zeolites, but could also occur for larger probe molecules in mesoporous materials; c diffusivities in different diffusion regimes [253, 254].

First, surface diffusion [240,241,242,243,244,245,246,247,248,249,250,251,252] of adsorbed species may play an important role in the overall transport of fluids through mesoporous materials, because of their high surface area, especially when there are strong interactions of adsorbed species with the wall. Surface diffusion bears similarities with diffusion in zeolites. Thus, it is typically activated, and shows similar dependencies on concentrations. In liquids, for fluids at high pressures and low temperatures, the contribution of surface diffusion, in parallel to that in the central body of the pore, may be significant. However, it can be difficult to isolate its contribution from other transport mechanisms experimentally. Hence, it is the least well understood and characterized form of diffusion in porous materials. We will return to it in Sect. 4.4.1.

Second, the available space for diffusion may have to be corrected for the fraction of the pore space near the wall that is inaccessible to diffusing probes. In the simplest case, this is equivalent to applying a reduced porosity and/or a reduced pore diameter in expressions for the effective diffusivity, which, for narrow pores, can be quite substantial.

Third, diffusing molecules interact with the walls—they collide with them, temporarily adsorb, then desorb and continue their trajectory in the bulk (see Fig. 12a); this effect strongly depends on the molecule–wall interactions, similar to what has been described in Sect. 3. The force fields will affect the molecular trajectory. Different from micropores, where species are constantly strongly affected by surrounding walls, a quickly decaying potential in the neighborhood of the wall means that there is often a central zone where diffusion is bulk-like—however, beyond the adsorbed layer, there may be an additional one or two molecular layers, where the walls exert a strong influence on the overall fluid transport (see Fig. 13). This phenomenon can be investigated using molecular dynamics, and there is considerable current interest in this, because of the importance in evaluating the performance of electrochemical devices and membrane separations, as well as for environmental, water and energy related applications.

Fig. 13
figure 13

a Subdivision of a pore in a wall region and an inner region; b fluid-wall interaction potentials in pores with different radii. Shown as an example is CH4 in a carbon nanotube (CNT). From Keil [255].

For gases, a unique nano-confinement phenomenon occurs that has been studied since the beginning of the 20th Century, and is called Knudsen diffusion [3, 256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271,272], after Martin Knudsen [273, 274], who first described it as “molecular flow” (German: Molekularströmung). When the mean free path of gas molecules, \(\lambda\), becomes larger or is of the same order of magnitude as the local channel diameter, \(d\), the frequency of the collisions between the molecules and the walls exceeds the intermolecular collision frequency. This is of enormous relevance in industrial practice, where this occurs in the majority of gas-phase processes using mesoporous materials. Therefore, we will discuss Knudsen diffusion in more detail in Sect. 4.4.2. For long cylindrical channels, the Knudsen diffusivity, \({D}_{K},\) is proportional to the channel diameter, \({D}_{K}\sim d\), and is lower than the bulk gas diffusivity, \({D}_{b}\sim \lambda\), when \(d<\lambda\) (Fig. 12c). In mesoporous materials, the pores are typically not all of the same diameter, which implies that the diffusivity changes throughout the material. Because the pore space of real mesoporous materials has a very complex geometry and the interactions with the walls are non-ideal—a far cry from the very long, smooth, cylindrical channels studied by Knudsen—this form of gas diffusion is challenging to model.

As illustrated in Fig. 12a, bulk, Knudsen and surface diffusion together combine to inform the overall diffusion behavior in single pores. Their respective contributions depend on the local, geometrical and physicochemical environment. In an equivalent resistance diagram, molecular bulk and Knudsen diffusion occur in series, and surface diffusion in parallel to those two. Bosanquet’s approximation [275] for the combined molecular bulk and Knudsen diffusion:

$$\frac{1}{D}=\frac{1}{{D}_{b}}+\frac{1}{{D}_{K}},$$
(7)

is remarkably accurate, as Pollard and Present [276] showed in 1948. Apart from molecular dynamics and kinetic Monte-Carlo simulations, there are approximate analytical theories to describe transport in pores. The best known of these is the dusty gas model (DGM) [250, 251, 277], an extension of the Stefan-Maxwell equations [278,279,280,281] for bulk fluids, which, apart from diffusion, is also able to include pressure-driven, convective flow and, as in the Onsager approach to irreversible thermodynamics [282], molecular transport induced by coupling with other than chemical potential gradients, such as temperature or electrical potential gradients. Kerkhof and Geboers [283] found some inconsistencies in the DGM and developed a new framework for multicomponent diffusion in fluids. Fortunately, two errors in the DGM almost cancel each other in most cases. Non-equilibrium adaptations of classical density functional theory (DFT) [284,285,286,287,288] are an alternative method to describe transport in nanopores. Classical DFT is based on a variational framework, in which the thermodynamic grand potential, related to the Helmholtz free energy, is minimized under constraints.

4.1.2 Accounting for disorder when modeling diffusion in mesoporous media

There is another substantial difference with the materials considered in Sect. 3, which fundamentally influence the way diffusion in mesoporous materials is modeled: they are inherently disordered at multiple levels. At the atomic level, mesoporous materials only rarely display crystallographic order, although exceptions occur in carbon nanotubes and graphene-based materials, or materials formed by the aggregation of nanocrystals. Almost all mesoporous materials used in industrial heterogeneous catalysis and separation processes—activated carbons, silica, alumina, zirconia or ceria—are amorphous. This literally means “shapeless” in Greek. Nevertheless, at larger scales, there are often various degrees of order that allow us to describe amorphous porous materials—from the structure of the pore walls to the organization of the pores or the particles that constitute the materials [289,290,291,292,293,294,295,296,297,298,299,300,301,302,303,304].

Fully atomistic descriptions are not often used for amorphous mesoporous materials, due to the absence of crystalline periodicity and, thus, the very large number of atoms required to represent the material. Nonetheless, Thyagarajan and Sholl have very recently published a database of over 200 amorphous porous materials that can be used for atomistic simulations. They remark, however, that the generation of new structure remains time-consuming [305]. There are three archetypical categories of models that lend themselves better to studying diffusion in disordered porous media: discrete particle models, pore network models and continuum models. Less used, at least for now, is a fourth archetype that involves the in silico growth of porous materials using computational techniques. Figure 14 illustrates the modeling techniques, which we briefly discuss here.

Fig. 14
figure 14

Adapted from Zalc et al. [266] with permission from Elsevier); b pore network models (Adapted from Ye et al. [322] with permission from Wiley); c continuum models (Adapted from Wang et al. [323] with permission from Elsevier); synthesis mimicking, atomistic models (Adapted with permission from Gelb and Gubbins [324]. Copyright 1999 American Chemical Society). Multiscale representations can combine these archetypes hierarchically at different scales

Archetypical model categories for disordered porous materials: a Discrete particle models (

4.1.2.1 Discrete particle models

Discrete particle models consider coarse-grained particles, which are used as building blocks to construct the material. The mesopores are the negative space, the voids in between these building blocks, which are typically simple shapes, like spheres, platelets or fibers, which do or do not overlap. This representation attempts to conform to the packing, aggregation, agglomeration, fusing or sintering of particles by which the porous material is synthesized experimentally [292,293,294,295,296,297,298,299,300,301, 306,307,308,309,310,311,312,313,314,315,316,317,318,319].

The particles can all be of the same size, or have a particle size distribution, and even include various shapes. The modeled material can be computationally realized by simulating synthesis steps, based on the elementary building blocks. Most commonly, however, the geometrical model results from a random packing that leads to overall properties in agreement with experiments, like porosity and surface area, and, possibly, a measured pore size distribution [320]. X-ray tomography [321] or electron microscopy might also inform such geometrical models. The creation of representations is an important subject all in itself, but which we will not discuss here in detail. Note that an atomistic representation for the particles can be embedded into discrete particle models, which is useful in molecular dynamics simulations. To simulate diffusion of fluids, molecules move through the open, accessible pore space in between the particles.

4.1.2.2 Pore network models

Pore network models [311,312,313, 325,326,327,328,329] are more frequently used to model diffusion, reaction and other phenomena in mesoporous materials. A representation of the pore space via a pore network model is a popular choice, because of the simpler way to directly represent the connected pore space, compared to discrete particle models, where this is indirect. However, the generation of this pore network is not an obvious step, because, in most synthesis methods, the pore network is a side result, formed as the “negative space” around the particle assembly that generates the porous material. Ignoring detailed features of the solid, these models directly represent the pore space as a network of pore bodies (channels), connected through nodes, which may be attributed a non-zero pore volume as well.

Simple pore shapes are generally assumed. By far the most common morphological approximation for the pores is that the pores are smooth cylinders, much longer than they are wide. Sometimes, depending on the material, connected spherical pores, slit-shaped channels or more complex shapes are used instead. Real mesopores in amorphous porous media, however, are never smooth. Their roughness is often deemed too complicated to account for, or a secondary effect. However, Avnir, Pfeifer and Farin [330, 331] demonstrated already in the 1980s that many amorphous porous materials have a self-similar (sometimes self-affine) fractal surface, meaning that similar features occur within a range of length scales, from a few chemical bonds to the pore diameter—a bit like a molecular-scale natural coastline, a classical example of a fractal [330, 332,333,334,335,336,337,338,339]. This implies a certain degree of “order” or symmetry (invariance under magnification), which allows to account for it in models. A fractal dimension, \({D}_{f},\) can be associated to the roughness and experimentally measured. For smooth surfaces, \({D}_{f}=2\), while \({D}_{f}=3\) for surfaces that are so convoluted that they are space filling. For many amorphous mesoporous materials, the fractal dimension of the surface is somewhere in between: \(2<{D}_{f}<3\). This was confirmed by many experimental studies [331, 333, 340], using adsorption measurements and small-angle X-ray scattering [341, 342]. The surface roughness that is inherent to amorphous materials can be accounted for in studies of diffusion and reaction in mesopores [336, 337]. In all these geometrical descriptions for the individual pores, the physicochemical nature of the solid itself is usually ignored, although it could be included through an atomistic description of the surface, especially when transport in individual pores is studied, a subject we will return to in Sect. 4.4.

Pore network models could account not only for the aforementioned morphological properties of the material, through the shape of the pore bodies and their spatial organization, but also for the topology of the pore space, through the connectivity of the pores [315, 343, 344]. The simplest, older representations for porous materials did not account for this, representing a porous material as a collection of individual pore channels, the connectivity of which was ignored. An effective diffusivity would then be obtained by simply averaging over diffusivities in single pores, using the pore size distribution, and then employing a correction factor to indirectly account for the pore network connectivity.

Whatever the nature of the pore network, the experimentally measured porosity, surface area and pore size distribution can be used to inform the parameters of the pore network. However, the connectivity is harder to obtain. Nitrogen adsorption and desorption measurements are the most common, accessible characterization technique to determine the textural properties of mesoporous materials. Because these techniques relate to pore network properties and pore size distributions, it again appears logical that pore network models are more widely applied than discrete particle models, despite their shortcomings. Nevertheless, these shortcomings can be considerable, as the interpretation of adsorption measurements using nitrogen or other probes, as well as (mercury) porosimetry, again requires a model, which involves modeling assumptions as well. The interpretation of porosimetry measurements is a mathematical “inverse problem”, thus inferring a pore size distribution and the pore connectivity is far from easy for a material with limited additional information on the pore shape [320, 345,346,347,348,349,350,351,352,353,354,355]. The danger in this is conflicting or circular argumentations. Ultimately, most studies fall back on the “familiar” long cylindrical pores with a circular cross-section, despite all evidence to the contrary. For example, when a porous material is synthesized by agglomeration of quasi-spherical particles, the pores will be anything but smooth cylinders, yet this is the common assumption. Underlying all this is the assumption that pore shape is of secondary importance—we will see that this is fundamentally incorrect.

Especially for mesoporous materials, the pore space representation is remarkably poor in most studies, while sophisticated modeling approaches are used for the thermodynamics and reaction kinetics. Even when using molecular dynamics simulations [4, 5], most studies try to infer general transport properties from single pores with simplified shapes. We advocate that more attention be paid to this topic, as we illustrate the significant effects of pore network properties on transport.

4.1.2.3 Continuum models

A third category of models for mesoporous media consists of continuum models [356,357,358]. These models remain the most commonly used, and seem the most straightforward, but their apparent simplicity can be misleading. Continuum models use averaging schemes to represent phenomena in porous media via sets of partial differential equations, such as the well-known diffusion–reaction equations, applied over macroscopic domains. The effect of the porous medium on the modelling parameters is accounted for by effective properties, often assumed constant throughout the medium, but not necessarily (as in the example shown in Fig. 14c)—spatial and temporal variations can be included.

To describe diffusion in mesoporous media using continuum models, the simplest assumption is to use effective diffusivities for each species i:

$${D}_{eff,i}=\frac{{\varepsilon }_{s}}{\tau }{D}_{i,}$$
(8)

where \({\varepsilon }_{s}\) is the porosity (accounting for the average reduction in diffusivity, as fluids cannot move through the solid itself), \({D}_{i}\) is the molecular diffusivity of species i, and \(\tau\) is the “tortuosity” [359,360,361,362]. Instead of a single diffusivity, a diffusion tensor is used for anisotropic media. The tortuosity lumps all other factors, apart from the porosity, that lead to a deviation of the diffusivity from that in a straight, smooth pore, whether deviations are due to pore constrictions, unaccounted morphological effects, and effects resulting from a finite pore network connectivity. Because species in the fluid cannot travel through the solid itself, there are fundamental limitations to this approximation, which are especially pronounced for strongly heterogeneous media and low connectivity. The latter could happen, e.g., when there are many dead-end pores or the porous medium approaches the percolation threshold, which could occur due to pore blockage (catalyst deactivation) or for gas/vapor–liquid processes. Clever volume averaging schemes [304] can account for a multitude of phenomena, on the basis of a discrete particle model, at the expense of more complicated equations. A major challenge remains to parameterize these equations, without which the tortuosity is a fitting parameter, and the model loses prediction power. Worse, due to the poor fundamental underpinning, there is no reason why the tortuosity would not depend on operation conditions. This is similar to representing the kinetics in heterogeneous catalysis by power law expressions that may only hold approximately within a range of concentrations.

These comments should not undermine the use of continuum equations. Every model is only a representation of reality. A model can only be as good as the available information, and, as discussed, some of the textural information is not so trivial to obtain from experiments, in particular a parameter like pore network connectivity, which is required to construct a representative pore network model. Solving continuum models is much quicker than solving pore network models. Furthermore, pore size distribution and morphological information at a pore level can be included in the expression for \({D}_{i}\), including nano-confinement effects, such as Knudsen diffusion or effects inferred from molecular dynamics simulations—implying a multiscale simulation approach.

For well-connected pore spaces that are not prone to change during operation, a continuum model, containing information from single-pore models, and with a purely network-related correction via the tortuosity (\(\tau \le 3\)), may be sufficient. Otherwise, one of the other types of models is required.

4.1.2.4 Synthesis mimicking, atomistic or coarse-grained models

Models for porous materials can also be created using computer simulations. A representative part of a porous medium could be grown in silico, using algorithms that emulate synthesis steps, like agglomeration, simulated bond breakage and formation steps in pyrolysis, simulated annealing and quenching [307, 314, 317, 363,364,365,366]. Already proposed by MacElroy [367], Gubbins and co-workers [314] in the 1990s, this method is likely to become more important, thanks to increased access to fast, parallel computing, better understanding of elementary synthesis steps, as well as high-resolution imaging of materials to compare the models with. Reactive dynamics [366] and reverse Monte-Carlo methods [301, 368] are used to construct materials varying from porous carbon to silicon and silica.

4.1.2.5 Multiscale models

Apart from these four broad categories of models, combinations are possible in multiscale representations. For example, a discrete particle model or pore network model can be used within a range of scales beyond the size of an individual pore, and up to the scale of a cell that is representative for the medium as a whole. Properties inferred within the “cell” (mesoscopic scales) can then be used in a macroscopic continuum model, typically by volume averaging, or periodic extension: unlike crystalline materials, the real cells would not be identical in an amorphous material, but they are statistically indistinguishable from each other.

Such techniques allow one to circumvent the intractability of a sample spanning pore network or a fully atomistic simulation at macroscopic scales. For materials very close to a percolation threshold or other materials with a fractal organization of the pores, like highly porous aerogels, discrete models are required, or models that explicitly account for self-similarity in the describing equations (e.g., via fractional diffusion equations) [369,370,371].

We now turn to mesoporous materials from a synthesis perspective, as a basis for modelling studies, discussed in later sections.

4.2 Mesoporous materials and their pore structures

4.2.1 Disordered and ordered mesoporous materials

Mesoporous materials come in a wide range of chemical compositions. They include, e.g., certain metal oxides (from silicon, aluminum, titanium, cerium, other metals and mixtures), sulfides, nitrides and phosphates, as well as carbons and carbon nitrides, non-oxide ceramics, and even polymers and protein crystals. Here, we will not elaborate on all of them, but only discuss a few general aspects and common examples (silica, ceria, titania, carbons and carbon nitrides) as background information, relevant to how to model them (i.e., to represent them properly), in particular for diffusion through their pores.

Until the late 1980s, almost all known mesoporous materials were structurally disordered, and most of them had a broad pore size distribution. Many oxides, in particular, were and are still produced by a sol–gel process [372, 373]. A “sol” is a colloidal solution of small particles, dispersed in a liquid and formed by hydrolysis and condensation reactions of metal alkoxide [Mn+(OR)n] monomers. From this sol, an integrated network or “gel” forms. To remove the fluid, gels are typically thermally treated and undergo various operations to yield the desired porous materials in particular morphologies, such as powders, pellets, fibers, films and monoliths. Most commonly, fluid is removed by drying (possibly after sedimentation and/or centrifugation), but in some cases by solvent exchange and extraction, to form a highly porous aerogel. Other high-temperature processes, like calcination, may follow. Through these forming and shaping steps, the materials, such as silicas and aluminas, employed as catalyst supports and adsorbents often have a bimodal pore structure, whereby the macroporous network (transport pores) ranges over the entire pellet or other structure. The mesopores branch off from the macropores and form mesoporous networks. Each step in the process is crucial in controlling the material’s texture, including porosity and pore size distribution, while maintaining other key properties for practical applications, such as mechanical stability. Despite years of research and the enormous practical relevance, this process is still very empirical, and as much an art as a science; more recently, computational research is starting to be used to understand and ultimately better control the various steps. For more details on sol–gel science, we refer to the outstanding books by Brinker and Scherer [372] and by Linsen [374] with specific examples on materials relevant to catalysis and adsorption.

In the 1990s, mesoporous materials started to attract more attention thanks to new synthesis methods that led to outstanding properties, such as a high surface areas (moving from the typical 100–300 m2/g to over 1000 m2/g), controllable pore sizes with narrow pore size distributions, alternative pore shapes, regular ordering of mesopores, and uniform nanosized frameworks, together with the possibility to functionalize the mesopore walls. Due to these combined properties, the potential applications of mesoporous materials expanded further into energy conversion and storage, heterogeneous, electro- and photocatalysis, advanced separations, drug delivery and biomedical engineering [373, 375,376,377].

These new synthesis strategies of mesoporous materials typically involve templating approaches, such as soft templating, hard templating, multiple templating methods and in situ templating pathways, but also include scaffolding and other template-free approaches [291, 376]. Unlike the structure-directing agents (SDAs) used to synthesize microporous zeolites, templating here involves larger, supramolecular assemblies of organic molecules (e.g., micelles) or solid particles, called, respectively, “soft” and “hard” templates. Around the templates, inorganic material nucleates and grows in a skin-tight manner, so that, upon the removal of the templating structure, its geometric characteristics are replicated by the inorganic materials, forming mesoporous and sometimes (especially for hard templates) macroporous materials. This is not quite the same as an SDA, where the formation of a specific framework is assisted by an organic compound, but where the voids in the inorganic structure, left after removal of the SDA, do not necessarily take the form of the organic molecule. In real templating there is an intimate relationship between the ultimate solid material’s lattice and the removed organic form, similar to a bronze sculpture and its original mold in the lost-wax process. Soft templating can, but need not proceed in a sequential process, e.g., via a preformed liquid crystal of micelles, which then templates an ordered silica; rather, surfactants and inorganic precursor species tend to co-assemble into ordered mesostructured composites.

Some of the first ordered mesoporous silicas and aluminosilicates that were made in this way, and also found industrial interest, were Mobil Corporation’s M41S family [378, 379], the hexagonal mesoporous sieve (HMS) and Michigan State University (MSU) materials [380, 381], and Zhao, Stucky and Chmelka’s Santa Barbara Amorphous no. 15 (SBA-15 [382]). The synthesis of these materials is based on sol–gel chemistry and soft-templating, but the first uses a cationic surfactant (typically, cetyl trimethyl ammonium bromide, CTAB), and the other ones a non-ionic one (for SBA-15, a tri-block copolymer of poly (ethylene oxide) heads and a poly (propylene oxide) core, called Pluronic P123). Different interaction mechanisms between surfactants and silica precursors, mediated by the solvent conditions, including the pH, lead to different ordered composites, and, eventually, mesoporous structures after removal of the micelles by calcination or extraction. The M41S family, illustrated in Fig. 15 consists of three types of mesophases, depending on the ratio of surfactant to silica precursor, namely lamellar MCM-50 (layered structure), cubic MCM-48 and hexagonal MCM-41 phases. MCM-41 possesses a highly regular hexagonal array of uniformly sized channels, whose diameters are in the range of ~2–10 nm, depending on the templates used, the reaction parameters, and the addition of auxiliary organic compounds to expand the pore size. SBA-15 [382] is synthesized using amphiphilic triblock copolymers in strong acidic media; its hexagonally ordered, cylindrical mesopores are sized between ~4.5–30 nm. It should be noted that both M41S and SBA-15 (and other members of the SBA family) are amorphous materials, so the crystalline order and corresponding space group, witnessed by X-ray diffraction or scattering, only relates to the organization of the pores. Therefore, the pore walls are rough, and, in SBA-15, which has thicker pores than MCM-41, there are often even micropores or narrow mesopores connecting the wider, main mesopores. Because of the thicker walls, SBA-15 tends to have higher thermal, mechanical and chemical resistance than MCM-41, which, in conjunction with the high specific surface area (several 100s to ~1000 m2/g) make it an excellent candidate as catalyst support. Further information on these and related mesoporous materials, including mesoporous cellular foams with even wider pores can be found in reviews by, for example, Kresge and Roth [378], Suib [375], Li et al. [376], Yang et al. [383] and Rahmat et al. [384]

Fig. 15
figure 15

(reproduced from Sun and Coppens [385] with permission from The Royal Society of Chemistry)

Ordered, amorphous mesoporous materials: M41S family (on the left, based on Kresge and Roth [378], with permission from The Royal Society of Chemistry) and the cubic MCM-48

Jansen et al. [386] introduced a templating method using small, cheap non-surfactant chemicals; this material, denoted as TUD-1, contains well-defined mesopores with an easily tunable size (2.5–25 nm in diameter), and three-dimensional connectivity in a foam-type structure, as well as surface areas up to about 1000 m2/g, and high thermal and hydrothermal stability. Template-free methods are also in use for the synthesis of mesoporous materials. Here, the mesoporous voids mostly stem from the aggregation of nanoscale building blocks.

Sun et al. [387,388,389] produced well-connected bimodal mesoporous materials using a dual-templating route, allowing them to independently control the pore size at two length scales: the smaller pores are nanopores ~3 nm in diameter, but possibly wider (using CTAB, as in MCM-41), while the wider pores are tunable in the wide mesopore range, from ~16 to >50 nm (using P123, as in SBA-15), so that they could serve to facilitate overall transport. Surface areas may exceed 1000 m2/g and pore volumes are up to 3.5 cm3/g.

Materials with a controlled mesopore and even macropore diameter, shape and pore fraction can also be synthesized by hard-templating routes, known as nanocasting [390]. Here, the mesoporous materials are made from preformed hard templates, such as mesoporous silica and carbon or from aggregates of nanoparticles. This approach avoids the need to control the hydrolysis and condensation of precursors and their co-assembly with surfactants. However, the hard templates are sacrificial, which could make the process expensive and wasteful. Full connectivity between the pores can be hard to achieve as well. On the other hand, the degree of control over, especially, wider pore size is very attractive to create secondary porosity. Indeed, more generally, various templating approaches (soft and hard) can be combined to produce hierarchical porous materials, such as hierarchically structured zeolites [383, 391, 392].

The inner and outer surfaces of mesoporous materials can be modified by introducing organic groups, metals or ions to affect the surface reactivity, make the surface hydrophobic by silylation to preclude water attack, protect the surface from chemical attack, increase the conductivity, etc. Applications of functionalized mesoporous materials in energy conversion and storage in solar cells, supercapacitors, rechargeable batteries and fuel cells are reviewed by Li et al. [376] Some catalytic applications are reviewed by Suib [375]. Size-controlled mesoporous nanoparticles for tunable drug release are described by Bouchoucha et al. [393] Applications of mesoporous silica for the controlled release of biological drugs has been reviewed by Siefker et al. [377] also, Moeller and Bein [394] present mesoporous nanoparticles for controlled intracellular release of bioactive substances. Clearly, efficient transport through the pore network is very important in these applications, and the mesopores form excellent hosts for larger guest molecules. Gunathilake and Jaroniec have developed mesoporous calcium oxide-silica and magnesium oxide-silica composites to adsorb and capture CO2 at ambient temperatures via physisorption [395]. Excellent catalytic performances could be achieved by employing a non-hydrolytic sol–gel (NHSG) process based on the reaction of chloride precursors with ether or alkoxide oxygen donors [396]. A decisive advantage of NHSG over conventional sol–gel and templated syntheses is that mesoporous mixed oxides are easily obtained, avoiding the use of expensive templating agents or high-pressure, supercritical drying. After calcination, the disordered mesoporous mixed oxides prepared by NHSG exhibit specific surface areas and pore volumes that are similar or in some cases much higher than in ordered mesoporous catalysts or aerogels. The lack of order in the porous structure is not a disadvantage for most catalytic applications. Furthermore, the authors state that NHSG mixed oxides are significantly more stable than ordered materials or aerogels [396].

4.2.2 Examples of non-siliceous mesoporous materials

4.2.2.1 Mesoporous titania

Due to its versatility, metals and metal oxides supported on TiO2 have been intensively studied in heterogeneous, thermal and photocatalysis, e.g., in a variety of mild oxidation reactions, including ethane to acetic acid, ethanol to acetaldehyde and oxidative dehydrogenation of propane to propylene. The properties and applications of TiO2 in catalysis are reviewed by Bagheri et al. [397], Niu et al. [398], and Zhang et al. [399], among others. TiO2 exists in three crystalline forms: anatase, rutile, and brookite. Anatase and rutile phases are the more important ones. Rutile is the thermally most stable phase. TiO2-based catalyst supports have outstanding resistance towards corrosion in different electrolytic media. Mesoporous TiO2, in particular anatase, interacts strongly with supported catalytic nanoparticles, which results in increased catalytic activity and stability. TiO2 mesostructures can be synthesized by template-free, soft-template, and hard-template routes. Multiple-template routes are also employed, in particular for hierarchical materials. The soft-template approach is based on the interaction between TiO2 precursors and surfactant molecules. By means of dual templates or post-treatment, such as calcination or ultrasonication, hierarchically porous structures may be created. The soft-template method has useful features, like controllable mesostructures, controllable pore structures, and tunable morphologies. The hard-template method leads to mesoporous materials of a high crystallinity. The template-free methods usually produce materials with randomly distributed mesopores. TiO2 nanofibers and nanosheets possess one-dimensional and two-dimensional porous structures, respectively. Many other structures have been created as well [399]. Various TiO2 supported catalysts, e.g., Au/TiO2, Pd/TiO2, Co/TiO2, and Pt/TiO2, have been developed for industrial applications, including the hydrodesulfurization of hydrocarbon oils, the epoxidation of propene, and the selective catalytic reduction of NO with NH3 [397]. For Fischer–Tropsch synthesis, Co/TiO2 turns out to be an effective catalyst, due to its high selectivity for long-chain linear paraffins, high resistance toward deactivation by water, and low activity for the competitive water gas shift reaction [400]. A further application of TiO2 worth mentioning is the V2O5/TiO2 system. Initially, V2O5-supported TiO2 was applied as a heterogeneous catalyst for the o-xylene oxidation to phthalic anhydride reaction. Then, its capability in oxidative dehydrogenation of propane with only a few byproducts was noted [401]. Furthermore, porous TiO2 materials have been developed for photocatalytic hydrogen generation, owing to their low cost [399]. The photocatalytic CO2 reduction into hydrocarbon fuels is a promising approach for the direct conversion of CO2 to products, like CO, methanol, etc. by sunlight.

4.2.2.2 Mesoporous ceria

Cerium occurs in two oxidation states: Ce (III) and Ce (IV). The extreme ceria compositions are CeO2 and Ce2O3. CeO2 crystallizes in the fluorite structure with a face-centered cubic (fcc) unit cell within the space group \(Fm\bar{3}m.\) Non-stoichiometric CeO2-y can be formed by oxygen release and reduction of Ce (IV) to Ce (III). CeO2 has many applications as a catalyst or as an active support for catalysts. CeO2 may be produced in one-, two- and three-dimensional structures [402, 403]. 1D nanostructured materials, such as nanowires, nanorods, and nanotubes, offer opportunities for fundamental research concerning the influence of size and dimensionality of a material on its physical and chemical properties [404]. 2D nanoplates and nanosheets have an easily controllable morphology by changing reaction parameters, such as precursor ratio, concentration, and reaction time. Compared to 3D CeO2 nanomaterials, the obtained CeO2 nanoplates have increased theoretical surface area to volume ratio and contain desirable (100) surfaces, which leads to a much higher oxygen storage capacity.

More generally, mesoporous ceria has shown great use as a versatile catalyst and catalyst support, owing to its elevated surface area and high dispersion of active secondary components. CeO2 is a very effective supporting material for catalysis by gold nanoparticles, because of its high oxygen storage and release capacity, facile oxygen vacancy formation, and the presence of a narrow Ce f-band. Experiments have shown that the surface of CeO2 can easily be enriched with oxygen vacancies, leading to the aforementioned CeO2-y. Meanwhile, there are some well-established applications of CeO2-based catalysts, the best known of which is as a promoter in the automotive “Three Way Catalysts” (TWCs). Here, doped Al2O3 acts as a thermally stable support, while the active phase consists of noble metals (Pt, Pd, Rh) and a CeO2-based promoter. Nowadays, a solid solution of CeO2-ZrO2 is mostly used, instead of pure CeO2. Ceria is also used in Solid Oxide Fuel Cells (SOFCs). Doped ceria is used as an electrolyte in some designs, or as a barrier layer for cathodes to prevent reaction with the electrolyte. Sometimes, CeO2 is added as a catalyst in both cathodes and anodes. A classic example is the Gd-doped ceria anode (CGO), which increases cell performance as a consequence of the increased oxygen vacancy concentration, caused by the doping. Polymer electrolyte-membrane fuel cells (PEMFCs) have also been explored in conjunction with ceria-based composites. Most studies focused on the half-reaction of the cell, typically implying a reforming process to produce H2, which is consumed in the PEMFC. Other applications of CeO2 are under development for reforming and water–gas shift reactions, oxidation of volatile organic compounds, dehalogenation, partial hydrogenation, photocatalysis, organic reactions, and thermochemical water splitting. Gangopadhyay et al. [405], Montini et al. [402] and Sun et al. [403] have written extensive reviews on ceria and its applications.

4.2.3 Mesoporous carbons and carbon nitrides

Mesoporous carbons [406,407,408,409] and carbon nitrides [410] are applied in adsorption, catalysis, electrochemistry and energy storage. These materials are usually prepared via a template carbonization route, using hard or soft-templating methods.

Pioneered by Ryoo et al., [411, 412] ordered mesoporous carbons (OMCs), such as the CMK (carbon mesostructured by KAIST) materials, are prepared by a hard-templating approach, whereby a mesoporous silica matrix with controlled architecture and a suitable carbon precursor are employed. This precursor is introduced into the mesopores by either wet impregnation or chemical vapor deposition (CVD). The resulting organic–inorganic composite is polymerized and carbonized at high temperature, followed by the removal of the silica mold by etching in HF or alkaline dissolution. The space once occupied by the hard template turns into the mesopores of the carbon material.

Currently, however, soft-templating routes dominate because the synthesis requires fewer steps. Self-assembly of a phenolic resin (monomers: phenol, resorcinol, phloroglucinol, and other ones) and block copolymer surfactants, as SDAs, creates three-dimensionally ordered mesostructures. When the surfactants are removed, mesoporous polymers form, with open pores of dimensions, shape and topology corresponding to the size and structure of the supramolecular aggregates. The last step is the carbonization of the mesoporous polymers to produce mesoporous carbons. Such directly synthesized carbon structures via soft-templating routes are more varied and mechanically stable, due to the continuous framework. Details of the synthesis approaches have been reviewed by, e.g., Xin and Song [406] and Ma et al. [409]. The mesopore size mainly depends on the hydrophobic groups of the surfactants. Pores can be enlarged by increasing the molecular weight of the hydrophobic blocks, rather than that of the copolymers. Adding organic swelling agents is another useful way to expand the pore size, which works by dissolving the hydrophobic organic species, like decane or hexadecane, inside the hydrophobic regions of the surfactant micelles, which causes them to swell; this is quite similar to pore swelling techniques used to expand the pore size of MCM-41.

The morphology of mesoporous carbons is important for industrial applications, such as carbon microwires as lightweight functional filler materials, films in catalysis and separation, monoliths in optics, and uniformly sized spheres in chromatography. There has been great progress in controlling shape and arrangement of mesoporous carbon materials on a macroscopic scale [409]. The synthesis temperature, stirring rate and template influence the crystal size and the shape of the materials considerably. Mesostructured carbon fibers can be obtained by growing them within a confined space, like in an anodic aluminum oxide (AAO) membrane. To obtain fibers with well-aligned mesopores, a shear-aligned block copolymer/polymeric matrix is usually needed. Aerosol-assisted co-assembly has been regarded as an efficient and productive route for simultaneously controlling the morphology and mesostructures of carbonaceous polydisperse spheres with diameters varying from 100 nm to 5 mm [409].

In order to optimize the behavior of the mesoporous carbons towards a given application, nonmetallic species, such as N, F, P, B, –SO3H, –COOH, –OH groups, and metallic species, such as Mg, Ca, Fe, No, Ti, have been successfully introduced into the mesoporous carbon matrix, giving them enhanced functionalities. Mesoporous carbons have shown high adsorption kinetics, and good stability. Their active surfaces provide the possibility of creating specific binding sites for any guest. The loading of magnesium oxides into ordered mesoporous carbons [407, 409] leads to enhanced adsorption capacity for CO2. N-doped OMCs can also be used for CO2 capture, and as an efficient oxygen reduction reaction catalyst, potentially useful in PEMFCs [409]. Mesoporous carbons are exceptionally suited to adsorb organic pollutants [408]. The presence of mesopores allows the adsorption of large adsorbates and often improves the adsorption rate in comparison to disordered microporous activated carbon. Industrial, agricultural and pharmaceutical components have been adsorbed onto mesoporous carbons. In electrode materials, mesoporous carbons work as conductive, stable materials for batteries or sensors.

Mesoporous carbon nitrides (MCNs) [410] with large surface areas and uniform pore diameters are semiconducting materials with highly versatile structural and excellent physicochemical properties. Applications as metal-free catalysts, photocatalytic water splitting, energy storage and convergence, gas adsorption, separation and sensing have been demonstrated. Mesoporous carbon nitrides are synthesized by the aforementioned templating techniques. Their surfaces can be functionalized metals, organic molecules, and inorganic semiconductors. The inherent basic functionality of MCNs is due to the presence of free –NH2 groups on the surface of their walls. Catalysis of Friedel Craft alkylation and acylation, cyclisation of nitriles and alkynes, coupling reactions, Knoevenagel condensation, and esterification have been demonstrated employing MCNs among other reactions. A CO2 adsorption capacity of 1.76 mol/kg at 25 °C and 1 atm has been shown for MCN. A special type of MCN with large pores and cages (MCN-7-T) could even adsorb 13.5 mol/kg at 0 °C and 30 bar.

4.3 Models of mesoporous materials as a basis for diffusion studies

Having surveyed mesoporous materials with some notable examples, focusing on their synthesis for use in catalysis and other diffusion affected processes, we return to modeling. In particular, we elaborate on the pore network models introduced in Sect. 4.1. This choice of presentation is guided by our aim to systematically discuss nano-confinement effects on diffusion in mesoporous materials. For the purpose of studying these nano-confinement effects, it is useful to formally distinguish individual pores, before turning to pore networks. Many of the fundamentals of nano-confinement can be studied in single pores, such as the effect of pore size and geometry (morphology), and even chemical structure. Thus, Sect. 4.3.1 introduces a few pore models, followed by a discussion on fractal models in Sect. 4.3.2. Then, in Sect. 4.3.3, the key characteristics of pore network models are discussed. Such models allow to specifically focus on effects of pore connectivity (topology) and to contrast results with continuum models.

It should be recalled from the Introduction that such division is not essential, but it is convenient for fundamental studies and also to design nano-confinement effects. The latter may be carried through to the particle, pellet and process level, after properly designing those scales as well, not to lose the benefits of nanoscale design. As tomography and parallel computing become cheaper and more accessible, we expect that discrete particle models (Sect. 4.1.2.1) and synthesis-mimicking, atomistic or coarse-grained models (Sect. 4.1.2.4) will become more commonplace, since they allow us to directly include aspects discussed in Sect. 4.2, integrating all effects (morphology, topology and chemical structure), possibly within multiscale representations (Sect. 4.1.2.5), and without artificial sub-divisions in pores and pore networks. For the purpose of design and understanding, however, pore and pore network models are extremely useful.

4.3.1 Single pore or capillary models

In chemical engineering, capillary models are in widespread use. Figure 16 shows a few examples, varying from smooth cylinders (Fig. 16a) to atomically resolved pore structures for carbon nanotubes (Fig. 16d) or silica mesopores, here representing the uniform pores of SBA-15, which are rough on molecular scales (Fig. 16e).

Fig. 16
figure 16

Pore models: a straight cylindrical pores; b tortuous pores [417]; c rough fractal pores [333]; d atomistic model of a carbon nanotube; e realistic, atomistic model of a rough, silica mesopore, e.g., as unit of SBA-15, with a lysozyme molecule adsorbed on its surface

By far the most frequently used model is the straight, smooth cylindrical pore model (Fig. 16a), either on its own, or as the constituting element of a pore network. Slightly more realistic are tortuous capillaries (Fig. 16b). The capillaries may have a pore diameter distribution. When applied to study transport (possibly with adsorption and reaction) in porous media, models differ from each other in how the molar fluxes in capillaries of different sizes are combined. The micro/macropore model (we would now say nanopore/wide pore) by Wakao and Smith [413], and the models by Johnson and Stewart [414], Feng and Stewart [415], and Wang and Smith [416] are well established in the chemical engineering community.

In the models by Feng and Stewart [415] and Wang and Smith [416], the capillary axes are randomly oriented and cross-linked (see Fig. 18a). They are not pore network models in a strict sense, because these models neglect the effect of the local fluctuating concentration field on the macroscopic flux. This approximation is called the “smooth-field approximation” (SFA). The SFA ignores the interaction of fluxes at the intersections of different capillaries; therefore, it may violate the material balance equations there. The SFA holds rigorously only for networks of infinitely long, straight, non-overlapping capillaries.

When applied to represent transport in porous media, these models use the aforementioned tortuosity, τ, as a fitting parameter. The SFA yields a theoretical tortuosity of 3. Aris [272] stated: “When models are made of the configuration of the pore structure τ can be related to some other geometrical parameter but in general it is a fudge factor of greater or less sophistication.”

In his seminal 1951 work on reaction rates and selectivities in catalyst pores, Wheeler [417] was one of the first who explicitly considered pore models, which included their roughness. He proposed that the mean radius \(\bar{\tau}\) and the length \(\bar{L}\) of pores in a pellet are determined in such a way that the sum of the surface areas of all the pores in a pellet is equal to the BET surface area [418] and that the sum of the pore volumes is equal to the experimental pore volume.

The average pore radius and the average pore length are given by:

$$\bar{r}=\frac{{2V}_{g}}{{S}_{g}}\sigma\left(1-{\varepsilon}_{s}\right)$$
(9)

and

$$\bar{L}=\sqrt{2}{V}_{p}/{S}_{x}$$
(10)

where Vg is the specific pore volume, Sg is the specific BET surface area, σ is a pore wall roughness factor, εs is the pellet porosity, Vp is the total volume of the catalyst particle and Sx is the external surface area of the catalyst particle. The factor \(\sqrt{2}\) in Eq. (10) was called “deviousness factor” by Wheeler, and he based this on the assumption that a molecule is, on average, deviating 45° from a straight path; thus, it corresponds to pore tortuosity. The model provides \(\bar{r}\) and \(\bar{L}\) in terms of measurable quantities, and Wheeler would study reaction rates and selectivity in catalyst pores, extrapolating to catalyst pellets via these relationships.

Remarkably, anno 2020, most representations of mesoporous media in chemical engineering do not go much further than the work from these pioneers. Effects beyond the single pore are all hidden in what Aris [272] called a “fudge factor”: a tortuosity that includes everything from actual pore tortuosity to constrictions, reduced connectivity, roughness and other structural effects, wrongly assuming that these effects are constant or insignificant. In the age of rational design and computational modeling, we feel that this position is undefendable and “devious”, indeed.

4.3.2 Fractal models

Of outstanding importance in adsorption, diffusion, and reaction is the structure and accessible area of the internal surfaces lining the pores. As a result of their synthesis, via sol–gel chemistry and other common routes described in Sect. 4.2, catalyst pellets and adsorbents are often made of microscopic, amorphous particles with rough pore walls. This structure influences molecular adsorption, surface diffusion, and Knudsen diffusion. Observed extents of reaction are shifted, owing to the interaction of reactants and products with the pore walls. The yardsticks to measure the accessible surface are the sizes, δ, of the diffusing and adsorbing molecules. As a consequence, small molecules “see” a much larger surface area, S, than larger ones (Fig. 17).

Fig. 17
figure 17

(reproduced from Feder [419] with permission from Springer), with its cascade of fjords (Df ~ 1.52) to b the internal surface of an amorphous porous catalyst (2 < Df < 3) (reproduced from Coppens [333] with permission from Elsevier)

Fractal, self-similar roughness, from a the Norwegian coastline

How S depends on δ is not arbitrary, however. This is because many natural objects with a seemingly disordered structure do, actually, possess scale symmetry: they look similar regardless of magnification. Such structures are called self-similar [332, 419]. An object is self-similar if it consists of N parts that are geometrically similar to the whole, but are (1/r) times smaller. Mandelbrot [332] discovered that self-similarity is a very common symmetry in nature, and part of a wider class of scale-invariance: he called such structures, fractals. Surprisingly, omnipresent natural fractals, like natural coastlines, clouds, lightning and trees, are formally akin to self-similar mathematical objects that were once deemed exotic exceptions, like the Cantor set (of zero length but containing an infinite number of points over parts of it, at all scales), the Koch curve (of infinite length within a finite area, with discontinuous derivative along the curve, at all scales), and the Menger sponge (fitting within a cube, having zero volume and infinite surface area, including self-similar porosity at all scales). A special characteristic of fractals is their fractal dimension, Df. If a fractal object within a 3D space can be completely covered by N(r) non-overlapping spheres of radius r (or disks in 2D, and segments in 1D), the following relation holds:

$$N(r)\sim {r}^{{-D}_{f}}$$
(11)

so that:

$${D}_{f}=\frac{\mathrm{ln}N}{\mathrm{ln}(1/r)}$$
(12)

Euclidean objects, such as straight lines, squares or spheres lead to N(r) ~ r−1, r−2, and r−3, respectively, so that the integer values for Df associated to their expected dimension of 1, 2 and 3, respectively, are recovered. However, the triadic Cantor set, the triadic Koch curve, and the classical Menger Sponge have a fractal dimension of Df = ln2/ln3 ~ 0.63, Df = ln4/ln3 ~ 1.26, and Df = ln20/ln3 ~ 2.73, respectively. These dimensions are non-integers. In general, the values for Df are smaller than the Euclidean dimension of the space in which they are embedded (1, 2 and 3, respectively), but higher than the topological dimension of the object (0 for points, 1 for a line, and 2 for a surface).

The examples of fractal objects we have just considered are referred to as exact or deterministic fractals, because their self-similarity is exact at every stage of their construction. Other examples, like percolation clusters and diffusion-limited aggregates (DLA), are statistically self-similar fractals [419]. These objects are self-similar in an average sense. If we look at one arrangement of a disordered system, it does not look self-similar on different scales, but if many arrangements of the system are superimposed on top of one another, the geometrical properties are self-similar on average.

In order to elucidate the details of how the surface geometry influences specific catalytic processes, the fractal dimension may be inadequate as the sole characteristic parameter [332, 420]. The concepts of lacunarity and multifractality with more parameters can eliminate this shortcoming [332, 335]. For further terms, like self-affinity, multi-affinity and detailed explanations of the mathematics of fractals, more detailed sources on the subject should be consulted, such as the books by Mandelbrot [332], Barnsley [421] and Feder [419].

Mandelbrot’s pioneering work [332] has inspired numerous books and papers on the fractal structure of disordered media. As mentioned earlier, Avnir et al. [330, 331] discovered that the surface of many amorphous porous media, including those relevant to catalysis and adsorption, are fractal and, in a sense, analogous to natural, self-similar coastlines (Fig. 17), but on a molecular scale, implying that the dependence S(δ) is not arbitrary, but a scaling law:

$$S\left(\delta \right)\sim {\delta }^{2-{D}_{f}}.$$
(13)

Investigations over the ensuing years have revealed that many phenomena in heterogeneous catalysis may be described using fractal concepts [311, 312, 330,331,332,333,334,335,336,337, 420, 421], including those related to diffusion, adsorption and reaction.

Fractal geometry has major repercussion on phenomena like steric hindrance and inaccessibility of inner parts of the surface (screening), which could contain active sites as a subset of all surface sites (chemical selectivity). There could also be a change in surface morphology in the course of the reaction (roughening and smoothing), as well as trapping of reactive molecules in cracks and narrow pores. These may be described in terms of scaling laws that depend on the fractal dimension, as discussed in the book edited by Avnir [331], and in the work by Sheintuch [422] and co-workers. For example, for catalytic reactions at a fractal interface, the initial rate of reaction with a surface may be described as follows:

$$rate \sim {{R}_{p}}^{{d}_{{re}^{-3}}}$$
(14)

where Rp is the particle radius and dre is the reaction dimension—a form of structure sensitivity. Details of the influence of the surface geometry on specific catalytic reactions has also been analyzed using multifractal scaling analysis (see, e.g., Havlin and Ben-Avraham [423]).

For the purpose of interpretation and application in design, it is essential to understand the fundamental basis for fractal scaling in amorphous porous materials—and whether it is the porous material that is a fractal object (like the Menger sponge), or only its (internal) surface lining the pores, or even the pores themselves. A typical origin for fractality is aggregation or agglomeration phenomena far from equilibrium, which occurs during common synthesis methods, such as sol–gel synthesis [424]. Some resulting porous materials are fractal themselves and are called mass fractals: their mass scales with the radius following a power law, with exponent Df. An example are highly porous aerogels, in which the tenuous structure is preserved thanks to careful solvent exchange [330, 331, 419, 424]. Normal thermal processing, including drying, does not preserve the fragile fractal architecture, which densifies to a Euclidean (3D) solid, porous object. Often, however, such materials still have a fractal surface. The synthesis process may lead to fractal surface roughness of the final porous materials at scales varying from a few bond lengths (few 0.1 nm) to the local mesopore diameter, in between the constituting particles [340]. Having a fractal surface, within a finite scaling regime (as opposed to infinitely self-similar, mathematical fractals, like a Koch curve) is the most common case of fractal geometry in mesoporous media. A model pore space with a fractal surface is shown in Fig. 17, and Fig. 16c illustrates a model fractal pore, with a molecule diffusing through it. In other material synthesis procedures, leaching or cracks may lead to pores with fractal tortuosity (Fig. 16c, top). Such materials are sometimes called pore fractals [331]. Unlike mathematical fractals, self-similar properties for the surface, mass or pore structure of real porous media in nature and technology only hold within a finite scaling regime [332, 335].

Unfortunately, the literature is rife with examples of “fractal models” that fit power laws of geometric and other properties (e.g., chemical reaction rates) to length scales (e.g., pore or even particle size), which rarely carry physical meaning beyond empirical curve fitting: those approximate power laws, within a narrow range of a characteristic yardstick (δ or Rp), happen to have a non-integer exponent, unrelated to an underpinning fractal geometry—similar to approximating non-integer power-law kinetics without molecular underpinning. For example, care should be taken not to confound experimental observation of a non-integer exponent in an expression like Eq. (14) to a fractal surface, without further evidence of the geometrical structure of the material. Rather than fractal crystallites or fractal distributions of active sites, such a power law could be the result of reactions occurring on sets of crystallite edges, corners and planes in proportions that vary with size, Rp. That said, where fractal properties have been experimentally confirmed, the repercussions are significant. For example, diffusion in disordered, fractal systems does not follow the classical laws describing transport in smooth or periodically ordered media, leading to many anomalous physical properties [331, 420]. In Sect. 4.4, we will return to diffusion in porous media with a fractal surface, given their common occurrence in industrial practice.

4.3.3 Pore network models

As mentioned in Sect. 4.1.2.2, individual pore models do not properly account for the connectivity between the pores (pore network topology), which affects transport and overall performance in, e.g., catalysis and separations. Progress was achieved from the 1980s on, by introducing increasingly sophisticated pore network models, first in two and then in three dimensions [313], as computational capabilities improved.

Some archetypical categories of pore networks are illustrated in Fig. 18; these models extend from purely mesoporous materials to hierarchically structured ones, with several levels of porosity: models as early as Wakao and Smith’s micro-macropore model [413] recognized the practical need to include at least two levels of porosity, corresponding to narrower and broader pores. What is called “micropores” in older articles really means nanopores in modern literature, and thus typically includes mesopores as well, while “macropores” in those models are wider pores to facilitate overall transport.

Fig. 18
figure 18

(adapted from Keil [425] with permission from Elsevier); b Bethe lattice, with pore tree (including connectivity, but no loops) (adapted from Sahimi et al. [313], with permission from Elsevier); c 2D and 3D pore networks, with random or correlated spatial organization of pores, accounting for pore size distribution, connectivity and loops (adapted from Keil [425] with permission from Elsevier); d Voronoi polyhedra [325] and two-level network of nano- and macropores (Creative Commons CC-BY license); e pore network construction on the basis of tomographic data (adapted from Larachi et al. [426] with permission from Elsevier); f 3D pore networks realized from pore network cutting algorithm, within particles or other structures of arbitrary shapes, with arbitrary pore network structural properties [427]. (Creative Commons CC-BY license)

Several generations of increasingly sophisticated and realistic classes of pore network models: a randomly directed straight cylindrical pores

Random three-dimensional networks are advantageous in several respects [312]:

  1. (1)

    Any type of network can be modeled (e.g., regular, irregular, and Voronoi tessellations).

  2. (2)

    The effect of pore connectivity can be taken into account.

  3. (3)

    Any pore size distribution and any pore geometry (e.g., cylindrical, slit-like) can be used. Long cylindrical pores are, by far, the most often applied. The nodes in the network can have a pore volume associated to them as well; most often, however, this volume is assumed to be negligible, which is questionable.

  4. (4)

    Local heterogeneities, for example, spatial variation in mean pore size, can be modeled. This is an important point, as Hollewand and Gladden [428, 429] detected via NMR studies that pore structures in catalysts are spatially quite heterogeneous up to the pellet scale.

  5. (5)

    Any distribution of catalytic active centers may be taken into account. Due to various impregnation and reduction methods, the distribution of active centers will be different. This causes different diffusional fluxes in the pore network.

  6. (6)

    Percolation phenomena [430, 431] can be described. This is particularly useful for modeling deactivation phenomena or gas/liquid multiphase transport, with capillary condensation of vapors or liquid transport.

  7. (7)

    The pore walls can be smooth, irregular, or fractal.

  8. (8)

    Of particular importance is that parameters like tortuosity can be avoided. Tortuosity is a fitting parameter that, in most cases, contains all model deficiencies.

There is a vast literature on pore network models that will not be reviewed here, although, in Sect. 4.5, we will return to specific aspects related to diffusion in these networks. Jerauld et al. [432] and Winterfeld et al. [433] have shown that, as long as the average coordination number of a topologically disordered system is equal to the coordination number of a regular network, the transport properties of the two systems are essentially identical. This is a significant finding, because it means that regular networks can typically be used in computational studies of diffusion and reaction in porous catalysts. However, Hollewand and Gladden demonstrated that this equivalence breaks down at low connectivity, C < 6, as the percolation threshold differs in random and regular networks, and transport properties differ quite significantly near the threshold (Fig. 19) [434].

Fig. 19
figure 19

Adapted from Hollewand and Gladden [434] with permission from Elsevier

Normalized conductivity of a network of conducting bonds (representative of diffusion in a pore network), constructed on random and regular lattices of connectivity C = 3 (hexagonal lattice) or C = 6 (cubic lattice). A randomly chosen fraction p of the bonds is open and conducting, effectively reducing the conductivity when p < 1. Conductivity is normalized to that at p = 1. Between brackets is the bond percolation threshold, pc, for each network; this is the minimal value of p where the overall conductivity (corresponding to the effective diffusivity) drops to zero.

Pore network models have been proven to be an effective research tool to upscale reactive transport processes from the local scale to the continuum scale. Pore network models allow for bridging the gap between laboratory simulations and larger scale transport. They suffer, as expected, also from some limitations. Firstly, one needs an adequate representation of the real porous medium. Pore network models usually simplify the geometry of the pore space, including the geometry of the pores and their spatial organization. Connectivity is an important parameter, but it is not easy to determine experimentally. Seaton [343], and Portsmouth and Gladden [347] have developed techniques to extract this pore connectivity from the hysteresis in nitrogen adsorption and desorption and mercury porosimetry measurements (see also Rigby and Daut [289]).

Furthermore, the exact molecular structures and intermolecular interactions cannot be represented explicitly in purely topological or simple geometrical pore network models. Sahimi and Tsotsis [316] generated pore networks by the Voronoi tessellation [298] of a solid material composed of hundreds of thousands of atoms, and by designating a fraction of the Voronoi polyhedra as the pores. They also used massively parallel molecular dynamics simulations, whereby a crystalline material has been melted at high temperatures, during which atomic bonds have been broken and molecules are removed from the system in order to generate nanopores. Then, the melt has been quenched to room temperature, hence generating a nanoporous material with an amorphous solid matrix. Raoof et al. [326, 328] have presented a method to generate a three-dimensional pore network representing porous media with a coordination number up to 26, with a pre-specified average value for the entire network. The method has been applied to reconstruct real sandstone and granular sand samples.

Several other algorithms for generating pore networks have been developed [316, 325, 327,328,329, 427]. Pore network modeling is conceptually scale indifferent, i.e. the approach can be applied to any length interval, where the structure of the pore space has been experimentally observed. Gostick et al. [327] published an OpenPNM package for pore network simulations. OpenPNM was designed with three overall objectives in mind: accessibility to a wide audience, generalizability to as many applications, and extensibility to simulate any type of physical process. OpenPNM was coded in Python. Great effort was expended creating detailed documentation, and the OpenPNM website (openpnm.org) has a continuously growing list of tutorials and examples. The package can model any cubic, fully random network, and any pore shape.

Ye et al. [427] introduced a pore network cutting technique to create regular or disordered pore networks within particles or other structures of any macroscopic shape, with arbitrary pore size distribution and pore network connectivity (Fig. 18f). Such models are very convenient to investigate effects of various network parameters on, e.g., diffusion and reaction in mesoporous materials, as illustrated for the catalytic hydrogenation of benzene on Pd/γ-alumina; it was shown that the random nature of the pore network as well as the particle shape significantly influence the results. Such insights are useful for catalyst particle design.

A remark about the numerical treatment of pore networks should be made. Instead of solving nth-order partial differential equations in pore networks, the pore space is treated as a network of connected channels. At the nodes, Kirchhoff’s laws [435, 436] should hold. Transport inside the channels of the network is modelled using analytical or numerical methods to calculate one-dimensional solutions of the relevant transport equations [435, 436]. A problem is to find a relation between the pore network model and the real material. Structural properties of the porous material can be obtained based on various experimental methods, such as gas adsorption (with nitrogen as the standard adsorptive) [351,352,353,354], mercury porosimetry [320, 345,346,347,348, 354, 355], small-angle X-ray scattering, electron and X-ray micro- and nanotomography [311,312,313,314,315,316,317,318,319, 321, 426, 437], and nuclear magnetic resonance [438,439,440,441,442]. The next problem is the pore network reconstruction from measurement data. There is a need for upscaling from a typically small imaged volume to larger domains. This leads to the need to construct statistically representative networks, based on the analysis of the image to extract size distributions of pores and throats and their connectivity. However, if the experimental data come from non-imaging techniques such as mercury porosimetry or gas adsorption, where not all pore space characteristics are readily available, regular pore network construction approaches are usually applied with assumed connectivity. Rigby and Chigada [354, 355] have used mean-field DFT [285, 286] to interpret data from integrated gas sorption and mercury porosimetry. The authors demonstrated that the experimental observations can be better understood in the light of mean-field DFT simulations of adsorption in representative pore models. This has led to a better description of the particular physical mechanisms underlying adsorption isotherms in disordered porous solids. In addition, the new method allows to obtain more details on the void space geometry, such as the ratio of pore neck length relative to pore body length. Other pore data evaluation approaches are presented in the review by Xiong et al. [325] More information on the description of disordered nanoporous media, including measurements and modelling of their structure, are presented in several books [307,308,309,310] and articles [311,312,313,314,315,316,317,318,319, 426].

4.4 Modeling diffusion in mesopores: surface diffusion and Knudsen diffusion

4.4.1 Surface diffusion

Diffusion phenomena inside mesoporous materials involve a number of specific features, introduced in Sect. 4.1. These arise in nanoporous media, in particular, as soon as a molecule enters a nanopore from a wider pore or the bulk fluid around the particle. For mesopores much wider than the diffusing molecules, such effects are typically less important.

Molecules in nanopores interact with the walls, so that surface diffusion occurs in conjunction with diffusion in the bulk of the pore (Fig. 12), which, for gases, may also include Knudsen diffusion, discussed in Sect. 4.4.2. Owing to an attractive potential, molecules move close to the pore wall. Here, many influences come into play, as lateral interactions between the molecules may be attractive or repulsive, depending on the adsorbing molecules, the wall material, the temperature and the coverage on the wall. At low temperatures, complex phases on the wall may occur, in particular on metal surfaces [443, 444]. The molecules coming from the gas phase into contact with the pore wall thermalize with the heat bath of the solid. The molecules adsorb on energetically preferred sites. There, they vibrate and hop after some time to a free adjacent site. These adsorption sites are separated by energy barriers, which are significantly lower than the energy barrier for desorption. The minimum energy difference between adjacent sites is the migration energy barrier, an activation energy, EA, of surface diffusion. As the hopping from one site to another requires the surmounting of a barrier, the temperature dependence of this process follows an Arrhenius-type exponential law. More details will be discussed below.

A complete microscopic picture of surface diffusion requires a detailed knowledge of the surface structure of the porous solid, and all interaction potentials involved. Surface diffusion has a many-body nature, i.e., it is a collective phenomenon. In nanopores, molecules move under continuous interaction with the pore wall. In general, the detailed molecular structure of walls of amorphous mesoporous materials can be very complex and is often not known precisely.

When the thermal energies are small (kbT << EA), the adsorbates are confined to the adsorption sites. For temperatures exceeding this value, surface migration is driven by the continuous energy exchange between the adsorbate and the pore wall. The corresponding energy fluctuations result in random jumps from one energy minimum to another, resulting in a stochastic hopping mechanism. Most of the time the adsorbates remain in the adsorption well, where they are vibrating, and only rarely is sufficient energy accumulated to overcome the diffusion barrier. Upon averaging over many events, a hopping rate can be defined.

If kbT ~ EA, the effect of the lateral surface corrugations on the adsorbate motion becomes smaller, until it eventually becomes negligible. Surface migration is thus less restricted and the adatoms transport rather freely, without confinement to specific sites. The hopping rate is now meaningless, and the diffusivity can be described as (non-activated) Brownian motion. At higher temperatures, a desorption from the pore walls occurs and, depending on the Knudsen number (Kn = λ/d, where λ is the mean free path and d is the pore diameter), Knudsen diffusion (Kn >> 1), molecular diffusion (Kn << 1), or both (Kn ~ 1) occur simultaneously (see Fig. 12).

Just like in zeolites, one has to distinguish between transport diffusivities and self-diffusivities. Transport diffusivity describes the transport of mass and the decay of density fluctuations in the system, while self-diffusion describes the diffusive motion of individual particles at equilibrium (see Sect. 2). The transport diffusivities are measured under non-equilibrium conditions in which concentration gradients (gradients in the chemical potentials) exist. Experimentally they can be determined by macroscopic measurements, for example, gravimetric, volumetric, or chromatographic approaches. Self-diffusivities are determined under equilibrium conditions by quasi elastic neutron scattering (QENS) and pulsed field gradient (PFG) NMR.

A problem is how to experimentally differentiate between surface diffusion and bulk (molecular) diffusion. There are many reviews on these phenomena in the literature [3, 241, 247, 250, 252, 265, 271, 279, 308, 311, 444, 445], among others. The macroscopic surface diffusion flux, Js, is equal to the molar amount of particles adsorbed on the surface that passes through a unit area during a unit time interval [241, 252]. It can be described as the product of the surface concentration, cs, and the average velocity of the moving molecules, vs: Js = csvs. According to irreversible thermodynamics [282], the driving force of surface diffusion is the gradient of the chemical potential, µs, of the adsorbed particles. The flux satisfies the following equation:

$${J}_{s}=-{L}_{s}\nabla \left(\frac{{\mu}_{s}}{T}\right),$$
(15)

where Ls is an Onsager coefficient, and T is the temperature. In engineering calculations, one often takes the surface concentration gradient, according to Fick’s first law:

$${J}_{s}=-{D}_{sur}\nabla {c}_{s}$$
(16)

where Dsur is the surface diffusion coefficient (we use this notation to distinguish with the unrelated self-diffusivity, Ds, used earlier). From (15) and (16), we find:

$${D}_{sur}={D}_{sur,0}\frac{1}{RT}{\left(\frac{\partial {\mu}_{s}}{\partial \ln{c}_{s}}\right)}_{T};\;\;\; {D}_{sur,0}=\frac{R{L}_{s}}{{c}_{s}}$$
(17)

Here, \({D}_{sur,0}\) is the surface diffusivity at zero loading. The surface diffusivity generally depends strongly on cs and T, as well as on the interaction between the adsorbate and the solid surface. Additionally, the structure of the solid surface may have a considerable influence on surface diffusion. Finding a functional form, Dsur (cs, T), is not easy. It is experimentally known that Dsur depends on the temperature, and the activated hopping process leads, as expected, to an Arrhenius-type of temperature dependence [446]. The activation energy is related to the adsorption energy, Ead but, as the molecules do not fully desorb in order to diffuse along the surface, the following relation should hold:

$$E_{A} \, = \,\alpha E_{{\alpha d}}$$
(18)

where α depends on the adsorbate-adsorbent system, and 0 < α ≤ 1. An Arrhenius-type relation like:

$${D}_{sur}{(\theta}_{s},T)={D}^{\infty}_{sur}{{(\theta}_{s})\mathrm{ exp}[-E}_{A}/(RT)]$$
(19)

where θs is the fractional surface coverage (a normalized surface concentration, cs), and \({D}^{\infty}_{sur}\) is the surface diffusivity at infinite (very high) temperature, has to be modified. If one plots Eq. (19) as (lnDsur) vs. (1/T), one sometimes finds a slightly convex curve, instead of a straight line. Therefore, Eq. (19) is modified to \({D}^{\infty}_{sur}={D}_{sur}^{*\infty}{T}^{n}\), with n between 0.5 and 1. Contrary to Eq. (18), it has been found experimentally in liquid–solid systems that sometimes EA > Ead, meaning that α > 1. Based on this observation, EA has been split into parts, namely one part that is related to the hopping process, Eh, and the other part, Evac, is related to a hole-making process (creation of a vacancy), such that EA = Eh + Evac. Thus, one obtains:

$${D}_{sur}{(\theta}_{s},T)={D}_{sur}^{*\infty}{T}^{n}{\mathrm{exp}[-(E}_{h}{+E}_{vac})/(RT)]$$
(20)

The hopping frequency can also be calculated by means of transition state theory (TST) [447]. The TST results in an upper bound for the true hopping rate [445]. A further requirement of the calculation of the time evolution of the particle’s movement can be achieved in the framework of Langevin dynamics [448,449,450]. Molecular Dynamics (MD) simulations of hopping adsorbates are widely employed [444, 451, 452] as well to explore potential energy surfaces [452, 453].

The dependency of \({D}_{sur}\) on the adsorbed amount (surface coverage) is based on the Darken equation [241]:

$${D}_{sur}={D}_{sur,0}{\left(\frac{\partial \ln c}{\partial \ln{c}_{s}}\right)}_{T}={D}_{sur,0}{\left(\frac{\partial \ln c}{\partial \ln{\theta}_{s}}\right)}_{T}$$
(21)

where \({D}_{sur,0}\) is the diffusivity at zero loading, and c is the concentration in the fluid phase. Limitations of the Darken equation are discussed by Skoulidas and Sholl [454] in the context of diffusion in zeolites, which bears some similarities to surface diffusion, but includes even stronger confinement effects [455]. Taking Eq. (17) into account, the Darken equation holds when:

$${\mu}_{s}={\mu}_{s0}+RT \ln c$$
(22)

with μs0 a concentration independent constant. If one applies, for example, the Langmuir isotherm to the Darken expression, one obtains:

$${D}_{sur}^{L}\left({\theta}_{s}\right)={D}_{sur,0}\frac{1}{1-{\theta}_{s}}$$
(23)

This relation cannot be correct for all values of the loading, as it would imply that Dsur → ∞ when the maximal loading, θs = 1, is reached. However, it is in remarkably good agreement with experimental data at least up to θs = 0.8 [455]. Kinetic Monte-Carlo simulations with TST provide a rationale for this, independent from a Darken-type theory, as shown by Higashi et al. (HIO model) [456]. Using a variation of the HIO model, Chen and Yang [455] proposed the following, more flexible single-parameter expression that agreed well with both surface diffusion on Vycor glass with 4.6 nm mesopores (using data from Gilliland et al. [446]) and diffusion of various probe molecules in zeolites:

$${D}_{sur}={D}_{sur,0}\frac{1-{\theta}_{s}+(\lambda/2){\theta}_{s}(2-{\theta}_{s})+{[H}_{s}(1-\lambda)](1-\lambda)(\lambda/2){{\theta}_{s}}^{2}}{[1-{\theta}_{s}+(\lambda/2){\theta}_{s}{]}^{2}}$$
(24)

Here, the parameter λ ≥ 0 is a measure of the blockade by other adsorbed particles and Hs is the Heaviside step function; λ = 0 reverts to the HIO model, Eq. (23). This is illustrated in Fig. 20. A list of surface diffusion models is given by Choi et al. [240], who distinguish monolayer and multilayer surface diffusion, and capillary condensate flow. Nevertheless, it remains difficult to describe surface diffusion in general, let alone unify this model with a theory for diffusion in zeolites, given the wide variety of interactions of molecules with heterogeneous surfaces. Both on surfaces and in zeolites, the diffusivity can level off or show a maximum as a function of loading, θs. A case in point was presented by Valiullin et al. [457], where PFG NMR diffusion and NMR adsorption data of acetone in mesoporous silicon with pores of 4 nm showed that the surface diffusivity first increased with loading, but then levels off at higher loadings, contrary to Fig. 20. For a sample with pores of 10 nm, a maximum was found. A diffusion model that included surface heterogeneity and a generalized Freundlich adsorption isotherm showed good agreement with experiments.

Fig. 20
figure 20

(adapted from Chen and Yang [455] with permission from Wiley). The experimental data for zeolites correspond to ethane on 4A at 50 °C (line 2), propane on 5A at 50 °C (line 3), benzene on ZSM-5 at 65 °C (line 4), and triethylamine on 13X at 190 °C (line 5) and 160 °C (line 6)

Unified modeling of surface diffusion on mesoporous materials (here, SO2 on Vycor glass at 15 °C, line 1; λ = 0) and diffusion in zeolites (other lines, λ > 0), using Eq. (24)

The heterogeneity has been modeled by several approaches, such as the parallel-path model [241] and the effective medium approximation (EMA) [252]. The two-dimensional EMA model best fits the experimental data. When the concentration of the adsorbate is near or beyond the region of monolayer coverage and the heat of adsorption is low, surface flow could be treated as a two-dimensional fluid slipping on the surface, and the laws of hydrodynamics can be applied to describe the surface flow [458]. Collective diffusion can be described by means of the hopping model by associating a coverage dependent, effective jump rate with the motion of individual adsorbates. Monte Carlo simulations have been employed to address a variety of questions in collective surface diffusion in an ensemble of adsorbates. Increased or reduced diffusivity was noted in the presence of next-neighbor repulsions or attractions, respectively [445]. Olivares and Reis [245] studied a model of random walks in the interstices of a simple cubic packing of solid spheres, to represent diffusion of a tracer interacting with the internal surface of that medium. A scaling approach showed three different regimes for diffusion in this medium: (1) dominant bulk residence, in which the tracer moves in the bulk most of the time and the diffusion coefficient is of the same order of the coefficient in free solution; (2) dominant surface residence with dominant bulk displacement, in which the tracer is adsorbed on the sphere walls most of the time, but with a very small mobility in that region, so that the average displacement is dominated by hops in the bulk; and (3) dominant surface residence with dominant surface displacement, in which the tracer is adsorbed on the sphere’s walls most of the time and executes most hops along those walls. The average residence times of tracers in the bulk and on the surface lead to a simple adsorption isotherm relating surface and bulk densities with kinetic parameters and the medium’s geometry. An extensive review of adsorption has been presented by Dabrowski [459]. For many examples, the relative contributions of surface diffusivity and bulk diffusivity to the total diffusivity have been discussed by Miyabe and Guiochon [247].

The effect of pore wall flexibility on surface diffusion in carbon nanotubes has been investigated by Jakobtorweihen et al. [81, 82, 460] At low loadings, diffusion is reduced in flexible pores, but, at higher loadings, the fluid–fluid interaction dominates. Complex quantum effects in the context of adsorption and diffusion in mesoporous media have been reviewed by Keil [461].

Experimental methods for measuring surface diffusion have been reviewed in a number of papers [241,242,243,244, 252, 439, 462] and books [3, 418], among other ones. Pivtsov et al. [463] have presented for the first time the precise localization of different guest molecules inside the mesoporous organosilica material UKON2a with a pore size of 6 nm. Applying a complementary set of different pulsed electron paramagnetic resonance (EPR) methods [249, 463], the authors obtained information about the dimensionality of the spatial distribution and local concentration via double electron–electron resonance experiments, orientation of the guest molecules and the pore walls via electron nuclear double resonance spectroscopy. This allowed localizing the guest molecules and shows that their spatial distribution in nanopores strongly depends on their polarity.

Krishna [246] has presented an extension of the Stefan–Maxwell approach to multicomponent surface diffusion of n adsorbed species. The approach treats the vacant sites as the (n + 1)th component in the diffusing mixture, drawing on the analogy with the dusty-gas model for diffusion in porous media, which we return to in Sect. 4.6. This treatment defines two types of diffusivities: the “intrinsic” diffusivities Ðiv, signifying the facility for diffuse exchange between species i and the vacant sites, and the binary diffusivities Ðij, describing the facility for counter-exchange between the adsorbed species i and j. As we have seen previously, the self-exchange Maxwell–Stefan diffusivity Ðii describes correlations among the same species. For negligible interactions between the adsorbed species, the intrinsic diffusivities Ðiv can be identified with the single-sorbate diffusivity. The generalized Maxwell–Stefan (GMS) counter-sorption diffusivity, in turn, is relatable to the coefficients Ðiv and Ðjv. The GMS approach clearly brings out the influence of the fractional surface coverage θi of the adsorbed species on the transfer behavior. This influence is given by the matrix of thermodynamic factors. For one-dimensional transport (coordinate z), the form of the equation describing (now scalar) surface fluxes by GMS is given by [246]

$${J}_{s}=-{n}_{i}^{s}{\left[B\right]}^{-1}\left[\Gamma \right]\frac{d \theta}{dz},$$
(25)

where \({n}_{i}^{s}\) is the total mixture surface molar concentration (in mol/m2), and [B] is the inverted matrix of GMS diffusion coefficients (in m2/s) with the following elements:

$${B}_{ii}=\frac{{\theta}_{i}}{{{{-}\!\!\!D}}_{i}}+\sum_{\small{\begin{array}{c}j=1\\j\ne i\end{array}}}^{n+1}\frac{{\theta}_{j}}{{{{-}\!\!\!D}}_{ij}} \;\;\; \left(i = 1, 2, \dots ., n\right)$$
(26a)
$${B}_{ij}=-{\theta}_{i}\left(\frac{1}{{{{-}\!\!\!D}}_{ij}}-\frac{1}{{{{-}\!\!\!D}}_{iv}}\right) \;\;\; (i, j = 1, 2, \dots, n;\, i \ne j)$$
(26b)

Here, θi is the fractional surface occupancy of species i. In addition, [Γ] is the matrix of thermodynamic correction factors with elements defined by:

$${\Gamma}_{ij}={\theta}_{i}\frac{{\partial \mathrm{ln}(f}_{i})}{{\partial \theta}_{j}} \;\;\; (i, j =\mathrm{ 1,2}, \dots , n)$$
(27)

For the Langmuir isotherm, one obtains:

$${\Gamma}_{ij}^{L}={\delta}_{ij}+\frac{{\theta}_{i}}{1-\sum_{i=1}^{n}{\theta}_{i}}\;\;\; (i,j = 1, 2 , \dots , n)$$
(28)

Furthermore:

$${\theta}_{v}=1-{\theta}_{t} ; \; \; \; {\theta}_{t}=\sum_{i=1}^{n}{\theta}_{i}=1-{\theta}_{v}$$
(29)

where θt is the fractional surface occupancy by the total mixture and θv is the fraction of uncovered sites. This approach can be combined with the Maxwell–Stefan model for Knudsen and bulk fluxes, as well as viscous flow, as discussed later on.

4.4.2 Knudsen diffusion

4.4.2.1 Knudsen diffusion in pore channels

If a molecule has gained enough energy, it desorbs from the surface and undergoes motion through the pore itself. As explained in Sect. 4.1, Knudsen diffusion arises for gases when their mean free path is longer than the (local) pore diameter. In that case, the molecules collide mainly with the walls, instead of with one another.

The value of the molecular diffusivity, Db, can be estimated from the kinetic theory of gases [464, 465]:

$${D}_{b}=\frac{\lambda}{3}\langle v \rangle$$
(30)

where λ is the mean free path, and \(\langle v \rangle\) is the mean molecular velocity:

$$\langle v \rangle =\sqrt{\frac{8RT}{\pi M}}$$
(31)

where \(R\) is the universal gas constant, \(T\) is the absolute temperature, and \(M\) is the molar mass of the gas species. There are formulae to calculate Db from measurable data as well [464, 466]. The mean free path [464] is given by:

$$\lambda =\frac{{k}_{b}T}{\sqrt{2}{\sigma}_{m}P}$$
(32)

where kb is Boltzmann’s constant, σm is the cross-sectional area of a molecule, and P is the pressure.

First formally studied by Knudsen [273] in 1909, the now so-called Knudsen diffusivity, \({D}_{K},\) in a very long (theoretically, infinitely long) cylindrical pore with a circular cross-section of diameter \(d\), is given by:

$${D}_{K}=\frac{1}{3}d\langle v \rangle =\frac{1}{3}d\sqrt{\frac{8RT}{\pi M}}$$
(33)

Thus, formally, the pore diameter replaces the mean free path in the expression for the bulk molecular diffusivity for an ideal gas, derived from the kinetic theory of gases. This expression uses Fick’s law of diffusion to describe the relation between molar flux and concentration gradient:

$${J}_{K}= -{D}_{K}\frac{dc}{dz}$$
(34)

here written for a one-dimensional channel, with its axis pointing in the z-direction.

Knudsen assumed that the collisions with the wall are diffuse, rather than specular; that is, upon colliding with the walls, molecules adsorb for a short time, after which energy is redistributed over their various degrees of freedom, and partially exchanged with the walls, so that, when molecules desorb back into the pore space and continue their trajectory, the angle of reflection is independent of the angle of incidence and follows a cosine distribution—similar to Lambert’s law of diffuse light reflection [3, 256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271,272]. Clausing [467] later showed that the cosine distribution is a direct result of the second Law of Thermodynamics, rather than requiring roughness. Nevertheless, different combinations of specular and diffuse reflection have been included in these molecule–wall interactions as early as Maxwell [468] and Gaede [469]. We will return to this later on, as it is becoming a hot research topic, due to the advent of tailored nanoporous media with atomically controlled pore surfaces, e.g., based on carbon and other nanotubes and graphene (oxide) sheets. Interest in the phenomenon of molecular flow during the early 20th century stemmed often from vacuum technology, when geometric confinement effects and wall interactions would dominate up to macroscopic scales. Clearly, the same information is key to model diffusion of gases in mesoporous materials, used in catalysis (as support or catalyst, or both) and in separations (e.g., in adsorbents, chromatography and membranes).

In 1910, Smoluchowski [470] derived an expression for non-circular cross-sections, showing appreciable effects of cross-sectional shape. In the 1930s, Clausing [471] proposed an elegant framework with expressions for the transmission probability, fT (the probability for a molecule to traverse the pore) as a function of pore channel length. Based on the kinetic theory of gases, the net flux through a channel (and leaving it, in the steady state) is:

$${J}_{K}={f}_{T}\frac{{C}_{0}\langle v \rangle}{4}$$
(35)

where C0 is the concentration at the entrance of the channel; a vacuum was assumed at the lower concentration end, which is possible for derivations, because all molecules move independently in the Knudsen regime. Under those same assumptions, using Eq. (34), the Knudsen flux is:

$${J}_{K}= {D}_{K}\frac{{C}_{0}}{L}$$
(36)

for a pore of length, L. Combining these two equations:

$${D}_{K}={f}_{T}\frac{L\langle v \rangle}{4}=f_T \frac{L}{4}\sqrt{\frac{8RT}{\pi M}}$$
(37)

which relates the Knudsen diffusivity to the transmission probability, for pores of arbitrary length.

For very long, cylindrical channels with circular cross-section, Clausing [471] confirmed Knudsen’s [273] result (with a correction from Smoluchowski [470]):

$${\mathrm{lim}}_{L/d\to \infty } f_{T}=\frac{4d}{3L}$$
(38)

However, he also derived approximate analytical expressions for cylindrical pores for the entire range of L/d, showing a marked influence of pore length on fT, and, thus, on DK. For example, it is trivial that:

$${\mathrm{lim}}_{L/d\to 0} f_{T}=1$$
(39)

and a fairly reasonable approximation for intermediate pore lengths is an interpolation between these results, already proposed by Dushman in the 1920s, so:

$${f}_{T}\sim \frac{4d}{3L+4d}$$
(40)

and:

$${D}_{K}= \frac{1}{1+(4d)/(3L)}\frac{d}{3}\sqrt{\frac{8RT}{\pi M}}=\frac{1}{1+(4d)/(3L)}{D}_{K,\infty }$$
(41)

where \({D}_{K,\infty }\) stands for the prediction of Knudsen’s originally (and most frequently used) formula for infinitely long channels. Note that, for L = d, a quite common situation, \(D_{K} = \left( \frac{3}{7} \right)D_{K,\infty } ,\) so the classical formula overpredicts the “real” Knudsen diffusivity by more than a factor of 2, further calling into question the meaning of “tortuosities”, Eq. (8), which are used to correct diffusivities in porous media! Following Clausing’s pioneering work, pore length effects on Knudsen diffusion have been investigated by a few others as well [256, 471,472,473]. Hesse [473] has demonstrated that pore lengths close to the average pore diameters lead to a higher transport resistance than obtained from the usual models, in agreement with Eq. (41). Despite this, Eq. (34) for the Knudsen diffusivity is almost invariably used, ignoring those factors, and lumping them in expressions for the effective diffusivity at macroscopic scales, where a “tortuosity factor” is used to include not only all these “imperfections”, but also the pore network topology. This is remarkable, and somewhat similar to lumping reaction kinetics by a single expression with an empirical rate constant and an effective reaction order. The effects on predictions of transport rates and overall rates and product distributions for diffusion limited reactions can be expected to be similarly significant.

In addition, there might be non-negligible pore entrance effects for finite-size pores, leading to a concentration jump at the pore entrance [267, 268]. Indeed, fluxes should be continuous, and chemical potentials identical across the pore entrance, yet this leads to an almost always ignored concentration jump, which can be calculated from the kinetic theory of gases, as shown by Coppens and Malek [474]. This jump has an effect on concentration profiles in nanopores for diffusion limited reactions. Such “surface barriers” are much better recognized for configurational diffusion in zeolites, but they can also occur for Knudsen diffusion in mesopores. Argönül and Keil [256] have developed a model that describes surface diffusion and Knudsen diffusion without making the assumption of adsorption–desorption equilibrium, and it incorporates automatically the variation of flow rates with pore length, as well as entrance and exit effects. The model also makes use of the outer surface area of the solid to calculate the boundary conditions for the surface flow. The results for the investigated conditions indicate that, if the surface flow rate is significant, then the independent flow and adsorption equilibrium assumption are not able to describe the behavior of the system adequately. The surface and gas flow rates, the impingement rate and the surface coverage distributions are considerably different from the predictions made based on such assumptions.

Dammers and Coppens [472] have studied the distribution of molecular hits on the wall of a finite cylindrical channel in the Knudsen regime. Particles entered the channel and either returned to the entrance or were transmitted to the opposite channel end. Using a first-passage approach, the authors have derived expressions for the spatial distributions of hitting probabilities. Monte Carlo simulations essentially confirmed the theoretical predictions, but significant boundary effects were observed. The authors have related these to the distribution of chord lengths {r}: g(r), characterizing chords connecting entrance positions and the locations of the first hit, and f(r), representing chords connecting two consecutive collision points. These distributions are very different. In particular, both distributions exhibit asymptotic power law behavior, but with different exponents, that is, g(r) ~ 1/r3 and f(r) ~ 1/r4, respectively. These results are only partially consistent with the results of Albo et al. [262] In particular, entrance effects are lacking in their hitting distributions, which is an artifact of their model. Albo et al. [260] have also performed Knudsen dynamics simulations in cylindrical pores with multiple sections of different diameters. Simulations performed in two-section pores have shown that the angle of transition between the different diameters has little effect on the transmission probability and the distribution of hits on the wall, except when the transition angle is very small, and the transition is extremely smooth.

Pollard and Present [276] have developed formulae that represent a rigorous solution to the long capillary diffusion problem, valid at all pressures and subject only to the limitations of the mean free path type of treatment. First, they found that Bosanquet’s simple formula—total resistance is equal to Knudsen + bulk molecular resistance, see Eq. (7)—is remarkably accurate (Fig. 21) [275]. In addition, they found that, starting with the Knudsen diffusivity at zero pressure, the specific flow in a long tube must initially decrease with pressure, pass through a minimum value, and at higher pressures increase toward the Poiseuille form. In irregular porous media, this minimum is not necessarily observed, masked by pore network and pore size distribution effects. Heydari-Gorji et al. [475] could show for MCM-41 and SBA-15 that adsorption performance is strongly dependent on the pore length. The shortest channels showed the highest capacity and fastest adsorption. These findings were associated with diminished diffusion resistance and enhanced accessibility of the pores inside.

Fig. 21
figure 21

Adapted from Pollard and Present [276] with permission from the American Physical Society

Analytically calculated diffusivity (lower full line) in a cylindrical channel of arbitrary radius, a = d/2, normalized by the mean free path, λ, as a function of a/λ, showing remarkable agreement with the Bosanquet interpolation formula (lower dashed line), Eq. (7). The Knudsen diffusivity (drawn from experiments by Knudsen) and the bulk diffusivity are also shown. Note that all diffusivities are normalized with respect to a < v > .

Krishna and van Baten [264] have carried out MD simulations to investigate the validity of the Bosanquet formula for calculation of the self-diffusivities of a variety of guest molecules within one-dimensional channels in the 2–10 nm range. For weakly adsorbing molecules, like H2, their simulations show that the Bosanquet formula is of reasonable accuracy for a wide range of concentrations inside the pores. Significant deviations from the Bosanquet formula were found with strongly adsorbing molecules at concentrations inside pores, where molecule–wall collisions are dominant. With increasing pore diameter, the bias in the molecular hops, introduced by adsorption has a decreasing influence in their simulations. The results obtained by the authors also cast doubts on the validity of the dusty gas model (Sect. 4.6).

Krishna and van Baten [263] have also investigated the validity of the Knudsen formula for diffusivities in mesoporous BTF-COF, which is a covalent organic framework that has one-dimensional 3.4 nm-sized channels. Their MD simulations showed that, with increasing binding energy of the guest molecules, the zero-loading diffusivity, Ði(0), falls increasingly below the values predicted by Eq. (33). The validity of this equation is restricted to cases where the binding energy of the molecules in negligibly small. For binary mixtures, for which species 2 has a stronger binding energy than species 1, the diffusion selectivity is significantly higher than the Knudsen selectivity, DK,1/DK,2 = (M2/M1)0.5, whereby M1,2 are the molecular weights of the molecular species 1 and 2, respectively. This ratio is called Graham’s law, and is the same for bulk, molecular diffusion. A similar effect has been observed for deviations form Graham’s law in carbon nanotubes [476].

On the other hand, Ruthven et al. [257] have carried out a detailed analysis of their experimental permeance data for several gases (He, Ar, N2, CH4, C3H8) in a mesoporous silica membrane, where they found that the experimental diffusivities are proportional to (T/M)0.5, in conformity with Knudsen’s model (Fig. 22). No obvious difference in behavior was found between the lighter and heavier species (Ar and C3H8). They concluded that, even under conditions of significant adsorption, the simple Knudsen model still provides a good representation of the permeance data. Of course, for porous media (as opposed to single channels), one has to know the tortuosity to obtain quantitative agreement. Petropoulos and Papadokostaki [258] reexamined previous results on Knudsen diffusion and they coined the term “quasi-Knudsen flow” in an attempt to eradicate confusion with genuine Knudsen flow behavior.

Fig. 22
figure 22

Adapted from Ruthven et al. [257] with permission from Elsevier

Experimental measurements of the effective diffusivity in a mesoporous silica membrane, Deff, normalized by the active layer thickness, L, as a function of Knudsen’s theoretical prediction for the dependence on temperature and molecular mass, DK ~ (T/M)0.5. The narrow pore size distribution of the membrane is shown in the inset.

As discussed in Sect. 4.3.2, the pores of amorphous materials are rarely straight and smooth. Pore constrictions understandably affect Knudsen diffusion, and Nakano et al. [261] are amongst those who demonstrated their effect via Monte-Carlo simulations. Zschiegner et al. [269] have studied molecular diffusion in linear nanopores with different types of roughness in the Knudsen regime. A complete compatibility with the laws of normal diffusion has been predicted, and equivalence between transport diffusion and self-diffusion was asserted on the basis of Fick’s and Einstein’s diffusion equations. In single pores, transport diffusion (due to a concentration or chemical potential gradient) and self-diffusion (single molecules, or equilibrium situation) are equal with respect to both their absolute values and their dependence on the surface roughness. A fundamental reason for this equivalence is that the molecules do not interact with one another in the Knudsen diffusion regime; therefore, a molecule does not “feel” potential gradients.

Coppens and co-workers have studied Knudsen diffusion in rough, self-similar fractal pores in detail, deriving an analytical formula for the effect of surface roughness, and validating it using Monte-Carlo simulations on a variety of fractal pore shapes (Fig. 23). The surface roughness acts as traps for diffusing molecules, where a higher fractal dimension increases the frequency of collisions with the walls, and complicates the return of molecules to the “main”, central pore space; this is especially so for smaller molecules, which see a larger surface area and are more easily trapped in the smallest “fjords” along the pores. Self-similarity allowed Coppens to obtain an analytical solution for the Knudsen diffusivity, as a function of molecular diameter, δ, and the fractal characteristics of the pore wall, by using a first-passage time approach, applied to a fractal perturbation on a smooth pore of the same diameter (where the Knudsen diffusivity would be DK0, e.g., given by Eq. (33) or Eq. (41) for cylindrical pores) [267, 268, 290, 477,478,479,480]:

Fig. 23
figure 23

(adapted from Coppens and Dammers [294] with permission from Elsevier); b Generation of a pore with statistically self-similar surface roughness (here, Df = ln13/ln3 ~ 2.33 and DC = ln8/ln3 ~ 1.89) (adapted from Malek and Coppens [290] with permission from Elsevier); c Normalized Knudsen diffusivity in two model pores with fractal walls, comparing first passage time based analytical calculations using Eq. (42) and kinetic Monte Carlo simulations (adapted from Malek and Coppens [268] with permission from the American Institute of Physics); d Generalized model for diffusion in a pore with a heterogeneous wall structure, leading to temporary trapping or adsorption. Adapted from Malek and Coppens [268] with permission from the American Institute of Physics

Effect of surface heterogeneity on Knudsen diffusion: a The movement through a rough pore can be represented as one through a smooth pore with wall perturbations and associated waiting or trapping times

$${D}_{K}=\frac{{D}_{K0}}{1+\alpha (1-{{\delta }^{\prime}}^{\beta })}$$
(42)

Here, the parameters α and β are analytically related to the fractal dimension, Df, and a shape dependent secondary parameter (related to another fractal dimension, DC); δ’ is the molecular diameter, δ, normalized by the outer cutoff of the fractal scaling regime, which could be approximated by the local pore diameter (or obtained from, e.g., SAXS experiments) [333, 335,336,337].

As a first-order approximation [336, 337], this expression can be simplified (mean-field approximation) to:

$${D}_{K}={\delta ^{\prime}}^{2-D_f}{D}_{K0}$$
(43)

This directly shows the inverse relation to the relative change in accessible surface area, see Eq. (13), which is more dramatic for higher fractal dimensions (rougher surface) or smaller molecules. Malek and Coppens [267, 268, 290, 474] carried out three-dimensional dynamic Monte Carlo simulations of Knudsen diffusion in model pores with random fractal surface roughness to probe the influence of surface roughness on transport and self-diffusion. Self-diffusion is strongly influenced by surface roughness, while transport diffusion was found to be roughness independent. Excellent agreement has been observed with analytical results, in particular Eq. (42), as shown in Fig. 23c. The independence on roughness for the Monte Carlo simulations for transport diffusion are derived from a calculation of the transmission probability, fT, through the pores (based on first-passages, similar to Clausing). The time for a molecule to traverse a rough pore is understandably longer (due to the many collisions) than that to cross a smooth pore of the same diameter. However, the probability to exit through the opposite end of the pore is not affected by roughness, because of diffuse reflections on the wall, unaffected by the random scattering through the fractal perturbations. Equation (35) and Eq. (37) suggest that the Knudsen transport diffusivity should then also be roughness independent; this paradox, contradicting the arguments stated earlier in the context of the work by Zschiegner et al. [269] could be explained by Coppens and Dammers [294] as due to a reduction in the projected (average) molecular velocity along the pore axis in a one-dimensional representation for rough pores, and, thus, a unique value for DK, corresponding to Eq. (42) or Eq. (43) is once again found.

Furthermore, Coppens and Dammers [294] have revealed a remarkable similarity between the effects of pore wall heterogeneity on diffusion in very disparate nanoporous media, from zeolites to rough amorphous mesoporous materials, and protein crystals. Their analysis has shown that the various types of diffusion behavior that can be explained via simple models in statistical physics are also observed in a broad variety of complex, heterogeneous nanoporous materials. This hints at a universal framework for transport in nanopores, based on a trapping and hopping model; the traps could be due to chemical or geometrical heterogeneity (Fig. 23d).

Very recently, Besser et al. [481] carried out experiments of Knudsen diffusion of a wide range of gases (similar to those by Ruthven et al. [257], plus CO, CO2, Ne and C2H4) in carefully synthesized membranes with a monomodal pore size distribution (23, 33, and 60 nm pore diameter), grafted with functional groups of broadly varied type and length. Their results are summarized in Fig. 24. Knudsen’s formula, in particular the proportionality to (T/M)0.5 and to the pore diameter, was confirmed. Surface diffusion did not play a significant role. Remarkably, the type of functional groups had no effect on the results, only the length of those groups. Non-binding interactions, similar to steric hindrance, but originating from a “quasi-liquid” layer at the pore walls, resulted in the scattering that is at the basis of classical Knudsen diffusion and the observed decrease in membrane flow. Thus, they found that the origin of the heterogeneity of the pore walls (geometric or chemical) causing delayed molecular retention and diffuse scattering is irrelevant, which they stated to be in agreement with the theoretical findings of Coppens and Dammers [294].

Fig. 24
figure 24

Copyright 2020 American Chemical Society

Experiments of the Knudsen permeance of a wide range of gases at several temperatures through mesoporous inorganic membranes with functional groups display results that are uniquely dependent on the group size, and not its chemical nature. Each result is averaged over multiple gases and temperatures. A reduction of the Knudsen diffusivity as a function of group size is noted, which depends on a dynamic trapping mechanism on the walls. Adapted with permission from Besser et al. [481].

Diffusion in mesopores is clearly a complex phenomenon with confinement effects that are theoretically still not fully understood. Careful experiments on tailored materials are required to investigate and properly account for surface heterogeneity, whatever its origin. This has become possible with the advent of nanostructured mesoporous materials. An important question is the effect of surface scattering, which, according to Clausing, is diffuse (as Knudsen assumed it to be in 1909) for purely thermodynamic reasons, essentially to satisfy microscopic reversibility. Experiments of gaseous diffusion through various inorganic materials tend to agree with this premise. However, as early as 1910, Smoluchowski [470], then Gaede [469], referring to Maxwell [468], included the possibility of partial specular reflection on the channel’s surface. When a fraction f of the collisions is diffuse (due to momentum accommodation) and a fraction 1−f specular, Smoluchowski argued that the formula for DK has to be multiplied by a factor (2/f−1), thus 50% specular reflection would already lead to a 3 times higher diffusivity. Different types of surface deflections and their significant influence on Knudsen diffusion are frequently brought up.

For materials where molecules interact very weakly with the surface, so that they do not temporarily adsorb on it, such deviations could indeed occur. This is unlikely for amorphous materials, such as silica and most materials discussed hitherto, but it could be the case for atomically smooth and ordered surfaces, such as those of carbon nanotubes (CNTs) and slits between two-dimensionally structured materials, such as graphene sheets, with a smooth energy landscape. There is a vast contemporary literature on transport through such materials for nanofluidic applications; we shall not discuss it here, but only mention some results in the context of the fundamental questions pertinent to Knudsen diffusion in mesopores.

Using MD simulations, Skoulidas et al. [482] predicted extraordinarily high transport fluxes for methane and hydrogen through CNTs, associated to similarly high transport diffusivities, almost independent of pressure; self-diffusivity was equally high at low loading, but rapidly dropped with pressure. Their predictions were orders of magnitude higher than those in zeolites with similar pore sizes, due to very high slip with the walls. Nevertheless, Jakobtorweihen et al. [81] showed the importance of including the influence of the flexibility of the walls, especially at low loadings, when this flexibility leads to a marked decrease in the predicted self-diffusivity. At high loadings, the results agree with those of Skoulidas et al., who used a rigid CNT model. Slip has also been invoked to explain high fluxes observed for water flowing through carbon nanotubes, following the experimental observations by Majumder et al. [483] For nitrogen transport through multi-walled CNT membranes, Hinds et al. [484] obtained excellent agreement with regular Knudsen diffusion; functionalization of the tubes further reduces transport. Striolo [485] found through molecular simulations that a few defects in a CNT would suffice to considerably drop the high flux predictions for water flowing through them and lose the effect of a high slip rate.

Worth mentioning for its new fundamental insights, aiming to bridge transport in micro- and mesoporous materials that include only partial momentum accommodation, especially in smooth channels, is the theory of Jepps, Bhatia and Searles [486], which includes an axial momentum gain between diffuse reflections for Lennard–Jones molecules diffusing at low densities through nanopores. Confirming their “Oscillator Model” via MD simulations, they showed that this effect—where molecules oscillate between diffuse wall reflections, as a result of a conservative 1D fluid–solid interaction potential—leads to a marked decrease of the transport diffusivity, below the value predicted by Knudsen’s formula (Fig. 25). Bhatia et al. [487] reviewed molecular transport in nanopores within a historical context, showing how this model can avoid arbitrary separation of transport in a gas-phase and surface diffusion component, as is the case in the dusty gas model and Maxwell–Stefan model (Sect. 4.6). In addition, at high gas–solid interactions in narrow pores, the trajectories are not straight, as is assumed in the classical models. They further reviewed their distributed friction model, which allows to treat multi-component diffusion at any Knudsen number (λ/d).

Fig. 25
figure 25

Adapted from Jepps et al. [486] with permission from the American Physical Society

Theoretical results for the transport diffusion coefficient, here denoted by D0, predicted by the Oscillator Model for low-density transport of methane through smooth, cylindrical silica pores. This model accounts for axial momentum gain between diffuse wall reflections. Results are compared with equilibrium (EMD) and non-equilibrium molecular dynamics (NEMD) simulations, as well as the classical expression for the Knudsen diffusivity and a corrected expression (using an effective diameter, d−2 × 0.92σsf).

Recently, Rhada et al. [488] investigated assemblies of atomically flat graphene sheets, separated by spacers to obtain channels with atomic-scale precision, having a height h = 0.35 N nm, where N varied from 1 to 30. For N = 1, the height of a single graphene layer, no water permeates. For N = 2–5, the flow rate increases vastly, which the authors attributed to large (but height independent) slip lengths and high capillary pressures, leading to a marked enhancement to Poiseuille flow. Helium, a permanent gas, was found to diffuse according to Knudsen’s formula in mesopores, but when N = 2–5 (micropores), an enhancement by up to two orders of magnitude was found, and the flow was even higher than that for N = 10, thus it appears that a breakdown of Knudsen’s formula occurs around the 2 nm transition between micro- and mesopores, both for CNTs and graphene, attributed to atomic smoothness. Water, however, never remained a vapor and was always transported as a liquid at much higher rates (up to 5 orders of magnitude increase) through capillary condensation in all pores.

In summary, well defined, long-range atomically structured materials and strong molecule-wall interactions can induce effects that cause deviation from Knudsen behavior for gases, while capillary flow occurs for liquids, which could also be enhanced due to slip.

4.4.2.2 Knudsen diffusion in mesoporous media

The fundamental dependencies of Knudsen diffusion on surface roughness, temperature, and molecular mass are readily extendable from single pores to porous media. However, most amorphous mesoporous media have a pore size distribution, \({p}_{d}(d)\), and pores are interconnected in complicated ways. This influences the overall, effective Knudsen diffusivity. In first instance, the effective Knudsen diffusivity is often estimated by combining Eqs. (8) and (33):

$${D}_{K,eff}=\frac{2}{3}\frac{{\varepsilon }_{s}}{\tau }\bar{r}\langle v \rangle =\frac{1}{3}\frac{{\varepsilon }_{s}}{\tau }\langle d\rangle \sqrt{\frac{8RT}{\pi M}}$$
(44)

which includes the average pore diameter, \(\langle d\rangle\), while pore network connectivity, constrictions and other geometric effects are contained in the tortuosity factor, τ. As the pore diameters vary inside porous pellets, Knudsen and molecular diffusion mechanisms could occur in series, at least in part of the pores. In this case, the Bosanquet equation [275] can be used, Eq. (7); note that this implies transport resistances in series within each pore, so that Eq. (44) is incorrect. Instead, the combined diffusivity associated to each pore diameter should first be estimated. Using a parallel pore model (Fig. 16a), this leads to the following expression for the effective diffusivity of a single component, i, in a porous medium with pore size distribution \({p}_{d}(d)\):

$${D}_{eff,i}=\frac{{\epsilon }_{s}}{\tau }\int \frac{{D}_{b,i}{D}_{K,i}\left(d\right)}{{D}_{b,i}+{D}_{K,i}\left(d\right)}{p}_{d}\left(d\right)d(d)$$
(45)

For multi-component mixtures, the full Stefan–Maxwell multi-component diffusion equation should be used. To explicitly account for the pore network topology, other methods, like the EMA should be employed (Sect. 4.5).

Whether the tortuosity, τ, is just a geometric factor or a function of the diffusion mechanism, reaction kinetics, and other parameters, has been extensively debated. Zalc et al. [362] have shown by tracer diffusion simulations within random porous structures that tortuosity is independent of the diffusion mechanism for all practical void fractions when an equivalent Knudsen diffusivity is correctly defined. Their model porous structures consisted of random-loose packings of spheres to accurately reflect the void space in practical porous solids. Effective diffusivities were estimated using tracer or flux-based Monte Carlo methods. Tortuosity estimated using the number-averaged distance between collisions, \(\langle {l}_{p}\rangle\), as the characteristic void length scale (instead of the single channel diameter) increased with increasing Knudsen number, Kn. A corrected length scale, first proposed by Derjaguin [489] in 1946, leads to a tortuosity independent of Kn. The expression for the corrected Knudsen diffusivity is as follows [490]:

$${{D}_{K}}^{\prime}=\frac{1}{3}\langle {l}_{p}\rangle \langle v\rangle \left(\frac{\langle {l}_{p}^{2}\rangle }{2{\langle {l}_{p}\rangle }^{2}}-{\beta}_{a}\right)$$
(46)

Derjaguin’s equation aims to describe diffusivities through porous solids with a disordered void structure. The classical equation for the Knudsen diffusivity appears in front of the brackets. It is corrected for non-exponential path distributions by the first term within the brackets and by βa for the nature of the wall reflections. For diffuse reflections (cosine distribution), βa = 4/13. Thus, the predicted (Knudsen) tortuosity is:

$$\tau ={\left(\frac{\langle {l}_{p}^{2}\rangle }{2{\langle {l}_{p}\rangle }^{2}}-\frac{4}{13}\right)}^{-1}$$
(47)

Lu and Torquato [491] showed that point tracers undergoing Knudsen diffusion within voids between impenetrable spheres follow free path distributions equivalent to the chord length distribution of the pore space, g(r). In this case, this distribution is exponential. The first ratio within the brackets of Eq. (47) becomes 1 for exponential distributions, and, thus, the tortuosity τ = 13/9 ≈ 1.44. This tortuosity, τ, is almost independent of the diffusion mechanism for many practical void fractions [362].

This picture changes, however, for slit pores [472] and fractal pore spaces, as investigated by Levitz [492] using a continuous time random walk (CTRW) formalism. In both cases, the second moment, \(\langle {l}_{p}^{2}\rangle\), diverges. For fractal pore spaces, the chord length distribution is a Lévy, power–law distribution, and Knudsen diffusion may become anomalous, time-dependent hyperdiffusion, at least up to the length scales where fractality holds (up to the outer cutoff of the fractal scaling regime, see Sect. 4.3.2).

Another situation is the one where the surface of the porous medium is fractal; the approach described for pore channels is readily generalized, through the fractal perturbation approach on a topologically equivalent medium with a smooth internal surface (Fig. 17b). In Eq. (42) or Eq. (43), DK0 is simply replaced by an expression like Eq. (46), or Eq. (45) is used for parallel pores, estimates for assemblies of discrete particles (e.g., spheres) or, of course, a pore network model is employed, accounting for roughness explicitly at pore level.

Sheintuch [422] and Coppens [333] have extensively reviewed the influence of fractal structures on reaction in porous media (see also the papers [330,331,332,333,334,335,336,337, 420, 421]). Adsorption experiments suggested that the internal surface of many porous media is fractal or (considering an adsorption site distribution) multi-fractal. As discussed in Sect. 4.3.2, the fractal surface morphology is a result of the way in which these catalysts are prepared. Fractal surfaces have an influence on diffusion [267, 268, 290, 477,478,479,480] and reaction [333, 335,336,337]. Using the principles of fractal geometry, the effects of roughness on diffusivities and reaction rates, and therefore also on conversions and product distributions, can be calculated.

Especially when Knudsen diffusion in mesopores plays a controlling role, the effect of the fractal surface roughness proves to be significant and noticeable on macroscopic, industrial reactor scales. A simple change in tortuosity as a fitting parameter is not enough, since it is species dependent. This has been demonstrated for the production of vinyl acetate [337] and catalytic reforming of naphta [333, 335, 479], where not only yields but also product distributions (selectivity) are considerably affected. For catalytic reforming, the mesoporous alumina support, which also acts as an acid catalyst, has a fractal internal surface. An increase in fractal dimension, Df, might lead to a higher active surface area, but its access is not the same for a small hydrogen molecule or a large hydrocarbon. The Knudsen diffusivity of hydrogen, which is more trapped by surface roughness, is disproportionally reduced over that of, say, benzene. The result is that the Knudsen diffusivities are closer than would have been estimated on the basis of the molar mass (where lighter hydrogen moves more rapidly than benzene). Coppens and Froment [333, 335, 479] showed that, when the fractal dimension of alumina was increased from the base case of Df = 2.3 (determined by SAXS/WAXS analysis) for an industrial PtRe/Al2O3 catalyst, the observed reaction rates of key reactions would generally increase, leading to an increase in the octane number of the gasoline produced by the catalytic reforming process. However, a too rough surface (Df > 2.6) would lead to severe Knudsen diffusion limitations and excessive hydrogenolysis and hydrocracking. Thus, there is an optimum fractal dimension to avoid undesired side reactions. This demonstrates the importance of accounting for the surface morphology and pore shape when analyzing or designing mesoporous catalysts.

Tjaden et al. [493, 494] have investigated the calculation of tortuosities in electrochemical devices. They have presented a systematic approach of calculating the tortuosity of different porous samples using image and diffusion cell experimental-based methods. Image-based analyses have included a selection of geometric and flux-based tortuosity calculation algorithms. Using diffusion-cell experiments, no clear dependency between tortuosity and porosity has been obvious over the whole range of analyzed temperatures and gas mixtures. The discrepancy between image and diffusion cell-based values might stem from a small pore diameter of the samples, which seemed to have a higher impact on the measured diffusion flux than the porosity alone, highlighting the importance of Knudsen diffusion effects in the analyzed samples.

As a final note that generalizes the conclusions of Tjaden et al. [493, 494], over the past century, numerous correlations for the tortuosity as a function of the porosity have been proposed, leading to some form of \(\tau =\tau ({\epsilon }_{s})\), based on empirical findings, theoretical calculations and Monte Carlo simulations over packings of spheres and other objects (e.g., to represent fibers), overlapping or not, with mono- or polydisperse distributions of sizes [313]. A popular correlation is \(\tau ={{\epsilon }_{s}}^{-a}\), sometimes called Archie’s Law (in petroleum engineering) or related to the Bruggeman’s [495] effective medium approximation in its simplest form (often used in electrochemical research [494]), where a = 0.5, while the model of Wakao and Smith [413] would predict a = 1. Put together, these lead to a range of tortuosity estimates for the same porosity. This is not surprising, as various pore networks and pore geometries could be conceived that lead to exactly the same porosity, but with a tortuosity that could vary from 1 (straight channels) to infinity (a single, meandering fractal channel, such as the Peano or Hilbert curve). Unless there is a clear indication to assume a particular geometry, based on synthesis (e.g., aggregation of spherical or cylindrical particles), such correlations are best avoided.

4.5 Diffusion in pore networks

Amorphous mesoporous media generally have a highly disordered pore space. The paths that they provide for fluid transport are almost never straight, but tortuous and meandering. To describe tortuous pathways, the terms “tortuosity” and “tortuosity factor” have been introduced [359,360,361,362]. As this is commonly defined as such in the diffusion literature, τ has been used here as a correction factor to relate the bulk diffusivity and the effective diffusivity, possibly correcting first by the porosity, εs, as well, viz. Eq. (8). By extension, the bulk diffusivity needs to be replaced to account for Knudsen diffusion as well. An alternative picture is that of Fig. 16(b) or (c), although this refers only to tortuous pores, and random intersections could be included as in Fig. 18(a). Then, a tortuosity [3, 271, 272] could be introduced to describe the longer connecting paths imposed by the tortuosity of the pores themselves or other obstacles within porous media relative to the straight-line length across the medium:

$${\tau}{^{\prime}}=\frac{\langle {L}_{p}\rangle }{{L}_{s}},$$
(48)

where < Lp > is the average path length through the real medium, and Ls the straight-line length. However, Epstein [360] discusses that, in fact, \(\tau={\tau^{\prime}}^{2}\); commonly, the square is omitted, which his derivations show to be a mistake. Strictly speaking, the term “tortuosity factor” should be used for τ and “tortuosity” just for τ′. This semantic difference is only important when using the parallel pore model and making adjustments to calculate the effective diffusivity.

Pore network models, introduced in Sect. 4.3.3, can properly account for the pore network topology, expanding from single pores to a porous particle, to study diffusion, adsorption and reaction. We are not concerned with surveying them here, as this is a vast subject area that has been reviewed in detail before, as in the extensive reviews by Sahimi, Gavalas and Tsotsis [313], and Keil [311, 312]. Nano-confinement effects can be studied fundamentally on single channels, but, typically, comparison with experiments requires a pore network. Thus, the question that concerns us here is: When are topological, network effects noted? Is there an effect of the presence of reaction on the effective diffusivity or, alternatively, tortuosity?

Sharratt and Mann [497] simulated diffusion and first-order reaction in a two-dimensional pore network model. The authors reported that the effective diffusivity varied with the reaction rate. On the other hand, Zhang and Seaton [498] have found that the effective diffusivity used in the continuum diffusion–reaction model, when defined appropriately, does not depend on the reaction rate. This conclusion applies only to catalysis with monodisperse structures and to microparticles within bidisperse, hierarchical catalysts (Sect. 5). Hollewand and Gladden [434] also warned about the importance of accounting for the hierarchical structure, rather than randomly distributing a bimodal pore size distribution over a pore network of defined connectivity—otherwise, unrealistic, extremely high tortuosities are found (Fig. 26). If a bidisperse catalyst is treated as a continuum, without taking separate account of diffusion in the macropore network and in the microparticles, the effective diffusivity fitted to the experimental data would be expected to vary with the reaction rate. This is because the relative contributions of the diffusion in the macropores and the microparticles in a catalyst pellet depend on the reaction rate. Close to the percolation threshold the continuum diffusion–reaction model breaks down. The results of Zhang and Seaton [498] suggest that the effective diffusivity is independent of the reaction rate and that deviations encountered in the analysis of experimental data may be due to an inadequate mathematical model.

Fig. 26
figure 26

Adapted from Hollewand and Gladden [434, 496]with permissions from Elsevier

Comparison of the tortuosity (factor), τ, predicted using a one-level, uncorrelated, random and b two-level random, bimodal (micro- or nano-/macro-)pore network models, when the connectivity C = 3 or C = 6, as a function of the fraction of the porosity associated to the nanopores (called micro-voidage here). Entirely different trends and values for τ are predicted for the (usually unrealistic) scenario a versus that in b.

Burganos and Sotirchos [499] proposed a general methodology to estimate the effective diffusivity of a pore network with arbitrary pore size distribution using the effective medium theory or approximation (EMA), in combination with the smooth field approximation (SFA). The SFA assumes that the overall, macroscopic and local, pore-based concentration gradients point in the same direction, which is inherent to some earlier pore-based flux models, such as those by Jackson and others [250, 251, 277], Johnson and Stewart [414], or Feng and Stewart [415]. However, this is usually not the case in networks of overlapping or interconnected pores, where the SFA might violate a material balance at the nodes. The SFA would predict a tortuosity τ = 3 for any random 3D network of long channels of the same conductivity (e.g., the same length and diameter). The EMA for networks of resistors, as generalized by Kirkpatrick [500, 501] from Bruggeman [495] and Landauer’s work for binary crystal mixtures, states that:

$$\underset{0}{\overset{\infty }{\int }}\frac{g-{g}_{e}}{g+\left(\frac{C}{2}-1\right){g}_{e}}f\left(g\right)dg=0$$
(49)

where \(f\left(g\right)\) is the conductivity distribution of the resistors in the network (\(f\left(g\right)dg\) is the number of pores per unit volume that have a conductivity between g and g + dg), and C is the connectivity or coordination number, e.g., C = 6 for a cubic network (recall that an average coordination number for disordered networks can usually be used for higher values of C, according to Jerauld et al. [432] and Winterfeld et al. [433], but care has to be taken for lower values of C, according to Hollewand and Gladden [496]). This equation is solved for ge, the EMA effective conductivity, which is the conductivity of a resistor (pore) in a uniform network of the same topology as the original network, but where all conductivities are identical. The EMA states that this uniform network has the same overall resistance to flow or diffusion as the original pore network [500, 501]. Note that g depends on the properties of the resistor, replaced by a pore, thus, at least, its diameter and length; for diffusion in a cylindrical pore, g(d) ~ D(d)d2, where D(d) is the local diffusivity of a component in the pore, which could include Knudsen diffusivity. Interestingly, by comparing the predictions with Monte-Carlo simulations, Burganos and Sotirchos found that the SFA could be used exactly on this uniform network, and that the EMA-SFA combination gave an excellent prediction for the effective diffusivity of the pore network:

$${D}_{eff}^{EMA-SFA}=\frac{{\varepsilon }_{s}}{{\tau }_{SFA}}\frac{{\langle D(d){d}^{2}\rangle }_{EMA}}{\langle {d}^{2}\rangle }$$
(50)

Here, \({\langle D(d){d}^{2}\rangle }_{EMA}\) denotes an average via Eq. (49), and τSFA = 3 for a 3D network or 2 for a 2D one. The average in the denominator is over the pore size distribution.

Sahimi et al. [503] introduced a renormalized effective medium approximation (REMA) to increase the accuracy of transport predictions close to a percolation threshold, where EMA fails. Such a situation might occur due to catalyst deactivation, leading to gradual pore blockage, and decreasing the effective connectivity of the network, as many pores become impenetrable, and the effective diffusivity drops to very low values. At the percolation threshold itself, the porous medium can no longer be traversed, and its network becomes fractal self-similar; thanks to this self-similarity, renormalization group theory can be used to describe transport at or very close to the percolation threshold. Building on the theories of Sahimi et al. [503] and Burganos and Sotirchos [499], Zhang and Seaton [502] found that direct simulations of both molecular and Knudsen diffusion in various pore networks were in excellent agreement with MC-REMA, even down to the percolation threshold; MC indicates the use of Monte Carlo sampling of the pore size distribution. EMA was also quite predictive, except near the threshold, in agreement with Sahimi’s results. They cited only 1 min simulation times on a modern workstation in 1992, which indicates that even such seemingly advanced methods should be easily implementable for more complex networks with more components today, perhaps in combination with the methods discussed in Sect. 4.6.

Dadvar and Sahimi [504] have carried out extensive numerical simulations of diffusion and non-linear reactions in porous materials, using a three-dimensional pore network model of the pore space. A second-order reaction and one according to the Michaelis–Menten kinetics were considered, and the pore size distribution, pore connectivity, and the rate coefficients were varied in order to study their effect on the effective diffusivities. The results of the simulations were also fitted to the classical continuum diffusion–reaction equation, typically used for phenomena at the macroscale. In all considered cases, the authors [504] found that the effective diffusivities under reactive and non-reactive conditions differed, sometimes greatly.

Despite the clear advantages offered by pore network models, they rely on parameters regarding the pore network topology that might not be readily available. Thus, it is comforting to know that a comparison of numerical simulations via continuum and pore network models by Zhang and Seaton [498] showed that the effective diffusivity was generally independent of the reaction conditions (Fig. 27d), and a continuum model could be used, except close to the percolation threshold (e.g., due to deactivation or otherwise severely hindered transport, as shown by, e.g., Ye et al. [322, 505], see Fig. 27b) or if the transport limitations are so significant that concentrations drop appreciably over a few pore lengths from the particle surface (including the surface of the particles within a pellet, for hierarchically structured porous media).

Fig. 27
figure 27

Adapted from Burganos and Sotirchos [499] with permission from Wiley); c MC-REMA effective Knudsen diffusivity predictions (full line) compared to EMA (dashed line) and direct simulations (points), as a function of the fraction of conducting bonds on a cubic lattice, p (Adapted from Zhang and Seaton [502] with permission from Elsevier); d Effective diffusivity predictions for pure diffusion (empty circles) almost coincide with those for diffusion and reaction (filled circles), according to Zhang and Seaton [498]. Reproduced from Zhang and Seaton [498] with permission from Elsevier

Pore network model predictions for the effective diffusivity, Deff: a SFA and EMA-SFA, normalized by exact (MC simulated) results for bulk and Knudsen diffusion, showing poor predictions for SFA, but very good to excellent predictions for EMA-SFA; b Comparison of effective Knudsen diffusivity predictions. All refer to a cubic lattice, with porosities ε1 and ε2 associated to randomly distributed wide and narrow pores with diameter ratio of 10 (

This means that, barring the just cited exceptions or if an explicit pore network topology is known from tomographic imaging and other reliable experimental techniques, a continuum modeling approach suffices for single-phase diffusion and reaction problems. This model should, however, properly include confinement effects induced by the material’s morphology and surface properties, as discussed in Sect. 4.4. For partially condensing vapors and gas–liquid reaction systems, pore network effects might again occur, due to capillary phenomena and partial flooding; these could even result in hysteresis effects, with different rates for increasing or decreasing process variables, like the temperature or pressure, as discussed by Ye et al. [505]

4.6 Fluid mixtures: the Maxwell–Stefan approach and the Dusty Gas Model

For binary mixtures of gases, there are several formulae available to calculate molecular diffusivities [466]. In many chemical engineering applications, one deals with multicomponent gas mixtures in porous media, in which case the Stefan-Maxwell equations [73, 251, 270, 277, 279,280,281, 414, 415, 506, 507] can be used with excellent accuracy, which can be rationalized and expanded using the Onsager formalism of irreversible thermodynamics. Taylor and Krishna [507] have presented the foundations and many applications of the Maxwell–Stefan approach in great detail. Their book includes simultaneous heat and mass transfer, turbulent transfer processes, and film theories developed by means of matrix methods. For multicomponent diffusion in liquids the Taylor diffusion method has proved to be accurate and reproducible.

Generalization of this approach to transport in porous media leads to the Dusty-Gas-Model (DGM) [73, 251, 270, 277, 279,280,281, 414, 415, 506, 507], introduced by Mason, Malinauskas and Evans, which has been the workhorse in multicomponent transport calculations, and can include surface diffusion and viscous flow, as well as thermal gradients. This approximation replaces the solid by an additional gaseous species of giant molecules (the “dust”) with infinite molar mass, so that it becomes effectively immobile. The kinetic theory of gases is then utilized for all species, including the dust, hereby accounting for Knudsen diffusion as well. The derivations of the DGM are complex, but the results are startingly simple, as they formally look like the Stefan-Maxwell equations, with an additional Knudsen diffusion flux, and effective transport parameters. One should keep in mind that the DGM is still an approximation and that there are some errors [283] that, fortunately, cancel each other out to a large extent.

Further simplification to a Fickian description, however, is rarely allowed. Indeed, in 1963 already, Rothfeld [508] showed that even for binary gas mixtures with non-equimolar counter-diffusion, the effective diffusivities in a porous medium depend on the local mole fraction, and a ratio of position dependent molar fluxes (arising from the Stefan–Maxwell or DGM description) must be included in the equations—so that the effective diffusivities are also position dependent. Only then could excellent agreement with experiments be obtained. This greatly affects observed membrane permeabilities and catalytic reaction rates. Even though more than half a century has passed, this is still largely ignored! There is no reason for this, with our current computer capabilities: the full Stefan–Maxwell equations can and should always be used [333, 479], barring situations with single components or equimolar counter-diffusion.

Application of various diffusion models in reactor calculations has been reviewed by Solsvik and Jakobsen [509], and Stegehake et al. [358], among others. Pinto and Graham [510] have extended the Stefan-Maxwell model to electrolyte solutions. Damköhler’s fundamental review article [511] is still of historical and scientific interest.

In the framework of the DGM, one describes the total flux, J, as the sum of diffusive (Knudsen and molecular) flux, JD, viscous flow, JV, and surface flux, JS, as follows:

$$\varvec{J}\, = \,\varvec{J}_{D}+ \varvec{J}_{V} \, + \,\varvec{J}_{S}$$
(51)

assuming that these processes occur in parallel (Fig. 28). The DGM is usually applied as a continuum model at particle or pellet level, assuming the smooth-field approximation (SFA); however, it can also be used to describe transport in single pores, as part of a pore network model. Including chemical reactions, one obtains the following set of component continuity equations for a single, cylindrical pore:

Fig. 28
figure 28

Adapted from Keil [425] with permission from Elsevier

Equivalent network of resistances used to evaluate the total flux in the Dusty Gas Model.

$$\frac{d{\varvec{J}}}{dz}-{\mathbf{\nu}}\frac{4}{d}\mathbf{r}\left(\mathbf{c},T\right)=0$$
(52)

with z the axial pore coordinate, d the pore diameter, c the concentrations of components, r the system of reaction rate expressions (per unit surface area), ν the matrix of stoichiometric coefficients, and T the temperature. The terms are written in matrix notation [425, 435, 436]. Clearly, similar equations can be derived for other pore geometries, e.g., with varying pore cross-section or with a fractal surface morphology. The single pores are then combined into a pore network. This network is solved numerically to obtain the concentration profiles of all components in the network. This is done either directly (Kirchhoff’s laws, see Sect. 4.3.3) or using the EMA or other modeling approaches, described in Sect. 4.5.

4.7 Summary and outlook

Despite more than a century of research on the effects of confinement on molecular transport in mesoporous materials—including some of the most common type of catalyst supports and separation media—fundamental questions remain. Furthermore, there is a deep gap between theoretical research by specialists in the field, and industrial practice, but, worryingly, also the vast majority of research employing those materials in chemical engineering and related fields, which is the research meant to relate to applications of mesoporous materials, either by analysis or materials design.

Analysis of diffusion in practical mesoporous materials is typically rudimentary, resorting to “tortuosity” as a cover-all correction factor for the various irregularities in amorphous materials. This tortuosity factor is assumed constant (popular choices are 2 or 3) or calculated as a function of porosity with a favorite correlation, say a power–law (Archie, Bruggeman, Wakao and Smith) or a linear relation (Maxwell). A positive trend is that multi-component diffusion is increasingly treated with the Stefan–Maxwell equations or the dusty gas model, rather than just Fick’s law; nevertheless, this typically amounts to a curve-fitting exercise with adaptable parameters and empirical correlations, ignoring the assumptions, and hereby risking to use those equations outside their applicable range. Surface diffusion is routinely overlooked. Knudsen diffusion is rarely handled correctly, assuming that Knudsen’s 1909 formula, which was derived for infinitely long, smooth, cylindrical pores with a circular cross-section also holds for rough, disordered mesoporous materials with short, interconnected pores; any deviation from that is covered by “tortuosity”. Fundamental theoretical work, validated by remarkably careful experimental research, some going back more than 50 years, has been reviewed here, next to more contemporary developments that seek to integrate various modes of molecular transport within a single consistent picture (such as the Oscillator Model) and complement molecular simulations. Unfortunately, these theories, as well as efforts from the later 20th Century on pore networks, are rarely applied, even though the methods presented would take little time to run on modern computers. The questions then are: Why is this so? Does it matter? And, where it does, what should be done?

First, the “why”. Amorphous materials are complex. They tend to involve disorder at multiple length scales (from the material itself to the organization of the pore space) and are harder to tackle than periodic, crystalline materials, including zeolites. Thus, theoreticians in the 20th Century had to resort to considerable simplifications. Computer simulations could not be carried out on materials with multi-scale complexity, both morphological and topological, which is the norm for amorphous materials. Only now are we starting to approach this possibility. Nevertheless, important insights were obtained on the dependence of diffusion on parameters at pore-level (pore length, shape, constrictions), fractal organization (pores or surface) or statistical models to account for disorder (related to, e.g., pore network organization, first-passage time calculations or the effect of chord length distributions on Knudsen diffusion). To investigate effects of topological disorder, pore network models can be used to test effects of network connectivity and pore size distribution using percolation theory and (renormalized) effective medium theory, as well as direct pore network simulations (Monte-Carlo based or using matrix solutions). While focusing on topology, most of these models tended to highly simplify pore shape (typically represented by a “resistor” or a long cylinder) and the pores were often spatially randomly distributed, even for bimodal pore size distributions, where the porosity is hierarchical. Thus, despite the vast amount of knowledge and fundamental insights, there has been a gap with real materials, where analysis techniques were rarely able to characterize the pore network architecture and surface properties to sufficient detail to parameterize these models. For tailormade materials, however, this becomes possible. Because of their narrow pore size distribution, Vycor or Controlled Pore Glasses and, then, ordered mesoporous materials, as well as some materials known to be fractal or fractally rough have helped to validate modeling approaches over the past decades.

Second, does it matter? Titze et al. have shown that even for complex nanoporous materials Fick’s first and second laws are applicable for deriving a diffusion coefficient, further validating that even if the precise details of a complex material are not known, there are situations where the overall diffusivity can be calculated using Fick’s laws [512]. Continuum models with a fitted tortuosity could be all one needs, unless there is a risk of very significant diffusion limitations, catalyst deactivation leading to (partial) pore blockage or dynamic multiphase vapor/liquid phenomena. Even so, the tortuosity should at least be reserved to the network effects, while morphology can be separately accounted for, given its different impact on distinct types and sizes of molecules, e.g., when the surface is geometrically or chemically heterogeneous. Furthermore, new classes of materials are emerging, such as those based on graphene or other ordered nanosheets, where reflections with the surface are no longer (purely) diffuse and slip might play a role. Pore entrance and exit effects might also play a significant role that, like for zeolites, is masked by effective diffusivities over samples that would be different if individual particles were observed—thus causing difficulties when extrapolating results for materials design. For hierarchical materials with a bimodal pore size distribution, a single tortuosity is unlikely to represent transport over the range of conditions the material is used for, or even properly represent where the transport limitations lie (something we return to in Sect. 5). Tortuosity may then depend on adsorption or reaction conditions, and change over time: then, “tortuosity” is merely a fudge factor [272], unsuited in many cases. Thus, yes, proper modeling and validation with experiments matters, both for analysis, and, even more, for design of hierarchically structured (meso-)porous materials.

Third, and finally, what should be done? More realistic representations of mesoporous materials can be obtained, thanks to tomographic imaging; however, even the most advanced electron and nano-X-ray tomographic tools are still unable to visualize amorphous materials at a resolution better than a few nanometers, and then only on a very small sample—thus typically ignoring surface roughness, and assuming macroscopic homogeneity and isotropy to extrapolate results. This calls for a merger of imaging with atomistic information and statistical modeling, aided by a combination of experimental techniques (porosimetry, scattering, etc.). Multiscale representations of amorphous materials are becoming possible, suitably integrated with pore network and morphological statistical representations that are parameterized on the basis of advanced materials characterization methods. In addition, synthesis-mimicking simulations can support efforts to accurately represent mesoporous materials. While fully atomistic materials representations are hard (or currently impossible) to achieve, this is not necessary either: one can combine statistical information, consistent with measurable quantities, such as pore volumes, (sub-)particle and pore size distribution, surface area, and even pore shape and connectivity—parameters hard to obtain a decade or more ago. Discrete particle models could be used, and pore networks could be derived from them, with “decorated” walls, informed by imaging, porosimetry and spectroscopy. Then, we advise to revisit those fundamental theoretical insights and approaches, with suitable modifications to include molecular-scale information (nature of the interactions between molecules, and between molecules and the walls). Such an approach would combine molecular dynamics simulations within representative sections of the material (not only single pores, but porous “voxels” or cells) with Monte Carlo simulations, pore network models, statistical volume averaging techniques or simple continuum models to scale up results from the nanoscale to the particle scale. For bimodal pore size distributions, the next level is again treated with a continuum or discrete modelling technique, depending on the degree of heterogeneity. Even macroscopic heterogeneity can be included in such a multiscale modeling approach.

5 Hierarchical porous materials

5.1 Introduction

The size and shape selective properties of zeolites that contribute to their success in industry can become one of their most prominent shortcomings when diffusion in the micropores becomes rate limiting. Size and shape selectivity can be regarded as the most useful properties that zeolites can impart on diffusing molecules by virtue of the diversity in micropore shapes and sizes. The separation of similar molecules by size implies that only the smaller molecule can enter a micropore, but it may also result in slow diffusion, due to the adsorbate kinetic diameter being similar to the pore size.

Strong diffusion limitations in nanoporous materials are the principal reason why many porous materials of practical interest are hierarchically structured (other reasons relate to process requirements, such as a sufficient pellet or macroscopic particle size in fixed bed reactors and separation processes to avoid too high pressure drops, with associated pumping or compressor costs). Particularly slow diffusion through their micropores led to the development of hierarchical zeolites with intercrossing micro- and mesopores to improve their transport properties [513]. Hierarchical MOFs have also been synthesized with a similar goal of improving the overall diffusivity. Bimodal meso/macroporous catalysts (e.g., those using mesoporous alumina as a support) are commonplace in industrial applications.

Among others, Schwieger et al. [391] and Yang et al. [383] have written comprehensive reviews on hierarchical materials. Schwieger et al. [391] define the following two properties as being constitutive for hierarchical materials: first, its structural elements (compartments) have to be characterized by more than one length scale; second, each of these compartments should have a distinct but complementary function. By organizing individual entities or structures in a hierarchical manner, the resulting hierarchical material can outperform the individual entities. Schwieger et al. distinguish three different forms of hierarchy, all related to the design of a desired property or function: structural hierarchy, transport hierarchy, and compositional hierarchy [391].

Structural hierarchy describes a specific, repetitive arrangement of materials with functional subunits, organized over different length scales. Transport hierarchy allows for fast, highly distributed flow. Compositional hierarchy relates to various, possibly self-organized arrangements of atoms, molecules, or larger entities. For porous materials, from the point of view of transport, one has to consider the interplay between the different pore levels, before one can name the overall pore system “hierarchical” [391]. Depending on the interconnectivity between differently sized pore systems (pore levels), they can be classified in two types (Fig. 29). In the first case, larger/wider pores subdivide into several smaller/narrower pores. In the second type, interconnected wide pores intersect the smaller-pore system, i.e., small pores branch off from a continuous large-pore system. Simulations to predict diffusivity in hierarchical materials can provide significant information on the unique diffusion mechanisms involved when different pore sizes intersect.

Fig. 29
figure 29

Reproduced from Schwieger et al. [391], with permission from The Royal Society of Chemistry

Schematic representation of transport hierarchy types, with hierarchical material examples in the bottom row. Distributed transport is required to call a material “hierarchical” from the point of view of transport.

5.2 Hierarchical zeolites

It could be surmised that decreasing the crystal size of zeolites may have the same effect of improving mass transport as in hierarchical zeolites. To some extent, this is true. Both the synthesis of hierarchical zeolites and the control of the synthesis parameters to tune the crystal size distribution and crystal morphology must be considered. USY, which has been employed in hydrocarbon cracking since the sixties, contains micro- and mesopores, with the presence of mesopores aimed to improve mass transport, in what may be the first commercial application of hierarchical zeolites. Mesopore formation affects the intrinsic kinetics, due to the way in which the mesopores are created, namely by dealumination. Desilication has also been employed to create mesopores and control the Si/Al ratio. Demetalation is a method to introduce mesopores or macropores into an already existing zeolite crystal. The principle behind this method is the selective removal or extraction of framework atoms (Al, Si, B or Ti) from crystalline zeolitic materials through acid or alkali treatment, steaming, radiation or a combination of them.

Much recent progress in synthesizing hierarchical zeolites stems from attempts to, first, delaminate zeolites and, second, to try to arrange the resulting layers, leading to hierarchical materials. Delamination is a procedure in which a layered zeolite precursor is treated in alkaline solutions to expand the interlayer distance with surfactants. The expanded layers are then completely separated by sonication to form a collection of poorly ordered zeolite layers, which are disconnected. Calcination of the delaminated materials then removes the surfactants. Materials obtained in this way usually show a high degree of intercrystalline (interlamellar) mesoporosity.

Layered zeolites (PREFER, ITQ-2) were initially synthesized by the groups of Marler [514, 515], Corma [516, 517], and Roth [518, 519], among others. Other prominent examples are ITQ-2, ITQ-6, and ITQ-18. These materials have been delaminated from their precursors, MCM-22, ferrierite, and Nu-6, respectively. The synthesis techniques to produce pillared materials, with polymers or aliphatic chains chemically bonded through silanization or sylilation, linking layers perpendicularly (such as MCM-36 [519]), have allowed to create unique hierarchical structures. Ryoo and coworkers [520] synthesized hierarchical zeolites using a specially designed polyquaternary ammonium surfactant that behaves both as surfactant (to generate mesopores) and as structure directing agent (to generate micropores), such as: C18H37–N+ (CH3)2–C6H12–N+ (CH3)2–C6H12–N+ (CH3)2–C18H37. This is an example of a “soft-templating” approach. Other procedures to generate additional porosity in zeolites include the so called “hard-templating” strategy (also called nanocasting), by using carbons, colloidal silica particles or polymers to imprint meso- or macroporosity during the synthesis.

The beneficial effects of introducing mesopores on zeolite mass transfer have been demonstrated by Christensen et al. [521], using meso/microporous ZSM-5 for the alkylation of benzene with ethene. This is a suitable test reaction, due to its relative simplicity, with a lot of available experimental and computational information, in addition to its large-scale industrial application, where diffusion limitations are present. Hansen et al. [522] used a combination of quantum chemical simulations, molecular simulations (Monte Carlo, Molecular Dynamics) and a continuum approach to simulate the ethylation of benzene over H-ZSM-5 particles. In that study, Maxwell–Stefan equations in combination with ideal adsorbed solution theory (IAST) were employed, whereby as many data as possible have been obtained from quantum chemical (reaction rates) and molecular simulations (adsorption isotherms, diffusivities) to make the model predictive. Integrating the microscopic simulation results from this work in a multiscale simulation of the ethylation of benzene in a catalytic reactor, Rao et al. [121] optimized the structure of macro/meso/microporous pellets, consisting of a composite of H-ZSM-5 and mesoporous silica. They determined the zeolite fraction that would minimize the effects of transport limitations on the pellet and the reactor level. Excellent agreement with reactor-scale experiments was found, but only when surface barriers across the external surface of the H-ZSM-5 crystals were accounted for, showing for the first time the importance of such barriers in practical processes.

Milina et al. [523] demonstrated the impact of different mesopore topologies on the catalytic properties of hierarchical ZSM-5 in the conversion of methanol and of propanal to hydrocarbons. The study was able to identify, through positron annihilation lifetime spectroscopy (PALS), the failure of constricted mesopores to enhance the connectivity of the intrinsic micropores with the crystal surface. A key aspect is that all synthesized samples contained the same micro- and mesoporosity, allowing to ascribe the results to the different mesopore connectivity.

A study by Gheorghiu and Coppens [524] tackled the aspects of different transport pore (broad meso- or macropore) connectivity and pore fraction for a given nanoporous catalytic material (e.g., a certain zeolite) in different models of diffusion and reaction in hierarchical porous catalysts. The study analyzed two kinds of bimodal pore size distributions, fractal-like for the wide pores, and bi-disperse pore structures (constant wide pore size). Molecular diffusion was studied numerically in multi-structured catalyst models in which the pore network geometry was optimized to maximize the yield of an isomerization reaction. The optimal structure that was obtained to maximize the effectiveness factor (per unit nanoporous material) had a fractal-like wide pore hierarchy. However, a uniform bi-disperse pore network maximized the yield per unit total volume. As bi-disperse structures are much easier to synthesize, the design of hierarchical zeolites and other hierarchically structured porous materials with nearly optimum performance is possible.

An atomistic simulation of hierarchical ZSM-5 was presented by Rezlerová et al. [525] who developed a model of ca. 45,000 atoms on ca. 8 × 16 × 8 nm unit cell, where a central cylindrical pore of 4 nm diameter (parallel to the straight channel direction) was produced by removing atoms and generating silanols. GCMC was used to calculate the adsorption isotherms, and the united-atom TraPPE force field was used for the alkane-alkane interactions. For the smallest alkanes tested (methane and ethane) microporous ZSM-5 adsorbed slightly more alkanes (in mmol/g) than the hierarchical (also called dual-porosity) ZSM-5, due to the greater stabilization of these molecules in the micropores. However, for propane and butane, their larger size reduces the effect of the micropores, and more adsorption is observed in the hierarchical than in the microporous ZSM-5. The trend is increasing at pressures close to the vapor pressure, due in part to the observation of a type II isotherm for propane and butane in hierarchical ZSM-5, in stark contrast with all the other cases, in which a type I isotherm was observed for both microporous and hierarchical ZSM-5.

Methods to understand the transport enhancement of hierarchical pore systems include the “Kärger equations” that describe the spatial and temporal changes in concentrations in mesopores (c1) and micropores (c2) [526]:

$$\frac{\partial {c}_{1}}{\partial t}={D}_{1}\frac{{\partial }^{2}{c}_{1}}{\partial {x}^{2}}-\frac{1}{{\tau }_{1}}{c}_{1}+\frac{1}{{\tau }_{2}}{c}_{2}$$
(53)
$$\frac{\partial {c}_{2}}{\partial t}={D}_{2}\frac{{\partial }^{2}{c}_{2}}{\partial {x}^{2}}-\frac{1}{{\tau }_{2}}{c}_{2}+\frac{1}{{\tau }_{2}}{c}_{1}$$
(54)

Here, D1 and D2 correspond to diffusivities in the mesopores and micropores, respectively, and \({\tau }_{1}\) and \({\tau }_{2}\) are the mean lifetimes of a molecule in either meso- or micropores. By non-dimensionalizing the model, a parameter γ is obtained that is the ratio of the mean lifetime in the micropores of a hierarchical material to the time constant of a purely microporous material. This has been called the exchange parameter, which is influenced by the meso/micropore interface. From this analysis, Kärger and coworkers have shown under what circumstances diffusion enhancement in a hierarchical material may be expected. Most interesting is the result that reducing the crystal size in a hierarchical material will only improve transport in the “fast-exchange” conditions where micro- and mesopores operate in parallel [526].

Simulations have also been used to study the diffusion enhancement of hierarchical materials. This has been addressed by Thomas and Subramanian [527] by creating models of different crystals of NaY zeolite with the composition 4 × 4 × 4 Na56Si136Al56O384. Using a force field for n-hexane-zeolite interactions taken from the Fuchs group [528], crystallite models with external surfaces in 1-D, 2-D and 3-D were compared with bulk crystals. The crystal size decreases from the 0-D to the 3-D case. Calculated self-diffusion coefficients included the effect of the intercrystalline diffusion. The molecular dynamics simulations (20 ns) of 128 n-hexane molecules (initially inside the crystal) at 200 and 800 K gave different occupation probabilities for intra- and intercrystalline regions. At 200 K, the stronger adsorption, relative to kbT, implies that less n-hexane molecules are located in the intercrystalline region in which the self-diffusivity is larger. The molecules do not exit the crystal, due to a certain energy barrier. Hence, the smaller crystal size from 0-D to 3-D implies a decreasing self-diffusivity as the external surface area increases. At 800 K, the opposite is observed. The n-hexane molecules now have sufficient kinetic energy to cross the surface barrier, and a considerable number of molecules desorb from the crystal, leading to a large population of the intercrystalline region in which the self-diffusivity is larger (ballistic). This leads to larger self-diffusion coefficients when going form 0-D to 3-D crystals.

A similar approach to simulating diffusion in hierarchical zeolites was taken by Bu et al. [529], who made two models of mesoporous ZSM-5 by considering large squared pores of 20 and 60 Å along the [001] direction and walls with 20 Å thickness. The molecular dynamics simulations extend over 4 ns and start with n-hexane molecules in the mesopore region, outside the crystal. At 300 K the molecules exhibit external surface diffusion, meaning that in 2 ns not so many molecules overcome the surface barrier to enter the micropores. This leads to a relatively small diffusion coefficient, which increases substantially with temperature. At 800 K, diffusion along the external surface becomes more prominent, and there is reduced uptake in the micropores. The Knudsen equation was compared with the hierarchical model, and better agreement (within one order of magnitude) between the two was found for low loading and large mesopore diameters.

Bai et al. [148] also studied n-hexane diffusion in microporous and mesoporous ZSM-5 with a mesoporous volume fraction of 0.55. The molecular dynamics simulations ran for 30 ns and (for the mesoporous case) started from the equilibrium conformations obtained from Monte Carlo simulations. At 363 K the authors find that self-diffusion is faster for the mesoporous material at loadings larger than 0.7 mmol/g (Fig. 30). Since these are overall loadings, a fraction of the loading will be in the mesopores, while for the purely microporous model, all molecules are located in the micropores. The authors also employ the equation Doverall = xmicro·Dmicro + (1 − xmicroDmeso for the diffusion of the mesoporous system, where xmicro is the average fraction of molecules in the microporous region. This allows the authors to compare both microporous and mesoporous contributions to the overall diffusivity. When the temperature is increased from 363 to 546 K, the relative diffusion coefficients (at those two temperatures) for the mesoporous system increase by factors of 12, 46 and 3.3 for loadings of 0.1, 0.39 and 1.55 mmol/g, respectively. Hence, as the loading increases there is a less pronounced effect of the temperature on the diffusion coefficient, due to the decreasing values of xmicro as the temperature increases. This indicates that the temperature has a larger effect on the micropore diffusion than on the mesopore diffusion.

Fig. 30
figure 30

Reproduced with permission from Bai et al. [148]. Copyright 2016 American Chemical Society

Trajectories of selected n-hexane molecules in mesoporous ZSM-5 at 363 K, whose fraction of time in the micropore region is close to the ensemble average, xmicro. Increasing loadings (from left to right) show that molecules tend to locate in the: microporous region (left), channels leading to a pore mouth (center), and mesoporous region (right).

Surface barriers may also form at external surfaces of the crystal and show a combined mechanism of surface diffusion and travelling across the mesopore. Both motions have very different effects on the overall diffusivity, with the latter (contributing to a ballistic regime) depending on the relative amount (size) of mesopores vs. the external surface. Finally, the access from the external surface to the micropores is affected by the presumed “degree of blocking” of the pores related to the fraction of impenetrable pore-mouths, which contribute to the surface barrier [530]. Also, the mesoporous volume fraction affects the overall diffusivity. If all these aspects are studied systematically, an improved understanding of how to tune hierarchy in order to maximize the diffusion of target molecules will emerge. The studies considered above have achieved great progress in this way towards mastering the complexity of hierarchical zeolites.

Zeolites are frequently combined with a mesoporous binder, whereby the mesopores act as transport pores for reactants. Fluid-catalytic cracking (FCC) is an industrially important example. Many other approaches to synthesize mesoporous hierarchical materials are presented by Yang et al. [383] To characterize the catalytic performance of hierarchical zeolite catalysts, Hartmann et al. [392] published a tutorial review of test reactions that are of industrial importance, like methanol-to-olefins, methanol-to-gasoline, alkylation, oxidation and cracking of model substances. The purpose of these test reactions is to compare the hierarchical porous material with standard zeolite catalysts.

A recent review by Coppens et al. [254] discusses the optimal properties of nature-inspired, hierarchically structured zeolites, guided by computer simulations. The authors note a considerable gap between most synthesized structures in academia and industrial design targets, a concern that was also raised by Chal et al. [531] This illustrates the importance of computation in guiding synthesis, to know which textural parameters are most important to control catalyst activity, selectivity and stability. Interconnected hierarchical pore systems, which are accessible from the outer surface of zeolite crystals are preferred due to having more prevalent mesoporous regions. Zeolites with narrow, intracrystalline mesopores may be useful to control the selectivity at microscopic scales, but this approach is not necessarily effective as a mechanism to reduce diffusion limitations. In fact, such nanopores may lead to additional transport limitations that affect the product distribution at macroscopic scales, in addition to barriers at the crystallite surface [254].

5.3 Hierarchical metal–organic frameworks

The plethora of organic linkers and variable connectivity of metal-oxide nodes described above allow some MOF topologies to possess a crystalline hierarchical pore structure. These materials could be excellent model materials for combined experimental and computational studies to further understand diffusion mechanisms involving exchange between mesopores and micropores. NU-1000 is a hierarchical MOF constructed from 8-connected Zr6 nodes linked with tetratopic 1,3,6,8-tetrakis(p-benzoic acid)pyrene and has 34 Å diameter mesopores with 8 Å diameter triangular micropores. The channels are connected by 8 Å diameter windows that allow molecules to exchange between the micro- and mesopores. A schematic of the NU-1000 pore structure that indicates the meso- and microchannels is shown in Fig. 31. Vargas and Snurr investigated linear alkane diffusion in NU-1000 at various adsorbate loadings, where self-diffusivities were shown to increase with loading up to a maximum and then sharply decrease [532]. Adsorbate trajectories were separated into contributions from the meso- and micropores, where micropore diffusivities were roughly an order of magnitude lower than mesopore diffusivities. Adsorbate density maps showed that molecules occupied mesopores for only brief periods of time and resided primarily in the micropores except when the micropores were saturated at high adsorbate loadings. The nonlinear dependence with adsorbate loading was attributed to the heterogeneous pore system, where the occupation of micro- and mesopores strongly influenced diffusion. Chen and Snurr further showed that, at low adsorbate loadings, diffusion in NU-1000 and analogous materials was controlled by the micropore size and not by the mesopores [533]. By varying the linker identity, the relative sizes of the micro- and mesopores could be modulated, and, by changing the MOF topology, different configurations of the micro- and mesopores could be constructed. The tunability of the relative pore sizes in this class of hierarchical MOFs provides an experimentally accessible system to further probe diffusion mechanisms in hierarchical materials.

Fig. 31
figure 31

Reproduced with permission from Vargas and Snurr [532]. Copyright 2015 American Chemical Society

a Structure of NU-1000 with hexagonal mesopore shaded. b Side view of NU-1000 with windows that connect micro/mesopores. c Connectivity of the micro- and mesopores.

5.4 Hierarchical mesoporous materials

In mesoporous catalysts and catalyst supports, the mesopore structure is a decisive property with respect to conversion and selectivity in catalytic reactions. Chemical reactions occur on active centers on the internal surface, or on (e.g., metallic) nanoparticles dispersed inside the supports. The actual composition of reactants and products at the active center is determined by diffusion processes to and from these centers. Reaction rates are, amongst other things, determined by the local composition of the reactants, and indirectly by the products, at the respective active centers. Here, porous structures come into play. If only rather narrow pores were present in the pellets, there would be many active sites because of the large surface area, but diffusion would be slow. If, on the other hand, only larger pores were present inside the pellets, there would be rather fast diffusion of the reactants, but there would be far fewer active sites per unit volume. Between these two extrema, there must be an optimum. Therefore, porous catalyst supports generally have a bimodal pore size distribution, which can be optimized with respect to conversion and selectivity by shifting the respective pore fraction and pore radii of the nanopores and of the wide pores [524, 534, 535]. The connectivity between these pores is important as well. While a smaller pellet diameter could reduce transport limitations, this diameter is often determined by reactor engineering constraints, e.g., the pressure drop due to the particles in a tubular reactor, which must be considered during catalyst design.

Wang and Coppens [536] optimized a washcoat catalyst consisting of a mesoporous catalyst for the selective catalytic reduction of NOx. They showed that the apparent reaction rate can be increased by a factor of 1.8–2.8 by introducing an optimal macroporosity of 0.2–0.4 for a washcoat thickness of 0.5–1.5 mm. In subsequent work, Wang and Coppens [537] also optimized a commercial, mesoporous Ni/Al2O3 catalyst for the autothermal reforming of methane, where they demonstrated that the selectivity toward hydrogen in the produced syngas mixture could be controlled by adjusting the macroporosity and the macropore size of the catalyst. Keil and Rieckmann [534] optimized a hydrodemetallation catalyst by a three-dimensional pore network model. Liu et al. [538] probed the short-term deactivation caused by sulfur condensation in alumina catalysts for the Claus reaction. They found that a Claus catalyst with reasonable hierarchical structure is more robust against this short-term deactivation. Ye et al. [322] have optimized a three-dimensional pore network in which deactivation by coke formation occurred. Propane dehydrogenation in a PtSn/Al2O3 catalyst particle was taken as a model reaction system. Again, a pore network model was necessary, for the reasons discussed in Sect. 4.5: access to the active sites changes due to reduced network connectivity and pore blockage over time—something that a simple “tortuosity” in a continuum model cannot capture. A snapshot of their simulations is shown in Fig. 14b. Catalyst particles with unimodal and bimodal pore size distributions were investigated. The porosity, connectivity, pore size, and their spatial distributions were optimized under two assumptions: constant intrinsic activity per unit catalyst weight and constant intrinsic activity per unit internal surface area. A significant improvement in the time-averaged propene formation rate, compared with a benchmark catalyst, could be achieved. Argönül and Keil [539] could show experimentally for ethene hydrogenation in pellets with different pore structure that the conversion is noticeably influenced by the pore structure.

More examples of catalyst optimization were recently discussed by Coppens et al. [254], who presented a general methodology with rules of thumb to determine the best pore size distribution, (micro-) particle size and porosity for given reaction kinetics on a nanoporous (micro- or mesoporous) catalyst to maximize yields, selectivity and stability against catalyst deactivation—all by controlling transport through the hierarchical pore network. This provides a systematic framework to guide catalyst synthesis, noting which parameters are of primary importance (namely, microparticle or crystal size, and porosity associated to the wide transport pores, as well as their minimal size), and which ones are of secondary importance (namely, the distribution of the wide pores).

5.5 Outlook

The large length scales required to include both nanopores and wide pore channels in many bimodal and other hierarchically structured materials make direct molecular simulations difficult, but modeling becomes possible with carefully chosen multiscale approaches. Over time, computational resources will become much faster, and so it is inevitable that larger and larger models will be attempted to study pore hierarchy. In addition, a better understanding of interconnected pore networks using high-resolution tomographic imaging allows for more realistic models to be constructed. Progress in the synthesis of fine-tuned hierarchically structured zeolites and the discovery of materials like NU-1000, a MOF with crystalline hierarchical pores, allow researchers to create model materials that can be used to probe the fundamental diffusion mechanisms in hierarchical pore networks. Hierarchically structured porous catalysts and separation media can also be designed by using multi-scale simulations that include microscopic details, in combination with a range of methods to simulate disordered meso- and macropore networks, discussed in Sect. 4. In summary, to create meaningful hierarchical models, it is necessary for theoreticians and experimentalists to collaborate, especially for characterization and model generation.

6 Conclusions

In this review we have covered both theoretical fundamentals and simulations at multiple scales to provide insight into the rich diffusion behavior induced by confinement effects in a wide range of nanoporous materials, from single nanopores to hierarchically structured porous particles. We have aimed to relate modeling efforts to experimental measurements, either to support or to promote experimental work. We have considered different pore shapes, sizes and pore network organizations, from crystalline microporous to amorphous mesoporous materials, including zeolites, MOFs, carbons, and oxides. We have shown that molecular simulations are powerful tools that can accurately predict experimental diffusion measurements in a plethora of different nanoporous materials if a suitable model is chosen—and we have discussed how to choose a proper modeling approach, depending on the material properties. In many cases, simulations were central in discovering the mechanism for diffusion by directly probing the free energy landscape in each material. Simulations may also be able to access conditions that are difficult or nearly impossible for experiments to probe, or to envision new opportunities for improved materials design, hereby guiding experimental efforts. To be most effective, close collaboration between theoreticians, simulation and experimental researchers are necessary to develop representative models that can capture the relevant characteristics of a real system and become increasingly predictive. Ultimately, any simulation is a representation of reality that must be experimentally tested and designed to intelligently capture the physicochemical processes of interest.

New developments in the synthesis of nanoporous solids have driven increased attention towards understanding and optimizing their diffusion properties for separations and catalysis. Computer simulation techniques are critical in understanding the diffusion mechanisms of these new materials, as well as predicting diffusion in as-yet non-existent materials. Newly parameterized force fields, integrated multiscale simulation approaches, and more efficient and faster computer resources increase the complexity of materials that can be simulated, while new experimental characterization techniques provide templates for more realistic models. With increasingly powerful characterization and synthesis discoveries for all the classes of nanoporous materials considered here, molecular and multiscale simulations will become even more accurate in predicting new and exciting materials to address challenging engineering problems.