Wanda Kellouai

Patrick Judeinstein

Marie Plazanet

Simon Baudoin

Martin Drobek

Anne Julbe

Benoit Coasne
email: benoit.coasne@univ-grenoble-alpes.fr

Gas Adsorption in Zeolite and Thin Zeolite Layers: Molecular Simulation, Experiment, and Adsorption Potential Theory

come

INTRODUCTION

Gas adsorption in nanoporous materials such as zeolites, metal organic frameworks and active carbons has been studied extensively using both experimental and theoretical methods. [START_REF] Raptopoulou | Metal-Organic Frameworks: Synthetic Methods and Potential Applications[END_REF][START_REF] Van Der Voort | Introduction to Porous Materials[END_REF][START_REF] Terranova | Carbon Nanomaterials for Gas Adsorption[END_REF][START_REF] Valencia | New Developments in Adsorption/Separation of Small Molecules by Zeolites[END_REF] From a theoretical viewpoint, such research efforts aim at elucidating the impact of confinement and surface forces on the thermodynamics and dynamics of nanoconfined fluids which drastically depart from their bulk counterpart. [START_REF] Coasne | Adsorption, intrusion and freezing in porous silica: the view from the nanoscale[END_REF][START_REF] Smit | Molecular Simulations of Zeolites: Adsorption, Diffusion, and Shape Selectivity[END_REF][START_REF] Van Speybroeck | Advances in theory and their application within the field of zeolite chemistry[END_REF] From a practical viewpoint, studies in this field are motivated by the increasing role played by such nanoporous solids in energy and environmental applications (catalysis and adsorption technologies, membrane processes for separation/filtration, energy storage/conversion, etc.). [START_REF] Van Der Voort | Introduction to Porous Materials[END_REF] Among available nanoporous materials, zeolites keep receiving a great deal of attention owing to their tunable large internal surface and nanometric pore size combined with their good thermal, mechanical and chemical stability. [START_REF] Davis | Design for sieving[END_REF][START_REF] Valtchev | Tailored crystalline microporous materials by post-synthesis modification[END_REF][START_REF] Weckhuysen | Recent advances in zeolite chemistry and catalysis[END_REF] Moreover, recent progress boost in the synthesis and shaping of zeolites have led to a new momentum in the field with the advent of novel materials such as zeolite films, [START_REF] Pina | Zeolite films and membranes. Emerging applications[END_REF] nanozeolites, [START_REF] Mintova | Advances in nanosized zeolites[END_REF] hollow zeolites, [START_REF] Pagis | Hollow Zeolite Structures: An Overview of Synthesis Methods[END_REF] hierachical zeolites, [START_REF] Valtchev | Hierarchical zeolites[END_REF][START_REF] Pérez-Ramírez | Zeolite Catalysts with Tunable Hierarchy Factor by Pore-Growth Moderators[END_REF] dendritic zeolites, [START_REF] Zhao | Engineering Fractal MTW Zeolite Mesocrystal: Particle-Based Dendritic Growth via Twinning-Plane Induced Crystallization[END_REF] etc.

While the confinement of fluids in the core porosity of zeolites is rather well-understood, the design of the novel zeolitic materials mentioned above raises additional basic and practical questions related to adsorption in their complex architecture. [START_REF] Mitchell | Structural analysis of hierarchically organized zeolites[END_REF][START_REF] Coasne | Molecular Simulation of Adsorption and Transport in Hierarchical Porous Materials[END_REF][START_REF] Cychosz | Recent advances in the textural characterization of hierarchically structured nanoporous materials[END_REF] Indeed, a common point to nanozeolites, hollow zeolites and hierarchical zeolites is their large external surface through which different porosities connect in a more or less topologically ordered fashion (we note that even regular zeolite powders already possess a large external surface area but to a lesser extent compared to the sample types listed here). The impact of this external surface is already documented with striking examples in adsorption uptakes measurements [START_REF] Saint Remi | The role of crystal diversity in understanding mass transfer in nanoporous materials[END_REF][START_REF] Brabec | Silicalite-1 Crystals Etched with Hydrofluoric Acid Dissolved in Water or Acetone[END_REF] as well as in separation/chromatography where it plays a key role in the balance between intra-particle diffusion and hydraulic transport [START_REF] Rybka | Molecular Dynamics Study of the Relation between Analyte Retention and Surface Diffusion in Reversed-Phase Liquid Chromatography[END_REF][START_REF] Melnikov | A Molecular Dynamics Study on the Partitioning Mechanism in Hydrophilic Interaction Chromatography[END_REF] (see also Ref. [START_REF] Bellat | Microporous volume and external surface of Y zeolites accessible to p-xylene and m-xylene[END_REF] for a discussion on the adsorption of xylene isomers at the external surface of zeolites). In catalysis, the effect of the external surface remains to be fully established but many papers report pore mouth catalysis mechanisms. [START_REF] Martens | Selective Isomerization of Hydrocarbon Chains on External Surfaces of Zeolite Crystals[END_REF][START_REF] Zecevic | Nanoscale intimacy in bifunctional catalysts for selective conversion of hydrocarbons[END_REF] In this context, the geometry and defects at the external surface was shown to affect both the catalysis of reacting molecules (e.g. Refs. [START_REF] Chizallet | Achievements and Expectations in the Field of Computational Heterogeneous Catalysis in an Innovation Context[END_REF][START_REF] Chizallet | Toward the Atomic Scale Simulation of Intricate Acidic Aluminosilicate Catalysts[END_REF] ) and the adsorption/dynamics of fluids (e.g. Refs. [START_REF] Coasne | Molecular Simulation of Adsorption and Transport in Hierarchical Porous Materials[END_REF][START_REF] Kortunov | Influence of Defects on the External Crystal Surface on Molecular Uptake into MFI-Type Zeolites[END_REF] ) in nanoporous solids.

Several authors have considered numerical approaches to assess the impact of the external surface on gas adsorption and diffusion (while a detailed review is out of the scope of the present paper, we provide in the following some important examples). Vlugt and coworkers [START_REF] García-Pérez | External Surface Adsorption on Silicalite-1 Zeolite Studied by Molecular Simulation[END_REF] computed surface adsorption isotherms for mixtures which were compared with adsorption in the zeolite core using surface excess concentrations. Using a simple in-silico model of a hierarchical zeolite, molecular dynamics and Monte Carlo simulations were also used to investigate adsorption and diffusion in systems displaying a large external surface area. It was shown that the overall adsorption isotherm can be decomposed as the sum of adsorption at the external surface and inside the zeolite core porosity, therefore providing a simple mean to rationalize adsorption in such samples. [START_REF] Coasne | Molecular Simulation of Adsorption and Transport in Hierarchical Porous Materials[END_REF][START_REF] Deliere | Adsorption in heterogeneous porous media: Hierarchical and composite solids[END_REF] Using a similar set-up, Inzoli et al. [START_REF] Inzoli | Surface Adsorption Isotherms and Surface Excess Densities of n -Butane in Silicalite-1[END_REF] showed by means of molecular dynamics simulations that intercrystalline diffusion in the zeolite porosity is a non isothermal process in which the driving force inducing transport from the external surface is governed by the adsorption enthalpy.

Despite the acknowledged impact of the external surface, a detailed picture of adsorption and transport across nanoporous materials displaying large external surface areas is still missing. In particular, very practical questions remain unanswered such as the influence on gas adsorption of the surface geometry (orientation, flexibility, etc.). Moreover, while theoretical frameworks are available to describe gas adsorption in zeolites, a unified formalism to model in a consistent fashion adsorption in the zeolite core and at its external surface is lacking. To gain fundamental insights into these issues, the present work reports a joint experimental and molecular simulation study on gas adsorption in zeolite materials. In more detail, using a prototypical system consisting of methane confined in silicalite-1, [START_REF]Database of zeolite structure -MFI[END_REF] we investigate gas adsorption at different temperatures in the zeolite core porosity as well as at the external surface of zeolite thin layers with different crystallographic orientations. Using our simulated and experimental data, we first extend the adsorption potential theory -originally developed by Polanyi to describe adsorption of condensable vapor onto surfaces [START_REF] Polanyi | Verhandlungen der Deutschen Physikalischen[END_REF][START_REF] Polanyi | The Potential Theory of Adsorption[END_REF][START_REF] Neimark | Adsorption-Induced Deformation of Microporous Solids: A New Insight from a Century-Old Theory[END_REF] -for fluids under supercritical conditions. Then, we show that this powerful adsorption theory can be used to predict adsorption at different temperatures both at the external surface and in the core porosity of zeolites.

The remainder of this paper is arranged as follows. In Section 2, we briefly introduce the experimental and computational methods used to investigate methane adsorption in silicalite-1 zeolite samples. In Section 3, both the simulated and experimental adsorption isotherms are presented. We first discuss methane adsorp-tion at different temperatures in a bulk zeolite. Then, using our simple model of zeolite layers, we discuss the role of the external surface on the overall adsorption phenomenon. In Section 4, using our experimental and simulated adsorption data, the adsorption potential theory is invoked to provide a formalism which allows predicting adsorption in the core of zeolites and at their external surface. In the final section, we summarize our findings and provide suggestions for future work.

COMPUTATIONAL AND EXPERIMENTAL METH-ODS

Molecular simulation

Molecular models. The silicalite-1 zeolite structure used in our study was taken from the International Zeolite Association database (IZA). [START_REF]Database of zeolite structure -MFI[END_REF] Silicalite-1, which possesses a crystalline network known as MFI (Mobil Five), is a pure silica structure -i.e. only composed of silicon and oxygen atoms. The experimental unit cell parameters for orthorhombic silicalite-1 structure are a = 20.09 Å, b = 19.738 Å and c = 13.142 Å [Figure 1(a)]. The MFI framework structure is composed of pentasil units which are linked to form pentasil chains connected via oxygen bridges to create corrugated sheets with 10-ring holes. Oxygen bridges between these sheets result in a 3-dimensional system with straight 10-ring channels parallel to the corrugations (along y) and sinusoidal 10-ring channels perpendicular to the sheets (along x).

The MFI 3D framework, consisting of straight and zigzag channels, is thus a rather complex zeolite framework type. The straight channels have circular openings of 0.54-0.56 nm and sinusoidal channels have elliptical openings of 0.51-0.55 nm. [START_REF] Julbe | In Encyclopedia of Membranes[END_REF] As shown in Figure S1(a), considering different silicalite-1 crystalline structures (e.g. obtained with different X-ray data refinement) leads to negligible changes in the simulated adsorption data. Two zeolite configurations were considered in our molecular simulation approach. The first zeolite configuration consists of an infinite zeolite crystal which is obtained by building a crystal supercell used with periodic boundary conditions. In more detail, a silicalite-1 supercell, consisting of 2 × 2 × 2 cells along the x, y and z axes, was obtained by duplicating the crystallographic unit cell as shown in Figure 1(b). This supercell was built to avoid finite size effects. The second zeolite configuration possesses an external surface as it corresponds to zeolite layers which are oriented along one of the crystalline axes. While periodic boundary conditions are applied in two directions of space, the third direction which is perpendicular to the external surface is in contact with bulk reservoirs on each side of the layer. The size of the gap (i.e. gas reservoir) between the two external surfaces of the zeolite film through the periodic boundary condition was set to 40.8 Å. Such a value was chosen as it exceeds by far the Lennard-Jones interaction cutoff (13 Å) used in our simulations; this ensures that fluid molecules at two opposite external surfaces do not interact as expected in real materials. Three layers were built from the MFI unit cell crystal by considering different crystallographic orientations for the external surface a, b, c (i.e. the crystallographic axes which are aligned with the x, y and z axes of the Cartesian reference frame). In practice, as illustrated in Figure 1(c), the layer thickness corresponds to 10 unit cells while the layer width corresponds to 2 unit cells in the two other directions. These layers were cut at well-defined positions along the normal direction to the zeolite surface to ensure that the two opposite external surfaces are equivalent (they possess the same silicon density and motif). Finally, whenever the surface has a terminating Si atom, an oxygen atom was added to build surfaces that are oxygen terminated. The use of an oxygen terminated external surface corresponds to a simplified version of real materials which are expected to display surfaces terminated with silanols (OH groups). Surface termination is expected to affect the behavior of many adsorbates at the zeolite external surface -especially polar molecules or molecules with complex chemical structure. On the other hand, with simple molecules such as methane considered here, the details of surface termination is assumed to be of limited impact (in particular, considering the very small polarizability of hydrogen atoms, the absence of H atoms at the external surface is believed to be a reasonable approximation).

To account for zeolite flexibility while maintaining the same overall crystalline structure, two different models were considered. In the rest of this article, "surface oxygen" denotes an oxygen atom which is only bonded to a silicon atom while "bulk oxygen" denotes an oxygen atom which is bonded to two silicon atoms. For the first surface flexibility model (s1), "surface oxygens" were maintained to their original position so that the surface is frozen while all other atoms are allowed to move. For the second surface flexibility model (s2), the "surface oxygens" were allowed to move during the simulation but their motion was restricted to the plane perpendicular to the layer axis. Grand Canonical Monte Carlo.

The methane adsorption isotherms at different temperatures in silicalite-1 were calculated using Grand Canonical Monte Carlo (GCMC) simulations. All Monte Carlo simulations were performed using the simulation package LAMMPS (version openmpi 1.8.1). [START_REF] Thompson | LAMMPS -a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales[END_REF] For both the zeolite crystal and the zeolite layers, the zeolite system is set in contact with a fictive bulk reservoir of methane imposing its temperature T and chemical potential µ (the latter can be directly converted into a pressure P through the equation of state as there is a unique relation between P and µ at constant T ). [START_REF] Frenkel | Understanding molecular simulation : from algorithms to applications[END_REF] At equilibrium, the methane adsorbed amount n a is calculated from the number of molecules N in the system that fluctuates around its equilibrium value N . In practice, to reach equilibrium, methane molecules are translated, inserted or removed randomly (Monte Carlo moves) with an acceptance probability defined using the Boltzmann factor in the Grand Canonical ensemble. The number of Monte Carlo moves needed to reach equilibrium depends on both T and P . In these calculations, the fluid/fluid and fluid/zeolite intermolecular potentials correspond to Lennard-Jones potentials with a cut-off r c ∼ 13 Å.

The cut-off value corresponds approximately to 3σ where σ is the Lennard-Jones parameter for methane. The United Atom model was used for CH 4 in which the methane molecule is treated as a single Lennard-Jones sphere. [START_REF] Ungerer | Optimization of the anisotropic united atoms intermolecular potential for n -alkanes[END_REF] The Lennard-Jones parameters used for the CH 4 -CH 4 interaction are σ = 3.73 Å, ε/k b = 147.9 K while those for the CH 4 -O interaction are σ = 3.214 Å, ε/k b = 133.2 K. These parameters were taken of Kar and Chakravarty. [START_REF] Kar | Diffusional Anisotropy of Simple Sorbates in Silicalite[END_REF] By treating methane as a single Lennard-Jones sphere, no intramolecular interactions have to be taken into account. This is known to be a good approximation for methane provided the temperature is high enough (like in the present work in which the minimum temperature considered is 200 K). As for the fluid/zeolite contribution, following previous works, [START_REF] Smit | Simulating the Adsorption Isotherms of Methane, Ethane, and Propane in the Zeolite Silicalite[END_REF][START_REF] Vlugt | Influence of Framework Flexibility on the Adsorption Properties of Hydrocarbons in the Zeolite Silicalite[END_REF] no Lennard-Jones potential between CH 4 and Si atoms was considered as dispersion interactions with Si can be neglected. Indeed, Si atoms have a small polarizability while O atoms have a large polarizability. We use the model by Vlugt and Schenk [START_REF] Vlugt | Influence of Framework Flexibility on the Adsorption Properties of Hydrocarbons in the Zeolite Silicalite[END_REF] to take into account the zeolite flexibility. With this model, simple harmonic potentials are used to describe the bonds between Si-O and O-O atoms connected to the same Si atom. The effect of flexibility on adsorption isotherms is shown Figure S1 (b). Adsorption isotherm under flexible conditions behaves in a similar way to rigid framework zeolite up to high pressure (∼ 10 4 bars).

Zeolites and Experiments

Sample Characterization. The experimental zeolite samples considered in our study were synthesized at the Institut Européen des Membranes in Montpellier, France according to the procedure described elsewhere. [START_REF] Motuzas | Rapid synthesis of silicalite-1 seeds by microwave assisted hydrothermal treatment[END_REF] For each sample, powder X-ray diffraction (XRD) and N 2 adsorption isotherm at 77 K were determined. There are four powder samples with different crystallite sizes (diameter) D g : 1 µm, 350 nm, 200 nm and 80 nm. More information can be found in the supporting information. The scanning electron microscope (SEM) images of the silicalite-1 powder samples are shown in Figure S2, accompanied by their XRD patterns and Periodic boundary conditions (PBC) are applied in each direction to mimic an infinitely large system and, hence, avoid finite size effects. (c) The MFI unit cell was duplicated 10 times along the chosen direction (the layer axis, here z) and 2 times along the other directions (perpendicular to the layer axis, here x and y) to form a zeolite layer with an external surface. Methane adsorption is simulated by setting the system in contact with a fictive bulk reservoir which imposes its temperature T and chemical potential µ.

nitrogen sorption isotherms in Figures S3 andS4, respectively. XRD analysis corroborates the presence of silicalite-1 diffraction lines at 2θ = 7.9°, 8.7°, 23.1°, and 24.4°, corresponding to the reflections of ( 101), (020), (501), and (303) crystallographic planes. Nitrogen sorption experiments reveal comparable microporous characteristics for the four prepared zeolite samples (typical for MFI zeolites), reaching all a specific surface area of S BET ∼ 400 m 2 /g. This measured specific surface area is mainly attributed to the internal zeolite network. The hysteresis loop observed at high p/p 0 values for samples A, B and C is attributed to the contribution of adsorption in inter-crystal spaces (i.e. mesopores) generated by crystal packing (adsorption outside zeolite crystals). The smallest mesopores were obtained with the smallest grains (sam-ple A), although no mesopores were detected for the sample with the largest crystals (sample D).

Methane adsorption.

Experimentally, methane adsorption in the zeolite samples was measured at the Institute Laue Langevin (ILL) in Grenoble, France. These experiments were carried out using a Hiden Isochema apparatus which relies on a volumetric method to measure adsorption at equilibrium. In short, the measurement principle is as follows: [START_REF] Rouquerol | Adsorption by Powders and Porous Solids[END_REF] methane is injected in a dosing cell prior to connection to the sample holder containing the zeolite matrix. The number of adsorbed molecules is determined from the difference between the number of injected molecules and the number of molecules left in the gas phase once equilibrium is reached (the latter being assessed from the final pressure and the cell volume corrected for the volume occupied by the zeolite sample). These experiments last between 10 minutes and few hours depending on the injected amount of methane and imposed temperature. Considering the experimental setup employed in this study, the adsorption data were recorded at low temperatures for which the operating conditions allowed us to obtain reliable data. As will be discussed later, using the adsorption potential theory, these data were found to be consistent with data obtained by other groups at larger temperatures.

ADSORPTION IN ZEO-LITE CRYSTALS AND THIN LAYERS

Bulk zeolite

Figure 2 shows the methane adsorption isotherms obtained at T = 300 K (red), 275 K (orange), 250 K (green), 225 K (blue), 200 K (black). The lines with crosses correspond to the simulation data while the symbols denote the experimental data. As expected for fluids confined in the very narrow pores of zeolites, adsorption is a reversible and continuous process which does not display any sharp increase in the adsorbed amount or hysteresis loop. [START_REF] Coasne | Molecular Simulation of Adsorption and Transport in Hierarchical Porous Materials[END_REF][START_REF] Coasne | Temperature Effect on Adsorption/Desorption Isotherms for a Simple Fluid Confined within Various Nanopores[END_REF] While the methane adsorbed amount n a increases upon increasing the pressure P at a given T , it increases upon decreasing T at a given P . For each T , the different regions of the adsorption isotherm correspond to different phenomena. Typically, at T = 300 K, the data in the pressure range between ∼ 0.4 and 11 bar corresponds to pore filling while the data in the pressure range from 11 to 10 3 bar corresponds to pore saturation. In the latter range, the zeolite pore volume is already filled so that the slope in the adsorption isotherm arises from the compressibility of confined methane which allows the insertion of additional methane molecules upon further increasing the pressure. Such a scenario (pore filling followed by saturation) was verified by careful inspection of the molecular configurations generated along the Grand Canonical Monte Carlo simulations at different chemical potentials. The effective methane compressibility inside the zeolite once pores are saturated was estimated from the slope of the adsorption isotherm in this pressure range: χ T = 1/ρ × ∂ρ/∂P . For example, at 300 K and P = 1200 bar, the compressibility from our simulating data is 3.52 Pa -1 for confined methane and 4.68 Pa -1 for bulk methane. The fact that the compressibility of confined methane is lower than that of bulk methane under similar thermodynamic conditions is in agreement with previous studies on nanoconfined fluids. [START_REF] Coasne | Effect of Pressure on the Freezing of Pure Fluids and Mixtures Confined in Nanopores[END_REF] In addition, Figure 2 shows the experimental adsorption isotherms for methane obtained at different temperatures for the four powder samples (open symbols). We also show experimental data obtained at 300 K from other groups taken from Ref. [START_REF] Smit | Simulating the Adsorption Isotherms of Methane, Ethane, and Propane in the Zeolite Silicalite[END_REF] (Abdul-Rehman et al., [START_REF] Abdul-Rehman | ternary, binary, and pure component sorption on zeolites. 1. Light alkanes on Linde S-115 silicalite at moderate to high pressures[END_REF] Ding et al., [START_REF] Ding | Adsorption equilibria of methane, nitrogen, and carbon dioxide on ZSM-5[END_REF] Golden and Sircar, [START_REF] Golden | Gas Adsorption on Silicalite[END_REF] Richard and Rees [START_REF] Rees | Sorption of N 2 , CH 4 and CO 2 in Silicalite-1[END_REF] ). Experimentally, the trend observed in our molecular simulation is also seen as adsorption increases upon decreasing the temperature (owing to the reduced thermal energy which promotes adsorption). The grain size D g in the zeolite powders does not strongly impact adsorption as can be observed from the data obtained at T = 250 K. This result was expected as the external surface area in the zeolite samples is small compared to the inter-nal surface in the zeolite porous structure. The simulation data at 300 K are found to be in reasonable agreement with the experimental data However, Abdul-Rehman et al. [START_REF] Abdul-Rehman | ternary, binary, and pure component sorption on zeolites. 1. Light alkanes on Linde S-115 silicalite at moderate to high pressures[END_REF] found larger adsorbed amounts compared to our simulation data but also with respect to other experimental data. Such discrepancies may stem from different material parameters including the type of silicalite-1 used by these authors; in particular, as explained in ref., [START_REF] Smit | Simulating the Adsorption Isotherms of Methane, Ethane, and Propane in the Zeolite Silicalite[END_REF] Linde S-115 sample corresponds to silicalite-1 mixed with a clay binder to form pellets so that its porosity and, hence, adsorption properties may differ from other (more conventional) samples. (circles), Golden and Sircar (squares), Richard and Rees (triangles). All these experimental data were reproduced from Smit et al. [START_REF] Smit | Simulating the Adsorption Isotherms of Methane, Ethane, and Propane in the Zeolite Silicalite[END_REF] The adsorbed amount n a is expressed in number of methane molecules per zeolite unit cell. The arrow indicates data obtained for increasing temperatures. The insert shows our experimental data obtained at 250 K for powder samples with different grain sizes: D g = 1 µm (squares), 350 nm (diamonds), 200 nm (circles), 80 nm (triangles).

Zeolite layers

Using molecular simulation, the adsorption isotherms were computed at T = 300 K, 275 K, 250 K, 225 K, and 200 K for the three layers aligned along the a, b and c directions, respectively. As explained in the Methods section, for each layer, the adsorption isotherms were computed for two surface flexibility models (s1 and s2). As will be discussed later, the adsorption isotherms for these two surface models were found to lead to similar adsorption isotherms so that only data for the model s1 will be discussed extensively. Moreover, for the sake of clarity, we will present adsorption isotherms for the layer a (data for the layers b and c at T = 300 K will be discussed later). To determine the impact of the zeolite external surface, the methane density profile was first computed at different methane pressures for each zeolite layer. An example obtained for the layer aligned along the c axis at T = 300 K and P = 4.6 bar is provided in Figure 3(a). The gray line corresponds to the local methane density in the system at a given position z -where z can lie in the bulk phase, at the external surface or in the core of the zeolite layer. These three regions are delimited by the two vertical black dashed lines. The first region, which is located on the left Figure 3(a), corresponds to the methane density in the external gaseous phase. The simulated density in this region was found to be close to the experimental density available from NIST [START_REF][END_REF] for bulk methane under the same thermodynamic conditions [see the blue horizontal dashed line in Figure 3(a)]. The second region, which is located between the two black vertical dashed lines in Figure 3(a), corresponds to the methane density at the external surface of the zeolite layer.Finally, the third region corresponds to methane adsorbed in the core of the zeolite layer as illustrated in Figure 3(a). The adsorbed density peaks in this region show marked variations corresponding to the varying methane concentration in the zeolite channels/cages. The average adsorbed methane density computed in this region -which does not include any contribution from the gaseous phase or methane adsorbed at external surface -was compared with the average adsorbed methane density estimated from the calculations for the infinite zeolite crystal [indicated by the horizontal red dashed line in Figure 3 (a)]. The two values are found to be close to each other; while the mean adsorbed methane density for the zeolite crystal is ρ inf = 44.14 kg.m -3 , the average methane density in the core of the zeolite layer is ρ c = 45.65 kg.m -3 . An important remark is in order here. Different choices can be made to define the boundaries used to estimate the methane density in the zeolite core, at the external surface, and in the bulk phase (leading to different excess quantities once volume boundaries and interfaces have been selected). However, to ensure consistency between volume and surface amounts, we made sure that our boundary definition in Fig 3 leads to bulk densities and adsorbed amounts in the zeolite core that match those obtained for bulk methane and the zeolite crystal without external surface (further discussion will be provided later when excess data for adsorption at the external surface are commented). In the same context, we highlight that Barrer and Robbins 53 developed a general approach that provides an equation of state for adsorbed films. In our case, while such a strategy could be used to define adsorption at the external surface, there is additional complexity as the surface on which adsorption occurs is nanoporous (zeolite layer beneath the external surface). Therefore, this raises the question of the boundary used to define surface adsorption versus adsorption in the core of the zeolite layer. Figure 3(b) shows the methane adsorbed amount as a function of pressure for the core of the zeolite layer a at different temperatures ranging from T = 200 K to 300 K. These data, which were extracted from the simulation data for the zeolite layer, correspond to the adsorbed amount in the innermost region of the zeolite layer (to avoid any surface effects from the zeolite external surface). In more detail, such data in molecules per unit cell correspond to the number of molecules adsorbed in a 2×2×2 supercell located at the center of the layer. As can be seen in Figure 3(b), these adsorption isotherms (open circles) match those obtained for the infinite zeolite crystal (full lines which correspond to data from Figure 2). Figure 3(c) shows the density of the bulk methane phase as assessed from the number of molecules in the region outside the zeolite layer. As expected, for all temperatures, methane behaves as an ideal gas (i.e. ρ e = P/RT ) provided the pressure does not exceed ∼ 100 bar. The density of bulk methane at P = 10 3 b and T = 300 K in our simulation (ρ e = 2.11 × 10 -2 mol.cm -3 ) is in good agreement with the experimental value from NIST (ρ N IST = 2.13 × 10 -2 mol.cm -3 ).

Figure 4 shows the adsorption isotherms obtained at different temperatures for methane at the external surface of the zeolite layer oriented with a normal vector along the crystallographic axis a. These data correspond to the absolute adsorbed amount obtained in the region between the two vertical dashed lines in Figure 3(a). As expected, the methane adsorbed amount ρ s in mol per unit of surface area increases upon increasing pressure but decreases upon increasing the temperature. Figure 4 also shows the data obtained for the different surface flexibility models (s1 or s2) and layer orientations (a, b or c) at 300 K. Considering the different layer orientations, in the lower pressure range [< 200 bar], the methane adsorbed amounts for a, b and c are found to be very similar. On the other hand, above 200 bar, at a given pressure, the methane adsorbed amount can be ranked as c > a > b. This result can be explained as follows. While the external surface in layers a and b possesses the same surface oxygen density, surface oxygens in layer b are distributed heterogeneously because of the large porosity corresponding to the MFI straight channels along the y-axis; layer a is thus believed to display a larger number of adsorption sites which, in turn, leads to larger adsorbed amounts. On the other hand, the layer c presents a larger oxygen surface density than the layers a and b, therefore leading to larger methane adsorbed amounts for c. As far as the effect of surface flexibility is concerned, the methane adsorbed amounts for the s1 and s2 models are found to be very close to each other over the entire pressure range. The insert in Figure 4 shows the contributions of the external phase, the adsorbed amount at the external surface and in the zeolite core to the adsorption isotherm for the zeolite layer a at T = 300 K. In more detail, for the different pressures considered in our study, these data show the ratios r = N e /N , N s /N , and N c /N -expressed in % -where N e , N s and N c are the number of methane molecules in the external phase, at the zeolite external surface and in the zeolite core while N is the total number of methane molecules. At low pressures (< 10 bar), due to strong confinement in the zeolite pores, the fraction of methane molecules in the zeolite core is the main contribution to the adsorption isotherm. For pressures larger than 10 bar, the zeolite pores are saturated with methane while the number of molecules at the external surface and in the gas phase keeps increasing. As a result, the contribution from these two regions increase while that from the zeolite core decreases. Following the work by Do et al. [START_REF] Do | Adsorption of argon from sub-to supercritical conditions on graphitized thermal carbon black and in graphitic slit pores: A grand canonical Monte Carlo simulation study[END_REF] on methane adsorption on surfaces at temperatures T > T c , we have also calculated surface excess densities which are illustrated in Figure S6 of the supporting information. Surface excess densities were estimated by subtracting the methane density of the external phase density from the methane surface density. As expected, in the low pressure range, the excess surface density is close to the absolute density as the external phase density is small. However, in the large pressure range, the external phase density is large so that the excess density becomes very small and even negative. These results are consistent with previous work by Vlugt and coworkers [START_REF] García-Pérez | External Surface Adsorption on Silicalite-1 Zeolite Studied by Molecular Simulation[END_REF] who showed that surface excess concentrations for pure ethane and propane increase and become positive when the adsorption isotherm reaches a plateau (therefore also supporting favorable adsorption at the surface). Similarly, using free energy considerations, Inzoli et al. [START_REF] Inzoli | Surface Adsorption Isotherms and Surface Excess Densities of n -Butane in Silicalite-1[END_REF] proposed that surface adsorption becomes more favorable as pores get saturated due to the increasing contribution from steric interactions inside the zeolite cages.

ADSORPTION POTEN-TIAL THEORY

In this section, we discuss a rational thermodynamic framework that allows describing in a consistent fashion adsorption in the zeolite core and at its external surface. With this goal, we invoke the adsorption potential theory which was introduced in his seminal work by Polanyi. [START_REF] Polanyi | Verhandlungen der Deutschen Physikalischen[END_REF][START_REF] Polanyi | The Potential Theory of Adsorption[END_REF][START_REF] Neimark | Adsorption-Induced Deformation of Microporous Solids: A New Insight from a Century-Old Theory[END_REF] Despite its theoretical nature, this model is mostly used in an empirical way with very practical applications through characterization equations such as the Dubinin-Radushkevich relation for instance. [START_REF] Rouquerol | Adsorption by Powders and Porous Solids[END_REF] In contrast, in the following section, we derive the adsorption potential theory by starting from the Hill's adsorption scheme which corresponds to a perturbation treatment of molecular adsorption. [START_REF] Hill | Statistical Thermodynamics of the Transition Region between Two Phases. II. One Component System with a Plane Interface[END_REF][START_REF] Hill | Advances in Catalysis[END_REF][START_REF] Hill | Extension of Fowler's Treatment of Surface Tension to Physical Adsorption[END_REF] Polanyi's adsorption potential theory can be derived by starting from a simple approximation which provides a physical model for gas adsorption on a solid surface. Let us consider a fluid phase adsorbed (a) at a temperature T on a flat surface in equilibrium with a bulk fluid phase (g) so that their chemical potential are equal µ a = µ g (Figure 5). Assuming the fluid phase behaves as an ideal gas, its chemical potential µ g can be written as µ g -µ 0 = RT ln[P/P 0 ] where µ 0 is the chemical potential at the saturating vapor pressure P 0 (T ) (we note that the same expression can be used for a nonideal phase but the pressure has to be replaced by the fugacity). Within this approximation, the chemical potential µ a of the adsorbed phase is described using a simple perturbation treatment by assuming that it corresponds to the chemical potential µ 0 of the bulk liquid at sat-ρ s (x10 -9 mol.cm (color online) (a) Adsorption isotherms for methane at the external surface of the zeolite layer a (these data are obtained for the surface model s1 in which surface atoms are frozen). Different colors correspond to the following temperatures: 300 K (red), 275 K (orange), 250 K (green), 225 K (blue), 200 K (black). The methane adsorbed amounts are expressed as surface densities ρ s in mol.cm -2 . The insert shows the following ratios at 300 K where N is the total number of methane in the system: methane at the surface uration corrected for the interaction U (z) with the solid surface:

µ a -µ 0 = U (z) (1)
This is the essence of this simple approximation which assumes that a molecule gets adsorbed at a position z from the surface when the interaction with the solid surface counterbalances the difference between the chemical potential of the fluid phase and that of the dense, cohesive liquid. In other words, at the chemical potential µ, the bulk liquid phase is not stable but the solid/fluid interaction stabilizes the adsorbed liquid. At very small z, the surface interaction potential U (z) is positive due to the repulsive interactions with the solid surface (hard core-like, repulsive part of the solid/fluid interaction). In contrast, at larger z, U (z) becomes negative with a minimum observed at a distance corresponding approximately to the position of the first adsorbed layer. At even larger z, U (z) goes to zero as interactions with the solid surface become vanishingly small. Note that Eq. ( 1) does not assume that the pressure or chemical potential depends on the position z; this would be unphysical from a thermodynamic viewpoint as chemical potential and pressure must be homogeneous because of chemical and mechanical equilibrium, respectively. In practice, the statistical film thickness t and surface potential U (z) can be related to the position z of the outermost adsorbed molecules as follows: t can be converted from the adsorbed volume N v as t = N v/S where N is the number of adsorbed molecules, v the molar volume assumed to correspond to that of the bulk liquid and S the surface area of the solid phase (Figure 5). For a given molecular configuration of adsorbed molecules, there is a unique function that links t and z; t = z + σ/2 where σ is the kinetic diameter of the adsorbate molecule. In other words, the surface interaction potential can be written as U (z) or equivalently U (t) -therefore offering a mean to relate the film thickness and surface potential using the perturbation treatment described in Eq. ( 1). While the treatment above shares similarities with the well-known Frenkel-Halsey-Hill approximation, it differs from this seminal treatment for the following reasons. First, in the Frenkel-Halsey-Hill approximation, the adsorption film is assumed to be a homogeneous phase whose properties are close to those of the bulk liquid. In contrast, considering that we treat the case of a highly heterogeneous adsorbed phase (above the bulk critical point and in a pressure range where the adsorbed phase is heterogeneous), our system departs from the conditions considered in the Frenkel-Halsey-Hill approximation. Second, within this approximation, one considers that the free energy difference be-tween the adsorbed film and the bulk liquid simplifies to the potential energy difference corresponding to the adsorbate/solid interactions (in other words, entropy terms and fluid/fluid contributions cancel out as the bulk and adsorbed phase are assumed identical). Here, despite the important differences between the adsorbed and bulk phase, we assume that the free energy difference can still be approximated as the interaction energy of an adsorbed methane molecule with the MFI zeolite channel or external surface. For these different reasons, our treatment should be considered only related to the Frenkel-Halsey-Hill approximation. In fact, in this respect, considering the treatment made and system conditions under study, this approach is closer to that introduced by Barrer and Robbins 53 and Hill [START_REF] Hill | Statistical Thermodynamics of the Transition Region between Two Phases. II. One Component System with a Plane Interface[END_REF][START_REF] Hill | Advances in Catalysis[END_REF][START_REF] Hill | Extension of Fowler's Treatment of Surface Tension to Physical Adsorption[END_REF] who proposed extension to non homogeneous adsorbed films (with local densities dependent on the distance from the surface).

Polanyi's adsorption potential theory states that adsorption equilibrium is governed by a simple bijective function f such as:

N v = f [µ -µ 0 ] ( 2
)
where N is the adsorbed amount in moles and v is the molar volume of the adsorbed phase. This function, which is unique for a given fluid/solid couple, has no explicit temperature dependence but we note that v, µ and µ 0 are temperature dependent. This equation is most often used as an empirical relation but its physical meaning can be understood as follows. Polanyi's model simply states that adsorption in sites having an adsorption energy E occurs as the chemical potential difference µ -µ 0 is such that µ -µ 0 = E. With this formalism, considering that intermolecular interactions are temperature independent since they derive from a Hamiltonian, it can be assumed that the latter condition is temperature independent and so is the function f . The characteristic function f does not have to be temperature independent. However, to allow practical use (i.e. to determine adsorbed amounts at a given temperature from available adsorption data at a reference temperature), it is often assumed that f is tem-perature independent -at least on the temperature range covering the reference temperature and that at which the data will be extrapolated.

In particular, as will be shown below, Polanyi's model is found to be relatively accurate by offering at least a semi-quantitative description of adsorption in zeolite materials. At equilibrium, the gas and adsorbed phase possess the same chemical potential, i.e. µ g = µ a = µ. Using the FHH approximation, we can write that molecules get adsorbed at a position z when the chemical potential of the adsorbed molecules is such that µ-µ 0 = U (z). Considering that U (z) is a bijective function of z in the range of z values where adsorption occurs, we can write z = U -1 (µ -µ 0 ) (where U -1 is the inverse function that transforms the chemical potential into the z position of the adsorbed molecule). By noting that t = N v/S ∼ z, this equation is equivalent to Polanyi's equation with this analogy leading to f ∼ U -1 . In practice, Polanyi's equation relies on the knowledge of an adsorption isotherm at a given temperature T to predict any adsorption isotherm at a different temperature. To illustrate Polanyi's model application, let us assume that the adsorbed amount at a state P , T is known. The chemical potential is also known through the appropriate equation of state at such pressure/temperature conditions (here, for the sake of simplicity, we assume that the ideal gas equation applies as it is relevant to the conditions considered in our study). Using the chemical potential/pressure relation for an ideal gas [i.e. µ -µ 0 ∼ k B T ln P/P 0 ], we can predict the pressure P at which the same chemical potential condition will be reached when the temperature is T : k B T ln P /P 0 = k B T ln P/P 0 where P 0 and P 0 are the saturating vapor pressure at T and T , respectively. Using Polanyi's model, i.e. N v ∼ constant, we obtain N v = N v when the same chemical potential conditions are met. In practice, with such a model, to predict data at a given temperature, one needs to shift the adsorbed amount along the y axis by considering that the molar volume v of the adsorbed film is identical to that of the bulk phase at the same temperature. Upon considering a different temperature, the data must also be shifted along the x axis to account for the shift in chemical potential (or equivalently in pressure) since Polanyi's model assumes that data can be described consistently through Eq.( 2) where µ 0 describes the saturating chemical potential. spheres) is µ g . The adsorbed phase (blue spheres) has a statistical film thickness t which can be converted from the adsorbed volume N v as t = N v/S where N is the number of adsorbed molecules, v the molar volume and S the surface area of the solid phase. For the adsorbed phase, the chemical potential µ a is described using a perturbation term by assuming that it corresponds to the chemical potential µ 0 of the bulk liquid at saturation corrected for the interaction U (z) with the solid surface: µ a -µ 0 = U (z). At very small z, the surface interaction potential U (z) is positive due to repulsive interactions with the solid surface. At larger z, U (z) becomes negative with a minimum located at a distance corresponding approximately to the position of the first adsorbed layer. At larger z, U (z) goes to zero as interactions with the solid surface become vanishingly small. Figure 6 compares the predictions from Polanyi's model with our simulated methane adsorption isotherms obtained at different temperatures. Both the theoretical and simulated data for adsorption at the external surface (a) and in the zeolite pores (b) are shown. For both adsorption types, the adsorption isotherm obtained at T = 250 K (green data) was used as reference data to predict the adsorption isotherms obtained at other temperatures T . Within the frame of the adsorption potential theory, we note that the accuracy of the inferred adsorption data from available data depends on the chosen reference temperature. In the present work, we chose to use data at 250 K as this temperature is intermediate between all temperatures under study (+/-50 K). While we think that this is the best option available, we checked that the absolute error bar remains of the same order -typically ∼ 1-2 molecules per unit cell -when choosing an extremum temperature (i.e. 200 K or 300 K). The reference data were smoothed using a nearest neighbor moving average (full green lines). In more details, data were first interpolated using splines and then 15 points were considered for the nearest neighbor smoothing -it was checked that, within numerical errors, the detailed treatment does not affect our conclusions below. Considering that all data shown in Figure 6 are obtained at temperatures above the critical point of methane (T c ∼ 190 K), Polanyi's model must be extended under such conditions. Indeed, for such temperatures, there is no saturating vapor pressure P 0 and, hence, no chemical potential at saturation µ 0 . Similarly, for a fluid phase under such supercritical conditions, one cannot define the liquid molar volume v as involved in the adsorption potential theory. To extend the validity of Polanyi's approach, we invoke the two following hypotheses to define pseudo-quantities equivalent to v(T ) and P 0 (T ) for T > T c . (1) P 0 (T ) is obtained by extrapolating Antoine's Law at temperatures above T c . Antoine law parameters were taken from NIST. [START_REF][END_REF] For convenience, the corresponding data showing the saturating vapor pressure P 0 (T ) as a function of temperature T is shown in Supplementary Figure S5. (2) v is obtained by fitting our simulated adsorption isotherms at the external surface once saturation is reached (Figure 4). In more detail, for each temperature T , the high pressure range of the adsorption isotherm at the external surface was fitted using the following function: ρ s (P, T ) = ρ s,0 (T ) + b exp[-cP ] where ρ s,0 (T ) = 1/v(T ) is the surface density at saturation while b and c are fitting parame-ters. In more detail, considering that the molar volume v of methane adsorbed at saturation (or equivalently its density since ρ ∼ 1/v) is not defined above the critical temperature T c , one has to estimate its extrapolated value from the saturation value observed in the adsorption isotherms. To do so, several functions were considered to fit the simulated adsorbed amounts obtained for the external surface. In practice, the simple mathematical function above was found to accurately describe the change in the methane adsorbed amount upon increasing the pressure. While purely mathematical, we note that this function does predict the right evolution of the methane adsorbed amount upon increasing the pressure, so that it can be used to predict the extrapolated molar volume at saturation as 1/ρ s,0 (T ) (we note that the parameters b and c were determined using the fit of the isotherm at 250 K and then were kept constant for the other temperatures as they relate to the same fluid/solid couple). The insert in Figure 6(a) shows the evolution of the ratio v(T )/v(T ) where, as mentioned above, T = 250 K is the reference temperature considered when applying Polanyi's model. Before discussing the predictions from the adsorption potential theory, we mention that the extrapolation of the concepts of the saturating vapor pressure and molar volume of the adsorbed phase is far from trivial. Such approximations assume that the adsorbed phase at temperatures above the critical point is reminiscent of the liquid phase under subcritical conditions. This is justified by the fact that, even under critical conditions, the dense adsorbed phase displays a density that is close to the random close packing density inherent to the liquid phase. Similarly, while the concept of saturating vapor pressure does not apply to fluids above their critical point, its extension at temperatures beyond T c characterizes the pressure (or chemical potential) at which surface adsorption increases sharply with pressure.

U(z) z U(∞) ~ 0 μ g z eq. t = Nv/S μ a
As shown in Figure 6, for both adsorption in the zeolite pores and at its external surface, the predictions from the adsorption potential theory are in good agreement with the simulated data. To further test the validity of the adsorp-tion potential theory, it was also applied to our experimental adsorption isotherms obtained for methane adsorption at different temperatures (two powder samples with grain sizes D g = 350 nm and D g = 1 µm were considered but as mentioned above the grain size was found to make almost no difference). Like with the simulated data, the experimental adsorption isotherm at 250 K for the powder with grain size D g = 350 nm was used as reference data to predict adsorption at a lower temperature for the other powder. The extrapolated values for P 0 and molar volumes at saturation v(T ) were taken identical to those obtained from the simulated data. As shown in Figure 6(c), a good agreement is observed between the predictions from the adsorption potential theory and the experimental data. This finding further confirms the applicability of Polanyi's theory to predict adsorption as a function of temperature from available reference data -even at conditions above the fluid critical temperature.

To assess the quantitative validity of perturbation treatment, we now consider explicitly the fluid/solid interaction potential U (z) to describe methane adsorption in the cylindrical channels aligned along the b-axis in silicalite-1. To do so, considering the cylindrical geometry of these channels and the Lennard-Jones potential employed in our molecular simulation to describe the methane/zeolite interactions, we use the interaction potential by Peterson et al. [START_REF] Peterson | Fluid behaviour in narrow pores[END_REF] which corresponds to an integrated version of the 12-6 Lennard-Jones interaction energy:

V (s 1 ) = 16 1 ρ z π × 7σ 12 1 512 K 9 (s 1 ) - σ 6 1 16 K 3 (s 1 ) (3)
where K 9 (s 1 ) and K 3 (s 1 ) are two mathematical functions that only depend on the distance of a methane molecule from the channel center s 1 . ρ z is the atom density of the cylindrical pore considered as homogeneous in this continuum-level equation while 1 and σ 1 are the fluid/solid Lennard-Jones parameters (for detailed calculations, the reader is referred to Ref. [START_REF] Peterson | Fluid behaviour in narrow pores[END_REF] ). Here, we use the same Lennard-Jones parameters as those used in the GCMC simu- lations for the CH 4 /zeolite interactions; However, while the solid/fluid interaction for an adsorbed molecule in GCMC simulations is given by a discrete sum over each O atom in the zeolite, the interaction potential by Peterson et al. assumes that the zeolite can be described as a continuum medium. By inserting such an integrated potential into the FHH approximation, we can predict the adsorption isotherm at 300 K for methane restricted to the straight channels in silicalite-1. For a straight cylindrical channel having a radius R and a length L, the radial position of the outermost adsorbed molecules z is related to the adsorbed amount

N (z) as N (z) = πρL[R 2 -(R -z) 2 ]
where ρ is the density of the adsorbed phase. The insert in Figure 6 (the FHH approximation model assumes the adsorbed phase to be incompressible so that adsorbed amount remains constant beyond pore filling). As can be seen from these data, the perturbation treatment used here in combination with the interaction potential by Peterson et al is not quantitatively accurate. We emphasize that this is not due to an intrinsic failure of the Frenkel-Halsey-Hill approximation but to the following points. (1) As already discussed, while the Frenkel-Halsey-Hill hypothesis relies on the free energy rather than the interaction energy, here we only consider the interaction energy of a single methane molecule so that it neglects the role of entropy contributions as well as of fluid/fluid contributions to the adsorbed phase free energy. (2) In our estimate of the interaction energy, we use Peterson's potential which relies on a continuum description of the Lennard-Jones interaction potential while the interaction energy landscape -even for a single molecule -can be more complex.

CONCLUSION

Using a combined experimental and molecular simulation study, we provide microscopic insights into the adsorption of simple gases in silicalite-1 zeolite crystals and thin layers. By considering both adsorption in the zeolite pores and at the external surface of zeolite layers with different crystalline orientations, we investigate the impact of the zeolite external surface on the overall thermodynamic behavior of the gas/solid couple. Due to strong confinement effects in the zeolite pores, we observe that methane first adsorbs inside the zeolite porous structure and then adsorbs at the external surface once the core porous volume gets filled.

In an attempt to provide a simple thermodynamic framework for adsorption in such zeolite samples, we employ the formalism of the adsorption potential theory to rationalize adsorption data at different temperatures. Considering that methane adsorption in our study occurs at temperatures above its critical point, we extend the adsorption potential theory for such supercritical conditions. Both our experimental and simulated data suggest that the adsorption potential theory is a robust thermodynamic modeling approach to capture adsorption at different temperatures in zeolite materials. While beyond the scope of the present paper, we believe that this general framework can be applied to any zeolite pore geometries (channels and/or cages). In particular, considering that the MFI zeolite displays a rather complex network structure compared to many other zeolites, we feel that such an extension should apply successfully. Moreover, our simulated data for the adsorption in the zeolite core porous volume and at its external surface suggest that the same parameters can be used to rationalize adsorbed amounts in these different regions. These findings provide a step toward the development of a consistent thermodynamic formalism to predict surface/volume adsorption in nanoporous materials with non-negligible external surface areas.

With the design of more and more advanced zeolite materials (e.g. hierarchical zeolites, zeolite films, nanozeolites, dendritic zeolites, etc.), we believe that the development of novel characterization techniques is needed. In this context, in addition to conventional techniques such as microscopy and tomography, adsorption-based approaches constitute methods of choice to probe the nanoporosity and large specific surface area in zeolite materials but also beyond with other nanoporous solids (Metal Organic Framework, carbonaceous materials, etc.). This method is also complementary of other techniques for complex porous samples with multiscale and/or disordered porosity such as NMR cryoporometry [START_REF] Mitchell | Nuclear magnetic resonance cryoporometry[END_REF] and hysteresis scanning [START_REF] Coasne | Domain theory for capillary condensation hysteresis[END_REF][START_REF] Schneider | Phase transitions in disordered mesoporous solids[END_REF] (which probe pore network connectivity and morphology). Such characterization efforts are needed as more and more applications in adsorption and catalysis science rely on the design 'on demand' of specific nanoporous samples with well-controlled porous volume and surface area. To assist such synthesis and design efforts, in parallel to the already available adsorption-based techniques cited above, a robust thermodynamic modeling framework must be developed to establish accurate structure-property relationships for gas adsorption. Beyond adsorption effects, such predictive thermodynamic tools is also needed to rationalize the impact of adsorption on diffusion and transport inside and at the external surface of nanoporous solids.

Graphical TOC Entry

Figure 1 :

 1 Figure 1: (color online) Silicalite-1 zeolite crystal. The Si and O atoms are represented as yellow and red spheres, respectively. a, b, c denote the crystallographic axes which are aligned with the x, y and z axes of the Cartesian reference frame. A single unit cell is represented in (a) while a 2×2×2 crystal is shown in (b). Periodic boundary conditions (PBC) are applied in each direction to mimic an infinitely large system and, hence, avoid finite size effects. (c)The MFI unit cell was duplicated 10 times along the chosen direction (the layer axis, here z) and 2 times along the other directions (perpendicular to the layer axis, here x and y) to form a zeolite layer with an external surface. Methane adsorption is simulated by setting the system in contact with a fictive bulk reservoir which imposes its temperature T and chemical potential µ.

4 TFigure 2 :

 42 Figure 2: (color online) Experimental and simulated adsorption isotherms n a (P ) for methane at different temperatures in silicalite-1 zeolite. The different colors correspond to the following temperatures: 200 K (black), 225 K (blue), 250 K (green), 275 K (orange), 300 K (red). The lines with crosses correspond to the results from our molecular simulations. The symbols without lines correspond to experimental data. The open symbols correspond to our experimental data. The closed symbols correspond to experimental data at 300 K from other groups: Abdul-Rehman et al. (diamonds), Ding et al.(circles), Golden and Sircar (squares), Richard and Rees (triangles). All these experimental data were reproduced from Smit et al.[START_REF] Smit | Simulating the Adsorption Isotherms of Methane, Ethane, and Propane in the Zeolite Silicalite[END_REF] The adsorbed amount n a is expressed in number of methane molecules per zeolite unit cell. The arrow indicates data obtained for increasing temperatures. The insert shows our experimental data obtained at 250 K for powder samples with different grain sizes: D g = 1 µm (squares), 350 nm (diamonds), 200 nm (circles), 80 nm (triangles).

Figure 3 :

 3 Figure 3: (color online) (a) Methane density for the silicalite-1 layer aligned along the c crystallographic axis. For visualization purpose, only the left side of the layer is shown (the system is symmetrical with respect to the position z = 0). The cavity size is about 1/3 of the Lennard-Jones size parameter for the methane molecule. (b) and (c) show the adsorption isotherms in the core of the zeolite layer and density in the external fluid phase, respectively. Different colors correspond to different temperatures T : 300 K (red), 275 K (orange), 250 K (green), 225 K (blue), 200 K (black). In (b), the lines show the adsorbed amounts obtained for the infinite zeolite while the open circles show the adsorbed amounts obtained by counting methane molecules adsorbed in the core of the layers.

Figure 4 :

 4 Figure 4:(color online) (a) Adsorption isotherms for methane at the external surface of the zeolite layer a (these data are obtained for the surface model s1 in which surface atoms are frozen). Different colors correspond to the following temperatures: 300 K (red), 275 K (orange), 250 K (green), 225 K (blue), 200 K (black). The methane adsorbed amounts are expressed as surface densities ρ s in mol.cm -2 . The insert shows the following ratios at 300 K where N is the total number of methane in the system: methane at the surface N s /N (line with open triangles), methane inside the core N c /N (line with open circles), methane in the external phase N e /N (line with open squares). (b) Adsorption isotherms for methane at the external surface for layers a, b, c which correspond to the lines with open circles, open squares, and open triangles, respectively. The symbols are obtained for the surface model s1 in which surface atoms are frozen (s1 model). The crosses correspond to data for the layer a with the flexible surface model (s2 model). Note that the data for the s1 and s2 models are nearly identical.

  Figure 4:(color online) (a) Adsorption isotherms for methane at the external surface of the zeolite layer a (these data are obtained for the surface model s1 in which surface atoms are frozen). Different colors correspond to the following temperatures: 300 K (red), 275 K (orange), 250 K (green), 225 K (blue), 200 K (black). The methane adsorbed amounts are expressed as surface densities ρ s in mol.cm -2 . The insert shows the following ratios at 300 K where N is the total number of methane in the system: methane at the surface N s /N (line with open triangles), methane inside the core N c /N (line with open circles), methane in the external phase N e /N (line with open squares). (b) Adsorption isotherms for methane at the external surface for layers a, b, c which correspond to the lines with open circles, open squares, and open triangles, respectively. The symbols are obtained for the surface model s1 in which surface atoms are frozen (s1 model). The crosses correspond to data for the layer a with the flexible surface model (s2 model). Note that the data for the s1 and s2 models are nearly identical.

Figure 5 :

 5 Figure 5: (color online) Schematic representation of the adsorption potential theory. The chemical potential of the fluid phase (red spheres) is µ g . The adsorbed phase (blue spheres) has a statistical film thickness t which can be converted from the adsorbed volume N v as t = N v/S where N is the number of adsorbed molecules, v the molar volume and S the surface area of the solid phase. For the adsorbed phase, the chemical potential µ a is described using a perturbation term by assuming that it corresponds to the chemical potential µ 0 of the bulk liquid at saturation corrected for the interaction U (z) with the solid surface: µ a -µ 0 = U (z). At very small z, the surface interaction potential U (z) is positive due to repulsive interactions with the solid surface. At larger z, U (z) becomes negative with a minimum located at a distance corresponding approximately to the position of the first adsorbed layer. At larger z, U (z) goes to zero as interactions with the solid surface become vanishingly small.

Figure 6 :

 6 Figure 6: (color online) Polanyi's adsorption potential theory applied to simulated methane adsorption isotherms in zeolite layer a (a) and in a zeolite crystal (b) at different temperatures T . In (c), Polanyi's model is applied to our experimental data. Different colors correspond to the following temperatures: 300 K (red), 275 K (orange), 250 K (green), 225 K (blue), 200 K (black). The open symbols correspond to the simulated/experimental adsorbed amounts as obtained using Monte Carlo simulations in the Grand Canonical ensemble/adsorption volumetric technique. For the three figures, the adsorption isotherm obtained at T = 250 K (green data) was used as reference data to predict the adsorption isotherms obtained at other temperatures. The dashed lines correspond to the predictions from Polanyi's model at different temperatures T . In (c) the adsorption isotherm for the powder of D g = 350 nm at 250 K was used as the reference to predict the adsorption of the powder of D g = 1 µm at 225 K. The adsorbed amount n a is expressed in mol.cm -2 (a), number of methane molecules per zeolite unit cell (b) and mol.cm -3 (c). The insert in figure (a) corresponds to the ratio of the molar volume at T over the molar volume at the reference temperature T = 250 K. The insert in figure (b) corresponds to the adsorption data predicted using the FHH approximation in combination with the interaction potential by Peterson et al.. These data correspond to the adsorbed amount in the straight channels of silicalite-1 in the pressure range from 1 to 100 bar. The black line corresponds to the adsorption isotherm for the zeolite channels obtained by means of molecular simulation. The dashed orange line corresponds to the adsorption isotherm calculated using the attractive part of the interaction potential by Peterson et al.58

  (b) compares the simulated adsorption isotherm for the straight channels (black line with crosses) with the predictions from the perturbation treatment based on the interaction potential by Peterson et al. (orange lines). The methane density ρ used to compute the methane adsorbed amount was adjusted to match the simulated adsorbed amount at saturation. The orange dashed line corresponds to data obtained using the perturbation treatment in combination with the attractive contribution of the interaction potential by Peterson et al. The data are plotted up to 20 bars only where the pore gets filled z = R

Acknowledgement W. K. is grateful to the Fédération Française de Diffusion Neutronique 2FDN (FR2004) for funding. This work was also supported by the French Research Agency (ANR CATCALL 19-CE07-0025).

Calculations were performed using the Froggy platform of the GRICAD infrastructure (https://gricad.univ-grenoblealpes.fr), which is supported by the Rhône-Alpes region (GRANT CPER07-13 CIRA) and the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investisse-ments d'Avenir supervised by the French Research Agency.

Supporting Information Available

The Supporting Information is available free of charge at DOI: XXX. Simulated adsorption curves for different zeolite models and characterization data for the zeolite samples (PDF Link)

Author Information

Corresponding Author * E-mail: benoit.coasne@univ-grenoble-alpes.fr Notes The authors declare no competing financial interest.


Léo Théodon
email: l.theodon@emse.fr

Jérôme Laurencin
email: jerome.laurencin@cea.fr

Peter Cloetens
email: cloetens@esrf.eu

Johan Debayle
email: debayle@emse.fr

A stochastic 3D model for three-phase heterogeneous microstructures in SOFC-electrodes

Keywords: 3D microstructure model, Image analysis, Simulation, Solid oxide fuel cell, Stochastic geometry

Introduction

The modeling of three-phase microstructures using stochastic geometry techniques is a major challenge in the context of developing digital twins for solid oxide fuel cell (SOFC) electrodes. The control of the microstructure properties and its optimization is crucial in order to improve the cell performance in operation [START_REF] Kong | Ni-YSZ gradient anodes for anode-supported SOFCs[END_REF]. Moreover, the agglomeration of the nickel phase, as well as its redistribution in the electrode after several hours of use have a significant impact on the cell life time [START_REF] Hubert | Impact of Nickel Agglomeration on Solid Oxide Cell Operated in Fuel Cell and Electrolysis Modes[END_REF] [START_REF] Zekri | Nickel Depletion and Agglomeration in SOFC Anodes During Long-Term Operation[END_REF].

In this context, many stochastic models have been developed to reflect the geometrical properties of three-phase microstructures. Such models can be based on hard spheres packing [START_REF] Westhoff | Generation of virtual lithium-ion battery electrode microstructures based on spatial stochastic modeling[END_REF] [START_REF] Moussaoui | Particle-based model for functional and diffusion layers of solid oxide cells electrodes[END_REF], Gaussian random fields [START_REF] Moussaoui | Stochastic geometrical modeling of solid oxide cells electrodes validated on 3D reconstructions[END_REF] [START_REF] Prifling | Stochastic 3D microstructure modeling of anodes in lithium-ion batteries with a particular focus on local heterogeneity[END_REF], or even random graphs [START_REF] Neumann | Stochastic 3D modeling of complex three-phase microstructures in SOFC-electrodes with completely connected phases[END_REF]. These models are effective in capturing the local heterogeneities in the typical SOFC electrodes with respect to properties such as specific surface area, tortuosity or constrictivity. All these models consider an homogeneous distribution of the phases in the microstructure. Nevertheless, Fig. 1 and Fig. 2 show that it is necessary to model structures with strong global heterogeneity due to phase gradients and/or the depletion of the metallic phase (Ni) at the interface with the electrolyte. For all these reasons, a model able to take into account this global phase heterogeneity is hereby proposed. The originality of the method comes from the fact that it enables to model SOFC-electrode interfaces using stochastic geometry tools.

In the following section, the characteristics of the proposed model are defined and examples of realizations are presented. Then, the optimization process is described, followed by the validation step and the obtained results. Finally, a conclusion summarizes the approach and opens new perspectives.

3D microstructures reconstruction and characterization 2.1 Image reconstruction

In the present study, the electrodes active layers of both cells are reconstructed by synchrotron X-ray holotomography [9] [10]. This method has the advantage of being able to represent 3D volumes with a large field of view (50 µm) and to capture details of the order of a tenth of a nanometer. The microstructure reconstructions at our disposal have a volume from 25 µm 3 to 50 µm 3 with a voxel size of 25 nm. Examples of 3D reconstructions of real microstructures with volumes of 17.5 µm 3 are shown in Fig. 1 and Fig. 2. The porous phase is in black, the ceramic phase in grey and the metallic phase in white. (Right) Evolution of the phase surface area density associated with each phase along the Z-axis.

Image characterization

The microstructure of an electrode can in principle be completely characterized by a set of properties that are directly measurable on the 3D volume. The computation of these properties and their meaning has already been described in detail in the literature, as for instance in articles [START_REF] Laurencin | Characterisation of Solid Oxide Fuel Cell Ni-8YSZ substrate by synchrotron X-ray nano-tomography: from 3D reconstruction to microstructure quantification[END_REF], [START_REF] Usseglio-Viretta | Quantitative microstructure characterization of a Ni-YSZ bi-layer coupled with simulated electrode polarisation[END_REF] and [START_REF] Neumann | Stochastic 3D modeling of complex three-phase microstructures in SOFC-electrodes with completely connected phases[END_REF].

• The volume fraction ε i represents the ratio of volume occupied by a phase i to the total volume of the microstructure.

• The specific surface area S i p that represents the normalized surface area of phase i. • The interfacial specific surface area S i,j p between two phases i and j allows to completely characterize the specific surface areas S k p for any phase k by linear combination. • The density of Triple Phase Boundary length or TPBl ξ TPBl is calculated by identifying each voxel edges at the triple contacts.

Throughout the rest of this article, the porous, YSZ and Nickel phases are referred to the label 0, 1 and 2 respectively. Table 1 shows some properties measured on the microstructure reconstructions in Fig. 1 and Fig. 2. The phase surface area density s i is defined on each slice of the microstructure along the Z-axis, which corresponds to the vertical axis on Fig. 1 and Fig. 2, and is the ratio of the surface area of phase i to the total area of the slice. Hence, we have the following result

ε i = 1 l z lz 0 s i (z)dz (1)
where l z is the length of the domain W along the Z-axis.

Table 1: Properties measured on real microstructures.

Properties Real µ-structure (Fig. 1) Real µ-structure (Fig. 2)

ε 0 0.253 0.250 ε 1 0.424 0.439 ε 2 0.323 0.310 S 0,1 0.878 µm -1 0.718 µm -1 S 0,2 p 0.236 µm -1 0.144 µm -1 S 1,2 p 0.728 µm -1 0.686 µm -1 ξ TPBl 1.200 µm -2 0.470 µm -2
Finally, to reduce the noise and fluctuations when computing s i (z), i = {1, 2} on the microstructures, the parameters of a function g are fitted using the least squares method so that g is an approximation of s i . The function g is defined by

g(z) = a + (1 -a) × 1 1 + exp(c × (z -b)) . ( 2
)
The intention of this article is to propose a model for generating microstructures with characteristics identical to those measured on the real ones, while sharing the same global heterogeneity.

Stochastic model

The proposed model to simulate the phase gradients located at the interfaces of Ni-YSZ anodes is based on spatial point processes and distance maps. The process of creating a synthetic microstructure is the following.

1. A set of points X i ⊊ R 3 is associated to each of the three phases. The sets X 0 and X 1 follow a Poisson point process depending only its density λ [START_REF] Kingman | Poisson Processes[END_REF] while the third set X 2 follows a Neyman-Scott point process [START_REF] Neyman | Statistical Approach to Problems of Cosmology[END_REF]. The Neyman-Scott point process requires a point density λ p for a parent set that follows a Poisson point process. To each parent point process is attached a subset of points randomly located from a normal distribution of standard deviation σ and whose number of elements depends on a Poisson distribution of intensity λ o (Fig. 3). Therefore, a minimum of three additional parameters is needed. Only the points of the child subsets are taken into account at the end of the process. Thus, three sets of points are generated using five parameters.

2. The three distance maps of the point sets are computed using the algorithm described in [START_REF] Felzenszwalb | Distance Transforms of Sampled Functions[END_REF]. Thus, a threephase microstructure can be built by attributing to each point of the space the set with which it has a minimum distance. Nevertheless, such a microstructure would be globally homogeneous. To introduce an heterogeneity in the microstructure, two functions f 1 (x, y, z) and f 2 (x, y, z) are used in order to modify the distance maps of the sets X 1 and X 2 .

3. Finally, a two-step smoothing process described in [START_REF] Neumann | Stochastic 3D modeling of complex three-phase microstructures in SOFC-electrodes with completely connected phases[END_REF] is applied. The parameter θ controlling the intensity of the smoothing allows adjusting the density of TPBl of the microstructure.

In the context of this paper, the functions f 1 and f 2 applied to the distance maps are based on sigmoids. In particular, the functions f i , i = {1, 2} are defined as follows.

f i (z) = 1 + a i × 1 - 1 1 + exp(c i × (z -b i )) (3)
where a i , b i and c i , i = {1, 2} are parameters of shape and position.

The reason for choosing sigmoid functions is purely based on the fact that the results provided are quite satisfactory, but it would be interesting to explore other classes of functions in future work.

Fig. 4 illustrates the process of generating a microstructure with the model described above while Fig. 5 shows some examples of realizations of synthetic phase gradient microstructures. The flexibility of the model comes from its large number of parameters: five for the point processes, two functions, and one smoothing parameter. The following section details the process of fitting the input model parameters using the real reconstruction and the method validation.

Model fitting and validation 4.1 Cost function definition

To fit the model on the real data, a cost function on R 12 must be defined since the model depends on 12 parameters. Indeed, as mentioned earlier, the model depends on five parameters for the point processes, six for the two functions acting on the distance maps and defined by (3), and one smoothing parameter. The proposed cost function is defined as follows.

F (ω) = 2 i=1 ∆ (∥s i (z)∥ L2 ) + 0≤i<j≤2 ∆(S i,j ) + ∆(ξ TPBl ) (4)
where ω ∈ R 12 is a set of model parameters, ∆ (∥s i (z)∥ L2 ) is the relative error on the L 2 norm of the phase surface area density function s i (z) of phase i, ∆(S i,j ) is the relative error on the interfacial specific surface area between phases i and j and ∆(ξ TPBl ) the relative error of TPBl density. The fitting process then consists of applying the CMA-ES algorithm [START_REF] Anne | A Restart CMA Evolution Strategy With Increasing Population Size[END_REF] to the cost function defined by equation ( 4) above.

Fitting process validation

In order to validate the method, it is necessary to apply the process to a synthetic structure generated by the model using a set of known parameters ω 0 and compare it to the set of parameters ω 0 obtained after minimizing the cost function. Table 2 and Fig. 6 illustrate the result obtained at the end of the fitting process. It can be seen that the relative error on each of the parameters of the model is relatively small, in the order of 3% to 4% on average. The parameters λ 0 and λ 1 correspond to the intensity of the two Poisson point processes of the two first set of points. The parameters λ p (parents), λ o (offspring) and σ control the Neyman-Scott point process. The parameter θ is the smoothing parameter and a i , b i and c i define the function f i as described in equation ( 3) for i ∈ {1, 2}.

Thus, Fig. 6 shows that the phase surface area density of each of the phases is very well approximated. The fitting process being validated, the following section presents the results obtained by applying it to real microstructures.

Model validation

The optimization process is performed on the two microstructures shown in Fig. 1 and Fig. 2 and the results are presented on Fig. 7 and Fig. 8. The fit on the YSZ phase is always good, whatever the microstructure considered. On the other hand, the fit on the Ni phase and the pore phase are more challenging.

Nevertheless, Table 3 and Table 4 show that the properties of the synthetic microstructures generated by the model after optimization are very close to those of the real microstructures. The higher errors on the interfacial specific surface areas are due to the simplicity of the model which, despite the use of a Neyman-Scott type aggregated point process, is not yet able to reproduce perfectly the Ni phase coarsening.

Fig. 6, 7 and 8 also show that the more locally homogeneous the microstructure, the better the results. This last hypothesis has been tested and confirmed on synthetic microstructures, in particular those shown in Fig. 5.

Conclusion

An original 3D stochastic model allowing to simulate real microstructures with strong gradient and global heterogeneity has been proposed. The flexibility of the model allows to generate microstructures of various local homogeneity, and with very diverse characteristics. Thus, it was possible to fit the model on real microstructure reconstructions with rather satisfactory results, both in terms of phase surface area density curves and global characteristics of the microstructure.

It became clear that the greater the local homogeneity, the better the results. Nevertheless, the primary objective of this work, namely the simulation of phase gradient globally heterogeneous microstructure has been achieved.

In the future, it would be worthwhile to extend the model to generate microstructures even closer to those shown in Fig. 2, where the Nickel phase forms very large clusters.

Figure 1 :

 1 Figure 1: (Left) 3D reconstruction of a real microstructure before use. Cube of 17.5 µm side (350 3 voxels). The porous phase is in black, the ceramic phase in grey and the metallic phase in white. (Right) Evolution of the phase surface area density associated with each phase along the Z-axis.

Figure 2 :

 2 Figure 2: (Left) 3D reconstruction of a real microstructure after 1000 h of use. Cube of 17.5 µm side (350 3 voxels). The porous phase is in black, the ceramic phase in grey and the metallic phase in white. (Right) Evolution of the phase surface area density associated with each phase along the Z-axis.

Figure 3 :

 3 Figure 3: Example of different point processes. From left to right: Poisson point process of intensity λ = 100, Neyman-Scott point process with λ p = λ o = 10 and σ = 0.075, and Neyman-Scott point process with λ p = λ o = 10 and σ = 0.025.

Figure 4 :

 4 Figure 4: Illustration of the microstructure generation process (2D slice equivalent). On the left, generation of the three point sets by two Poisson processes (red and green) and a Neyman-Scott process (blue). In the center, distance map of the ceramic (gray) phase and resulting microstructure. On the right, application of a sigmoid function to the distance map of the ceramic (gray) phase and resulting microstructure.

Figure 5 :

 5 Figure 5: Example of different realizations of synthetic microstructures with phase gradient from the proposed model on a domain W = [0; 350] 3 .

Figure 6 :

 6 Figure 6: Comparison of phase surface area density for the original synthetic microstructure (solid lines) and the synthetic microstructure obtained after optimization of the cost function (dashed lines). The area around those corresponds to the standard deviation. The corresponding microstructure is shown on the right in Fig. 5.

Figure 7 :

 7 Figure 7: (Left) 3D reconstruction of a synthetic microstructure with the same properties as the microstructure in Fig. 1. (Right) Evolution of the phase surface area density. Real data (solid lines) and simulations (dashed lines). The area around those corresponds to the standard deviation. Average over a set of 8 realizations.

Figure 8 :

 8 Figure 8: (Left) 3D reconstruction of a synthetic microstructure with the same properties as the microstructure in Fig. 2. (Right) Evolution of the phase surface area density. Real data (solid lines) and simulations (dashed lines). The area around those corresponds to the standard deviation. Average over a set of 8 realizations.

Table 2 :

 2 Comparison of the values used to generate a synthetic microstructure with those returned by the fitting process.

	Parameters Original values (ω 0 ) Fitted values ( ω 0 ) Error
	λ 0	12.67 µm -3	12.24 µm -3	-3%
	λ 1	15.70 µm -3	14.98 µm -3	-5%
	λ p	0.95 µm -3	0.98 µm -3	4%
	λ o	19.51	18.9	-3%
	σ	0.182 µm	0.189 µm	4%
	θ	3.31	3.25	-2%
	a 1	-0.61	-0.59	3%
	b 1	0.16	0.16	-1%
	c 1	-55.44	-57.24	3%
	a 2	2.04	1.96	-4%
	b 2	0.22	0.21	-4%
	c 2	-42.56	-39.84	-6%

Table 3 :

 3 Comparison of the properties of the real microstructure and the synthetic microstructure in Fig.7.

	Properties Real values	Simulation	Error
	ε 0	0.253	0.250	-1%
	ε 1	0.424	0.433	2%
	ε 2	0.323	0.317	-1.8%
	S 0,1	0.878 µm -1 0.890 µm -1	1.4%
	S 0,2 p S 1,2 p ξ TPBl	0.236 µm -1 0.230 µm -1 -2.5% 0.728 µm -1 0.686 µm -1 -5.7% 1.200 µm -2 1.162 µm -2 -3.3%

Table 4 :

 4 Comparison of the properties of the real microstructure and the synthetic microstructure in Fig.8.

	Properties Real values	Simulation	Error
	ε 0	0.250	0.248	-0.8%
	ε 1	0.439	0.439	0%
	ε 2	0.310	0.312	0.5%
	S 0,1 p S 0,2 p S 1,2 p ξ TPBl	0.718 µm -1 0.842 µm -1 0.144 µm -1 0.138 µm -1 -4.2% 17% 0.470 µm -1 0.442 µm -1 -5.9% 0.640 µm -2 0.640 µm -2 0%

Acknowledgment

The author(s) acknowledge(s) the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-18-CE05-0036 (project ECOREVE).