Kévin Potier

Kristina Ariskina

Amaël Obliger

Jean-Marc Leyssale
email: jean-marc.leyssale@u-bordeaux.fr

Molecular Simulation of Argon Adsorption and Diffusion in a Microporous Carbon with Poroelastic Couplings

come

Introduction

Gas adsorption in microporous carbons is a fundamental thermodynamic process at the core of numerous technological applications in the fields of energy supply and environment protection like gas storage 1 and separation, 2 mitigation of pollutants from wastewater, [START_REF] Bandosz | Activated Carbon Surfaces in Environmental Remediation[END_REF] or gas (and oil) recovery from coal bed or shale reservoirs. [START_REF] Omari | The progress of enhanced gas recovery (EGR) in shale gas reservoirs: A review of theory, experiments, and simulations[END_REF] Furthermore, gas adsorption is also one of the most suited characterization technique to investigate the pore network of microporous carbons, allowing to determine useful quantities such as pore surface, volume and size distributions. [START_REF] Stoeckli | Microporous carbons and their characterization: The present state of the art[END_REF][START_REF] Rodríguez-Reinoso | Textural and chemical characterization of microporous carbons[END_REF] Microporous carbons have been long-known for being prone to significant adsorptioninduced deformations [START_REF] Bangham | The Expansion of Charcoal accompanying Sorption of Gases and Vapours[END_REF] even though only few experimental reports are available so far. [START_REF] Gor | Adsorption-induced deformation of nanoporous materials -A review[END_REF] These effects can be particularly important in the case of compliant natural porous carbons like coal, for which swelling of a few percent in volume have been reported due to high pressure CO 2 adsorption. [START_REF] Ottiger | Competitive adsorption equilibria of CO2 and CH4 on a dry coal[END_REF] Although numerical simulations like density functional theory (DFT) and grand canonical Monte Carlo (GCMC) simulations have become increasingly popular for predicting and rationalizing gas adsorption, adsorption-induced deformations are generally neglected in such simulations. Kowalczyk et al. proposed a first thermodynamic model based on GCMC simulations of Ar adsorption in graphitic slit pores at 243 K, to account for adsorption induced deformations. [START_REF] Kowalczyk | Adsorption-Induced Deformation of Microporous Carbons: Pore Size Distribution effect[END_REF] In the same publication, the model was applied to the experimental results of Yakovlev et al. [START_REF] Yakovlev | Adsorption-stimulated deformation of microporous carbon adsorbent[END_REF] to predict the pore size distribution (PSD) and bulk modulus of a car-bide derived microporous carbon based on adsorption and swelling data. Few years later, Brochard et al. [START_REF] Brochard | Poromechanics of microporous media[END_REF] proposed a linearized microporomechanical model that allows predicting the adsorption-induced volumetric strain from the sorption isotherm computed with GCMC in the non-deformable matrix. It was then applied to predict coal swelling under competitive CO 2 /CH 4 adsorption [START_REF] Brochard | Adsorption-Induced Deformation of Microporous Materials: Coal Swelling Induced by CO 2 -CH 4 Competitive Adsorption[END_REF] and compared to the experimental data of Ottiger et al. [START_REF] Ottiger | Competitive adsorption equilibria of CO2 and CH4 on a dry coal[END_REF] This model, which requires knowledge of the material's bulk modulus K and poroelastic coupling constant C (see below) was used in other studies of coal [START_REF] Zhang | Molecular simulation of CO 2 -CH 4 competitive adsorption and induced coal swelling[END_REF] and kerogen [START_REF] Huang | Molecular Insights into Kerogen Deformation Induced by CO 2 /CH 4 Sorption: Effect of Maturity and Moisture[END_REF] swelling, yet assuming the same values of K (2.65 GPa) and C -6.05 and 7.60 for CH 4 and CO 2 , respectivelythan in Brochard et al. [START_REF] Brochard | Poromechanics of microporous media[END_REF][START_REF] Brochard | Adsorption-Induced Deformation of Microporous Materials: Coal Swelling Induced by CO 2 -CH 4 Competitive Adsorption[END_REF] However, while the models by Kowalczyk et al. and Brochard et al. allow predicting the volumetric strain induced by adsorption, they both create a disconnection between volumetric and adsorption data as the associated change in pore volume on the sorption isotherms is not considered.

Molecular simulation of gas adsorption in microporous carbon explicitly accounting for adsorption-induced deformations comes with some significant additional difficulty and computational cost compared to the usual GCMC calculations under the rigid (or frozen) matrix approximation in which the atoms constituting the porous carbon are frozen. [START_REF] Obliger | Development of Atomistic Kerogen Models and Their Applications for Gas Adsorption and Diffusion: A Mini-Review[END_REF] In particular, it is crucial that the atomistic model of the host material captures well its pore space (PSD), composition (atomic H/C, O/C, etc.) and elastic properties (especially its bulk modulus K for an isotropic material). This disqualifies oversimplified models based on graphitic slit pores or disordered models based on functionalized aromatic flakes, either flat or curved. [START_REF] Di Biase | Systematic development of predictive molecular models of high surface area activated carbons for adsorption applications[END_REF] Also of importance, the model, as described by the adopted force-field, has to be stable so that no major chemical or structural change can occur during the simulation.

Using Molecular Dynamics (MD) simulations, Obliger et al. [START_REF] Obliger | Poroelasticity of Methane-Loaded Mature and Immature Kerogen from Molecular Simulations[END_REF] have investigated the swelling of immature (aliphatic matrix with high H/C ratio) and overmature (aromatic with H/C ∼ 0) kerogen as a function of methane loading under explicit shale reservoir conditions (T and P). They have observed significant differences between the two matrices, showing that immature kerogen can swell up to a 20 % increase in volume while the swelling of the overmature kerogen in the same conditions remains far below 1 %. The large swelling experienced by the immature kerogen was subsequently shown to lead to an increase in diffusion coefficient with increasing loading at constant reservoir (P,T) conditions, [START_REF] Obliger | Methane Diffusion in a Flexible Kerogen Matrix[END_REF][START_REF] Ariskina | Free Volume Model for Transport in Flexible Kerogen of Source Rock's Organic Matter[END_REF] which is the opposite behavior to the commonly observed one for rigid hosts. More recently, a few studies have been reported in which adsorption-induced deformations in kerogen are accounted for via an hybrid method alternating GCMC cycles and short MD simulations in the isothermal-isobaric (NPT) ensemble. [START_REF] Yu | Determination of CH 4 , C 2 H 6 and CO 2 adsorption in shale kerogens coupling sorption-induced swelling[END_REF][START_REF] Chong | Molecular characterization of carbon dioxide, methane, and water adsorption in micropore space of kerogen matrix[END_REF][START_REF] Wu | Molecular investigation on CO 2 -CH 4 displacement and kerogen deformation in enhanced shale gas recovery[END_REF] Doing this, the authors simulate an osmotic ensemble [START_REF] Coudert | Thermodynamics of Guest-Induced Structural Transitions in Hybrid Organic-Inorganic Frameworks[END_REF] which corresponds to unjacketed adsorption conditions where the fluid and host are simulateneously put under thermal, chemical and mechanical equilibrium. In all these works, significant swelling, of 1-10 % under typical geological pressures were observed, even for supposedly mature kerogen. [START_REF] Yu | Determination of CH 4 , C 2 H 6 and CO 2 adsorption in shale kerogens coupling sorption-induced swelling[END_REF][START_REF] Wu | Molecular investigation on CO 2 -CH 4 displacement and kerogen deformation in enhanced shale gas recovery[END_REF] While these recent investigations constitute a significant body of progress, many questions remain unanswered regarding adsorption-induced deformations in microporous carbon. First, as pointed out in ref. 18, viscoelastic effects can be important for such materials, which clearly question the sampling time used in simulations and the error bars, rarely reported, regarding obtained modulus, adsorption, swelling and diffusion results. Second, the effect of poromechanics on sorption isotherms (i.e. the amount of adsorbed fluid at given T and P) has not been reported either. Third, the validity of the poromechanical model proposed in ref. [START_REF] Brochard | Poromechanics of microporous media[END_REF] has not really been discussed in literature, aside from a brief discussion in ref., [START_REF] Obliger | Poroelasticity of Methane-Loaded Mature and Immature Kerogen from Molecular Simulations[END_REF] and the same applies to the evolution of pore space induced by adsorption.

Furthermore, common to all the recent hybrid GCMC/MD studies [START_REF] Yu | Determination of CH 4 , C 2 H 6 and CO 2 adsorption in shale kerogens coupling sorption-induced swelling[END_REF][START_REF] Chong | Molecular characterization of carbon dioxide, methane, and water adsorption in micropore space of kerogen matrix[END_REF][START_REF] Wu | Molecular investigation on CO 2 -CH 4 displacement and kerogen deformation in enhanced shale gas recovery[END_REF] is the use of the so-called "molecular" kerogen models of Ungerer et al. [START_REF] Ungerer | Molecular Modeling of the Volumetric and Thermodynamic Properties of Kerogen: Influence of Organic Type and Maturity[END_REF] These models, constructed by packing small and identical molecules, accurately reproduce the chemistry of the material, and thus the strength of host-guest interactions. However, due to the limited size of their constituents, it is unclear whether they can accurately capture poromechanics, which is governed by the nature and connectivity of the carbon skeleton of the host. As discussed in a recent review, [START_REF] Obliger | Development of Atomistic Kerogen Models and Their Applications for Gas Adsorption and Diffusion: A Mini-Review[END_REF] models constructed without making any assumption on the size of the constituents, as is the case for models obtained using hybrid reverse Monte-Carlo (HRMC) [START_REF] Jain | Molecular modeling and adsorption properties of porous carbons[END_REF][START_REF] Bousige | Realistic molecular model of kerogen's nanostructure[END_REF] or the replica exchange molecular dynamics (REMD) method [START_REF] Atmani | From cellulose to kerogen: molecular simulation of a geological process[END_REF][START_REF] Atmani | Simulating the Geological Fate of Terrestrial Organic Matter: Lignin vs Cellulose[END_REF][START_REF] Leyssale | Replica Exchange Molecular Dynamics Simulation of Organic Matter Evolution: From Lignin to Overmature Type III Kerogen[END_REF] appear as better alternatives, even though the chemistry can be less accurately pictured.

In this work we present an exhaustive GCMC/MD investigation of the unjacketed adsorption of Ar at 243 K in a mature kerogen model, as a typical microporous aromatic carbon with moderate H/C ratio (0.5), obtained by simulating the geological evolution of lignin with the REMD method. The chosen conditions coincides with those adopted in the seminal experimental report by Yakovlev et al. [START_REF] Yakovlev | Adsorption-stimulated deformation of microporous carbon adsorbent[END_REF] and theoretical analysis by Kowalczyk et al. [START_REF] Kowalczyk | Adsorption-Induced Deformation of Microporous Carbons: Pore Size Distribution effect[END_REF] Also Ar adsorption at such temperature is a common case study for adsorption studies in carbon [START_REF] Malbrunot | Adsorption measurements of argon, neon, krypton, nitrogen, and methane on activated carbon up to 650 MPa[END_REF][START_REF] Fan | Computer Simulation of Argon Adsorption on Graphite Surface from Subcritical to Supercritical Conditions: The Behavior of Differential and Integral Molar Enthalpies of Adsorption[END_REF] and has the advantage, with respect to other gas like CO 2 , N 2 , etc., of not presenting any specific electrostatic interactions with functional groups in the adsorbent, which presence is neglected in the model. After a brief presentation of the porous carbon model and simulation methods, we present in sequence results regarding (i) the viscoelastic properties and PSD of the host, (ii) the effect of swelling and matrix flexibility at constant volume on the sorption isotherm, (iii) the adsorption-induced volumetric strain and its analysis through Brochard's poromechanical model, (iv) the evolution with gas pressure of the host pore space and pore occupancy, and (v) the effect of adsorption-induced deformations on the gas self-diffusion coefficient.

Methods

Porous carbon model

The microporous carbon model used in this work is shown in Fig. 1. It was obtained by simulating the geological evolution of lignin at 423 K and 25 MPa using Replica Exchange Molecular Dynamics (REMD) [START_REF] Atmani | From cellulose to kerogen: molecular simulation of a geological process[END_REF][START_REF] Atmani | Simulating the Geological Fate of Terrestrial Organic Matter: Lignin vs Cellulose[END_REF][START_REF] Leyssale | Replica Exchange Molecular Dynamics Simulation of Organic Matter Evolution: From Lignin to Overmature Type III Kerogen[END_REF] and a ReaxFF forcefield. [START_REF] Van Duin | A reactive force field for hydrocarbons[END_REF] Prior to this work, the final ex-lignin model in ref. 29 was slightly modified. First, the 15 oxygen-containing groups, mostly hydroxyls, were replaced by equivalent oxygen-free groups (H for OH and CH 2 for ether O) as the force field used in the simulations (see below) does not have parameters for O atoms. This should not affect the results because i) 15 O atoms are very little compared to the 1377 C atoms and 682 H atoms in the model and ii) Ar only interacts with kerogen via van der Waals interactions which are of similar magnitude with C and O atoms. Second, the model was duplicated in the shorter periodic direction (X in Fig. 1) to increase its volume.

The resulting model contains 2754 carbon atoms and 1364 hydrogen atoms enclosed in an orthorhombic cell of ∼ 48 × 39 × 34 Å3 , with periodic boundary conditions in all directions.

Its H/C ratio, 0.495, locates it in the mature coal/kerogen category according to the definition by Van Krevelen, [START_REF] Van Krevelen | Coal 3 rd edition: Typology -Physics -Chemistry -Constitution[END_REF][START_REF] Vandenbroucke | Kerogen origin, evolution and structure[END_REF] implying a dominant aromatic structure. Detailed structure analysis reveals that carbon atom hybridizations, determined from the numbers of closest neighbors, are 87.4 % sp 2 , 7.9 % sp 3 and 4.7 % sp (or radical sp 2 ). Texture analysis shows that 83.4 % of the carbon atoms are involved in ring clusters, with an average cluster size of 33 C atoms (9 rings) while the remaining atoms form small chain structures (2 atoms on average). Rings are predominantly hexagons (51 %), pentagons (34 %) and heptagons (9 %). The large amount of pentagons with respect to heptagons indicates a significant (positive) curvature of the carbon layers, ensuring an isotropic texture.

As is common practice in molecular simulation, [START_REF] Obliger | Poroelasticity of Methane-Loaded Mature and Immature Kerogen from Molecular Simulations[END_REF][START_REF] Falk | Subcontinuum mass transport of condensed hydrocarbons in nanoporous media[END_REF][START_REF] Falk | Effect of Chain Length and Pore Accessibility on Alkane Adsorption in Kerogen[END_REF] the kerogen model considered in this work is used as an archetypal model of disordered microporous aromatic carbon. As discussed in the introduction, and in more details in ref. 16, the main factors governing the adsorption and poromechanical behavior of a porous carbon with respect to a non-polar fluid adsorbate are the H/C ratio and the porosity, which are both strongly related to the bulk modulus. With a bulk modulus of 1-2 GPa, a PSD extending from 3.5 to 10.5 Å and a porous volume at rest of 0.37 cm 3 /g (see below for the details), this model is highly relevant for a comparison with the carbide derived carbon studied by Kowalczyk et al. [START_REF] Kowalczyk | Adsorption-Induced Deformation of Microporous Carbons: Pore Size Distribution effect[END_REF] who predicted values of 7 GPa for the bulk modulus, 3.5-8 Å for the PSD and 0.47 cm 3 /g for the pore volume -note that the error bars on these data, not given in the publication, are certainly quite significant considering the underlying assumptions. The model properties also show intermediate values between those of the popular CS400, CS1000 and CS1000a microporous carbons prepared by annealing (and activating) saccharose cokes [START_REF] Pikunic | Structural Modeling of Porous Carbons: Constrained Reverse Monte Carlo Method[END_REF][START_REF] Jain | Effects of Activation on the Structure and Adsorption Properties of a Nanoporous Carbon Using Molecular Simulation[END_REF][START_REF] Jain | Molecular Modeling of Porous Carbons Using the Hybrid Reverse Monte Carlo Method[END_REF] and commonly used as proxies for coal [START_REF] Brochard | Adsorption-Induced Deformation of Microporous Materials: Coal Swelling Induced by CO 2 -CH 4 Competitive Adsorption[END_REF] or kerogen. [START_REF] Falk | Subcontinuum mass transport of condensed hydrocarbons in nanoporous media[END_REF][START_REF] Falk | Effect of Chain Length and Pore Accessibility on Alkane Adsorption in Kerogen[END_REF][START_REF] Lee | Activated desorption at heterogeneous interfaces and long-time kinetics of hydrocarbon recovery from nanoporous media[END_REF]

Molecular simulations

Interatomic interactions within the carbon matrix are described with the AIREBO potential [START_REF] Stuart | A reactive potential for hydrocarbons with intermolecular interactions[END_REF] with the reparameterization of the Lennard-Jones terms for kerogen/hydrocarbon systems proposed by Obliger et al. [START_REF] Obliger | Poroelasticity of Methane-Loaded Mature and Immature Kerogen from Molecular Simulations[END_REF] This potential has been shown to be well suited to investigate poromechanics and transport in kerogen. [START_REF] Obliger | Poroelasticity of Methane-Loaded Mature and Immature Kerogen from Molecular Simulations[END_REF][START_REF] Obliger | Methane Diffusion in a Flexible Kerogen Matrix[END_REF] Argon is described using a single Lennard-Jones site and the Lorentz-Berthelot combining rules are applied for cross terms. [START_REF] Allen | Computer Simulation of Liquids[END_REF] A cutoff of 3 σ ij , where i and j refers to the types of atoms, is applied to van der Waals interactions. The Lennard-Jones parameters are summarized in table 1.

A combination of molecular dynamics simulations in the canonical ensemble (i. e. at constant number of atoms N , volume V and temperature T ), isobaric-isothermal ensemble (i.e. at constant N , pressure P and T ) and Monte Carlo simulations in the grand canonical ensemble (GCMC) at constant V , T and Ar chemical potential (µ Ar ) is used to investigate the behavior of Ar at 243 K in the carbon matrix. These simulations are all performed with the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package. [START_REF] Plimpton | Fast Parallel Algorithms for Short-Range Molecular Dynamics[END_REF][START_REF] Thompson | LAMMPS -a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales[END_REF] Molecular dynamics simulations use an integration timestep of 0.25 fs. A three-member Nosé-Hoover chain method, with a time constant of 100 fs is used to fix the temperature (NVT and NPT) and an anisotropic Nosé-Hoover barostat with a time constant of 500 fs is used to fix independently the three diagonal components of the stress tensor, [START_REF] Martyna | Explicit reversible integrators for extended systems dynamics[END_REF] ensuring a perfectly isotropic pressure within the system. GCMC simulations are used to simulate Ar adsorption in the carbon matrix. They consist of random attempt to insert or suppress an Ar atom in the volume of the matrix, with equal probability of insertion and deletion attempts, accepted or rejected with the usual GCMC acceptance criterion. [START_REF] Frenkel | Understanding Molecular Simulation[END_REF] For simplicity, the fugacity f = exp -µ Ar k B T , which corresponds to the pressure of the ideal gas in equilibrium with the adsorbed Ar gas, is fixed in the GCMC simulations instead of the chemical potential. The relationship between f and the pressure P of the real Ar fluid was computed via GCMC simulation of the bulk fluid (Fig. S1) and fitted with polynomials, allowing to plot sorption isotherms as a function of P .

Prior to the unjacketed adsorption simulations, it is important to determine the volumetric properties of the matrix -volume and porosity -which will serve as a reference state to characterize adsorption-induced deformations at the suited temperature. The carbon matrix is thus extensively relaxed using NPT MD at 243 K and zero pressure for 600 ns. The initial 100 ns are discarded from the analysis of the average properties of the empty matrix and ten independent configurations are stored for porosity analysis. the matrix is both at chemical and mechanical equilibrium with a bulk fluid at pressure P and temperature T , fully capturing adsorption-induced deformations. For this, the simulation breaks down to an alternance of a short GCMC run (2.5 × 10 5 adsorption or deletion attempts at fixed f, V and T ) during which only the gas degrees of freedom are sampled and a short MD run (at 4 × 10 5 MD steps at fixed N, P and T ) during which both the gas and matrix degrees of freedom are considered, defining a cycle. The NPT simulation is run for a total of 3000 cycles, amongst which 1000 are used to bring the system to both chemical (via GCMC) and mechanical (via NPT MD) equilibrium, where both loading and volume fluctuate around constant average values. From the statistical mechanics perspective, this hybrid GCMC/MD simulation is an accurate approximation of the osmotic ensemble. [START_REF] Coudert | Thermodynamics of Guest-Induced Structural Transitions in Hybrid Organic-Inorganic Frameworks[END_REF][START_REF] Zang | Osmotic ensemble methods for predicting adsorptioninduced structural transitions in nanoporous materials using molecular simulations[END_REF] The two additional types of simulation are merely aimed at evaluating the effect of commonly neglected matrix degrees of freedom in the adsorption simulation literature. The NVT scheme is almost identical to the NPT case, except that the MD runs are performed in the NVT ensemble, the volume being fixed at its equilibrium value for the (empty) matrix at zero pressure. The NVT scheme thus allows for internal rearrangements of the matrix along gas adsorption yet with no net volume change. Finally, the RIGID scheme is a pure GCMC simulation in which no matricial mobility is accounted for. Equilibration being much faster in the NVT and RIGID schemes than in the NPT scheme, as no volume change is taking place, 500 equilibration and run cycles are sufficient to obtain well converged properties.

Ar adsorption is investigated under three different approximations (see table 2). The

The convergence of properties in the equilibration stage was carefully monitored for all the considered conditions. Examples of the convergence of the number of adsorbed Ar atoms and volume obtained at 22.2 bars with the NPT scheme are shown in Fig. S2 andS3, respectively.

Adsorption simulations are performed at fixed pressure values (exponentially spaced) in the 0.05-2116 bar range. The extended pressure range, and especially the largest, somehow unrealistic, values, is considered to better capture trends in the poromechanical behavior (swelling) which is subject to large relative uncertainties at low and moderate pressures (see below). Sorption isotherms are then fitted using the Tóth isotherm model: [START_REF] Tóth | Uniform interpretation of gas/solid adsorption[END_REF][START_REF] Terzyk | Developing the solution analogue of the Toth adsorption isotherm equation[END_REF] ω(P ) = ω max κ T P

1 + (κ T P ) t 1/t (1)
where ω is the adsorbed amount (in mmol/g) and ω max , κ T and t are adjustable parameters.

Porosity characterization

Geometric algorithms similar to those describe in refs. 18,19,29 are used to characterize the porous structure of the matrix and its evolution with Ar pressure. First, a regular orthorhombic mesh is constructed over the matrix volume and all the lattice vertices which can possibly be part of an Ar atom, without overlapping with the van der Waals spheres of the matrix atoms, are considered as part of the accessible porous volume V p . The latter is computed as the sum of the "porous" lattice elements multiplied by the volume of lattice elements.

The accessible porous surface S sp is calculated using the Marching Cube algorithm, [START_REF] Lorensen | Marching Cubes: A High Resolution 3D Surface Construction Algorithm[END_REF] implemented in the scikit-image library. [START_REF] Van Der Walt | the scikit-image contributors, scikit-image: image processing in Python[END_REF] This algorithm makes use of the previously determined porosity grid to determine the surface of division S sp between the pores and the matrix elements.

We note that while V p converges to the exact accessible volume when the voxel size tends towards zero, S p only converges in the specific case of non-curved surfaces. Therefore, the calculated surfaces will be overestimated by at most 10 %, which corresponds to the deviation obtained in the limiting case of maximum curvature (sphere).

Other interesting descriptors of the porous volume are the accessible porosity defined as the ratio of the accessible porous volume over total volume ϕ a = Vp V , the accessible free volume defined as

ϕ f = Vp-Vo V
, where V o is the volume occupied by the fluid (sum of the volume having their centers within [START_REF] Serafin | Highly microporous activated carbons from biomass for CO2 capture and effective micropores at different conditions[END_REF] 2 σ Ar of the center of an Ar atom) and the pore occupancy ratio χ = Vo Vp . Finally, the pore size distribution (PSD) is computed exactly as in the work by Atmani et al., [START_REF] Atmani | Simulating the Geological Fate of Terrestrial Organic Matter: Lignin vs Cellulose[END_REF] by first locating all the mesh centers than can accomodate the center of an Ar atom without van der Waals overlap, and for each of them performing a random search of the largest spherical pore than can enclose it. The separation between mesh centers is fixed to about 0.1 Å for the fast converging PSD calculations and to about 0.025 Å for the slow converging porous volume and surface calculations.

Error bars on all properties are computed as the standard error on the mean from five independent configurations (porosity) or block averages (other quantities), given with a 99 % confidence interval (3 σ).

Results

Fig. 2 one, operating on a subnanosecond scale, correspond to purely elastic vibrations around an almost constant average volume value, over several tens of nanoseconds. The slowest one shows a few abrupt changes in the average volume value and is associated to internal reorganizations of some aromatic fragments within the matrix (see below). For simplicity, we associate these motions, involving the crossing of energy barriers to sample different potential energy wells, to viscous deformations as opposed to the elastic deformations in which the system remains in the same potential well, although acknowledging that experimentally, such motions would probably not be considered as viscous due to their short timescale (∼ 100 ns).

Furthermore, the matrix can still be considered as rigid despite these limited reorganizations that remain mostly local, as opposed to a truly viscoelastic material. Nevertheless, This timescale separation explains the need for extensive simulation time (here 500 ns) to properly capture average volumetric properties and their fluctuations. Fig. 3 shows the time-dependent bulk modulus K τ computed using the volume fluctuation

equation 1 Kτ = V 2 τ -V 2 τ k B T V τ
, where k B is the Boltzmann's constant and ... increasing time, down to a plateau at which all phonons are accounted for. In the present case, this plateau seems to be located around 10 ns, at which a pure elastic bulk modulus K e ∼ 1.7 ± 0.1 GPa is observed. However, integrating to larger times, we observe a second decrease in K τ , albeit large error bars, indicating the occurrence of viscoelastic deformations.

While longer simulations would be required to really converge to the viscoelastic plateau, we here adopt the value at the largest averaging time (167 ns) as the viscoelastic modulus K ve = 1.2 ± 0.7 GPa. We note that, although this value is probably an upper estimate, the error is certainly compensated by the large error bar (0.7 GPa).

To further illustrate the equilibrium viscoelastic deformations taking place in the porous carbon matrix, we show in Fig. 4(a) the distribution of non-affine atomic displacements u na between the configurations at 80 and 180 ns (i. e., the configurations of maximum and minimum porosity, respectively, in Fig. 2). The non-affine displacement of an atom between configuration A and B is defined as

u na α = r B α -r A α -
α = r A α × L B α -L A α L A α , L B
α and L A α being the lengths in the α direction of the deformed and undeformed cell, respectively. As can be seen in Fig. 4, while the vast majority of atoms show only small displacements, with more than 95 % being within 2 Å, some atoms show significantly larger displacements, up to about 1 nm. They correspond to an internal reorganization of some groups within the matrix, as shown in Fig. 4(b,c), that allows for a significant increase in the local density, explaining the important decrease in the matrix volume observed in Fig. 2.

Overall, the average volume of the matrix, V = 63.9 ± 0.6 nm 3 (Fig. 2), leads to the apparent density of the matrix d a = 0.89 ± 0.02 g/cm [START_REF] Bandosz | Activated Carbon Surfaces in Environmental Remediation[END_REF] . Considering now the average accessible porosity, ϕ a = 0.330 ± 0.007 (Fig. 2), the skeleton (or wall) density reads d s = d a /(1 -ϕ a ) = 1.33 ± 0.08 g/cm 3 for an argon probe. Fig. 5 shows that the matrix present a rather homogeneous pore size distribution (PSD), from the Ar probe size (3.4 Å) to about 11 Å, confirming the microporous nature of the material. The average pore size φ = 7.09 ± 0.09 Å.

Fig. 6 shows the sorption isotherms computed using the NPT, NVT and RIGID approaches. As can be seen in the inset, neglecting adsorption-induced deformations, under Normalized PSD φ(Å) either the RIGID or NVT frameworks, does not have a significant effect on low pressure adsorption. The obtained Henry's constants, 1.10, 1.04 and 1.09 mmol/g/bar for the NPT, NVT and RIGID schemes, respectively, are very similar. While the three data sets are well fitted with the Tóth isotherms in the whole pressure range, significant quantitative differences are observed amongst them at large pressures, typically above 10 bar. Neglecting matrix swelling (NVT) and relaxation (RIGID) indeed lead to significant underestimation of the high pressure adsorption capacity, up to about 19 (NVT) and 28 (RIGID) % at saturation (see Fig. S5). An important consequence of performing the adsorption simulations at constant volume is that, while the chemical potentials of the free (reservoir) and adsorbed gas are equal, defining a chemical equilibrium, the matrix and the reservoir deviate from mechanical equilibrium.

Figure 7 shows the evolution of the virial pressure of the simulated system P v with respect to the reservoir pressure P computed for the NVT and NPT simulations. In the NVT case, while P v ≈ P at low P, a significant deviation becomes apparent when P exceeds a few bar, then progressively increases up to a positive deviation of about 1000 bar at P ≈ 1000 bar. This is clear evidence of the unappropriateness of constant volume conditions (either NVT ot RIGID) in adsorption simulations at large pressures. Conversely, as shown in Fig. 7, within the NPT scheme, P v and P remains almost perfectly equal, whatever P, ensuring both perfect chemical and mechanical equilibria and, thus, perfectly simulating unjacketed adsorption. The volumetric strain of the matrix induced by adsorption ε

V = V -V 0 V 0
where V and V 0 are the average matrix volumes computed from the NPT simulations at pressure P and at 0 pressure (i.e. in the absence of adsorbed fluid), respectively, is given as a function of P in Fig. 8. A first important point to notice is that comparing the swelling observed under NPT conditons to the increase in P v observed in the NVT case (Fig. 7) clearly shows that matrix swelling is responsible for the pressure equilibration in the NPT scheme.

Going into further details, we observe that no significant swelling strain can be detected at P below ∼ 5 bar, which agrees well with the results of Yakovlev et al. to methane adsorption at room temperature were observed, seem to confirm this tendency.

Yet, definitive proof of these observations would require much larger simulation times to reduce error bars.

Above 5 bar, as in refs. 11 and 10, a clear increase in strain is observed with increasing pressure, up to a value of about 8 % at P = 1000 bar. This evolution is significantly sublinear at large pressures. Hence, two functions are proposed to fit the evolution of ε V with P . The first one is based on the observation that ε V evolves mostly linearly with Ar loading in the whole pressure range (see the linear fit of the evolution of ε V with ω N P T in Fig. S6) justifying the use of a Tóth isotherm function (see Eq. 1). We note that a similar linear correlation between volumetric strain and loading was recently observed in the case of methane adsorption on Barnett and Eagle Ford shale samples. [START_REF] Pang | Evaluation of Matrix Swelling Behavior in Shale Induced by Methane Sorption under Triaxial Stress and Strain Conditions[END_REF] The second form used is a stretched exponential function : ε V = c × 1 -exp -a * P b where a, b and c are adjustable parameters. The main difference between the two models is that in the latter swelling plateaus at high pressure. However, considering the size of the error bars, it is not possible to conclude on whether it is, or not, more realistic.

Some years ago, Brochard et al. have proposed a microporomechanical model relating

the adsorption isotherm of a deformable matrix to the corresponding (virtual) isotherm of the same matrix, yet underformable, via a first order expansion in volumetric strain. [START_REF] Brochard | Poromechanics of microporous media[END_REF] In this model, valid at small deformations (ε V << 1), the unjacketed isotherm n(P ) reads

n(P ) ≈ n 0 (P ) (1 + C(P )ε V (P )) (2)
where n(P ) and n 0 (P ) are the molar conentrations of fluid in the matrix at P with respect to the volume V 0 of the undeformed matrix, and C(P) is a pressure dependent coupling coefficient. Noticing that n(P ) = ω(P )m 1000V 0 and n 0 (P ) = ω 0 (P )m 1000V 0 where ω and ω 0 are the isotherms computed for the deformable (ω = ω N P T ) and non-deformable (ω 0 = ω N V T or ω RIGID ), respectively, and m the mass of the matrix, the poromechanical coupling coefficient can be computed as Eq 3/Sim. Eq 3/Tóth Eq 3/Exp. Eq 4/Tóth Eq 4/Exp.

C(P ) = ω(P ) -ω 0 (P ) ω 0 (P )ε V (P ) (3)

2.75

Figure 9: Poromechanical coupling constant C(P) computed from simulation data using Eq. 3 (squares) with ω 0 = ω N V T . Straight lines are predictions based on Eq. 3 using the fitted isotherms (Fig. 6) and the (dark-red) Tóth or (blue) exponential fits of ε V (P ) (Fig. 8); Dashed lines are predictions based on Eq. 4 using the fitted NPT isotherm and swelling (same color code as for the straight lines); The dotted (green) line correspond to a constant fit at P > 10 bar.

Fig. 9 shows the evolution of the coupling coefficient with pressure obtained using the simulation data obtained within the NPT (ω and ε V ) and NVT (ω 0 ) frameworks. Because of the small values and large error bars for ε V , the predicted data are rather scattered and subject to huge error bars at low pressures. However they tend to converge towards a constant value of about 2.75 when P > 10 bar. Some analytical predictions of C(P ) are also proposed in Fig. 9. A first series of predictions is also based on Eq. 3 and makes use of the fitted isotherms (see Fig. 6) and of the two fitted expressions for ε V (see Fig. 8). As can be seen, the two predicted curves show very similar behaviors and reaches values at large P that are close to the predicted constant. One also observes that, with both ε V functions, C(P ) tends toward large negative values at low P . This result is clearly an artifact of the fitting process. Indeed, while, by definition of the analytical expressions, ε V is strictly positive, it turns out that at low pressures, the fitted isotherms using Tóth model predicts that ω N P T is slightly lower than ω N V T , which is not supported by the data (see the linear model in the inset to Fig. 6). This explains the negative values predicted for C(P ) using Eq. 3.

The second model used for C(P ), also proposed by Brochard et al., [START_REF] Brochard | Poromechanics of microporous media[END_REF] reads:

C(P ) = 1 + K dε V dP (P ) n(P )V b (P ) -ε V (P ) × 1 + K dε V dP (P ) (4)
where K is the bulk modulus, here taken as the lower (viscoelastic) limit (K = 1.2 GPa ± 0.7) and V b is the molar volume of bulk Ar at the considered pressure and temperature. V b was computed from NPT MD simulation of bulk Ar in a large pressure range (see Fig. S7).

A major difference with Eq. 2 & 3 is that in this model, the expression C(P ) does not refer to an undeformable matrix reference and can be directly computed from the NPT calculations and the bulk fluid and matrix properties. As can be seen in Fig. 9, the predictions of C(P ) using Eq. 4 show much less P dependence than those using Eq. 3, even at low P values. Building on all these results, the assumption of a constant value of C = 2.75 seems reasonable.

Fig. 10 is equivalent to Fig. 9 with the only difference that data obtained using the RIGID scheme are used as the underformable reference (ω 0 ) in Eq. 3. Although very similar, some noticeable differences between Fig. 9 and Fig. 10 can be observed. The most significant difference is that when the RIGID reference is used in Eq. 3, C(P ) converges to a slightly larger value, 4.47, instead of 2.75 when the NVT reference or Eq. 4 are used. Eq.3/Sim. Eq 3/Tóth Eq 3/Exp. Eq 4/Tóth Eq 4/Exp. 4.47

Figure 10: Same as Fig. 9 with ω 0 = RIGID in Eq. 3. Data for Eq. 4, identical to those in Fig. 9 are reploted for comparison puropose.

This result demonstrates the importance of accounting for internal (i.e. at constant volume) matrix deformations in addition to macroscopic volume expansion, to accurately describe poromechanics in mature kerogen. We also note that the predicted expressions for C(P ) using Eq. 3 lead to large positive values at low P . This is due to exactly the same artifact that was discussed for Fig. 9, yet with the fitted ω N P T being larger than the fitted ω RIGID .

Fig. 11 shows the evolution with applied pressure of the geometric specific surface area S sp and porous volume V p . At zero pressure, S sp and V p values reads 2034 ± 44 m 2 /g and 0.37 ± 0.01 cm 3 /g, respectively. Upon pressure increase, both S sp and V p increase, as expected regarding the observed swelling of the matrix. Values of 2705 ± 47 m 2 /g and 0.49 ± 0.02 cm 3 /g are obtained at 2116 bar, corresponding to increases of about 33 % for both quantities, with respect to their initial values. While we could expect that swelling would lead to an increase in pore size, and hence to a decrease in the surface over volume ratio of the pores, Fig. 11 shows that it is not the case (see the inset showing an almost constant S sp /V p ). This is also confirmed by looking at the evolution with pressure of the PSD (Fig. S8) or the average pore size Φ , which remains within the error bars of the value at 0 pressure, 7. 0.1 (Fig. S9). This demonstrates that in this specific case, the adsorption induced increase in porosity is achieved via the creation of additional pores of the same size as those present in the matrix at rest.

Regarding the absolute values of V p and S sp , it is interesting to note that values computed for V p are typical of the microporous volumes measured for typical activated carbons, 10,55 yet the specific surface areas are slightly larger (about twice as large) as those usually measured for carbons having this range of V p . This discrepancy between computed and measured S sp may results either from: (i) the small average size of the pores (7 Å) which tend to increase the surface over volume ratio; (ii) some intrinsic error in the definition of the geometric surface using the marching cube algorithm, known to overestimate the area of curved surfaces (up to 10 % for a sphere); the different nature between the computed geometric surface and the measured specific surface which are determined from adsorption isotherms under the Brunauer-Emmett-Teller (BET) theory. [START_REF] Brunauer | Adsorption of Gases in Multimolecular Layers[END_REF] Finally, it is also worth noting that experimental values of S sp are determined from a theory (BET) considering a constant surface during the adsorption process. While, in our case, this may be valid at low pressures, Fig. 11 shows that this approximation is not justified at larger pressures. Figure 12: Evolution of the accessible porosity (ϕ a ), free porosity (ϕ f ) and pore occupancy ratio (χ) with Ar pressure using the hybrid NPT/GCMC (NPT) method. ϕ a at zero pressure is indicated with horizontal lines (thick line: mean; thin lines: standard error).

Fig. 12 shows the evolutions with applied Ar pressure of the accessible porosity ϕ a , free volume ϕ f and pore occupancy χ ratios. While ϕ a increases from about 0.33 at zero pressure to about 0.41 at 2116 bar, ϕ f , identical to ϕ a at zero pressure, by definition, decreases down to about 0.22, due to the increase in the volume occupied by the fluid. As shown by the snapshots in Fig. 13, adsorption essentially takes place on the pore surface up to a few bar of Ar pressure (Fig. 13(b)). Above this pressure, pore filling takes place (Fig. 13(d)) and the pore occupancy ratio χ increases significantly, up to a value of about 0.45. The latter is relatively close to the theoretical maximum packing ratio for disordered homogeneous spheres, 0.634. [START_REF] Song | A phase diagram for jammed matter[END_REF] It is also very close, yet slightly lower than the jammed packing fraction, 0.494, at which hard spheres cannot diffuse anymore. [START_REF] Torquato | Toward the jamming threshold of sphere packings: Tunneled crystals[END_REF] It is also important to add that our geometric description of pore space is only semi-quantitative, as about 50 % of the Ar atoms, whatever the actual value of pressure, are actually slightly overlapping with the van der Waals sphere of a matrix atom, hence, by definition, are not entirely located in the accessible pore volume. Some Ar atoms, see for instance the one at the bottom right in Fig. 13(b) can even be located in an area where no accessible porosity is observed. Fig. 14 shows the evolution of the computed Ar self-diffusivity as a function of applied pressure. As expected due to the small size of the matrix, D s is not isotropic in our system and we overall observe that D z s > D x s > D y s where the superscript refers to the Cartesian direction. Nevertheless, some trends can be distinguished in Fig. Figure 15: Evolution of the self-diffusivity coefficient D s computed using the NPT scheme with the free volume ratio ϕ f . Squares: x direction; circles: y direction; triangles: z triangles. Dashed lines are linear fits to the data.

Relating the anisotropy of D s to the details of the pore space is a complex problem. The determination of the limiting pore diameter (LPD) [START_REF] Vasileiadis | Modeling of Bulk Kerogen Porosity: Methods for Control and Characterization[END_REF][START_REF] Vasileiadis | Transport Properties of Shale Gas in Relation to Kerogen Porosity[END_REF] along percolation paths in the three cartesian coordinates could be helpful in that aim, even though both LPD and percolation paths may be ill-defined for a flexible matrix as the one used here. Instead, recognizing that the self-diffusivity D s (Fig. 14) and the free volume ratio ϕ f (Fig. 12) seem to show similar evolutions with pressure, we show in Fig. 15 the evolution of the three components of D s computed using the NPT scheme (i.e. in unjacketed conditions) with ϕ f . We indeed observe a linear evolution (increase) of D s with ϕ f in all directions. This correlation between diffussion and free volume was also observed in former works [START_REF] Obliger | Methane Diffusion in a Flexible Kerogen Matrix[END_REF][START_REF] Ariskina | Free Volume Model for Transport in Flexible Kerogen of Source Rock's Organic Matter[END_REF] where adsorption induced swelling was significant enough to compensate for the fluid steric hindrance.

Discussion and Conclusion

In this work we have presented a thorough investigation of unjacketed Ar adsorption in a mature kerogen model, as a proxy for a generic microporous carbon with H/C ∼ 0.5, using a coupled GCMC/MD(NPT) strategy. At the chosen temperature, 243 K, the matrix exhibits a viscoelastic behavior, with two distinct timescales for volume relaxation: the purely elastic behavior, associated to phonons, converges within 1 to 10 ns and is followed by viscoelastic deformations characteristic of molecular reorganizations within the matrix. The latter renders extremely difficult the determination of volumetric properties, including density, porosity and adsorption-induced volumetric strains. Here MD relaxation of a few hundred nanoseconds at every investigated pressure was required to obtain reasonable predictions of such properties.

Adsorption isotherms and volumetric strains were computed up to a very large pressure, 2116 bar, in order to accurately characterize the poromechanical behavior of the material and its effects on the absolute adsorbed amounts. Adsorption isotherms were also computed using hybrid GCMC/MD(NVT) and pure GCMC strategies, for comparison. We observed that neglecting macroscopic volume change leads to an underestimation of fluid loading of about 19 % at the maximum pressure. Hindering any matrix reorganization leads to an even stronger underestimation of up to 28 %.

The observed swelling was then rationalized with respect to the microporomechanical model of Brochard et al. [START_REF] Brochard | Poromechanics of microporous media[END_REF] and we have shown that the linearized model n(P ) = n 0 (P ) ×

(1 + Cε(P )) applies with C = 2.75 or 4.47 when n 0 (P ) is computed within the NVT or RIGID frameworks, respectively. However, an important limitation of the microporomechanical model, even in the case where C can be considered constant with respect to P, is that this constant is clearly specific to the adsorbate-adsorbent pair, and may also depend on temperature. This means that this model, although being an excellent tool in interpreting adsorption-induced deformations as done here, is not predictive, as C has to be determined via explicit swelling calculations. The swelling predictions reported for coal [START_REF] Zhang | Molecular simulation of CO 2 -CH 4 competitive adsorption and induced coal swelling[END_REF] and kerogen [START_REF] Huang | Molecular Insights into Kerogen Deformation Induced by CO 2 /CH 4 Sorption: Effect of Maturity and Moisture[END_REF] based on C constants determined for another porous carbon, should thus be reconsidered.

A clear strength of the explicit poromechanical predictions obtained using GCMC/MD simulations is that they allow characterizing the change in pore space accompanying adsorptioninduced deformations. We have indeed observed that the accessible pore volume increases by ∼ 30% during the simulation, which, in a purely elastic matrix would necessarily come with an increase in the average pore size. Our results show that it is not the case for the considered matrix and gas conditions where the PSD is only marginally affected by the swelling, due to viscoelastic deformations.

We have also shown that the self-diffusivity of Ar decreases with increasing pressure (hence increasing loading). We have attributed this decrease in D s with the decrease in free volume ϕ f observed with increasing pressure, despite significant swelling. The positive correlation between self-diffusivity and free volume agrees well with former studies, [START_REF] Obliger | Methane Diffusion in a Flexible Kerogen Matrix[END_REF][START_REF] Collell | Transport of Multicomponent Hydrocarbon Mixtures in Shale Organic Matter by Molecular Simulations[END_REF][START_REF] Obliger | Free Volume Theory of Hydrocarbon Mixture Transport in Nanoporous Materials[END_REF] the linear dependence observed here is typical of relatively small variations of free volume. We have also observed that neglecting swelling and internal matrix flexibility leads to significant underestimations of the self-diffusivity coefficient, yet not as dramatic as those reported in ref. 19 at low values of ϕ f .

Overall, we have shown in this work that adsorption-induced deformations can significantly affect both the thermodynamics and dynamics of the adsorbed fluid. However, only one set of carbon matrix and gas conditions was adopted here due to the long simulation times required for the accurate determination of volumetric properties. Further work is definitively required to to generalize these findings to other carbon matrices, showing different (visco-)elastic behaviors, and gases. In this regard, the problematic of enhanced methane recovery from kerogen, via CO 2 injection and sequestration, is of particular interest.

Supporting information

Figure 1 :

 1 Figure 1: Snapshot of the microporous carbon model equilibrated at 243 K and zero pressure. cyan: C; white: H.

  The self-diffusion coefficient (D s ) is obtained from the slope of the time evolution of the Ar mean square displacement (MSD) over three directions (x, y, z) using MD simulations in the N V T ensemble. For this, five independent configurations from each of the adsorption simulation schemes (NPT, NVT and RIGID) are taken at six different gas pressure values (from 1.44 to 2116 bar) as initial configurations. The simulations are run for 6 ns and the MSD computed using time origins spaced of 25 fs. MSDs were found to reach the linear regime after approximately 100-500 ps. Examples of the MSDs obtained at 22.2 bars with the NPT scheme are shown in Fig. S4.

Figure 2 :

 2 Fig. 2 shows the volume fluctuations along an equilibrium NPT MD of the empty porous carbon at 243 K and zero pressure. Two distinct timescales are clearly visible. The fastest

Figure 3 :

 3 Figure 3: Evolution of the bulk modulus K τ computed with the volume fluctuation method as a function of the averaging time τ . Error bars are given as the standard error (3σ) on the mean obtained from a number of blocks n b = 500 τ (ns) . Thick blue and red lines indicate estimates of the elastic (K e ) and viscoelastic (K ve ) moduli, respectively.

Figure 4 :

 4 Figure 4: (a) Distribution of non-affine displacements u na AB between C atoms of configurations A (80 ns) and B (180 ns). Close-ups of configurations (b) A and (c) B, color-coded according to the values of u na AB . Atoms with large u na AB correspond to fragments located in a low density area. Their reorganization via rotations of the fragments is associated to the densification of the area.

Figure 5 :

 5 Figure 5: Pore size distribution (Ar probe) of the porous carbon at 243 K and zero pressure.

Figure 6 :

 6 Figure 6: Simulated sorption isotherms obtained with the hybrid NPT/GCMC (NPT), hybrid NVT/GCMC (NVT) and pure GCMC (RIGID) schemes. Inset is an enlargement on the low pressure (Henry) regime. Lines in the main panel and inset are Tóth isotherm and linear fits, respectively.

Figure 7 :

 7 Figure7: Evolution of the deviation between the computed virial pressure P v and applied reservoir pressure P with respect to the latter. Squares: hybrid NPT/GCMC simulations (NPT); circles: hybrid NVT/GCMC simulations (NVT).

Figure 8 :

 8 Figure 8: Swelling isotherm computed with the hybrid NPT/GCMC (NPT) method. The average volume of the empty matrix at 243 K and zero pressure is used as reference V 0 . Straight line (dark-red): Tóth isotherm fit; dashed (blue) line: stretched exponential fit; dashed (orange) line: linear fit; dotted (grey) line: ε V = 0.

  11 and Kowalczyk et al. 10. Instead, below 5 bar our results rather suggest a possible contraction (negative strain) of the kerogen, again in agreement with the contraction of ∼ 0.02 % observed in ref. 11 and theorized in ref. 10, even though our results are not fully conclusive due to large error bars. It has to be noted that adsorption-induced contraction of microporous carbons has been observed since the early work of Haines and McIntosh, 53 even though not yet clearly characterized in molecular simulations. This work, and the recent study by Obliger et al., 18 in which signs of the contraction of an overmature kerogen matrix subjected

Figure 11 :

 11 Figure 11: Evolution of the specific surface S sp (circles) area and porous volume V p (squares) with applied Ar pressure. Inset (triangles) is the evolution of the surface over volume ratio.

Figure 13 :

 13 Figure 13: Snapshots of equilibrium configurations at 243 K and P = 0 bar (a), 1.44 bar (b), 22.2 bar (c) and 591 bar (d). Matrix atoms: cyan (C) and white (H) sticks; adsorbed Ar atoms: large pink spheres; surface grid points: small yellow spheres.

Figure 14 :

 14 Figure 14: Evolution of Ar self-diffusivity D s as a of pressure computed using the NPT (squares), NVT (circles) and RIGID (triangles) schemes. One dimensional D s contributions are given along the three Cartesian directions: x (top); y (middle); z (bottom).

Figures

  Figures S1-S9; a relaxed configuration in xyz format of the matrix at 243 K and zero pressure ; and the modified AIREBO force-field file in LAMMPS format.

Table 1 :

 1 Lennard-Jones parameters.

	Atom	C	H	Ar
	σ( Å)	3.4 2.52 3.405
	ε(meV ) 3.41 1.80 10.32

Table 2 :

 2 Adsorption simulation parameters for the NPT, NVT and RIGID schemes. × 10 5 2.5 × 10 5 2.5 × 10 5

		NPT	NVT	RIGID
	MC step/cycle 2.5 MD ensemble NPT	NVT	-
	MD cycle duration 100 ps	100 ps	-
	Equilibration cycles 1000	500	500
	Production cycles	2000	500	500
	most accurate one, noted NPT in what follows, simulates unjacketed adsorption in which

  14. First, despite large error bars, we see that obtained D s values are overall larger with the NPT scheme and lower with the RIGID scheme, the NVT scheme giving intermediate values. This shows that both swelling and internal mobility have a positive effect on diffusion. Also, a clear trend appears in the evolution of D s with pressure, whatever the direction considered: (i) D s is relatively independent of P up to about 50 bar; then (ii) D s significantly decreases with increasing P at larger pressures.

		7
	/s)	5 6	x y z
	D s (10 -9 m 2	2 3 4
		1
		0	0.2	0.22	0.24	0.26	0.28	0.3	0.32	0.34
						ϕ f

Acknowledgements

Access to the computing facilities of the Mésocentre de Calcul Intensif en Aquitaine (MCIA) are gratefully acknowledged.

Graphical TOC Entry