Nicolas Sobecki

Carlos Nieto-Draghi

Angela Di Lella

Yu Didier

Ding

Phase behavior of hydrocarbons in nano-pores

Keywords: Confined fluid, critical shift, vapour-liquid equilibria, molecular simulation, shale gas

Inside nanopores, solid-fluid interactions are of the same order of magnitude as intermolecular interactions of fluid molecules. This fact strongly modifies the thermodynamic properties of confined fluids with respect to bulk phases. Tight oil and shale gas reservoirs, where the proportion of micro (below 2 nm) and mesopores (between 2 and 50 nm) can reach more than 20% of the volume distribution, represent an environment with such problems and industrial challenges to hydrocarbon fluid pressure/volume/temperature (PVT) modeling. This study provides a detailed understanding of the thermodynamic behavior of confined fluid and reference data of the thermodynamic properties of pure components (methane, ethane, n-pentane, n-decane) and mixtures (methane/ethane, ethane/n-pentane) confined in graphite slit pores. Furthermore, a detailed explanation of the different pressures considered in a porous medium with nano-pores is given. The Gibbs Ensemble Monte Carlo (GEMC) NVT simulation is used for pure components instead of the more traditional Grand Canonical Monte Carlo ensemble (GCMC) simulation to get more precise results of liquid and vapor confined pressure avoiding the phase change location determination problem. The evolution of critical temperature and pressure versus pore radius is compared to literature correlations and confined vapor and liquid densities are calculated. A new and robust method in the GEMC ensemble called GEMC NPT Bubble point Monte Carlo (BPMC) completed with GEMC NVT simulations has been developed to get thermodynamic properties including pressures at equilibrium of confined mixtures. Pressure versus density diagrams, pressure versus molar fraction isotherms and examples of pressure versus temperature diagram for a specific composition are built. The phase envelope of the confined fluid is shifted and closes with respect to phase envelope of bulk fluid. The critical temperature and pressure are shifted from the bulk value to a lower value and the bubble point pressure is decreased as the dew point pressure is increased. With regards to the selectivity of the confined system compared to the bulk fluid, for the methane/ethane and ethane/n-pentane mixtures, the heavier component is preferentially adsorbed in the vapor phase and the lighter component is preferentially adsorbed in the liquid phase. All these results for pure components and mixtures provide relevant information concerning the understanding of the phase behavior in confined systems such as shale gas and tight oil reservoirs, emphasizing the difference from the bulk fluid. Furthermore all these data may be used as references for pore radius dependent equation of state (EOS) calibration.

Introduction

Growth in the world economy requires more energy and demand is projected to increase by 30% by 2035 [1]. Fossil fuels will still account for more than three-quarters of world energy consumption through 2040. Among oil and gas production, unconventional resources such as shale gas and tight oil have received much attention in the past decade and become the focus of the petroleum industry for the development of energy resources worldwide. Indeed, according to the U.S. Energy Information Administration [START_REF]Energy Information Administration, International Energy Outlook[END_REF] world tight oil production will more than double from 2015 to 2040 and will represent almost 10% of world oil production. Shale production will account for around 60% of the increase in gas supplies to 2035 according to BP [1] and will represent 30% of world natural gas production in 2040 [START_REF]Energy Information Administration, International Energy Outlook[END_REF].

Email address: carlos.nieto@ifpen.fr (Nicolas Sobecki, Carlos Nieto-Draghi*, Angela Di Lella and Didier Yu Ding) Unlike conventional reservoirs where pore size distribution has a micrometer scale, tight oil and shale gas reservoirs are made up of a very heterogeneous pore size distribution ranging from several nano-meters to micro-meters. The proportions of micro (below 2 nm) and mesopores (between 2 and 50 nm) are generally much smaller than that of the macropores (more than 50 nm) but can reach more than 20% of the volume distribution [START_REF] Kuila | Understanding Pore-Structure And Permeability In Shales[END_REF]. They are mainly associated with clay minerals and kerogen [START_REF] Pommer | Quantitative assessment of pore types and pore size distribution across thermal maturity, Eagle Ford Formation[END_REF]. As hydrocarbon molecules range between 0.5 and 10 nm [START_REF] Nelson | Pore-throat sizes in sandstones, tight sandstones, and shales[END_REF], interaction forces between confined fluid and pore wall molecules become as significant as inter molecular interactions within the confined fluid, which can dramatically change the fluid phase behaviour. Understanding and studying the thermodynamic behavior of such fluids in such confined space become therefore crucial for optimal production forecasts and accurate estimates of fluid in place, which are currently very challenging especially for the production Gas-Oil-Ratio (GOR) prediction [START_REF] Kumar | Upper and Lower Bakken Shale Production Contribution to the Middle Bakken Reservoir[END_REF][START_REF] Nojabaei | Effect of Capillary Pressure on Phase Behavior in Tight Rocks and Shales[END_REF][START_REF] Tian | The Eagle Ford Play, South Texas: Regional Variations in Fluid Types, Hydrocarbon Production and Reservoir Properties[END_REF]. This is also of great interest for chemical, food and pharmaceutical industries where micro and meso-porous mate-rials are widely used for pollution control, mixture separation and catalysts. Several approaches have been used to study the confinement effect in the phase behavior of confined fluids.

Firstly, the most reliable methods to measure fluid properties are still the experimental techniques. However getting all fluid properties in confined media through experiments is very challenging because of many limitations. Some nanofluidic experiments are showing that the bubble point temperature is increased in nano-channels compared to bulk [START_REF] Alfi | Experimental investigation of confinement effect on phase behavior of hexane, heptane and octane using lab-on-a-chip technology[END_REF][START_REF] Wang | Experimental Study and Modeling of the Effect of Nanoconfinement on Hydrocarbon Phase Behavior in Unconventional Reservoirs[END_REF]. Luo et al. [START_REF] Luo | Effect of confinement on the bubble points of hydrocarbons in nanoporous media[END_REF], Cho et al. [START_REF] Cho | Bubble Point Measurements of Hydrocarbon Mixtures in Mesoporous Media[END_REF] conducted the same kind of experiments on pore glasses and mesoporous materials respectively and measured also an increase in bubble point compared to bulk. Al Ismail and Horne [START_REF] Al Ismail | An Investigation of Gas-Condensate Flow in Liquid-Rich Shales[END_REF] conducted core-flooding experiments with gas condensate sample. The core-flooding was conducted through Marcellus shale and conventional Berea sandstone. The results show that the magnitude of variation in the gas composition along the direction of flow during depletion in the Marcellus shale was less than in the Berea sandstone core.

Secondly, a lot of work has been done on extension of equations of state (EOS) for confined fluid behavior modeling. The first main method used includes capillary pressure in the flash calculation considering two phases in equilibrium having two different pressures [START_REF] Alharthy | Multiphase Compositional Modeling in Small-Scale Pores of Unconventional Shale Reservoirs[END_REF][START_REF] Firincioglu | Thermodynamics of Multiphase Flow in Unconventional Liquids-Rich Reservoirs[END_REF][START_REF] Haider | Impact of Capillary Pressure and Critical Property Shift due to Confinement on Hydrocarbon Production in Shale Reservoirs[END_REF][START_REF] Rezaveisi | Compositional Simulation Including Effect of Capillary Pressure on Phase Behavior[END_REF][START_REF] Sandoval | The Phase Envelope of Multicomponent Mixtures in the Presence of a Capillary Pressure Difference[END_REF][START_REF] Stimpson | Effects of Confined Space on Production from Tight Reservoirs[END_REF][START_REF] Teklu | Phase Behavior and Minimum Miscibility Pressure in Nanopores[END_REF][START_REF] Xiong | Effect of Large Capillary Pressure on Fluid Flow and Transport in Stress-sensitive Tight Oil Reservoirs[END_REF][START_REF] Zhang | Capillary Pressure Effect on Hydrocarbon Phase Behavior in Unconventional Reservoirs[END_REF]. Some authors improve this method by considering the thickness of the adsorbed layer [START_REF] Dong | Phase Behavior of Hydrocarbon Mixtures in the Organic Nanopores of Unconventional Gas Condensate Reservoirs[END_REF][START_REF] Liu | Phase Behavior of Fluid Mixtures in a Partially Confined Space[END_REF]. The second main method takes into account a unique pressure for both phases and adds shift of critical temperature and pressure in the flash calculation [START_REF] Alharthy | Multiphase Compositional Modeling in Small-Scale Pores of Unconventional Shale Reservoirs[END_REF][START_REF] Haider | Impact of Capillary Pressure and Critical Property Shift due to Confinement on Hydrocarbon Production in Shale Reservoirs[END_REF][START_REF] Teklu | Phase Behavior and Minimum Miscibility Pressure in Nanopores[END_REF][START_REF] Dong | Phase Behavior of Hydrocarbon Mixtures in the Organic Nanopores of Unconventional Gas Condensate Reservoirs[END_REF][START_REF] Devegowda | Phase Behavior of Gas Condensates in Shales Due to Pore Proximity Effects: Implications for Transport, Reserves and Well Productivity[END_REF][START_REF] Jin | Investigating The Effect of Pore Proximity on Phase Behavior And Fluid Properties in Shale Formations[END_REF][START_REF] Sanaei | Effect of Pore Size Distribution and Connectivity on Phase Behavior and Gas Condensate Production From Unconventional Resources[END_REF]. The correlations used for critical temperature and pressure versus pore radius are either analytical [START_REF] Meyra | Thermodynamic equations for a confined fluid at nanometric scale[END_REF] or built from molecular simulation results [START_REF] Jin | Investigating The Effect of Pore Proximity on Phase Behavior And Fluid Properties in Shale Formations[END_REF][START_REF] Sanaei | Effect of Pore Size Distribution and Connectivity on Phase Behavior and Gas Condensate Production From Unconventional Resources[END_REF]. The two methods of capillary pressure consideration and shift of critical properties are also applied together [START_REF] Teklu | Phase Behavior and Minimum Miscibility Pressure in Nanopores[END_REF][START_REF] Dong | Phase Equilibria of Confined Fluids in Nanopores of Tight and Shale Rocks Considering the Effect of Capillary Pressure and Adsorption Film[END_REF][START_REF] Lopez Jimenez | Effects of Thermodynamic and Rock Properties on the Performance of Liquids-Rich Nano-Porous Shale Reservoirs[END_REF][START_REF] Zuo | Impact of Capillary Pressure and Nanopore Confinement on Phase Behaviors of Shale Gas and Oil[END_REF]. Other extensions of EOS have been proposed by including the pore/fluid interaction effect. However such input parameters must be fitted for each component with either experimental or molecular simulation results. [START_REF] Barbosa | Extending an equation of state to confined fluids with basis on molecular simulations[END_REF][START_REF] Islam | A theory-based simple extension of Peng-Robinson equation of state for nanopore confined fluids[END_REF][START_REF] Travalloni | Phase equilibrium of fluids confined in porous media from an extended Peng-Robinson equation of state[END_REF][START_REF] Travalloni | Critical behavior of pure confined fluids from an extension of the van der Waals equation of state[END_REF].

Finally, molecular simulation appears to be the best way to approach the reality of the thermophysics of confined fluids. Traditionally the study of confined fluid using Monte Carlo molecular simulation is performed by the Grand Canonical Monte Carlo ensemble (GCMC), which considers one box of simulation with constant chemical potential, volume and temperature (µVT). It was first introduced by van Megen and Snook [START_REF] Van Megen | Physical adsorption of gases at high pressure[END_REF] for gas adsorption isotherm in slit like pores and has been widely used for confined pure components studies [START_REF] Balbuena | Theoretical interpretation of adsorption behavior of simple fluids in slit pores[END_REF][START_REF] Binder | Capillary condensation in the lattice gas model: A Monte Carlo study[END_REF][START_REF] Page | Phase equilibrium in a molecular model of a fluid confined in a disordered porous material[END_REF][START_REF] Pitakbunkate | Effect of Confinement on Pressure/Volume/Temperature Properties of Hydrocarbons in Shale Reservoirs[END_REF][START_REF] Porcheron | Monte Carlo simulations of nanoconfined n-decane films[END_REF][START_REF] Radhakrishnan | Effect of the fluid-wall interaction on freezing of confined fluids: Toward the development of a global phase diagram[END_REF][START_REF] De Hemptine | Application to Petroleum Engineering of Statistical Thermodynamics -Based Equations of State[END_REF]. Pitakbunkate et al. [START_REF] Pitakbunkate | Phase Behavior of Methane-Ethane Mixtures in Nanopores[END_REF] have applied GCMC to methane/ethane mixture phase behavior calculation. Jin et al. [START_REF] Jin | Molecular simulation of the pore size distribution effect on phase behavior of methane confined in nanopores[END_REF] have modified the GCMC ensemble to create the gauge-GCMC method and studied pure system, binary and ternary system. The precise determination of the chemical potential corresponding to the condensation still remains a challenge and therefore impacts the precision of the liquid/vapor thermodynamic properties. Peterson and Gubbins [START_REF] Peterson | Phase transitions in a cylindrical pore[END_REF] have proposed a method using integration of the grand free energy of pure compounds, Pitakbunkate et al. [START_REF] Pitakbunkate | Effect of Confinement on Pressure/Volume/Temperature Properties of Hydrocarbons in Shale Reservoirs[END_REF] have identified the phase change by plotting density versus fugacity and have observed the gap in density value. These meth-ods may cause some errors in the location of the phase equilibrium also affecting the estimation of the critical point. The Gibbs Ensemble Monte Carlo (GEMC) developed by Panagiotopoulos [START_REF] Panagiotopoulos | Direct determination of phase coexistence properties of fluids by Monte Carlo simulation in a new ensemble[END_REF] and improved by Panagiotopoulos et al. [START_REF] Panagiotopoulos | Phase equilibria by simulation in the Gibbs ensemble[END_REF] does not have this problem as it considers one vapor box and one liquid box at equilibrium. It has been widely used for confined pure fluid equilibria [START_REF] Gelb | Phase separation in confined systems[END_REF][START_REF] Panagiotopoulos | Adsorption and capillary condensation of fluids in cylindrical pores by Monte Carlo simulation in the Gibbs ensemble[END_REF][START_REF] Rahmani Didar | Confinement Effects on Hydrocarbon Mixture Phase Behavior in Organic Nanopore[END_REF][START_REF] Jiang | Vapour-liquid equilibria in two-dimensional Lennard-Jones fluids: unperturbed and substrate-mediated films[END_REF] but very few studies have been done for mixtures [START_REF] Pathak | Experimental and Molecular Modeling Study of Bubble Points of Hydrocarbon Mixtures in Nanoporous Media[END_REF][START_REF] Li | Vapor-liquid equilibrium properties for confined binary mixtures involving CO2, CH4, and N2 from Gibbs ensemble Monte Carlo simulations[END_REF]. Other authors have used Grand Canonical Transition Matrix Monte Carlo (GCTMMC) [START_REF] Siderius | Extension of the Steele 10-4-3 potential for adsorption calculations in cylindrical, spherical, and other pore geometries[END_REF] and Configurational-biais grand-canonical transition-matrix Monte Carlo simulations [START_REF] Singh | Vapor-Liquid Phase Co-existence, Critical Properties, and Surface Tension of Confined Alkanes[END_REF][START_REF] Singh | Effect of pore morphology on vapor-liquid phase transition and crossover behavior of critical properties from 3D to 2D[END_REF] but only for confined pure component property calculations. Whatever the ensemble used, all authors cited above agree on the fact that with confinement critical temperature and pressure of the fluid are shifted from their bulk value and vapor density increases while liquid density decreases. These observations are independent of pore shape and composition [START_REF] Singh | Effect of pore morphology on vapor-liquid phase transition and crossover behavior of critical properties from 3D to 2D[END_REF][START_REF] Jin | Phase Behavior in Shale Organic/Inorganic Nanopores From Molecular Simulation[END_REF] and pore size distribution [START_REF] Jin | Molecular simulation of the pore size distribution effect on phase behavior of methane confined in nanopores[END_REF]. However, to the best of our knowledge, analyses related to confined mixture properties and liquid and vapor pressures using GEMC ensemble are not available.

In this work, we propose a new and robust method in the GEMC ensemble aiming to get thermodynamic properties of mixtures and confined liquid and vapor pressures at equilibrium. Thermodynamic properties of pure fluids and mixture will be determined in confined spaces, emulating those in kerogen pores in shale reservoirs using GEMC ensemble. The pores will be modelled by slit pores with graphite walls with different sizes. These data may be used as references for radius dependent EOS calibration. This article is organized as follows: section 2 is devoted to the clarification of the different pressures considered in a porous medium with nano-pores. Section 3 is dedicated to the methodology and the workflow used to model the thermodynamic behavior of confined fluids. In section 4 the main results devoted to pure fluids and mixtures are presented. The main conclusions and perspectives are drawn in section 5.

Considerations of the reference pressure in confinement

Saturated bulk pressure and saturated confined pressure are often compared in the literature in studies using EOS modification by critical point shift [START_REF] Alharthy | Multiphase Compositional Modeling in Small-Scale Pores of Unconventional Shale Reservoirs[END_REF][START_REF] Haider | Impact of Capillary Pressure and Critical Property Shift due to Confinement on Hydrocarbon Production in Shale Reservoirs[END_REF][START_REF] Teklu | Phase Behavior and Minimum Miscibility Pressure in Nanopores[END_REF][START_REF] Dong | Phase Behavior of Hydrocarbon Mixtures in the Organic Nanopores of Unconventional Gas Condensate Reservoirs[END_REF][START_REF] Devegowda | Phase Behavior of Gas Condensates in Shales Due to Pore Proximity Effects: Implications for Transport, Reserves and Well Productivity[END_REF][START_REF] Jin | Investigating The Effect of Pore Proximity on Phase Behavior And Fluid Properties in Shale Formations[END_REF][START_REF] Sanaei | Effect of Pore Size Distribution and Connectivity on Phase Behavior and Gas Condensate Production From Unconventional Resources[END_REF] or capillary pressure method [START_REF] Alharthy | Multiphase Compositional Modeling in Small-Scale Pores of Unconventional Shale Reservoirs[END_REF][START_REF] Firincioglu | Thermodynamics of Multiphase Flow in Unconventional Liquids-Rich Reservoirs[END_REF][START_REF] Haider | Impact of Capillary Pressure and Critical Property Shift due to Confinement on Hydrocarbon Production in Shale Reservoirs[END_REF][START_REF] Rezaveisi | Compositional Simulation Including Effect of Capillary Pressure on Phase Behavior[END_REF][START_REF] Sandoval | The Phase Envelope of Multicomponent Mixtures in the Presence of a Capillary Pressure Difference[END_REF][START_REF] Stimpson | Effects of Confined Space on Production from Tight Reservoirs[END_REF][START_REF] Teklu | Phase Behavior and Minimum Miscibility Pressure in Nanopores[END_REF][START_REF] Xiong | Effect of Large Capillary Pressure on Fluid Flow and Transport in Stress-sensitive Tight Oil Reservoirs[END_REF][START_REF] Zhang | Capillary Pressure Effect on Hydrocarbon Phase Behavior in Unconventional Reservoirs[END_REF]. Comparison between classic EOS and modified EOS for confined fluid modeling showed that the bubble point pressure decreases with confinement. Pitakbunkate et al. [START_REF] Pitakbunkate | Effect of Confinement on Pressure/Volume/Temperature Properties of Hydrocarbons in Shale Reservoirs[END_REF][START_REF] Pitakbunkate | Phase Behavior of Methane-Ethane Mixtures in Nanopores[END_REF] built phase diagram of methane-ethane mixture in nanopores and compared it to a bulk fluid. According to the phase diagrams, the bubble point pressure decreases with confinement. All the authors cited above only used a unique pressure as reference to describe the system without giving further details. However the pressure definition for a confined system in a pore network is not trivial and need some explanations. It is important to clarify the definition of the link between pore network topology and the different pressures that can be observed in confinement. To compare the bulk fluid to the confined fluid, we compare in reality two systems: a bulk fluid and a bulk fluid in thermodynamic equilibrium with a confined fluid. The two systems have a bulk pressure but the oil is represented in green and gas is represented in red. a) Large pore connected with pressure between saturation pressure of a single bulk fluid (P bulk sat ) and saturation pressure of bulk fluid connected to confined fluid where first bubble appears in confined fluid (P bulk sat con f ). Gas appears in every 3 pores. b) Large pore connected with pressure below saturation pressure P bulk sat con f , volume of gas is increasing homogeneously in the 3 pores. c) Large pore (1, 2, 3) connected with nanopore 4 with pressure between P bulk sat and P bulk sat . Gas appears in large pores but not in nanopore. d) Large pore (1, 2, 3) connected with nanopore 4 with pressure below P bulk satcon f . Gas volume in large pores is increasing and gas begins to appear in nanopore. thermodynamic properties of the fluid at equilibrium are not the same: the first one has no interactions with the graphite slit when the second does. In the case of a bulk fluid in thermodynamic equilibrium with a confined fluid, three pressures are observed for a given equilibrium state: the bulk pressure of the fluid and the pressuresof the vapor and liquid phases of the confined fluid. Figure 1 illustrates a schematic representation of two different pore networks in tight oil reservoirs. The first system considered is constituted by three large pores completely connected (1,[START_REF]Energy Information Administration, International Energy Outlook[END_REF][START_REF] Kuila | Understanding Pore-Structure And Permeability In Shales[END_REF] in figure 1.a and b; the liquid and vapor pressures are the same in each pore as the saturation pressure, therefore during the depletion gas phase appears at the same time in each pores during depletion. The second system considered is constituted by three large pores (1, 2, 3) and one nanopore (4) in figure 1.c and b. In this system, for bulk pressure between the bulk fluid saturation pressure and the confined fluid saturation pressure (figure 1.c), gas phase appears only in large pores and not in the nanopores. In contrast, for bulk pressure under the confined fluid saturation pressure, gas also appears in the nanopores and three pressures are present in the system (figure 1.d): the bulk pressure and the liquid phase and the vapor phase of confined fluid pressures are all different. In conclusion in order to compare a bulk fluid to a confined fluid two different systems are considered: a bulk fluid and a bulk fluid pressure in thermodynamic equilibrium with a confined fluid. The bulk pressure is considered as reference. As GCMC simulation only considers one pressure, Pitakbunkate et al. [START_REF] Pitakbunkate | Effect of Confinement on Pressure/Volume/Temperature Properties of Hydrocarbons in Shale Reservoirs[END_REF][START_REF] Pitakbunkate | Phase Behavior of Methane-Ethane Mixtures in Nanopores[END_REF] naturally choose to use the pressure corresponding to the chemical potential of the GEMC simulation.

Besides giving thermodynamic properties of confined pure components and mixtures, this paper will give the link between the different pressures observed in a porous medium with nanopores.

Methodology

From now on, GCMC ensemble is used for the calculation of thermodynamic properties of mixtures at equilibrium in confined space by molecular simulation. However the precise determination of the chemical potential of each component of the mixture at liquid/vapor thermodynamic equilibrium remains a challenge and leads to a lack of precision in the results. Furthermore no ensemble is now able to calculate confined liquid and vapor pressure at thermodynamic equilibrium of pure components or mixtures. We therefore propose in this article a new methodology in the GEMC ensemble allowing the calculation of precise thermodynamic properties of mixtures in confined space at equilibrium and confined liquid and vapor pressures. In this section we explain the methodology and the workflow used to model confined fluid thermodynamic behavior and how to get its properties such as density, phase composition and pressures. We give the details of the different systems studied, the statistical ensemble used, the molecular simulation parameters and the data post-processing. The Gibbs code from IFPEN and the Laboratoire de Chimie Physique (LCP) at University Paris-Sud has been used for all the molecular simulations cited in this paper [START_REF] Ungerer | Applications of molecular simulation in the oil and gas industry: Monte Carlo methods[END_REF].

Case studies

The case study is hydrocarbon fluid confined in a nanometric pore, such as kerogen pores present in tight oil and shale gas reservoirs. The simplified model is a slit pore with graphite walls. Two infinite parallel slices are in the orthogonal directions to the slit pore (i.e. x and y) and the slit pore has a width of length H in the z direction. In order to first validate the intermolecular potential models with experimental or analytical results, all simulations were first performed for bulk fluids. The chosen pure components are methane, ethane, n-pentane and n-decane and the chosen mixtures are methane/ethane and ethane/n-pentane. Two different workflows are used for pure components and mixtures. These workflows are presented in the schematics figure 2 and they are detailed in section 3.2 and 3.3.

Pure components workflow

As mentioned earlier, confined fluid studies using Monte Carlo molecular simulation are usually performed by the Grand Canonical Monte Carlo ensemble (GCMC) [START_REF] Balbuena | Theoretical interpretation of adsorption behavior of simple fluids in slit pores[END_REF][START_REF] Binder | Capillary condensation in the lattice gas model: A Monte Carlo study[END_REF][START_REF] Page | Phase equilibrium in a molecular model of a fluid confined in a disordered porous material[END_REF][START_REF] Pitakbunkate | Effect of Confinement on Pressure/Volume/Temperature Properties of Hydrocarbons in Shale Reservoirs[END_REF][START_REF] Porcheron | Monte Carlo simulations of nanoconfined n-decane films[END_REF][START_REF] Radhakrishnan | Effect of the fluid-wall interaction on freezing of confined fluids: Toward the development of a global phase diagram[END_REF][START_REF] Pitakbunkate | Phase Behavior of Methane-Ethane Mixtures in Nanopores[END_REF][START_REF] Jin | Molecular simulation of the pore size distribution effect on phase behavior of methane confined in nanopores[END_REF]. However this method has several drawbacks. Firstly, the identification of phase change chemical potential remains challenging and complex [START_REF] Peterson | Phase transitions in a cylindrical pore[END_REF]. Secondly, the pressure of vapor and liquid inside the pore at thermodynamic equilibrium cannot be obtained easily. The pressure inside a simulation box is estimated by the virial equation (equation 15) which corresponds to the average of the diagonal terms of the pressure tensor. In addition to the lack of precision of the transition phase The results will therefore not be usable. That is why we propose to use the Gibbs Ensemble Monte Carlo (GEMC) to study such fluids using the total number of particles (N), the total volume (V) and the temperature (T) as constants. Two boxes of simulations are in thermodynamical equilibrium, one for the liquid phase and one for the vapor phase. Different trial moves are randomly proposed and accepted according to their appropriate probability. The Monte Carlo moves considered in this study are translation and rotation of particles inside a box, transfer of a randomly selected particle from one box to the other and the volume change in such way that the total volume remains constant. In our case of confinement in a slit pore, modifications have been done to allow volume changes only in the directions orthogonal to the slit pore (i.e. x and y). It is important to mention that solid-fluid interactions are modeled by an analytical intermolecular potential which is function of z (see details in 3.4). This means that solids are considered as a continuum media, therefore, it is possible to modify the volume of the pores in x and y directions without any problems. For longer chain molecules such as n-pentane and n-decane, internal translation and rotation moves are added. GEMC NVT simulation allows us to get accurate readings of the liquid/vapor equilibria properties but it needs good initial estimation of densities of each phase to converge. The confined GCMC µVT simulations will therefore be used to get initial vapor and liquid densities for pure components. These densities will then be used to calculate the number of particles in each phase in order to initialize the vapor and liquid boxes of the GEMC NVT simulation for pure components. µVT simulation is performed at constant volume and temperature for several chemical potentials. Therefore the fluid is described from ideal gas to compressed liquid. The gap of density value in function of the fugacity allows us to determine roughly the liquid/vapor equilibrium. The initial values of the input fugacities come from an EOS calculation of the bulk phase. All these simulations have been performed for several isotherms and pore widths. More details about GCMC and GEMC simulations can be found in a reference book [START_REF] Frenkel | Understanding Molecular Simulation: From Algorithms to Applications[END_REF].

Mixtures workflow

As for pure components, the application of GEMC NVT simulation for mixtures needs a good initial estimation of densities and molar fraction of each component in each phase to converge. Because of the importance of convergence issues and the complex use of GCMC ensemble for mixtures, a new ensemble called GEMC NPT BPMC have been used in this study in order to initialize the Monte Carlo simulations for mixtures. The description of the classic GEMC NPT version can be found in Panagiotopoulos et al. [START_REF] Panagiotopoulos | Phase equilibria by simulation in the Gibbs ensemble[END_REF]. The Gibbs ensemble GEMC-NPT-BPMC has a constant number of particles, temperature and pressure with a bubble point Monte Carlo (BPMC) pseudoensemble suggested by Ungerer et al. [START_REF] Ungerer | Direct calculation of bubble points by Monte Carlo simulation[END_REF] and improved by Ferrando et al. [START_REF] Ferrando | Ethanoled gasoline bubble pressure determination: Experimental and Monte Carlo modeling[END_REF] for liquid mixture bubble pressure calculation. In Ferrando et al. [START_REF] Ferrando | Ethanoled gasoline bubble pressure determination: Experimental and Monte Carlo modeling[END_REF] pseudo ensemble, two different boxes of simulation for vapor and liquid phases are used. Only the number of particles of the liquid phase is imposed, which allows us to set the liquid composition and chemical potential. In addition the total volume and the temperature are imposed. In order to guarantee the equality between liquid and vapor chemical potential, insertion and deletion of particles are performed exclusively in the vapor phase. After the convergence of the simulation, saturated vapor composition and density are known and the bubble point can be calculated. The partition function of this pseudo ensemble and the probability density of a configuration can be found in Ferrando [START_REF] Ferrando | Potentiels intermoléculaires et algorithmes de monte carlo: Application à l'étude des composés oxygénés[END_REF]. The acceptance probability of the insertion and deletion moves are equivalent to the acceptance probability of the transfer move in the standard Gibbs ensemble. We have adapted this method to confined fluid mixture thermodynamic equilibrium property calculation. As illustrated in figure 3, two boxes are in thermodynamic equilibrium, one box for the bulk saturated liquid fluid and one box for the confined fluid. The pressure in the non-confined box is constant as is the system's temperature. The BPMC is characterized by a Monte Carlo move consisting of keeping constant the number of particles of the bulk saturated liquid mixture and inserting and deleting particles in the confined fluid box. In addition, the volume of the bulk box changes whereas the volume of the confined fluid box remains constant. Translation and rotation of particles inside the boxes are performed according to the acceptation probability. The bulk liquid box is initialized thanks to a classic GEMC NPT simulation for a bulk fluid and the confined fluid box is initially empty of molecules. As represented in figure 4 for an ethane/n-pentane mixture, the number of particles of the different components in the confined box will fluctuate between the liquid and the gas phase during the simulation. Average values of the number of particles in each phase can be obtained by plotting histograms of the number of particles for each component which can be fitted to a sum of two Gaussian curves (equation 1). α 1 i , α 2 i , β 1 i and β 2 i are the constants to fit for one component i and N i , N v i and N l i correspond to the number of particles of component i, the average number of particles in the vapor and the liquid phase for a component i at equilibrium respectively. These number of particles in each phase can be used as an initialization for a confined GEMC NVT simulation for accurate results on thermodynamic properties of the confined fluid and vapor and liquid pressures determination.

f (N i ) = α 1 i e - (N i -N v i ) 2 2α 2 i 2 + β 1 i e - (N i -N l i ) 2 2β 2 i 2 (1)
The summary of all the workflows and statistical ensembles used for molecular simulation in this paper are given in figure 2. Concerning the pure components, GEMC NVT simulations have been performed for the bulk fluid and confined GCMC µVT has been used as an initialization for the confined GEMC NVT simluations. Pure components simulations will be performed for several isotherms and pore widths. Concerning the mixtures, GEMC NPT simulations will be done for the bulk fluid. The resulting liquid properties will be used to initialize the bulk liquid box of the GEMC NPT BPMC simulations, finally confined liquid and vapor properties will be used for the confined GEMC NVT simulations. These simulations will be performed for different isotherms and for a pore width of 3 nm.

Force field and simulation parameters

The main forces taken into account in the system are the interaction potentials between particles. There are three interaction types: fluid-fluid interaction, solid-fluid interaction and solid-solid interaction. The Steele potential [START_REF] Gelb | Phase separation in confined systems[END_REF][START_REF] Steele | The physical interaction of gases with crystalline solids: I. Gas-solid energies and properties of isolated adsorbed atoms[END_REF] is used to model solid-fluid interactions. This potential considers only the first layer of the graphite pore wall; the remaining layers are considered as continuum solid [START_REF] Siderius | Extension of the Steele 10-4-3 potential for adsorption calculations in cylindrical, spherical, and other pore geometries[END_REF], and solid-solid interaction is fully neglected.

The Steele 9-3 potential is written as (equation 2):

u s f (z) = 2 3 πρ s ε s f σ s f 3 2 15
σ s f z 9 - σ s f z 3 ( 2
)
where z corresponds to the distance between solid and fluid particles, ρ s is the atomic density of solid, σ s f represents the distance between two atoms and the surface where attractive and repulsive forces are canceled, ε s f is an energy and represents the depth of the potential well of fluid-solid interactions at the minimum of the function. The values of these parameters can be found in Porcheron et al. [START_REF] Porcheron | Monte Carlo simulations of nanoconfined n-decane films[END_REF].

The fluid-fluid interaction is modeled by the Lennard-Jones 12-6 potential (equation 3):

U i j LJ (r i j ) = 4ε i j σ i j r i j 12 - σ i j r i j 6 ( 3
)
where r i j is the distance between two particles i and j, ε i j is the Lennard-Jones well depth and σ i j is the van der Waals radius. Lorentz-Berthelot combining rules are used for solid/liquid and liquid/liquid interactions (equation 4):

σ i j = σ ii +σ j j 2 ε i j = √ ε ii ε j j (4)
Periodic boundary conditions are applied in all directions for the bulk system and in x and y directions for the confined fluid. In order to increase the calculation performance, the Lennard-Jones interactions are only calculated in a sphere of radius r c (cut-off radius), which is generally the half length of the simulation box. Standard long-range correction for the energy and the pressure is applied for bulk fluid calculations only.

Long-range correction is not used for the slit pore as there are no periodic boundary conditions in the slit width direction. All studied molecules are described using the AUA4 (Anisotropic United Atoms model) optimized parameters [START_REF] Bourasseau | Monte Carlo simulation of branched alkanes and long chain n -alkanes with anisotropic united atoms intermolecular potential[END_REF][START_REF] Ungerer | Optimization of the anisotropic united atoms intermolecular potential for n -alkanes[END_REF]. The CH 2 and CH 3 molecules are represented with a single center of force located near the geometric center of the atoms of each molecule. The AUA model consists of a displacement of the Lennard Jones centers of force toward the hydrogen atoms. The magnitude of the shift between the carbon center and the interaction site is the adjustable parameter δ .

All bond lengths are kept fixed. For long-chain n-alkanes, intramolecular interactions are considered by means of additional energy terms including bending and torsion terms. In addition, LJ interactions are applied between atoms of the same molecule separated by four bonds. Atoms separated by two bonds interact via a harmonic bending potential (equation 5):

U bend k B = 1 2 k bend (cosθ -cosθ 0 ) 2 (5)
where k B is the Boltzmann constant, k bend is the bending constant, and θ and θ 0 are the bending angle and the equilibrium bending angle, respectively. Atoms separated by three bonds are governed by a torsional potential of the following form (equation 6):

U tors k B = 8 ∑ n=0 a n (cosφ ) n ( 6
)
where φ is the torsional angle and a n are constants. The non-flexible and flexible molecule parameters describing the force field are given in [START_REF] Porcheron | Monte Carlo simulations of nanoconfined n-decane films[END_REF][START_REF] Bourasseau | Monte Carlo simulation of branched alkanes and long chain n -alkanes with anisotropic united atoms intermolecular potential[END_REF][START_REF] Ungerer | Optimization of the anisotropic united atoms intermolecular potential for n -alkanes[END_REF] and summarized in tables 1 and 2. The simulation parameters are summarized in table 3 where L x , L y , L z corresponds to the box lengths in the x, y, z directions respectively. The L x and L y lengths were adapted in function of the critical point proximity for all simulations ensembles. The different move probabilities and Monte Carlo steps used for non-flexible and flexible molecules are summarized in tables 4, 5, and 6 of Appendix 6.1. The simulation data post-processing are given in Appendix 6.2. Bond length

r 0 (nm) C-C 0.1535 Bend CH 2 θ ( • ) k bend (K) CH 3 -CH 2 -CH 2 120.15 3.6072 CH 3 -CH 2 -CH 2 Torsion a i (K) CH 3 -CH 2 -CH 2 -CH 2 a 0 = 1001.35 a 5 = 1965.93 CH 2 -CH 2 -CH 2 -CH 2 a 1 =

Results

This section will present the results of the molecular simulation workflows described in 3.2 and 3.3 with simulation parameters and data post-processing detailed in [START_REF] Kumar | Upper and Lower Bakken Shale Production Contribution to the Middle Bakken Reservoir[END_REF]

Pure fluids

We start the results section by showing the effect of confinement on phase density of pure fluid. Figure 5 shows the density phase diagram of C 10 H 22 at bulk condition and with pore confinement. The results for all the pure components studied are in Appendix 6.3, figure 16. Bulk simulations well match the reference data. The confined fluid vapor density increases and the confined fluid liquid density decreases with confinement for all the studied molecules. Fluid/pore interaction attracts particles near the pore wall and creates an adsorption layer. Consequently the vapor density, which is the average density inside the entire pore, will be larger than the bulk density because of this adsorption layer. Pore walls participate in particle cohesion: close to the walls, the molecules are highly structured imposing a translational order from the surface wall to the pore. This imposed order creates a strong layering slightly increasing the inter-particle distance in z direction, which is larger than the one observed in disordered bulk phases. Consequently confined liquid density is less than the bulk one [START_REF] Kuchin | Influence of the Disjoining Pressure on the Equilibrium Interfacial Profile in Transition Zone Between a Thin Film and a Capillary Meniscus[END_REF]. Another conse- The pressure is estimated by the virial equation (equation 15), which takes into account the diagonal elements of the pressure tensor. The anisotropic effect induced by the presence of the walls are therefore take into consideration in our simulations. The critical pressure is calculated using confined vapor pressure and equation 19. The evolution of critical temperature versus pore length differs from the evolution of critical pressure (see figure 6 and Appendix 6.3, figure 17). The curve is different for each component and naturally tends towards bulk value at high pore length. It is worth remembering that pressure calculation in confined pores does not include long-range correction of pressure and energy as it is usually done for bulk calculations. That is why critical pressure for large pores slightly differs from the bulk value.

The comparison with data from the literature (figure 6) shows that results for critical temperature are close to those from Singh and Singh [START_REF] Singh | Effect of pore morphology on vapor-liquid phase transition and crossover behavior of critical properties from 3D to 2D[END_REF], Pitakbunkate et al. [START_REF] Pitakbunkate | Effect of Confinement on Pressure/Volume/Temperature Properties of Hydrocarbons in Shale Reservoirs[END_REF]. Several correlations are used in the literature to describe the evolution of critical temperature and pressure with confinement [START_REF] Meyra | Thermodynamic equations for a confined fluid at nanometric scale[END_REF][START_REF] Jin | Investigating The Effect of Pore Proximity on Phase Behavior And Fluid Properties in Shale Formations[END_REF]. Concerning the critical temperature, the correlation from Jin et al. [START_REF] Jin | Investigating The Effect of Pore Proximity on Phase Behavior And Fluid Properties in Shale Formations[END_REF] is well adapted and matches the observed results. Our results for the critical pressure are not as close to the literature data [START_REF] Singh | Effect of pore morphology on vapor-liquid phase transition and crossover behavior of critical properties from 3D to 2D[END_REF][START_REF] Pitakbunkate | Effect of Confinement on Pressure/Volume/Temperature Properties of Hydrocarbons in Shale Reservoirs[END_REF] as those from critical temperature because the methods used by the previous authors are different. Pitakbunkate et al. [START_REF] Pitakbunkate | Effect of Confinement on Pressure/Volume/Temperature Properties of Hydrocarbons in Shale Reservoirs[END_REF] performed GCMC simulations and used the pressure of the bulk fluid in equilibrium with the confined fluid and monitored the end of density discontinuity versus pressure, which is very challenging. Singh and Singh [START_REF] Singh | Effect of pore morphology on vapor-liquid phase transition and crossover behavior of critical properties from 3D to 2D[END_REF] ran simulations in the Grand Canonical Transition Matrix Monte Carlo (GCTMMC) ensemble to get saturation pressure for different temperatures. Then the critical pressure was obtained by fitting equation [START_REF] Stimpson | Effects of Confined Space on Production from Tight Reservoirs[END_REF]. That is probably one of the reasons the difference between our results and those from the literature. Concerning the correlation for those critical pressure, the analytical solution from Meyra et al. [START_REF] Meyra | Thermodynamic equations for a confined fluid at nanometric scale[END_REF] shows a better match with molecular simulation results than the one of Jin et al. [START_REF] Jin | Investigating The Effect of Pore Proximity on Phase Behavior And Fluid Properties in Shale Formations[END_REF], which leads to negative values of P c con f P c bulk when applied to longer alkane chains inside small pores.

Unlike the bulk fluid where vapor and liquid pressure are equal at equilibrium, confined fluids have different values of vapor and liquid pressures. The pressure can no longer be considered as a reference for thermodynamic equilibrium but chemical potential or fugacity must be used instead. Comparing bulk saturation pressure with confined vapor and liquid pressure in a P-T phase diagram, for example, is of no practical interest, as the bulk and the confined fluids are not at the same thermodynamic equilibrium state. An alternative of using chemical potential as reference is to the use the bulk pressure of both systems: bulk fluid and bulk fluid in thermodynamic equilibrium with a confined fluid (see section 2). An example of the difference in pressure between the vapor and liquid phase for different pore confinement for the n-pentane is shown in figure 7. The pressures are calculated using the Virial in the two boxes of the NVT simulation. The vapor pressure is positive for all pore widths and increases with confinement at a constant temperature as can be seen in figure 7. As the gas molecules [START_REF] Pitakbunkate | Effect of Confinement on Pressure/Volume/Temperature Properties of Hydrocarbons in Shale Reservoirs[END_REF]. The '+' symbols correspond to the results of Singh and Singh [START_REF] Singh | Effect of pore morphology on vapor-liquid phase transition and crossover behavior of critical properties from 3D to 2D[END_REF]. Finally the line correspond to the correlations of Jin et al. [START_REF] Jin | Investigating The Effect of Pore Proximity on Phase Behavior And Fluid Properties in Shale Formations[END_REF] for critical temperature and Meyra et al. [START_REF] Meyra | Thermodynamic equations for a confined fluid at nanometric scale[END_REF] for critical pressure.

become closer due to the smaller slit pore length constraints, the interaction between molecules is enhanced and the Virial pressure increases. The liquid pressure increases also with confinement, but its value goes from negative to positive. The slit pore is by assumption infinitely rigid, therefore the two graphite sheets cannot deform themselves under the action of capillary pressure or van der Waals solid-fluid-solid interactions as in reality [START_REF] Schoen | Pore-lattice deformations in ordered mesoporous matrices: Experimental studies and theoretical analysis[END_REF]. For bigger pores (5 nm and 6 nm) the liquid tends as to densify itself but this is prevented by the pore constraints.

That is why the virial pressure is negative. For smaller pores, the molecules are much more closer, so repulsions could occur and the sum of the total forces gives a positive virial pressure.

The mechanical constraint of the slit length and the Steele potential between wall and fluid causes inhomogeneity of pressure inside the fluid. Disjoining pressure occurs in the perpendicular direction of the slit surface in function of the length of the adsorbed layer [START_REF] Kuchin | Influence of the Disjoining Pressure on the Equilibrium Interfacial Profile in Transition Zone Between a Thin Film and a Capillary Meniscus[END_REF][START_REF] Benet | Disjoining Pressure, Healing Distance, and Film Height Dependent Surface Tension of Thin Wetting Films[END_REF]. We have also analyzed the behavior of the capillary pressure as can be seen in figure 7. For the bigger pores (5 nm and 6 nm) the capillary pressure is positive and decreases with temperature until reaching the critical point where the capillary pressure must be zero. The positivity of the capillary pressure shows the wettability of the surface to the liquid phase, which is common in conventional reservoirs between oil and gas. The values of the capillary pressure for the 5 nm slit pore are higher than the ones obtained for the 6 nm slit pore, which is consistent with the Laplace equation. However the behavior is totally different for the smaller pores (2 nm and 3 nm) where the capillary pressure is negative and increases with temperature towards zero at the critical point. At such small scales the Laplace-Young equation is no longer applicable [START_REF] Benet | Disjoining Pressure, Healing Distance, and Film Height Dependent Surface Tension of Thin Wetting Films[END_REF][START_REF] Churaev | The shape of the transition zone between a thin film and bulk liquid and the line tension[END_REF].

The confined NVT simulation has allowed us to precisely calculate thermodynamic properties of several pure components. The calculation of critical temperature and pressure evolution for several pore sizes has provided reference data to validate the more convenient correlations proposed in the literature. Furthermore confined vapor and liquid pressure has been calcu-lated for one pure components which is a total novelty.

Mixtures

Liquid/vapor thermodynamic equilibrium calculations are performed for two mixtures: methane/ethane and ethane/npentane using the workflow described in figure 2 in section 3.1. A pore length of 3 nm is used for both mixtures.

Methane-Ethane

The mixtures workflow of figure 2 has been performed for the mixture CH 4 -C 2 H 6 at five isotherms: 200 K, 220 K, 230 K, 240 K, 260 K. Post-processing explained in Appendix 6.2 allowed us to obtain the different equilibrium thermodynamic properties of confined fluids.

In order to understand why a confined NVT simulation is needed after a GEMC NPT BPMC simulation, an example of results obtained for a specific isotherm for these two simulations is shown in Appendix 6.3 figure 18. The GEMC NPT BPMC and confined NVT results are quite similar until they get close to the critical point where we observe strong fluctuations of the number of particles in the GEMC NPT BPMC simulations. Indeed, close to the critical point, the two modes of the number of particle histogram are very difficult to detect (see figure 4C). Furthermore the statistical uncertainty of the GEMC NPT BPMC results (around 10 molecules) is much higher than the NVT results (around one molecule) at high temperatures close to the critical point. Finally GEMC NPT BPMC can not give confined liquid and vapor pressures. That is why we have adopted the strategy described in figure 2, where the GEMC NPT BPMC results are used to initialize the confined Gibbs NVT simulations to get better results.

As observed for pure components, the vapor density increases and the liquid density decreases with confinement for all isotherms (see Appendix 6.3 figure 19). The critical pressure decreases with a value below the bulk value for all isotherms (see figure 8). for different isotherms obtained in Appendix 6.3 figure 20, a pressure-temperature diagram can be built for a specific concentration of methane/ethane. Indeed, one value of ethane molar fraction corresponds to one value of dew-point pressure and bubble-point pressure for each isotherm. An example of this diagram for a mixture of 34.9% methane and 65.1% ethane is shown in figure 9. It is observed that the phase envelope of the confined fluid has shifted and closed itself from its bulk value. The critical temperature and pressure have shifted from the bulk value to a lower value and the bubble-point pressure decreases as the dew-point pressure increases. This observation is valid regardless of the proportion of methane and ethane in the mixture.

In the pressure-molar fraction diagrams for the confined fluid (see figure 8), the bulk pressure corresponds to the pressure in the pores 1, 2 and 3 in the schematic figure 1, when the first gas bubble appears in the nanopore 4. The molar fraction of ethane presented in the pressure-molar fraction diagram (see figure 8) corresponds to the molar fraction of ethane in the fluid inside the nanopore 4.

The critical pressures and temperatures for the different mixtures of methane/ethane are shown in figure 10. Regardless of the mixture, the confined fluid always has a smaller value of the critical pressure and critical temperature than the bulk. The results of GEMC NVT simulations have been compared to critical-point calculations for methane/ethane mixtures described by Peng-Robinson EOS using the PVTFlow TM software (IFPEN-Beicip-KAPPA partnership [72]). Critical pressure and temperature of methane and ethane for a 3 nm slit width calculated in section 4.1 have been used to model numerically, from EOS, the mixture critical point of the confined fluid. We identify a critical point whenever a point of phase envelope has an equilibrium constant equal to unity. The bulk values of critical pressure and temperature between GEMC NVT simulation and EOS calculation are quite close as can be seen in figure 10. Concerning the confined fluid, the numerical results obtained from EOS show the same trend as the results obtained from molecular simulation. The method used by several authors [START_REF] Alharthy | Multiphase Compositional Modeling in Small-Scale Pores of Unconventional Shale Reservoirs[END_REF][START_REF] Haider | Impact of Capillary Pressure and Critical Property Shift due to Confinement on Hydrocarbon Production in Shale Reservoirs[END_REF][START_REF] Teklu | Phase Behavior and Minimum Miscibility Pressure in Nanopores[END_REF][START_REF] Dong | Phase Behavior of Hydrocarbon Mixtures in the Organic Nanopores of Unconventional Gas Condensate Reservoirs[END_REF][START_REF] Devegowda | Phase Behavior of Gas Condensates in Shales Due to Pore Proximity Effects: Implications for Transport, Reserves and Well Productivity[END_REF][START_REF] Jin | Investigating The Effect of Pore Proximity on Phase Behavior And Fluid Properties in Shale Formations[END_REF][START_REF] Sanaei | Effect of Pore Size Distribution and Connectivity on Phase Behavior and Gas Condensate Production From Unconventional Resources[END_REF] to model the confined fluid by adding shift of critical pressure and temperature of pure components in the flash calculation is therefore consistent with molecular simulation results for mixtures. Unlike the bulk fluid, the spacial distribution of ethane and methane molecules inside the silt pore is not homogeneous due to the solid-fluid interactions. An example of the density profile in z direction for methane/ethane mixture at 220 K and 2 MPa is shown in figure 11. The initial composition is made of 39.5% methane molar fraction and 60.5% ethane molar fraction, which gives at equilibrium for 220 K and 20 MPa the liquid molar fractions x CH 4 = 0.265, x C 2 H 6 = 0.735 and the vapor molar fractions y CH 4 = 0.642, y C 2 H 6 = 0.358. The density profile corresponds to the number of particles of methane and ethane per nm 3 in each slices of 0.003 nm width along z direction for the slit pore of 3 nm width. The simulation parameters and post-processing details are given in Appendix 6.2. For both vapor and liquid the density of molecules is higher close to the wall than at the center of the pore. The selectivity of the confined system towards methane or ethane in comparison with the bulk one can be calculated with equation 7.

S C 2 H 6 /CH 4 = x con f ined,C 2 H 6 /x con f ined,CH 4 x bulk,C 2 H 6 /x bulk,CH 4 ( 7
)
This selectivity is calculated in the adsorbed layers of the vapor and the liquid phase. As expected we obtain a symmetric profile where particles accumulate close to the wall surface. The adsorbed layer is defined as the first minima or the last value before the plateau of the density profile. The layer length is 0.6 nm and 0.5 nm for the vapor and liquid phase respectively. Our results give a selectivity of S C 2 H 6 /CH 4 =1.17 for the vapor phase and S C 2 H 6 /CH 4 =0.81 for the liquid phase. It means that compared to bulk, ethane is preferentially adsorbed in the vapor phase, whereas methane is preferentially adsorbed in the liquid phase of the confined system.

The confined NVT simulation initialized by the NPT BPMC ensemble has allowed us to calculate thermodynamic properties of confined methane/ethane mixture and build a phase diagram with the bulk pressure as reference for the thermodynamic equilibrium. In order to validate the method, the same workflow will be applied to another mixture and confined pressures will be calculated.

Ethane/n-Pentane

We applied the same approach as in section 4.2.1 to simulate the behavior of the mixture C 2 H 6 -C 5 H 12 for five isotherms: 320 K, 330 K, 340 K, 350 K, 360 K and 370 K. The results of the pressure versus density of the mixture are shown in Appendix 6.3 figure 21. The details of the simulation parameters and data post treatment are provided in section 6 and as explained previously in section 4.2.1, all the results presented come from the NVT simulations.

Similar to the previous case, the vapor density increases and the liquid density decreases with confinement for all isotherms with respect to the non-confined bulk phase (see Appendix 6.3 figure 21). The critical pressure decreases with a value below the bulk value for all isotherms (see figure 12). The pressuremolar fraction diagrams for different isotherms shown in figure 12 allow the construction of a pressure-temperature diagram for a specific ethane/n-pentane concentration (59.7% ethane and 40.3% n-pentane) as can be seen in figure 13. Regardless of the proportion of ethane and n-pentane in the mixture, the phase envelope of the confined fluid is shifted inwards and closes itself from its bulk value. The critical temperature and pressure are shifted from the bulk value to a lower value and the bubblepoint pressure decreases as the dew-point pressure increases.

Critical temperature and pressure calculated from GEMC NVT simulation results have been compared to critical-point calculations for ethane/n-pentane mixtures described by Peng-Robinson EOS using the PVTFlow TM software in Appendix 6.3 figure 22. The critical pressure and temperature of ethane and n-pentane obtained in section 4.1 from our confined simulations in a slit pore of 3 nm width have been used as an input for an EOS to model numerically the mixture critical point of the confined fluid mixture. The bulk values of critical temperature and pressure obtained from GEMC NVT results and EOS calculation are very close. Concerning the results for the confined fluid, we observe a better agreement of both approaches for the critical temperature than the critical pressure.

An example of the density profile in the z direction for the ethane/n-pentane mixture at 340 K and 2 MPa is shown in figure 14. The initial composition is made up of 50.4% ethane molar fraction and 49.6% n-pentane molar fraction, which gives at equilibrium the liquid molar fractions x C 2 H 6 = 0.351, x C 5 H 12 = 0.649 and the vapor molar fractions y C 2 H 6 = 0.758, y C 5 H 12 = 0.242. The vapor and liquid molar fraction values are coherent with figure 14 where there is more methane in the vapor phase and more ethane in the liquid phase. Details of parameters and post-processing are given in Appendix 6. Here again, the densities of molecules in the vapor and liquid phases are higher close to the walls due to fluid/wall interactions. In this case, the selectivity of the adsorbed layers of the confined system towards ethane and n-pentane compared to bulk can be calculated using equation 8. The adsorbed layer length is 0.7 nm and 0.6 nm for the vapor phase and the liquid phase respectively. The results give a selectivity of S Figure 10: Evolution of critical pressure (left) and temperature (right) versus ethane molar fraction for bulk and confined fluid. The blue line corresponds to numerical results obtained from Peng Robinson EOS for the bulk fluid. The red line corresponds to numerical results obtained from Peng Robinson EOS with critical temperature and pressure shift from Jin et al. [START_REF] Jin | Investigating The Effect of Pore Proximity on Phase Behavior And Fluid Properties in Shale Formations[END_REF] and Meyra et al. [START_REF] Meyra | Thermodynamic equations for a confined fluid at nanometric scale[END_REF] correlations respectively. The blue '+' symbols correspond to critical point NPT results for bulk and the red ones correspond to critical point NVT results for confined fluid in H=3 nm pore length. that compared to bulk, n-pentane is preferentially adsorbed in the vapor phase of the confined system, whereas ethane is preferentially adsorbed in the liquid phase.

S C 5 H 12 /C 2 H 6 = x con f ined,C 5 H 12 /x con f ined,C 2 H 6 x bulk,C 5 H 12 /x bulk,C 2 H 6 (8)
A pressure-molar fraction diagram of ethane/n-pentane mixture at 320K for bulk and confined fluid is given in figure 15. The green curve corresponds to the pressure of the bulk fluid (P bulk sat ) versus n-pentane molar fraction, the blue curves give the pressures of a bulk fluid connected to a confined fluid when the first bubble appears in the confined fluid (P bulk sat con f ) and the red and the orange curves give respectively the vapor (P con f v ) and liquid pressure (P con f l ) of the confined fluid connected to a bulk fluid as explained in figure 1. All these pressures are different but only P bulk sat con f , P con f v and P con f l are at the same thermodynamic equilibrium state, i.e. the same chemical potential (dash line in figure 15). In order to compare bulk fluid and confined fluid (see section 2), bulk pressure can be used as the reference in both systems: a single bulk fluid and a bulk fluid, thermodynamically connected to a confined fluid. Then the comparison between P bulk sat and P bulk sat con f enables to say that the bubble point is decreased with confinement. This conclusion cannot be applied into a system at the same thermodynamic equilibrium where confined liquid and vapor pressures are higher than the bulk pressure (figure 15).

The application of the workflow using GEMC NPT BPMC ensemble as an initialization of the confined NVT simulation has been applied to two mixtures and showing very good results. A Phase diagram of confined methane/ethane and ethane/n-pentane has been built with bulk pressure as reference. Furthermore confined vapor and liquid pressures for ethane/npentane have been calculated and compared to bulk pressure which is a total novelty.

Conclusions

The thermodynamic equilibrium properties of several pure components and mixtures in slit graphite pores have been calculated and analyzed. Furthermore, a detailed explanation and a calculation of the different pressures considered in a porous medium with nano-pores have been performed. The pure components studied were methane, ethane, n-pentane and n-decane. The mixtures studied were methane/ethane and ethane/n-pentane.

In this work, the confined GEMC NVT simulation has been used instead of the more traditional GCMC simulation because of its limitations. The GCMC simulation could lead indeed to some errors in the liquid/vapor equilibrium correspond to the bulk pressure of a bulk fluid in thermodynamic equilibrium with a confined fluid. Its critical point is 'CP conf'. P con f v and P con f l are vapor and liquid pressures of a confined fluid in thermodynamic equilibrium with a bulk fluid. P bulk sat is the bulk pressure of a single bulk fluid with critical point 'CP bulk'.

properties because of the challenge in identificate chemical potential phase change. Furthermore the GCMC simulation is not able to give accurate values of confined liquid and vapor pressure. As the GEMC NVT simulation considers two confined boxes for each phase, the accuracy of the results is considerably improved. Since the GEMC NVT needs a good initialization in order to converge, two different workflows for pure components and mixtures have therefore been proposed to overcome this problem. The confined GCMC µVT simulation is used getting approximate initial vapor and liquid densities of pure components for confined GEMC NVT simulation. A new ensemble: the Gibbs ensemble constant temperature and pressure Monte Carlo with bubble-point movement (GEMC NPT BPMC) has been used for the initialization of the confined GEMC NVT simulation for mixtures. The mixture workflow is divided into three steps. The first step uses the standard Gibbs ensemble Monte Carlo at constant pressure and temperature (GEMC NPT). It enables to establish the bulk saturated liquid equilibrium used in the bulk liquid box of the GEMC NPT BPMC simulation. In the GEMC NPT BPMC simulation, two boxes are in thermodynamic equilibrium, the bulk box and the confined box. This simulation allows to get approximate values of densities and molar fraction of each components in each phase. These values are finally used to initialize the confined GEMC NVT simulation. Besides giving the thermodynamic properties of a confined fluid at equilibrium, this workflow allows us to make the link between the bulk fluid pressure and the confined vapor and liquid pressures of a confined fluid.

For all pure components considered in this study, the vapor density increases and the liquid density decreases with confinement compared to the bulk. The critical temperature and pressure are reduced as the pore width decreases and they approach to the bulk values at large pore sizes. However the evolution of critical temperature and pressure versus pore size follow different trends and should be treated with different correlations. Unlike the bulk fluid, where vapor and liquid pressure are equal at equilibrium, confined fluids have different vapor and liquid pressures. The change of sign of capillary pressure as the pore size decreases attests to a more complex behavior where the standard Laplace equation is no longer valid and disjoining pressure may occur [START_REF] Benet | Disjoining Pressure, Healing Distance, and Film Height Dependent Surface Tension of Thin Wetting Films[END_REF][START_REF] Churaev | The shape of the transition zone between a thin film and bulk liquid and the line tension[END_REF].

For the two mixtures studied regardless of the composition, the observations for density are the same as for pure components. The vapor density increases and the liquid density decreases with confinement compared to the bulk. The pressure versus molar fraction diagrams allowed us to build an example of a pressure versus temperature diagram for a specific composition. The phase envelope of the confined fluid is shifted inwards and closes with respect to phase envelope of bulk fluid. The critical temperature and pressure are shifted from the bulk value to a lower value. Finally the bubble-point pressure decreases as the dew-point pressure increases. For a given pore size, the critical temperature and pressure values of pure components have been used to model numerically the critical point of mixtures by the mean of a EOS. These calculated results have been compared to the critical point of mixtures obtained through the use of GEMC NVT simulation. The numerical model using EOS results is consistent with the molecular simulation results. The selectivity of the confined system compared to bulk for the mixtures methane/ethane and ethane/n-pentane in given thermodynamic conditions has been studied. The observations are that the heavier component is preferentially adsorbed in the vapor phase, whereas the lighter component is preferentially adsorbed in the liquid phase. That is to say that compared to bulk fluid, the composition at liquid/vapor equilibrium of the confined fluid has heavier components in the vapor phase and lighter components in the liquid phase than the bulk fluid composition. Finally the different pressures of ethane/n-pentane in a 3nm pore size have been calculated and compared to bulk. As mentioned in section 2 two systems were in reality compared: a single bulk fluid and a bulk fluid in thermodynamic equilibrium with a confined fluid. Only one pressure is present in the first system: the bulk pressure of single bulk fluid. Three different pressures are present in the second system: the bulk pressure of a bulk fluid in thermodynamic equilibrium with a confined fluid, the confined vapor pressure and the confined liquid pressure. The calculation and comparison of the different pressures observed in the confined systems in one of the main contribution of the present work and it is, to the best of our knowledge, a novelty in this field.

All these results for pure components and mixtures provide relevant information concerning the understanding of the phase behavior in confined systems such as shale gas and tight oil reservoirs. This behavior is completely different compared to the bulk fluid. Furthermore, all these data may be used as reference values for the development of radius dependent EOS calibration.

Appendices

Simulation parameters

The different move probabilities and Monte Carlo steps used for non-flexible and flexible molecules are summarized in table 4 for GCMC µVT simulations, in table 5 for GEMC NVT simulations and in table 6 for GEMC NPT BPMC simulations.

Simulation data post-processing

When the stationary state is reached in the simulations, the number of particles inside the simulation box at each step fluctuates around its average value and then the system is at equilibrium. The average of a macroscopic propriety X is then calculated using equation 9.

< X >= 1 n n ∑ i=1 X(r n i ) (9)
where n is the number of configurations and r n i the positions of the particles in configuration i (or sampling i).

The density of component i at equilibrium is calculated using equation 10.

ρ i = < N i > M i V N av (10)
with N i the number of particles i, M i the molar mass of the particle i, V the volume of the simulation box and N av , the avogadro number.

For mixtures, phase densities and molar fractions are calculated using equations 11 to 14.

ρ l = ∑ i < N l i > M i V N av (11)
ρ v = ∑ i < N v i > M i V N av (12)
x i = < N l i > ∑ i < N l i > (13)
y i = < N v i > ∑ i < N v i > (14)
where l subscription refers to liquid phase and v refers to vapor phase.

The pressure is estimated using the virial equation (equations 15)

< p >= <N>k B T V + 1 3V ∑ i ∑ j>i r i j . f i j = <N>k B T V + 1 3V ∑ i ∑ j>i ((r i j f i j ) xx + (r i j f i j ) yy + (r i j f i j ) zz ) (15)
f i j = -∇(U( r i j )) (16)
where N is the average number of particles inside the simulation box, k B is the boltzmann constant, T is the temperature and f i j is the intermolecular force (equation 16) with r i j the distance between particle i and j and U the interaction potential. In this calculation we have included the solid-fluid contribution in z direction of the Steele potential. This fact makes the different elements of the virial (xx, yy and zz) anisotropic, contrary to it is observed in bulk phases when all elements are equal in average. This explains why we have two different pressures inside the liquid and the vapor phases. For the calculation of the critical pressure, we have used the value of the vapor phase.

The critical parameters are estimated with the least square fit of the following scaling law (equation 17) [START_REF] Vega | Phase equilibria and critical behavior of square-well fluids of variable width by Gibbs ensemble Monte Carlo simulation[END_REF].

ρ l -ρ v = B(1 - T T c ) β (17)
where β =0.325 and B is the constant to fit. The critical temperature estimated is then used to calculate the critical density from the least square fit of the following equation (equation 18). [START_REF] Lambert | Critical temperature and density from liquid-vapor coexistence data: Application to refrigerants R32, R124, and R152a[END_REF]

ρ l -ρ v 2 = ρ c + A(1 - T T c ) (18)
where ρ c is the critical density and A is the constant to fit. Finally the critical pressure is obtained from the Antoine equation (equation [START_REF] Stimpson | Effects of Confined Space on Production from Tight Reservoirs[END_REF], which is derived from the Clausius Clapeyron equation, which results from the fitting of the vapor pressure values obtained at each temperature.

ln(P

c ) = C + D T c (19)
where C and D are constants to adjust. C and D correspond respectively to the intercept and the slope of ln(P v ) versus 1/RT, where P v is the vapor pressure.

In the case of a binary mixture, an initial estimate of the critical pressure P c is found, assuming the following scaling law (equation 20) with α=0.325:

ρ l -ρ v = γ(P -P c ) α (20)
The procedure is identical to T c determination for a pure compound. The estimate of critical pressure is used to perform the regression of λ and µ, assuming the second following scaling law (equation 21) with alpha=0.325:

x v -x l = λ 1 (P c -P) + µ(P c -P) α (21)
where x l and x v are the liquid and the vapor composition. The critical pressure P c corresponds to the minimum dimensionless error on (equation 20) and (equation 21)

The regression of critical composition xc is on the basis of equation 22 in the same way as ρ c for a pure compound.

x v + x l 2 = x c + λ 2 (P c -P) (22)
The coexistence densities and compositions can be calculated according to equation 22 to 26

x l = x c + (λ 2 - λ 1 2 )(P c -P) - µ 2 (P c -P) α (23)
x v = x c + (λ 2 + λ 1 2 )(P c -P) + µ 2 (P c -P) α (24)
ρ l = ρ c + γ(P -P c ) + λ 2 (P c -P) α (25)
ρ v = ρ c + γ(P -P c ) - λ 2 (P c -P) α (26)
Local z density of the confined fluid has been built by averaging the number of particles per slice. 100 000 000 Monte Carlo iterations have been performed on a equilibrated GEMC NVT simulation with an output of the particles coordinated every 5000 steps, then there are in the end 20 000 configurations which have been averaged per slice. The space between the two walls has been divided into 100 slices. As the slit width is 3 nm, each slice represents 0.003 nm. The liquid and the vapor boxes used are of lengths 5*5*3 nm for the CH 4 -C 2 H 6 mixture. The liquid box length are 9*9*3 nm and the vapor boxes are 7*7*3 nm for the C 2 H 6 -C 5 H 12 mixture.

Results

The different results for pure components and mixtures are given in figure 17, figure 16 Figure 22: Evolution of critical pressure (left) and temperature (right) versus n-pentane molar fraction for bulk and confined fluid. The blue line corresponds to numerical results obtained from Peng Robinson EOS for the bulk fluid. The red line corresponds to numerical results obtained from Peng Robinson EOS with critical temperature and pressure shift from Jin et al. [START_REF] Jin | Investigating The Effect of Pore Proximity on Phase Behavior And Fluid Properties in Shale Formations[END_REF] and Meyra et al. [START_REF] Meyra | Thermodynamic equations for a confined fluid at nanometric scale[END_REF] correlations respectively. The blue symbols '+' correspond to critical point NPT results for bulk and the red ones correspond to critical point NVT results for confined fluid in H=3 nm pore length.

C

Figure 1 :

 1 Figure1: Schematic of large pore and nanopore network in tight oil reservoirs; oil is represented in green and gas is represented in red. a) Large pore connected with pressure between saturation pressure of a single bulk fluid (P bulk sat ) and saturation pressure of bulk fluid connected to confined fluid where first bubble appears in confined fluid (P bulk sat con f ). Gas appears in every 3 pores. b) Large pore connected with pressure below saturation pressure P bulk sat con f , volume of gas is increasing homogeneously in the 3 pores. c) Large pore (1, 2, 3) connected with nanopore 4 with pressure between P bulk sat and P bulk sat . Gas appears in large pores but not in nanopore. d) Large pore (1, 2, 3) connected with nanopore 4 with pressure below P bulk satcon f . Gas volume in large pores is increasing and gas begins to appear in nanopore.

Figure 2 :

 2 Figure 2: Schematic of the simulation workflow used in this study for pure components and mixtures.

Figure 3 :

 3 Figure 3: Schematic of the GCMC NPT BPMC method for an example of ethane/n-pentane mixture. µ corresponds to chemical potential and the arrows correspond to the different trial moves.

Figure 4 :

 4 Figure 4: GEMC NPT BPMC post processing, example of ethane/n-pentane. A) Confined fluid box with constant volume. B) fluctuation of the number of particles of C 2 H 6 and C 5 H 12 inside the confined slit pore. C) Histogram of the number of particles of each species showing the bimodal probability of particles in the vapor and liquid phase.

Figure 5 :

 5 Figure 5: Results of the liquid vapor NVT simulations for C 10 H 22 for different pore widths. The black curves are reference bulk values from the National Institute of Standards and Technology database website [68].

Figure 7 :

 7 Figure 7: Pressures (left) and capillary pressures (right) of C 5 H 12 . The 'x' symbols corresponds to confined vapor pressure NVT results. The '+' symbols correspond to confined liquid pressure NVT results and the circles correspond to critical point 'CP'. The star symbols on the right plot correspond to the difference between confined vapor and liquid pressure.

Figure 9 :

 9 Figure 9: Example of a pressure-temperature diagram for a mixture of 34.9% methane and 65.1% ethane. The blue line corresponds to EOS Peng Robinson results with the critical point in black diamond. The 'x' symbols correspond to bulk NPT results and the red circles correspond to confined NVT results for H=3 nm with the critical point in black.

  C 5 H 12 /C 2 H 6 = 1.16 for the vapor phase and S C 5 H 12 /C 2 H 6 = 0.7 for the liquid phase. It means 0

Figure 11 :Figure 12 :

 1112 Figure 11: Local z density profile of methane/ethane in a 3 nm slit pore at 220 K and 2 MPa for vapor (left) and liquid (right) for an initial composition of 39.5% methane molar fraction and 60.5% ethane molar fraction.

Figure 13 :

 13 Figure 13: Example of a Pressure-Temperature diagram for a mixture of 59.7% of ethane and 40.3% of n-pentane. The blue line correspond to EOS Peng Robinson results with the critical point in black diamond. The 'x' symbols correspond to bulk NPT results and the red circles correspond to confined NVT results for H=3nm with the critical point in black.

Figure 14 :

 14 Figure 14: Local z density profile of ethane/n-pentane in a 3 nm slit pore at 340 K and 2 MPa for vapor (left) and liquid (right) for an initial composition of 50.4% ethane molar fraction and 49.6% n-pentane molar fraction.

  , figure 18, figure 19, figure 20, figure 10, figure 21 and figure 22.

Figure 18 :CH 4 C 2 H 6 (CH 4 C 2 H 6 (CH 4 C 2 H 6 (Figure 19 :Figure 20 :C 2 H 6 C 5 H 12 C 2 H 6 C 5 H 12 C 2 H 6 C 5 H 12 C 2 H 6 C 5 H 12 Figure 21 :

 1842642642619202651226512265122651221 Figure 18: Pressure molar fraction (left) and pressure density diagrams (right) for methane/ethane at 240 K. The black line corresponds to EOS Peng Robinson results. The red diamonds correspond to bulk NPT results with the critical point (CP) in black circle. The blue symbols 'x' correspond to NPT BPMC results for H=3 nm. The blue symbols '+' correspond to NVT results for H=3 nm with the critical point in blue circle.

  Pure components: 𝐶𝐻 4 , 𝐶 2 𝐻 6 , 𝐶 5 𝐻 12 , 𝐶 10 𝐻 22

				Mixtures:
			𝐶𝐻 4 + 𝐶 2 𝐻 6 , 𝐶 2 𝐻 6 + 𝐶 5 𝐻 12
	Bulk	Confined fluid	Bulk		Confined fluid
	GEMC	Confined	GEMC	ρ 𝑖 𝐿 , 𝑥 𝑖	GEMC NPT
	NVT	GCMC µVT	NPT		BPMC
		ρ 𝐿 , ρ 𝑉		ρ 𝑖 𝐿 , ρ 𝑖 𝑉 , 𝑥 𝑖 , 𝑦 𝑖
		Confined GEMC			Confined GEMC
		NVT			NVT

Table 1 :

 1 Fluid-fluid and solid-fluid parameters.

	Atom	ε	σ	δ	ρ s	M
		(K)	(nm)	(nm)
	CH 4	149.92 0.37372			16.043
	AUA CH 2 86.291 0.34612 0.038405		14.03
	AUA CH 3 120.15 0.36072 0.021584		15.03
	Steele	47.0651 0.38663		3.3

(10 -5

nm -3 ) (10 -3 kg/mol)

Table 2 :

 2 Intramolecular force field parameters.

  .1 and 6.2 respectively. The pure components studied are CH 4 ,C 2 H 6 ,C 5 H 12 ,C 10 H 22 and the mixtures are CH 4 -C 2 H 6 and C 2 H 6 -C 5 H 12 .

Table 3 :

 3 Simulation parameters.

	Simulation		box length (nm)	Monte Carlo steps
	GCMC µVT	L x	L y		L z
		3	3		pore width	5*10 7
	GEMC NVT	Liquid box		Vapor box
	L x L y	L z	L x L y	L z
	5	5 pore width 9	9 pore width	Appendix 6.1
	NPT BPMC	bulk liquid box	confined fluid box
	L x L y	L z	L x L y	L z
	6	6 pore width 7	7 pore width	Appendix 6.1
	quence of this behavior is the reduction of the critical tempera-
	ture with respect to the bulk. The critical temperature is indeed
	estimated with the least square fit of the scaling law function of
	liquid and vapor densities (equation 18).

  Comparison of critical temperature (left) and pressure (right) versus pore length (H) evolution with literature. The stars correspond to the NVT results for the studied pure components. The 'x' symbols correspond to the results of Pitakbunkate et al.

	1.0				1.2
	0.9				1.0
	0 0.3 0.4 0.5 0.6 0.7 0.8 conf /T T c c bulk	2	4	6 H (nm) CH4 NVT 8 10 12 C2H6 NVT C5H12 NVT C10H22 NVT Pitakbunkate CH4 Pitakbunkate C2H6 Singh CH4 Singh C4H8 Jin CH4 Jin C5H12 Jin C10H22	0.8 0 0.0 0.2 0.4 0.6 conf /P P c c bulk	2	4	6 H (nm)	8 10 12 CH4 NVT C2H6 NVT C5H12 NVT C10H22 NVT Pitakbunkate CH4 Pitakbunkate C2H6 Singh C4H8 Meyra CH4 Meyra C5H12 Meyra C10H22
	Figure 6:

  Using the pressure-molar fraction diagrams

	Pressure (MPa)	0.2 0.4 0.6		Pv H=2nm Pl H=2nm CP H=2nm Pv H=3nm Pl H=3nm CP H=3nm Pv H=5nm Pl H=5nm CP H=5nm Pv H=6nm Pl H=6nm CP H=6nm	P v -P l (MPa)	0.4 0.2 0.0 0.2
		0.2 0.0	300	350 Temperature (K) 400	450	0.6	Temperature (K) 300 325 350 375 400 H=2nm H=3nm H=5nm H=6nm

Table 4 :

 4 Simulation parameters for GCMC simulations.

		Move	Probability
	Non flexible molecules Translation	0.25
		Rotation	0.25
		Insertion/deletion	0.5
	Flexible molecules	Translation	0.1
		Rotation	0.1
		Internal regrowth	0.1
		Internal rotation	0.1
		Insertion/deletion	0.6

Table 5 :

 5 Simulation parameters for GEMC NVT simulations.

	Type of molecules	Period Monte Carlo steps Moves	Probability
	Non flexible molecules	1	10 6	Translation	0.3
				Rotation	0.3
				Volume change	0.4
		2	10 8	Translation	0.2
				Rotation	0.2
				Volume change	0.05
				Particle transfer	0.55
	Flexible molecules	1	10 6	Translation	0.25
				Rotation	0.25
				Internal regrowth	0.25
				Internal rotation	0.25
		2	10 8	Translation	0.175
				Rotation	0.175
				Internal regrowth	0.15
				Internal rotation	0.145
				Volume change	0.005
				Particle transfer	0.35

Table 6 :

 6 Simulation parameters for GEMC NPT BPMC simulations.

	Type of molecules	Period Monte Carlo steps Moves	Probability
	Non flexible molecules	1	10 6	Translation	0.3
				Rotation	0.3
				Volume change	0.4
		2	10 8	Translation	0.2
				Rotation	0.2
				Volume change	0.05
				BPMC	0.55
	Flexible molecules	1	10 6	Translation	0.1
				Rotation	0.1
				Volume change	0.6
				Internal regrowth	0.1
				Internal rotation	0.1
		2	10 8	Translation	0.1
				Rotation	0.1
				Internal regrowth	0.1
				Internal rotation	0.1
				Volume change	0.05
				BPMC	0.55

  10 H 22 Results of the liquid vapor NVT simulations for CH 4 ,C 2 H 6 ,C 5 H 12 and C 10 H 22 for different pore widths. The black curves are reference bulk values from the National Institute of Standards and Technology database website [68]. Evolution of critical temperature (left) and pressure (right) versus pore diameter (H) for CH 4 ,C 2 H 6 ,C 5 H 12 and C 10 H 22 .

	0.2 0.3 0.4 0.5 0.6 0.7 Bulk pressure (MPa)		EOS bulk NPT bulk CP NPT bulk NPT BPMC H=3nm NVT H=3nm CP NVT H=3nm	0.2 0.3 0.4 0.5 0.6 Bulk pressure (MPa)	EOS bulk CP NPT bulk NPT bulk NPT BPMC H=3nm NVT H=3nm CP NVT H=3nm
	Figure 16: 0 0.70 0.75 0.80 0.85 0.90 0.95 1.00 T c conf /T c bulk 17: 0.2 0.1	0.4	2 x C 2 H 6 , y C 2 H 6 4 H (nm) 0.6 0.8	6 1.0 CH4 C2H6 C5H12 C10H22	8	H=2nm H=3nm H=5nm Bulk 4 H (nm) 100 200 300 400 500 0 2 6 0.4 0.5 0.6 0.7 0.8 0.9 conf /P P c c bulk CH4 C2H6 C5H12 C10H22 0 CH 4 C 2 H 6 (kg/m 3 )	8

Acknowledgments

The authors would like to thank Pascal Mougin for the fruitful discussions of the analysis of the different thermodynamic data.


Ari Bq

Tau D V1370

Asgarian

A Poro

New ephemeris of two W Uma-type binaries

Keywords: Variable star, W UMa, Photometry

The new light curve and ephemeris of two eclipsing binary systems BQ Ari and V1370 Tau are presented. The observation was made with CCD through C filter and we used the Czech astronomical society's database to extract raw data, to obtain the time of minimums we used Table Curve software and then calculate O-C with the new ephemeris.

Introduction

We study on BQ Ari (00646-00946) and V1370 Tau (01305-01430) are two W UMa-type eclipsing binary. The W UMa-type is very significant to the study about the formation and dynamical evolution of contact binary [START_REF] Haghi | Journal of Occultation and Eclipse (JOE)[END_REF]. The periods of this group of variable stars are typically short and range between 0.25 days to around 1.0 days. Each of the components in the binary are of nearly the same spectral type, from around middle A to early K, with the bulk concentrated in middle F through early G. W UMa has exhibited some slight period changes since its discovery in 1903. One theory to describe this phenomenon is mass transfer between the components of the system, which is viewed to change the period on the basis of angular momentum redistribution [START_REF] Rucinski | Contact Binaries of the W UMa Type, The Realm of Interacting Binary Stars[END_REF]. The visual magnitude of BQ Ari changes from 10.28 to 10.55 and for V1370 Tau is 10.9 to 11.31 magnitude. Both of these binary systems have short period; for BQ Ari is 0.282336 day and for V1370 is 0.295523 day (https://www.aavso.org/vsx/index.php?). The spectral type of these stars are unclear.

In this paper, the new light curve and ephemeris are presented. The data which used had inserted into the variable star section of Czech astronomical society's database (http://var2.astro.cz). We used these data to obtain the time of minimums, then calculate O-C with the new ephemeris for each stars.

Observation

The observations for BQ Ari were made on 16 and 17 October 2017 at Broumov NMO station and for V1370 Tau were made on 15 and 16 November 2017 at Vígľaš station. These data were observed by CCD photometry through the filter clear. The information of stars and observations is given in Table 1. These data transforming from JD to HJD and changing to phase, based on the last ephemeris (Anton Paschke, 2011) the new light curves are presented in Figures 1 and2.

Times of minimum

To determine the exact time of minimums through filter clear, the data were imported into the Table Curve software. Then by fitting a Gaussian curve, we get the minimum. The results are given in Table 2. 5 and6. The diagrams are presented in Figure 3. Also, the new period is given in Table 7 for each stars.

Figure 1 .Figure 2 .

 12 Figure 1. The light curve of BQ Ari

Table 1 . Information of observers and instruments Star Star Coordination (j2000)

 1

	Observer	Instrument	Coordinates

Table 2 . The new Times of minimum in Clear filter

 2 To calculation the new period, at first we should calculate Epoch for each new minimum, so we have used equation (1). In this equation, P (Period) and T0 (ToM) is related to previous ephemeris (Table3) and T is the new Time of Minimum. The results of this equation present in table 4.

	Star	MIN Ⅰ (HJD)	MIN Ⅱ (HJD)
	BQ Ari	2458043.453995±0.00062529	2458043.596128±0.00095971
	V1370 Tau	2458073.547877±0.00054767	2458073.402151±0.00045486

Table 3 . Previous ephemeris of BQ Ari and V1370 Tau

 3

	Star	ToM	Period	Reference
	BQ Ari	2451382.9000 (HJD)	0.2823	(Paschke, 2011)
	V1370 Tau	2456664.23601 (HJD)	0.2955	Bob Nelson's Database

Table 4 . New Epoch and ToM for BQ Ari and V1370 Tau

 4 To obtain a new period, we use the gradient of the diagram. The diagram contains previous ephemeris and calculations. They are given in Tables

	Star	Epoch	ToM	Type
	BQ Ari	23590.96107 23591.46449	58043.45400 58043.59613	Ⅰ Ⅱ
	V1370 Tau	22352.59270 22352.09958	58073.54788 58073.40215	Ⅰ Ⅱ

Table 5 . Previous calculation of BQ Ari

 5

	Epoch	ToM	Type	O-C	Reference
	0	51382.902	-	0	G C V S 4
	12104	54800.2853	Ⅰ	-0.011644	IBVS 5918
	14913	55593.361	Ⅰ	-0.017768	OEJV 0142
	15880	55866.378	Ⅰ	-0.01968	BAVM 225
	18349	56563.4596	Ⅰ	-0.025664	IBVS 6118
	18349.5	56563.6016	Ⅰ	-0.024832	IBVS 6118
	18384.5	56573.48279	Ⅱ	-0.025402	IBVS 6114
	18385	56573.62509	Ⅰ	-0.02427	IBVS 6114
	19572.5	56908.8959	Ⅱ	-0.02746	JAVSO..42..426

Table 6 . Previous calculation of V1370 Tau

 6

	Epoch	ToM	Type	O-C	Reference
	0	51467.79345	-	0	IBVS 5570
	14804.5	55842.85	-	-0.04627088	IBVS 6018
	15179.5	55953.6685	Ⅱ	-0.04972088	IBVS 6029
	15192	55957.3634	Ⅰ	-0.04888588	IBVS 6084
	16431.5	56323.6548	Ⅱ	-0.06097128	IBVS 6063
	17346	56593.90371	-	-0.069858199	IBVS 6092
	17584.5	56664.38308	Ⅰ	-0.07324688	OEJV 0168

Diagram of BQ Ari and V1370 Tau Table 7. The new period of stars

	Star	Period (d)
	BQ Ari	0.2823
	V1370 Tau	0.2955
	O-C = (E -E') × P0	(2)

Finally, to calculating O-C we used the equation 2 (Zahabi, 2014) and the result is given in Table

8

.

Table 8 . O-C of each minimum

 8 ConclusionThis paper's goal is to present new light curve and ephemeris of two W Uma-type eclipsing binaries BQ Ari and V1370 Tau. The ephemeris is calculated in the following order for each stars: -C parameter is calculated based on Equation 2 and the result shown in Table8.

	Star	O-C	MIN
	BQ Ari	0.010990 0.010024	Ⅰ Ⅱ
	V1370 Tau	0.001873 -0.000159	Ⅰ Ⅱ
	BQ Ari MINⅠ(HJD)= 2458043.453995 (±0.00062529) + 0.2823	(3)
	V1370 Tau MINⅠ(HJD)= 2458073.547877 (±0.00054767) + 0.2955	(4)

O

Acknowledgment

This project was prepared by the International Occultation Timing Association Middle East section during IOTA/ME educational and research course 2019.