Kinetic Monte Carlo Modelling to Study Diffusion in Zeolite Understanding the Impact of Dual Site Isotherm on the Loading Dependence of n-Hexane and n-Heptane Diffusivities in MFI Zeolite , as Revealed by QENS Experiments

Kinetic Monte Carlo Modelling to Study Diffusion in Zeolite: Understanding the Impact of Dual Site Isotherm on the Loading Dependence of n-Hexane and n-Heptane Diffusivities in MFI Zeolite, as Revealed by QENS Experiments — This study concerns the diffusion of single-component molecules in zeolites, characterised by an isotherm represented by a dual-site Langmuir model with a point of inflection. The systems investigated are n-hexane and n-heptane in MFI zeolite at 300 K. Experiments conducted using the Quasi-Elastic Neutron Scattering (QENS) technique have demonstrated that this inflection has an impact on the loading dependence of the transport Dt and corrected DC diffusion coefficients of these systems. The results of these experiments are described here. A Kinetic Monte Carlo study is then conducted, showing how the energy levels of the molecule adsorption sites in a zeolite affect the loading dependence of the diffusion coefficients of these molecules. Oil & Gas Science and Technology – Rev. IFP, Vol. 64 (2009), No. 6, pp. 773-793 Copyright © 2009, Institut français du pétrole DOI: 10.2516/ogst/2009065 Catalysts and Adsorbents: from Molecular Insight to Industrial Optimization Catalyseurs et adsorbants : de la compréhension moléculaire à l'optimisation industrielle D o s s i e r 10_ogst09016 26/11/09 15:32 Page 773


INTRODUCTION
Zeolites, due to their crystalline structure, possess welldefined pore dimensions, determined by the number of SiO 4or AlO 4 tetrahedra in the ring describing the pore opening.Molecules diffusing within the zeolite pores are always in close contact with the walls, leading to a specific type of diffusion mechanism, the so-called "configurational diffusion" (Weisz and Frilette, 1960), comparable to surface diffusion.In such cases, diffusion is an activated process, and the motion of the molecules is close to a discontinuous site-tosite hopping as described for the first time by Barrer in 1941.A molecule adsorbed on a site corresponding to potential well, vibrates, and its energy fluctuates until it becomes high enough to pass the energy barrier to move to another site.As a result, strongly-binding or tight-fitting guest-zeolite systems are characterised by residence times in adsorption sites that are much longer than travel times between sites.In such cases, Molecular Dynamics (MD) simulation techniques require to integrate the Newton's law of motion on very small time steps, and are not adapted as the molecule cannot explore a representative volume of the pore space due to excessive computational time.
Specific techniques are required to filter vibrations and handle separately fast and slow events by using an intrinsically discontinuous mechanism.These techniques are Dynamical Monte Carlo (DMC) also called Kinetic Monte Carlo (KMC) simulations and Transition State Theory (TST), and are related to the Markovian master equation that describes the evolution over time and space of the system as jumps between successive spatial configurations.The pore space of the zeolite through which the molecule diffuses can be mapped as a network or a lattice of adsorption sites connected by bonds, and diffusion can be viewed as successive random jumps from an adsorption site to another by crossing free energy barriers along bonds connecting neighbouring sites.As reported by Keil et al. (2000), combining Kinetic Monte Carlo (KMC) algorithm and lattice model is the more adapted technique to study diffusion at finite loadings of poorly connected systems.By its flexibility, it allows to consider different types of sites and to account for their local environment.The jump frequencies of molecules are estimated by Transition State Theory (TST) (June et al., 1991;Dubbeldam et al., 2005) or by MD (Smit et al., 1997) and are used as input data in KMC algorithms.KMC simulations has been used to study single component loading dependence of diffusivities of benzene in NaY (Saravanan and Auerbach, 1997), methane and perfluoromethane in MFI (Paschek and Krishna, 2001a) and 2-methylhexane in MFI (Paschek and Krishna, 2000).It has been extended to study binary mixtures diffusion of real system such as CH 4 /CF 4 in MFI (Paschek and Krishna, 2001a), or hypothetical ones by imposing arbitrary transition rates (Maceiras and Sholl, 2002;Paschek and Krishna, 2001b).KMC has also been used to model surface growth (Metiu et al., 1992) and adsorption kinetics on surfaces (Fichthorn and Weinberg, 1992).
This study concerns the diffusion of single-component molecules in zeolites, characterised by an isotherm represented by a dual-site Langmuir model with a point of inflection.Molecular Modelling studies have demonstrated that this inflection has an impact on the loading dependence of the transport D t and corrected D C diffusion coefficients of this type of system (Chong et al., 2005;Krishna et al., 2004Krishna et al., , 2005a, b), b).This effect has been checked experimentally for the first time recently (Jobic et al., 2006), on nC6/MFI and nC7/MFI systems at 300 K, using the Quasi-Elastic Neutron Scattering (QENS) technique, the only experimental method which can be used to measure both self-diffusivity D self and transport diffusivity D t .The corrected diffusivity D C can then be calculated using the thermodynamic correction factor, if the adsorption isotherm is known.
The purpose of this paper is to examine using the KMC technique how the energy levels of the molecule adsorption sites in a zeolite affect the loading dependence of the diffusion coefficients of these molecules.The systems studied are nC6 and nC7 in MFI zeolite, in order to compare the simulation results with the QENS experimental results published by Jobic et al. (2006).In their paper, the authors show that the experimental Self-and Fick diffusivities are in good agreement with those simulated by Molecular Dynamics simulations, even though simulated values are one order of magnitude higher.The general loading dependence on both diffusivities is reproduced as well as the impact of the strong isotherm inflexion of nC7 that involves a sharp maximum of the nC7 Fick diffusivity.

QENS EXPERIMENT: APPARATUS, PRINCIPLE AND PREPARATION OF SAMPLES
The objective of the neutron scattering experiments is to study the influence of the loading θ on the diffusion coefficients D self , D t and D C of hexane and heptane in silicalite-1 at 300 K. Deuterated molecules C 6 D 14 and C 7 D 16 as well as C 6 H 14 were used.In addition, the interchange coefficient, -D 11 (θ), derived from the single-component Maxwell-Stefan formulation, was determined.The silicalite-1 samples studied were activated by heating up to 773 K under oxygen flow.Then, after cooling, the zeolite samples are heated again to 773 K under vacuum to 10 -4 Pa.The alkanes are introduced into the zeolite by pressure difference.All the samples are then left at 400 K for 4 h in order to homogenise the distribution of molecules in the solid.For nhexane, hydrogenated molecules and deuterated molecules were studied, whereas for n-heptane, only deuterated molecules were used.The loadings studied for each molecule are given in Table 1.Since the quantity of product used for each sample is in the region of 4 g to 8 g, the accuracy on θ is estimated to be 10%.A sample containing only the same quantity of dehydrated silicalite-1 was also prepared.
The experiments were conducted at the Institut Laue Langevin, Grenoble, on the IN16 spectrometer, schematised in Figure 1.The resolution of this machine is Δω ≈ 1 μeV, enabling scattering phenomena to be measured at nanosecond scale..9, 1.8, 2.7, 3.6, 4.5, 5.4, 6.3 The continuous beam from the reactor is chopped into neutron "pulses" by the choppers, the monochromator being used to select an incident wavelength of 6.27 Å.The neutrons striking the sample are deflected towards the detectors positioned at predefined angles corresponding to momentum transfer values Q comprised between 0.19 and 1.8 Å -1 .
In a system containing different types of isotope, the intensity measured by neutron scattering can be divided into two separate contributions, one coherent and the other incoherent: (1) S coh (Q,ω) and S inc (Q,ω) are called coherent and incoherent dynamic structure factors.These functions are Fourrier transforms of the correlation functions describing the movements of the atoms in time and space, σ inc and σ coh are the incoherent and coherent cross-sections respectively.When a neutron collides with a nucleus, it is either absorbed or scattered.hωand -hQ, the energy and momentum transfers of the neutrons which have collided with a nucleus (Fig. 2), are written: (2) (3) QENS principle, deflection of the neutron on colliding with the nucleus.Q → is the wave-vector transfer, k 0 → the wave vector of the incident neutron, k 1 → that of the deflected neutron, m the mass of the neutron andhω the energy transfer.
S inc (Q,ω) is the incoherent part of the signal describing the individual movements of the atom in time and space, to obtain the self-diffusion coefficient D self .The molecular movements can be separated into translational, rotational and vibrational movements and are assumed to be dynamically uncoupled.The vibrational movements, occurring at much smaller time scales, have a limited contribution in the region of the quasi-elastic peak.
S coh (Q,ω) is the coherent contribution of the signal.It describes how the collective movements of the nuclei in the structure are correlated, in order to access the global properties such as D t , by studying the changes in local concentration gradients.One characteristic of the coherent signal is that it is highly dependent on Q: (4) When the molecules studied contain hydrogen atoms, the scattering is mainly incoherent since the incoherent crosssection σ inc of the hydrogen atom is high (82 barns).The use of hydrogenated compounds therefore provides a means of examining the dynamics of a proton or an individual molecule and accessing the self-diffusion coefficient With deuterated molecules, the situation is more complex.The cross-section of carbon is totally coherent (5.554 barns), in fact, while that of the deuterium atom is both coherent and incoherent (respectively 5.597 and 2.04 barns).It is more difficult to process the spectra in this case, since the coherent and incoherent contributions must be separated.This is the situation in particular for high loadings, when the incoherent contribution must be taken into account, in order to determine the coefficients D self and D t by analysing the spectra.At low loadings, the incoherent contribution remains negligible and only D t can be extracted.
The coherent structure factor S coh (Q) which governs the intensity due to coherent scattering, is equal to the reciprocal of the thermodynamic factor, Γ, when the momentum transfer Q tends to zero: (5) The extrapolations of S(Q) and Γ can be related since they correspond to measurement of the fluctuations in the number of molecules in a given volume.According to the kinetic theory of gases, particle fluctuations are related to isothermal compressibility.Consequently, we can expect a higher diffusion power at low loadings (where the adsorbed phase is highly compressible, as in a gas) rather than at saturation (similar case to a liquid which is virtually incompressible, giving low QENS intensities).
Neutron scattering experiments on deuterated compounds therefore provide a means of estimating the transport coefficient D t and the thermodynamic factor 1/Γ for each loading θ.The loading dependence of the corrected diffusion coefficient D C can therefore be calculated from the Darken relation ( 6): (6)

MODELLING THE DIFFUSION BY KMC: METHODOLOGY
In the context of this study, an approach by lattice model combined with a Kinetic Monte Carlo (KMC in the remainder of this document) was used to determine the influence of the microscopic phenomena on the overall behaviour of molecules of n-hexane or n-heptane diffusing in the MFI structure zeolite.

Principle of Kinetic Monte Carlo Simulations
Unlike Grand-Canonical Monte Carlo simulations, extensively used to study adsorption in zeolites, Kinetic Monte Carlo simulations, also called Dynamic Monte Carlo simulations, investigate the trajectories of molecules and study properties over time.In KMC simulations, the pore space through which molecules diffuse is modelled by a lattice of adsorption sites connected by bonds.Molecules occupy a certain fraction θ of the sites, and are assumed to hop from one site to the neighbouring sites along connecting bonds.The evolution over time of the system between successive configurations is described by a Markovian master equation: where i and j refer to the configuration of the molecules in the lattice, P's are the probabilities of the configurations, and W 's are the transition probabilities per unit time.The Master Equation can be solved using Kinetic Monte Carlo algorithms (see Binder (1986) for a detailed derivation).The KMC algorithm used in this work is the residence time algorithm (Bortz et al., 1975).For a given configuration of molecules on the lattice, a list of all possible events is generated, i.e. a list of all possible jumps characterised by a jump rate k h .The probability distribution of all elementary jumps is calculated P h = k h /Σ j k j and the ith event is selected, such as P i-1 < υ < P i where υ ∈ [0,1] is a random deviate.The actual time step of the selected event is randomly chosen from a Poisson distribution: Δt = In(u)/Σ j k j where u ∈ [0,1] is another random deviate.The selected event is performed leading to a new configuration, and the list of all possible events is updated.Such Monte Carlo step need to be repeated for a sufficiently large number of configurations and molecular trajectories to be followed.After the system has equilibrated, statistics for the calculation of diffusion properties are collected.

Definition of the Pore Lattice Model
In a lattice model approach, the MFI zeolite is represented by a lattice of adsorption sites connected by bonds.Each node in the lattice, modelling an adsorption site, represents a position of minimum energy, i.e. a position in the zeolite structure where there is a high probability of finding a molecule.The number and the positions of the adsorption sites depend on the system studied (Fig. 3).Special attention must be paid to the local connectivity which varies for each type of site.The connectivity represents the number of connections of a given site with its nearest neighbours.It can be interpreted as being the number of possible directions that the molecule can take to move.In this study, the lattice is considered as independent of the temperature and of the loading θ.

Initial Configuration
The system is studied at thermodynamic equilibrium.Consequently, depending on the required temperature T and loading θ, a fixed number N of molecules is inserted in the lattice.As with molecular dynamics, the KMC simulations are carried out in the statistical assembly 〈N,V,T〉 or 〈N,N sites ,T〉.Since the zeolite volume is constant, a fixed number N sites of adsorption sites is defined.N, the number of molecules, is fixed by the adsorption isotherm.
Since the systems studied are considered as Langmuirian, each site contains only one molecule.The molecules are positioned in the vacant sites either randomly or according to the knowledge acquired concerning their preferential positions in certain adsorption sites.This type of information is provided, in particular, by the Grand Canonical Monte Carlo (GCMC) type simulations.For example, nC6, at high Θ, preferentially position itself in the zigzag channels of MFI silicalite (Vlugt et al., 1999;Pascual et al., 2003) whereas 22DMB only positions itself in the intersections, irrespective of the loading (Pascual et al., 2003).

Jump Rates
The crucial parameter of the KMC algorithm is the jump rate, or transition rate, written k i → j , for the transition from site i to site j.It is calculated using a numerical technique at atomic scale, either by Molecular Dynamics (MD) or by the Transition State Theory (TST).The jump rates resulting from MD and TST simulations therefore contain the influence of microscopic phenomena.Since this type of simulation was not used in the context of this study, the jump rate values are adjusted so as to reproduce the results of the QENS experiment at infinite dilution.Consequently, these jump rates do not contain the influence of the interactions between the molecules.By definition, the KMC algorithm is a stochastic algorithm generating a Markov chain whose events are independent.There is therefore no deterministic relation between the diffusion coefficients and the jump rates, especially since each jump rate, combined with a particular event, is affected locally by the surrounding molecules.Consequently, the parametric estimation of the jump rates is not based on minimising an error function, but uses a Robbins-Monro approximation algorithm implementing a stochastic gradient descent (Bartoli and Del Moral, 2001).

Taking into Account the Interactions Between Molecules
As mentioned above, the jump rates implemented are determined at infinite dilution.Consequently, they do not take into account the interactions between molecules.When the loading increases, however, these interactions must be taken into account.The approach chosen to model these interactions consists in modifying the amplitude of the energy barrier to be crossed (see Fig. 4) depending on the type of interaction with the surrounding molecules, according to Equation (8) (Krishna et al., 2004):  Diagram of the energy reduction to take into account the influence of the nearest neighbours.δE AB represents the reduction of the energy barrier generated by the presence of a molecule in a neighbouring site B (Krishna et al., 2004). ( where n is the number of closest neighbours, ƒ is the Reed-Ehrlich parameter.k A→C is the jump rate from site A to the vacant site C, B represents the neighbouring sites of A which are occupied.
with δ Several profiles of ƒ can be implemented to model the forces of attraction (f < 1) or repulsion (f > 1) between the molecules.ƒ can be constant or variable depending on the loading θ, according to Equation ( 9) (Krishna et al., 2004): In the work presented in this paper, the Reed-Ehrlich parameter (ƒ) is estimated by matching the QENS experimental self-diffusivities (D self ) at the lowest loading available using the Robbins-Monro algorithm.Estimated values of ƒ are determined for both nC6 and nC7.As ƒ represents the intrinsic pair-pair interactions between two molecules, it will be considered constant whatever the loading.Nevertherless the loading effect is taken into account by the number of closest neighbours, n, representing the local environment of each molecule.Implementing a loading dependence of the Reed-Ehrlich parameter (Krishna et al., 2004) can improve simulation-experiment matching but introduces long distance correlations between molecule-molecule interactions.

Calculating the Diffusion Coefficients
The self-diffusion and corrected diffusivity coefficients can be determined using the following relations, based on the Mean Square Displacement calculation: -relative to the individual paths of the molecules: (10) where d represents the movement dimensionality, N the number of molecules, r i (t) and r i (0) the respective position vectors of molecule i at instant t and t = 0; -relative to the path of the centre of mass of the set of molecules: (11) where R(t) is the position vector of the centre of mass of the set of molecules at time t.By recording the positions of each molecule during the simulation, we can calculate the mean square displacements as a function of time.The calculations of the diffusion coefficients using Equations ( 10) and ( 11), obtained from Einstein's law (asymptotic law), can only be made in the range where the mean square displacement is proportional to time.The difficulty therefore consists in identifying the minimum observation time, t obs,min , from which the MSD is proportional to time.The molecules then evolve in the "diffusive" regime and the diffusion coefficients are calculated using Equations ( 12) and ( 13), for t > t obs,min : In addition, Equations ( 12) and ( 13) are only valid if the mean distances travelled by the molecules are greater than the size of a unit cell, in each crystallographic direction.For each simulation, therefore, the mean distances travelled by each molecule are determined in order to validate the assumptions made to calculate the diffusion coefficients.

Adsorption of n-Hexane and n-Heptane
in Silicalite-1

Adsorption Isotherms
Figures 5a and 5b show a comparison of adsorption isotherms respectively for nC6 and nC7, published by various authors.These isotherms were obtained from Grand Canonical Monte Carlo (GCMC) calculations made by the group of Krishna and Van Baten (Jobic et al., 2006) and by Pascual et al. (2003) and, for nC6, also come from Tapered Element Oscillating Microbalance (TOEM) experiments (Zhu et al., 2001).These isotherms are then represented by a Dual Site Langmuir (DSL) model as demonstrated by Millot et al. (1998Millot et al. ( , 1999)), whose mathematical expression is given by Equation ( 14): (14 where indices A and B correspond to the two types of adsorption site of the MFI structure, of different capacity and adsorption force.b A and b B are the model parameters, ratios of the adsorption and desorption rates on the site considered.Θ sat,A and Θ sat,B represent the number of molecules at saturation respectively of sites A and B, Θ(p) is the total number of molecules on all sites A and B at pressure p.
The parameters of the DSL models described are shown in Table 2.
The adsorption isotherm for nC6 exhibits a slight point of inflection, whereas that of nC7 shows a much more pronounced inflection.The point of inflection is the consequence of separate sites of different energies ΔG.The greater the difference in energies, the more pronounced the inflection (Millot et al., 1998).In silicalite-1, some authors explain this by the fact that the size (i.e. the length) of these molecules is close to that of the zigzag channels in the zeolite (the length of nC6 is similar to that of these channels) (Vlugt et al., 1999).(Zhu et al., 2001) The capacities at saturation Θ sat are expressed in molecules/u.c., b A and b B in Pa -1 .

Estimation of the Thermodynamic Correction Factor Γ
The 1/Γ profile against the loading can be calculated from the adsorption isotherm in several ways: -from the DSL model, according to expression ( 15): (15) where Θ A and Θ B represent the number of molecules per unit cell located in sites A and B respectively, and Θ is the total number of molecules per unit cell in silicalite-1; -from the CBMC simulations using the formula derived by Reed and Ehrlich based on the fluctuation in the number of molecules N in silicalite-1: -1/Γ can also be determined experimentally by QENS according to relation (5).Figures 6a and 6b show the profile of 1/Γ against the loading determined by these various methods, for nC6 and nC7 respectively in silicalite-1 at 300 K.
With nC6, the inflection is more pronounced on the profile of 1/Γ against loading Θ (see Fig. 6a) than on the adsorption isotherm calculated by CBMC (see Fig. 5a), because of the derivative involved when calculating 1/Γ.The minimum obtained on the adsorption isotherm of nC7 for Θ = 4 molecules/u.c.results in a virtually zero value for S(Q) and 1/Γ, indicating that the molecules in the silicalite-1 are in stable configuration (see Fig. 6b).
The data resulting from the QENS experiments are not exactly comparable to the inverse of the thermodynamic correction factor.The coherent structure factor, S coh (Q), cannot be measured at Q = 0 and is therefore extrapolated.In addition, integration of S coh (Q) is carried out in this case over the entire accessible transfer energy range ω, and is therefore only partial.This is a major problem at high loadings when the signal, which is more spread out, leaves the measurement window.However, taking into account the experimental error bars, the trends obtained are similar, which confirms that 1/Γ can be obtained by QENS from the coherent structure factor S coh (Q) for Q → 0.
Nonetheless, to guarantee a rigorous approach, in order to estimate the corrected diffusion coefficient D C from the transport coefficient D t according to the Darken Equation ( 6), the calculations are carried out with 1/Γ estimated by Equation ( 15), taking the parameters of the DSL model obtained by Jobic et al. (2006).

Location of n-Hexane and n-Heptane Molecules
in Silicalite-1 Location of the molecule adsorption sites is essential data in KMC simulations.This information can be obtained by analysing the adsorption isotherm and the results of Grand Canonical Monte Carlo (GCMC) simulations.
All publications in the literature, whether experimental or numerical, indicate that about 8 nC6 molecules/u.care adsorbed at saturation (see Tab. 2).They also demonstrate that the nC6 molecule can be positioned at three types of site, intersections, straight channels and zigzag channels.In addition, considering the size of the nC6 molecule (6.4 Å) with respect to the length of the straight channels (4.9 Å) and zigzag channels (7.7 Å) (Pascual et al., 2003), it is assumed that each site can only contain a single molecule.Note, however, that it is difficult to estimate the size of the channels, especially the length of the zigzag channels which are not really straight.Consequently, the lattice chosen consists of unit cells containing 12 adsorption sites (see Fig. 3), 4 intersection sites, 4 sites in the straight channels and 4 sites in the zigzag channels.
However, although the nC6 molecules can be positioned in the 3 types of site defined in Figure 3, it appears that there is no consensus over the loading dependence of the distribution of molecules in the solid.
As mentioned previously, the adsorption isotherm of nC6 in silicalite-1 at 300 K is generally represented by a Dual Site Langmuir (DSL) model, revealing the existence of two types of quite distinct sites (see Fig. 5a and the associated DSL parameters in Tab. 2).The difference in adsorption between these two types of site is related to the ratio of the DSL model Langmuir parameters, b A /b B , (see Eq. 14).The molecules are preferentially positioned in the first (the more stable) then in the second type of site.
The studies conducted by Vlugt et al. (1999), Zhu et al. (2001) and Pascual et al. (2003) demonstrate that the nC6 molecules are uniformly distributed in all sites for Θ < 4 molecules/u.c.Above this value, they are preferentially positioned in the channels through a loss of entropy.The molecules are in fact stabilised in the channels by favourable contact energies with the walls.Although there is more space in the intersections, they are less favourable from the energy point of view since the nC6 molecule has to bend (June et al., 1992).According to the probability density calculations on the presence of nC6 molecules close to saturation carried out by Vlugt et al. (1999), the 8 molecules would be preferentially positioned in the channels.This observation is corroborated by Pascual et al. (2003), who postulate that the size of the nC6 molecules is 6.4 Å, while the length of the straight and zigzag channels is 4.9 Å and 7.7 Å respectively.A molecule positioned in a straight channel would therefore block an intersection site (in terms of adsorption).Since they are roughly the same length as the nC6 molecules, the zigzag channels are preferential adsorption sites compared with straight channels (Vlugt et al., 1999;Zhu et al., 2001).These authors speak of "commensurate freezing" in the zigzag channels.Jacobs et al. (1981), who have a different point of view, demonstrated that the molecules are preferentially positioned in the straight channels.
The authors therefore agree on loading dependence of the distribution of molecules on the adsorption sites, mainly due to a loss of entropy of the molecules.Around 4 molecules/u.c., Profile of the inverse of the thermodynamic factor against loading for nC6 a) and nC7 b) in silicalite-1 at 300 K.The results obtained from the QENS experiments (Eq.5) are compared with those of the CBMC simulations (Eq.16) and those obtained from the DSL model calculations (Eq.15).From Jobic et al. (2006) they observe a phase transition corresponding to a rearrangement of the molecules in the solid.However, none of them has really identified the two types of site defined by the DSL model.
Very few studies concerning adsorption of n-heptane (nC7) in silicalite-1 are available in the literature.As with nC6, the GCMC simulations provide information concerning the distribution of molecules in the various types of site.However, since the size of the nC7 molecule is 7.8 Å, greater than the length of both types of channel, its position is less precise.Vlugt et al. (1999) put forward the same evolution in the distribution of positions of the nC7 molecules as for nC6: the molecules are uniformly distributed for Θ < 4 molecules/u.c. then positioned preferentially in the zigzag channels for very high loadings.Pascual et al. (2003) adopt a different point of view.According to these authors, the nC7 molecules first fill the zigzag channels up to Θ = 4 molecules/u.c.Since the size of the molecule is greater than that of the zigzag channel, the molecule also occupies a large volume of the intersection, blocking the insertion of another molecule in the structure, which accounts for the plateau on the adsorption isotherm.The adsorbed molecules must undergo a rearrangement (a loss of entropy is required) before other nC7 molecules can be inserted.

Approaches Tested in KMC Concerning the Adsorption Sites
In the KMC simulations, MFI zeolite is represented by a lattice of adsorption sites connected by bonds.To define a lattice representative of the system, it is essential to know which of these sites correspond to positions of minimum energy, i.e. locate the zones where there is a high probability of finding a molecule (see Fig. 3).
Three types of molecule distribution were tested, in order to study the influence of this distribution on the loading dependence of the diffusion of the molecules: -Approach 1: All adsorption sites have the same energy level.This approach consists in ignoring the existence of the two types of sites identified by the nC6 and nC7 isotherms.
This assumption is nevertheless in agreement with the studies conducted by Vlugt et al. (1999), Zhu et al. (2001) and Pascual et al. (2003) demonstrating that the nC6 molecules are uniformly distributed in all sites for Θ < 4 molecules/u.c.A distinction is made, however, concerning the rates along the straight channels and along the zigzag channels of the MFI; -Approach 2: The energy of the sites in the channels is lower than that of the intersection site.This approach considers the channels as being the most stable sites (site A of the DSL model), with no distinction between the channel types, and the intersection site is site B identified in the DSL model; -Approach 3: The sites in the straight channels and in the zigzag channels have different energy levels, identified at sites A and B of the DSL model, and the intersection site is a site of higher energy.This approach considers the channels as being sites A and B of the DSL models and the intersection site is a passage site.The cases identifying site A as being the straight channels firstly, and the zigzag channels secondly, have both been tested.

Diffusion of n-Hexane and n-Heptane in Silicalite-1
The influence of loading on the transport coefficient for nC6 and nC7 in silicalite-1 at 300 K was determined from the spectra obtained by QENS using the deuterated compounds C 6 D 14 and C 7 D 16 .The values of D self calculated for nC6 are derived from experiments on the hydrogenated compounds as well as on the deuterated molecules at high loading.Under these conditions, in fact, the incoherent contribution of the signal is no longer negligible and can be taken into consideration.Furthermore, the influence of loading on the corrected diffusion coefficient D C was determined from the Darken relation (6) using the thermodynamic correction factor Γ obtained from (5).
The results are given in Tables 3 and 4 and in Figures 7  and 8.
Special attention is required when extracting the experimental diffusivities from QENS spectra.When S(Q) is high, the QENS spectra measured using deuterated molecules are dominated by coherent scattering.The data can then be adjusted with a single Lorentzian, corresponding to the transport diffusion.When S(Q) is low, the spectra can no longer be adjusted with a single Lorentzian (Jobic et al., 2006).To take into account the fact that the incoherent signal is now no longer negligible, the entire spectrum must be modelled with 2 Lorentzians corresponding to the coherent and incoherent contributions.D self can be determined by interpreting the incoherent signal.The error bars on the measurements of D t at high loadings are therefore larger due to the use of two Lorentzians.The slight inflection in the nC6 isotherm, evidenced by the inflection in 1/Γ (see Fig. 6a), is not strong enough to influence the behaviours of D t , D C and D self (see Fig. 7).The experimental values of D self for nC6 were determined using the incoherent scattering generated by the hydrogenated molecules and from the incoherent contribution observed at high loading for the deuterated compounds.(molecules/u.c.) (×10 -10 m 2 s -1 ) (×10 -10 m 2 s -1 ) (×10 -10 m 2 s -1 ) (×10 -10 m 2 s -1 ) 1 8.47 ± 1.7 6.9 --2 8.9 ± 1.8 The inflection observed in the nC7 isotherm for Θ = 4 molecules/u.c.(see Fig. 5b, 6b) generates a marked maximum for the transport coefficient in the same loading range (see Fig. 8).
This maximum is similar to that observed by Shah et al. (1995) for the transport coefficient of benzene in silicalite-1 at Θ = 4 molecules/u.c.The isotherm inflection also generates an inflection in the corrected diffusivity D C , but no inflection is observed on D self for Θ = 4 molecules/u.c.especially since only the values at high loadings were obtained.
The behaviours of D C measured for nC6 and nC7, corroborated by the behaviour of ethane in silicalite-1 at 300 K (Chong et al., 2005), confirm that the approximation made on the Darken Equation ( 6) according to which D C is independent of the loading is not valid, unlike the behaviour of CH 4 simulated by molecular dynamics, which exhibits a constant D C at 300 K (Maginn et al., 1993).
The interchange coefficients D 11 show in both cases trends similar to those of the corresponding values of D self (see Fig. 7, 8).Consequently, the assumption sometimes made that D 11 = D C is valid neither for nC6 nor for nC7 in silicalite.Note, however, the negative value for nC7 at 4.5 molecules/u.c,which is not shown on the graph for reasons of clarity.It is generated by the corresponding high value of D C .
Experiments have been conducted at various temperatures (from 230 K to 350 K) to determine the activation energy of n-C6.The diffusivities vary exponentially with the activation energy E a .By plotting the log of the diffusion coefficient against 1/T, therefore, E a can be calculated from the slope of the linear regression line (see Fig. 9).Activation energies of about 5.8 kJmol -1 and 5.1 kJmol -1 are obtained, respectively for D t (1 molecule/u.c.) and for D self (2 molecules/u.c.), in agreement with the molecular dynamics results (Maginn et al., 1996).

Diffusion Coefficients D t , D C and D self Obtained by KMC
Several approaches based on the differences in energy levels of the adsorption sites were developed in order to conduct KMC studies on the impact of the adsorption isotherm on the loading dependence of nC6 and nC7 diffusivities in zeolite.
For each approach tested, expressions to relate the jump rates to the DSL model parameters were developed and the jump rates were estimated at infinite dilution in order to reproduce the experimental results.
Approach 1: all adsorption sites have the same energy level This approach amounts to considering that the energy levels of all sites are identical.Consequently, the jump rates are not related to the DSL model parameters.Nevertheless, bearing in mind that the molecules circulate more easily in the straight channels than in the zigzag channels, two different jump rates are implemented, depending on the direction of movement of the molecules.The only way to discriminate the movement of the molecules towards either the straight or zigzag channels is through the use of the energy barriers to be crossed (see Fig. 10).A ratio K dir is introduced, expressing this direction anisotropy: where k i ↔ j represents the jump rate from site i to site j and vice versa.
An approximation to the ratio k str /k zz can be obtained using Equation ( 19), provided that the diffusion coefficients in the crystallographic directions x and y are known: where a and b represent the true dimension of the unit cell, x and y are the crystallographic axes.

Site str
Site inter E a,inter ↔ str E a,inter ↔ zz E(q) q Site zz Figure 10 Diagram of the approach implementing anisotropic rates as a function of the type of direction.E a,i↔j represents the activation energy corresponding to the jumps between sites i and j. inter, str and zz correspond to the sites in the intersections, the straight channels and the zigzag channels respectively.not depend only on the anisotropy of the solid, but also on the correlations between the molecules.Consequently, the ratios determined at infinite dilution using the MD data of Runnebaum and Maginn (1997) were chosen: K dir = k str /k zz = 6.0 for nC6 and K dir = 9.5 for nC7.
Firstly, in order to obtain a basis for comparison, simulations were carried out without interaction between the molecules (ƒ = 1.0).The profiles of D self and D C both decrease as a function of Θ, the number of molecules per unit cell.The curve of the normalised D self (see Fig. 11b) is convex whereas the profile of the normalised D C indicates a linear decrease as a function of Θ: where θ represents the coverage rate, i.e. the ratio between the number of molecules inserted, Θ, and the number of sites in the lattice.The linear decrease observed on the normalised D C can be represented by an equation similar to Equation (20), indicating conventional behaviour.However, the occupancy fraction θ involved in this equation is no longer calculated with respect to the capacity at saturation of nC6 but with respect to the number of adsorption sites per unit cell.Extrapolation of the curve shows that D C tends to zero for Θ = 12 molecules/u.c.However, the trends simulated do not agree with the experimental trends measured by QENS.
To obtain different behaviour, in particular different from that of linear decrease as a function of Θ for D C , interactions between the molecules were introduced in the model.To obtain an idea of the value of parameter ƒ, it was estimated using the Robbins-Monro algorithm.This estimation was carried out for Θ = 2 molecules/u.c., i.e. for the lowest possible Θ for which experimental data of D self is available.The estimated value of ƒ is 0.62.This value was implemented for simulations over the entire range of Θ.The profiles of the normalised diffusion coefficients, D self and D C , obtained are convex curves, in agreement with the behaviour observed for CH 4 for ƒ < 1 in silicalite (Krishna et al., 2004).The trends determined agree qualitatively with the experimental trends for low loadings, but a difference still remains for high loadings.The sharp decrease observed experimentally when approaching saturation, especially with D self , is not reproduced.Although parameter ƒ, modelling attraction between the molecules, tends to reduce the diffusion coefficients by increasing the residence time of the molecules in the adsorption sites, the lattice considered has 12 adsorption sites.Consequently, even under conditions close to saturation, at least a third of the sites are always vacant.The molecules therefore retain the possibility of moving to a large number of sites, with non-negligible probabilities.
For nC7, as when studying nC6, simulations were initially carried out without taking into account the interactions between the molecules.
Constant interactions were then introduced.These interactions were estimated using the Robbins-Monro algorithm for Θ = 3.6 molecules/u.c., i.e. for the lowest possible Θ for which experimental data of D self is available.The estimated value of f is 0.45.This value was implemented for simulations over the entire range of Θ.The results of Figure 12 show a linear decrease as a function of (1 -θ).The marked experimental decrease at high loadings and the minimum value of D C for Θ = 4 molecules/u.c. are not reproduced.The deviations with respect to the experimental values are larger for nC7 than for nC6.The explanations provided for the nC6 simulations are also valid for nC7.
Approach 2: the energy of the sites in the channels is lower than that of the intersection site This approach assumes that the two types of site defined by the DSL model consist firstly of the intersections and secondly of the sites in the undifferentiated channels, i.e. no distinction between the two types of channel (see Fig. 13).It is corroborated by the MD simulations of June et al. (1992) which demonstrate that the channels are energetically more favourable sites than the intersections.
In order to take into account the preferential position of the molecules in the channels at high loadings, the jump rates from the channels are reduced with respect to those of the inverse jump, introducing the ratio of Langmuir parameters of the DSL model b A /b B according to Equation ( 21): (21) where k i → j represents the jump rate from site i to site j.
In addition, the ratio between the jump rates along the straight channels and along the zigzag channels is still K dir : Influence of parameter f on the loading dependence of the self-diffusion coefficient, D self a), and the corrected diffusivity coefficient, D C b), for nC7 in silicalite-1, implementing anisotropic jump rates depending on the directions of movement (approach 1).

Site
Diagram of the approach implementing anisotropic rates as a function of the type of site, considering the two sites defined by the DSL model as being the undifferentiated channels and the intersections.E a,i↔j corresponds to the activation energy to move from site i to site j, related to the jump rate k i→j by the Arrhenius type relation.
inter, A and B correspond to the sites respectively in the intersections, and in the two sites A and B defined by the DSL model.E* inter↔i corresponds to the transition energy between the intersection and site i.
For nC6, in this approach, the ratio b A /b B = 16.7 obtained by Jobic et al. (2006) was considered.The set of rates is adjusted to reproduce the value of D self at infinite dilution supplied by QENS: D self = 8.0 × 10 -10 m 2 s -1 .The impact of interactions with the neighbouring molecules was studied by performing a first series of simulations with no interaction (f = 1.0) and a second series with an estimated factor f for Θ = 2 molecules/u.c.(f = 0.2).
The simulation results shown in Figure 14 indicate a general decrease in D self and D C as a function of the loading Θ.A linear decrease as a function of (1 -θ), where θ represents an occupancy fraction, is observed for the profiles of the normalised coefficients D self and D C during the simulations performed with no interaction.The decrease relative to D self is more marked than that relative to D C , in accordance with the results obtained from the QENS experimental measurements.The simulated values of the two coefficients nevertheless overestimate the experimental values.
During these simulations, the mean residence time of a molecule in the channels is higher than that in the intersections.This observation agrees with the jump rates implemented which, being unfavourable to the jump rates from the channels to the intersections, favour the mean residence time in the channels.However, the greater the number of molecules inserted in the lattice, the greater the probability of finding molecules in the intersections, necessary point of passage when moving between the various channels.The mean residence time of the molecules in the intersections therefore increases.Note that this increase is more marked for Θ > 4 molecules/u.c., the loading at which half of the channels would be occupied if all the molecules were in the channels.
If a constant value of ƒ less than 1 is implemented (ƒ = 0.2 estimated at Θ = 2 molecules/u.c.), the trends observed for the normalised diffusion coefficients deviate from the linear decrease observed for ƒ = 1.0.At low values of Θ, the simulated values are closer to those measured experimentally, especially for D self .Since attraction between molecules is taken into account, however, the underestimation of the two coefficients for Θ > 2 molecules/u.c.becomes significant.
At low loadings, the molecules position themselves preferentially in the channels, the jump rates undergoing no substantial modification caused by their local environment.As the loading increases, the mean residence time of a molecule in the intersections increases for two main reasons.Firstly, the greater the number of molecules inserted, the higher the possibility of finding a molecule in an intersection site (fewer vacant sites).Secondly, at high loading, as soon as a molecule positions itself in an intersection, it is subjected to the channel-intersection attractive forces implemented via parameter ƒ.The jump rates of the events associated with this molecule are decreased, which increases its residence time in the intersection.The attraction phenomena implemented compensate for the preferential location of the molecules in the channels dictated by the jump rate ratios.Consequently, the resulting distribution of molecules is not in agreement with the MD observations.
Values of f varying with the loading were then tested, based on the work of Krishna et al. (2004).The profile of f varies exponentially and corroborates ƒ = 0.2 at Θ = 2 molecules/u.c.. Since this profile generates a change in the attractive and repulsive interactions at higher loading, the transition observed can be approximated at around Θ = 4 molecules/u.c.Influence of parameter f on the loading dependence of D self , a), as well as of D C , b), for single-component nC7 in silicalite-1, using approach 2 and ratio b A /b B = 2647.The KMC results are compared with the QENS results of Jobic et al. (2006).
Diagram of the approach implementing anisotropic rates as a function of the type of site, considering the two sites defined by the DSL model as being the two types of channel.The intersection corresponds to a "passage" site.E a,i↔j corresponds to the activation energy to move from site i to site j, related to the jump rate k i→j by the Arrhenius type relation.inter, A and B correspond to the sites respectively in the intersections, and in the two sites A and B defined by the DSL model.E* inter↔i corresponds to the transition energy between the intersection and site i. for D C .In addition, a sharp decrease in D self is observed for loadings close to saturation, which was not obtained on the previous simulation results, related to very high (even unrealistic: ≥ 53) values of parameter ƒ.Nonetheless, implementing a profile of interaction between paraffin molecules changing from attractive to repulsive phenomena does not seem realistic and is probably a way of making up for a phenomenon which was not explicitly taken into account in the modelling.
With nC7, estimating the interactions at 3.6 molecules/u.c.results in a value much less than 1: ƒ = 0.05 necessary to compensate for the very high rates from the intersections to the channels, associated with the very high ratio b A /b B (b A /b B = 2647 (Jobic et al., 2006)).The trends obtained partially agree with the experimental trends (see Fig. 15).However, the marked drop of D self and the minimum value of D C for Θ = 4 molecules/u.c. are not reproduced.
Approach 3: sites A and B of the DSL model are the channels and the intersection site is a site of higher energy This time, sites A and B evidenced by the DSL model are assumed to be the straight channels and the zigzag channels.The intersection site, being energetically unfavourable, is considered as a "passage" site.Each case, considering either the straight channels or the zigzag channels as being site A (the most stable site), is tested.Equations ( 23) and ( 24) detail the implemented relations between jump rates.In this approach, ratio b A /b B favors the moves from an intersection to site A rather than moves to site B (see Fig. 16).The already defined parameter K dir (Eq. 18), expressing the direction anisotropy between straight and zigzag channels is also used.Moreover, another parameter, k int needs to be introduced in order to determine all jump rates.This parameter, defining a relation between the energy levels of the intersection sites and sites B, allows to tune the relative affinity of a molecule for these two types of sites.
Site A = zigzag channels and site B = straight channels:  Loading dependence of normalised D self and D C for nC7 in silicalite-1, using approach 3, and assuming that the straight channels are the most stable site (site A).Jump rates defined by Equation ( 23 Figure 18a details the loading dependence of the positions of the molecules.It can be compared with the distribution obtained at equilibrium between sites A and B (Fig. 18b) using the DSL model with parameters b A and b B obtained from Jobic et al. (2006).
The results of Figure 17 clearly demonstrate that the model implementing relations (23) and considering the zigzag channels as stable site is unsuitable, since the profiles of D self and D C diverge for Θ ≥ 4 molecules/u.c.In contrast, the simulations assuming the straight channels as stable site give results which follow a good trend.Loading dependence of D self , a) and b), and of D C , c) and d), for nC7 in silicalite-1, using a combination of the two approaches."Direction" and "Site" correspond to the implementation of anisotropic jump rates depending respectively on the type of direction (approach 1) and the type of site (approach 2).a), c) show the absolute values and; b), d) the values normalised by the value at infinite dilution.The KMC results are compared with the QENS results of Jobic et al. (2006).
For nC7, only the case considering the straight channels as being the most stable site (site A) was taken into account, in view of the results obtained on nC6.These simulations were performed by implementing jump rates according to Equation ( 24), with b A /b B = 2647 and k int = 2650.The trends observed of coefficients D self and D C simulated by KMC do not quite agree with those determined experimentally (see Fig. 19).Although the values concerning D C agree qualitatively, those concerning self-diffusion underestimate the experimental values in particular.In this approach, the mean residence times of the molecules demonstrate that they occupy only the straight channels for Θ ≤ 4 molecules/u.c.Above this value, molecules are transferred to the zigzag channels, the intersections being occupied for only a very short period of time.
Combined approach: for nC7 exhibiting a strong inflection in the isotherm For nC7, characterised by a very strong inflection in the adsorption isotherm (and therefore a very high ratio b A /b B ), an approach combining approaches 1 and 2 was tested.
This approach was developed so as to reproduce the probability density distributions of the presence of molecules in each type of site.As we have seen previously, therefore, approach 1 where the sites are assumed to have the same energy level, generates a virtually uniform residence time distribution on the various types of site, whereas heterogeneous distributions are obtained when the sites have different energy levels.
Consequently, a model consisting of a combination of the two approaches was tested: -approach 1: all sites with the same energy level for Θ ≤ 3.6 molecules/u.c. with implementation of parameter ƒ = 0.45, giving the best results for this loading range; -approach 2: for Θ ≥ 3.6 molecules/u.c,the channels are assumed to be sites (site A) which are more stable than the intersections (site B).Since the jump rates of approach 1 have already been estimated at infinite dilution, a set of rates for approach 2 guaranteeing the continuity of the two approaches must be implemented.The new set of rates in approach 2 was therefore estimated using the Robbins-Monro algorithm in order to reproduce the mean residence time of a configuration at Θ = 3.6 molecules/u.c.generated by approach 1.Since the estimated jump rates of approach 2 depend on the presumed value of parameter ƒ at the presumed point of continuity, two estimations were carried out for ƒ = 0.45 (value estimated previously) and ƒ = 1.0.
The results are given in Figure 20.
The continuity of the two approaches for Θ = 3.6 molecules/u.c. is respected, to within the uncertainty values.The simulations for Θ ≤ 3.6 molecules/u.c.agree qualitatively with the experimental results, with a more marked deviation for D C .Those for Θ ≥ 3.6 molecules/u.c. also agree qualitatively with the general experimental trend for the two coefficients.
The marked decrease at high Θ for D self and the minimum value of D C , however, are not reproduced.Moreover, the interaction forces have little influence on the two coefficients.In approach 2, in fact, since the jump rates from the channels to the intersections are very small compared with those in the opposite direction, the molecules are trapped in the channels.Consequently, the interaction forces, which apply only between molecules located in a channel and an intersection, have only very little influence on the movements of the molecules.

CONCLUSIONS
This study concerned two systems, n-hexane and n-heptane in silicalite-1, MFI structure zeolite, at 300 K, characterised by Dual Site type isotherm.This work is based on an experimental study using Quasi-Elastic Neutron Scattering and on a Kinetic Monte Carlo modelling study.
Several points were identified from the experiments conducted using QENS: -the profile of the thermodynamic correction factor 1/Γ against Θ obtained by QENS from the coherent structure factor S coh (Q) for Q → 0 is comparable, to within experimental errors, with the profiles obtained by MD modelling or from the DSL model parameters; -the study of the profile of 1/Γ against Θ identifies the presence of Dual Sites more clearly than the isotherm alone (especially for nC6 in silicalite-1); -the influence of the loading θ on the diffusion coefficients D self , D t and D C of hexane and heptane in silicalite-1 at 300 K was studied, using deuterated molecules C 6 D 14 and C 7 D 16 as well as C 6 H 14 .In addition, the interchange coefficient, -D 11 (θ), derived from the single-component Maxwell-Stefan formulation, was determined; -the slight inflection in the nC6 isotherm at 300 K, evidenced by the inflection in 1/Γ, is not strong enough to influence the behaviours of D t , D C and D self ; -the inflection observed in the nC7 isotherm for Θ = 4 molecules/u.c.generates a marked maximum for the transport coefficient in the same loading range.The inflection in the isotherm also generates an inflection in the corrected diffusivity D C , but no inflection is observed on D self for Θ = 4 molecules/u.c;-the behaviours of D C measured for nC6 and for nC7 in silicalite-1 at 300 K confirm that D C cannot be considered as independent of the loading.
The KMC modelling study showed how the energy levels of the molecule adsorption sites in a zeolite affect the loading dependence of the diffusion coefficients of these molecules: -since two different types of sites A and B were detected by the DSL model of the adsorption isotherm, these sites must be identified amongst the possible adsorption sites of nC6 and nC7 in an MFI, i.e. the straight channels, the zigzag channels and the intersections; -if all the adsorption sites have the same energy levels, the molecules are distributed uniformly on each possible site and neither the marked decrease in D self close to saturation, nor the possible inflection in D C (nC7) can be reproduced; -if the straight and zigzag channels are identified as having lower energy (stable site A) and the intersections as being site B, the simulated diffusion coefficients are overestimated.The trend can nevertheless be improved by introducing locally attractive interaction forces or changing from attractive to repulsive as a function of the loading, but this seems to be mainly a way of making up for a phenomenon which was not explicitly taken into account in the modelling; -for nC6 in silicalite-1, characterised by a slight inflection in the isotherm, the evolution in the simulated diffusion coefficients gives correct trends if sites A and B of the DSL model are identified as being the straight channels and the zigzag channels, the intersections being considered as a passage site.In contrast, if we assume that the zigzag channels are the most stable site, the simulated diffusion coefficients diverge; -for nC7 in silicalite-1, characterised by a strong inflection in the isotherm, if the straight channels as identified as being site A and the zigzag channels as site B, the simulated diffusion coefficients are underestimated.In this case, another approach must be considered, assuming that the energy levels of the sites depend on the loading.At low loading (Θ ≤ 3.6 molecules/u.c.), all sites have the same energy level, so the molecules can be uniformly distributed in the various possible sites.At Θ = 3.6 molecules/u.c., the molecules undergo a rearrangement in the zeolite since the profile 1/Γ disappears.Above Θ = 3.6 molecules/u.c,it is assumed that the channels become more stable than the intersections.The results obtained agree qualitatively with the general experimental trend for the two coefficients.The marked decrease at high Θ for D self and the minimum value of D C , however, are not reproduced.
Figure 3Representation of the cell unit defined for the KMC simulations of nC6 in MFI structure zeolite.The straight channels (dark grey circles) and the zigzag channels (light grey circles) are interconnected, defining intersections (black circles).The unit cell parameters are a = 20.1 Å, b = 19.9Å and c = 13.4Å.

N
Figure 5Adsorption isotherms estimated by the DSL model of nC6 a) and nC7 b) in silicalite-1 at 300 K.The DSL parameters were estimated using CBMC simulations conducted byPascual et al. (2003) andJobic et al. (2006) or TOEM experiments conducted byZhu et al. (2001).The loading is expressed in molecules per unit cell and the pressure in Pa.
Figure 6 Figures 7 and 8 represent the influences of the loading for D t , D C and D self , for nC6 and nC7 respectively, calculated from the QENS experimental results.From data on D self and D t , the interchange coefficient D 11 obtained from the Maxwell-Stefan formulation was determined according to equation (17):

N
Laloué et al. / Kinetic Monte Carlo Modelling to Study Diffusion in Zeolite: Understanding the Impact of Dual Site Isotherm on the Loading Dependence of n-Hexane and n-Heptane Diffusivities in MFI Zeolite, as Revealed by QENS Experiments Figure 7Influence of the loading on coefficients D t , D C and D self and the interchange coefficient D 11 of nC6 in silicalite-1 at 300 K obtained byQENS (Jobic et al., 2006).
Figure 8Influence of the loading on coefficients D t , D C and D self and the interchange coefficient D 11 of nC7 in silicalite-1 at 300 K obtained byQENS (Jobic et al., 2006).

N
Laloué et al. / Kinetic Monte Carlo Modelling to Study Diffusion in Zeolite: Understanding the Impact of Dual Site Isotherm on the Loading Dependence of n-Hexane and n-Heptane Diffusivities in MFI Zeolite, as Revealed by QENS Experiments 783 Figure9 Figure 11Influence of parameter f on the loading dependence of D self , a) and b), as well as of D C , c) and d), for nC6 in silicalite-1, using approach 1. a) and c) represent the absolute values and; b) and d) the values normalised by the value at infinite dilution.The KMC results are compared with the QENS results ofJobic et al. (2006).

N
Figure12 Figure 14 Influence of f on the profiles of D self , a) and b), and of D C , c) and d), as a function of Θ for nC6 in silicalite-1 using approach 2 and ratio b A /b B = 16.7.a) and c) represent the absolute values and; b) and d) the values normalised by the value at infinite dilution.The KMC results are compared with the QENS results of Jobic et al. (2006).

N
Figure15 Figure 17 Loading dependence of D self , a) and b), and of D C , c) and d), for nC6 in silicalite-1, using approach 3. The most stable site A is indicated.Ratio b A /b B and parameter k int are kept respectively at values of 16.7 and 17.0.Parameter f varies.a) and c) represent the absolute values and; b) and d) the values normalised by the value at infinite dilution.The KMC results are compared with the QENS results of Jobic et al. (2006).
Figure 18 a) Mean residence time spent in each site (%) corresponding to the simulations shown in Figure 17; b) relative occupation on site A and site B using DSL model for nC6 isotherm.
Figure 19 ) with b A /b B = 2647 and k int = 2650.The KMC results are compared with the QENS results of Jobic et al. (2006).Site A = straight channels and site B = zigzag channels: (24) All the simulations described for nC6 are based on several parameters: definition of the most stable site A, k int , ƒ, and ratio b A /b B .The simulations were performed to test relations (23) and (24) considering either the sites in the straight channels or those in the zigzag channels as stable.They were performed with ratio b A /b B = 16.7(Jobic et al., 2006) and parameter k int = 17.The influence of parameter ƒ was studied for As previously, Figure17shows the loading dependence of coefficients D self and D C for the various cases tested. Figure20

TABLE 1 Loadings
Laloué et al. / Kinetic Monte Carlo Modelling to Study Diffusion in Zeolite: Understanding the Impact of Dual Site Isotherm on the Loading Dependence of n-Hexane and n-Heptane Diffusivities in MFI Zeolite, as Revealed by QENS Experiments

TABLE 4
Self-diffusion coefficients D self , corrected D C and transport D t diffusivities and interchange coefficient D 11 of C 7 D 16 at 300 K in silicalite-1 measured by QENS N Laloué et al. / Kinetic Monte Carlo Modelling to Study Diffusion in Zeolite: Understanding the Impact of Dual Site Isotherm on the Loading Dependence of n-Hexane and n-Heptane Diffusivities in MFI Zeolite, as Revealed by QENS Experiments