domain: phys sdu

text:
Antoine Patt

Sylvain Picaud
email: sylvain.picaud@univ-fcomte.fr

Molecular selectivity of CH 4 -C 2 H 6 mixed hydrates : A GCMC study

Keywords: Clathrate, Ethane, Methane, Grand Canonical Monte Carlo, Simulations, Planetary Science

In this paper, we report the first Grand Canonical Monte Carlo simulation study aiming at characterizing the competitive trapping of CH 4 and C 2 H 6 molecules into clathrate hydrates in temperature conditions typical of those encountered at the surface of Titan. Various compositions of the fluid in contact with the clathrate phase have been considered in the simulations, including pure methane, pure ethane and mixed fluids made of various methane:ethane ratios. The trapping isotherms obtained from the simulations clearly show that ethane molecules can be enclathrated at lower pressures than methane molecules. In addition, they evidence that the methane molecules can occupy both small and large cages of the clathrate lattice, whereas the ethane molecules have a strong preference for the large cages, in accordance with experimental conclusions. However, increasing the pressure may also lead to the trapping of ethane in the small cages of the clathrates, leading to a possible competition between methane and ethane molecules for these small cages at high pressure, if both molecules are concomitantly present in the fluid phase. The above mentioned features could strongly influence the composition of a mixed methane:ethane fluid phase in contact with the clathrate phase, which might be thus first impoverished in ethane before methane starts to be 1 trapped into the clathrate. However, this conclusion strongly depends on the clathrate structure considered in the simulations.

Introduction

Clathrate hydrates are ice-like inclusion compounds made of water molecules forming cages stabilized by the presence of entrapped small guest molecules. [START_REF] Sloan | Clathrate Hydrates of Natural Gases[END_REF] Their structure depends largely on the size of the guest species and are most often either made of small pentagonal dodecahedral cages, denoted 5 [START_REF] Mousis | Equilibrium composition between liquid and clathrate reservoirs on Titan[END_REF] (12 pentagonal faces in the cage), and large tetrakaidecahedral cages, denoted 5 [START_REF] Mousis | Equilibrium composition between liquid and clathrate reservoirs on Titan[END_REF] 6 2 (12 pentagonal faces and 2 hexagonal faces in the cage) for the structure I (sI) or of 5 [START_REF] Mousis | Equilibrium composition between liquid and clathrate reservoirs on Titan[END_REF] cages combined with large hexakaidecahedral cages, 5 [START_REF] Mousis | Equilibrium composition between liquid and clathrate reservoirs on Titan[END_REF] 6 4 (12 pentagonal faces and 4 hexagonal faces in the cage) for the structure II (sII). Small hydrocarbons such as methane and ethane induce formation of sI structures while larger molecules such as propane and isobutane rather result in sII structure. [START_REF] Sloan | Clathrate Hydrates of Natural Gases[END_REF] These solids are thought to be ubiquitous on Earth, especially within the surface layer of the frozen lands near the poles (the so-called permafrost) and in oceanic sediments [START_REF] Sloan | Clathrate Hydrates of Natural Gases[END_REF] where they contain not only primarily methane (CH 4 ), but also ethane (C 2 H 6 ), propane (C 3 H 8 ) and, possibly, other gases in much smaller amounts. [START_REF] Takeya | Coexistence of structure I and II hydrates formed from a mixture of methane and ethane gases[END_REF] Clathrate hydrates are also of great concern for oil industry because their formation may completely block flows in pipelines.

Meanwhile, the hydrocarbons that natural clathrates likely contain in huge amount is also viewed as accessible gas resources thanks to recent technological advances allowing their industrial exploitation. [START_REF] Sloan | Clathrate Hydrates of Natural Gases[END_REF][START_REF] Rajnauth | Potential Industry Applications Using Gas Hydrate Technology[END_REF][START_REF] Hassanpouryouzband | Gas hydrates in sustainable chemistry[END_REF] Natural clathrate hydrates are also conjectured to be present in various extraterrestrial environments, especially on other planets and/or their satellites in the Solar System. [START_REF] Mousis | Volatile Inventories in Clathrate Hydrates Formed in the Primordial Nebula[END_REF][START_REF] Mousis | Methane clathrates in the Solar System[END_REF][START_REF] Dartois | Ethane clathrate hydrate infrared signatures for solar system remote sensing[END_REF] While their presence has been invoked to infer, for instance, the possible composition of the interior oceans of Enceladus and Europa, [START_REF] Mousis | Methane clathrates in the Solar System[END_REF][START_REF] Bouquet | Possible evidence for a methane source in Enceladus' ocean[END_REF] one on the most favorable places for the formation of clathrate remains Titan, a Saturn's satellite, where both the presence of light hydrocarbons and the temperature and pressure conditions permit the existence of stable clathrate hydrate phases both on the surface and in its interior. [START_REF] Mousis | Methane clathrates in the Solar System[END_REF][START_REF] Mousis | The fate of ethane in Titan's hydrocarbon lakes and seas[END_REF][START_REF] Vu | Rapid formation of clathrate hydrate from liquid ethane and water ice on Titan[END_REF] More specially, clathrates are likely to form when liquid hydrocarbons enter in contact with the exposed icy crust of Titan. [START_REF] Mousis | Sequestration of ethane in the cryovolcanic subsurface of Titan[END_REF] Their formation is then suspected to strongly influence the compositions of lakes and seas present on Titan, when these liquid reservoirs interact with the clathrate layer. [START_REF] Mousis | The fate of ethane in Titan's hydrocarbon lakes and seas[END_REF] Thus, a model study predicted that, if Titan's lakes interacted with clathrates, they could be strongly depleted in methane and, as a consequence, dominated by ethane and/or propane. [START_REF] Mousis | Equilibrium composition between liquid and clathrate reservoirs on Titan[END_REF] However, this conclusion has been shown to strongly depend on the initial lake composition introduced in the calculations, and re-examination of the same model in different conditions recently led to the demonstration that a sea could become methane-dominated for any initial ethane mole fraction below 0.75. [START_REF] Mousis | The fate of ethane in Titan's hydrocarbon lakes and seas[END_REF] The question of the composition of clathrates from the entrapment of methane and ethane is thus of crucial importance to better quantify the exact composition of the lakes and seas present on Titan.

From an experimental point of view, it has been long shown that both methane and ethane single-guest clathrates usually exhibit sI structure, whereas the clathrate formed from mixtures of methane and ethane molecules may have sI or sII structures (see Figure 1) depending on the mixture composition. [START_REF] Takeya | Coexistence of structure I and II hydrates formed from a mixture of methane and ethane gases[END_REF][START_REF] Murshed | Kinetic studies of methane-ethane mixed gas hydrates by neutron diffractions Raman spectroscopy[END_REF][START_REF] Ohno | Raman studies of methane-ethane hydrate metastability[END_REF][START_REF] Kida | Dissociation behavior of methane-ethane mixed gas hydrate coexisting structures I and II[END_REF] Moreover, both sI and sII structures may also coexist upon certain conditions, because of different metastability. [START_REF] Takeya | Coexistence of structure I and II hydrates formed from a mixture of methane and ethane gases[END_REF][START_REF] Murshed | Kinetic studies of methane-ethane mixed gas hydrates by neutron diffractions Raman spectroscopy[END_REF] However, it is noteworthy that most of these experiments have been performed in thermodynamic conditions that do not correspond to those of Titan and extrapolating their conclusions to the low temperature range of Titan's surface (typically around 90 K) is thus not obvious.

On the other hand, the composition of such mixed clathrate hydrates at any given temperature and pressure can be theoretically estimated by using the statistical approach of van der Waals and Platteeuw that relies on the Langmuir adsorption model. 16 However, the use of this method also implies some assumptions and the knowledge of parameters that are generally issued from experimental results. Again, this raises the question of the transferability of some parameters to thermodynamical conditions different from those where they have been fitted. [START_REF] Sloan | Clathrate Hydrates of Natural Gases[END_REF] More recently, computer simulations based on an atomic-scale description of the systems under consideration have appeared as a promising tool to characterize the composition and the stability of clathrates. 17,[START_REF] Tsimpanogiannis | Monte Carlo simulation studies of clathrate hydrates: A review[END_REF] In this respect, calculations performed on the grand canonical (µ, V, T ) ensemble, where the number of trapped molecules can vary in the simulations, is a particularly suitable tool for characterizing the fraction of the enclathrated gases in the different types of cages, which is one of the most important data sought. Indeed, in the grand canonical Monte Carlo (GCMC) simulations, [START_REF] Adams | Grand Canonical Ensemble Monte Carlo for a Lennard-Jones Fluid[END_REF] the chemical potential rather than the number of the molecules is fixed and thus, the number of enclathrated molecules can be calculated, as a function of their chemical potential or of their partial pressure. Notice that this numerical method can be used for single-guest clathrate as well as to characterize the fractionnal occupancy of multi-guest clathrates. [START_REF] Tsimpanogiannis | Monte Carlo simulation studies of clathrate hydrates: A review[END_REF] Thus, we have recently performed GCMC simulations to characterize the ammonia enclathration in conditions relevant to astrophysical environments, [START_REF] Fàbiàn | Ammonia Clathrate Hydrate As Seen from Grand Canonical Monte Carlo Simulations[END_REF] and to thoroughly investigate the composition of the multiple-guest clathrate formed in contact with a gas mixture of CO and N 2 , at very low temperatures (typically around 50 K) in the context of the comet 67P/Churyumov-Gerasimenko. [START_REF] Lectez | A 32-70 K Formation Temperature Range for the Ice Grains Agglomerated By Comet 67 P/Churyumov-Gerasimenko[END_REF] In addition, the results of GCMC simulations performed at higher temperatures up to 150 K, have been shown to agree fairly well with recently available measured data, especially when regarding the clathrate selectivity which appeared to strongly favor CO at the expense of N 2 trapping. [START_REF] Patt | A grand canonical Monte Carlo study of the N 2 , CO, and mixed N 2 -CO clathrate hydrates[END_REF][START_REF] Petuya | Molecular Selectivity of CO-N 2 Mixed Hydrates: Raman Spectroscopy and GCMC Studies[END_REF] Here, we thus used the GCMC method to investigate extensively, for the first time, the selectivity in the CH 4 -C 2 H 6 mixed clathrate hydrates at various temperatures relevant for the conditions of Titan, and for different compositions of the initial methane-ethane mixture in contact with the forming clathrate. We consider both sI and sII structures because of the possible co-existence of these two phases for mixed clathrates, as evidenced in some experiments (although performed at higher temperatures). [START_REF] Takeya | Coexistence of structure I and II hydrates formed from a mixture of methane and ethane gases[END_REF][START_REF] Murshed | Kinetic studies of methane-ethane mixed gas hydrates by neutron diffractions Raman spectroscopy[END_REF][START_REF] Ohno | Raman studies of methane-ethane hydrate metastability[END_REF][START_REF] Kida | Dissociation behavior of methane-ethane mixed gas hydrate coexisting structures I and II[END_REF] The present paper is organized as follows. The details of the GCMC simulations are provided in Section 2 and the corresponding results for both single-and multiple-guest clathrate phases are discussed in detail in Section 3. Finally, the main conclusions of this study are summarized in section 4.

Methodology

Simulation details

The trapping of methane and ethane molecules in clathrate hydrate lattices of both structures I and II has been investigated by performing a set of GCMC simulations, at the temperature of 91 K typical of the surface of Titan. [START_REF] Mousis | The fate of ethane in Titan's hydrocarbon lakes and seas[END_REF][START_REF] Mousis | Equilibrium composition between liquid and clathrate reservoirs on Titan[END_REF] In addition, two other, lower (70 K) and higher (120 K), temperatures have also been considered for investigating the temperature effects on the results. Notice that, although single-guest methane and ethane clathrates have sI structure, we have also considered here the sII structure for the single-guest clathrates, as a case study to allow a better understanding of the situation that is experimentally observed for mixed clathrates where both sI and sII structures have been evidenced. [START_REF] Takeya | Coexistence of structure I and II hydrates formed from a mixture of methane and ethane gases[END_REF][START_REF] Murshed | Kinetic studies of methane-ethane mixed gas hydrates by neutron diffractions Raman spectroscopy[END_REF][START_REF] Ohno | Raman studies of methane-ethane hydrate metastability[END_REF][START_REF] Kida | Dissociation behavior of methane-ethane mixed gas hydrate coexisting structures I and II[END_REF] The simulations have been performed by considering a flexible clathrate, i.e., a system where rotational and translational degrees of freedom of water molecules are allowed. As in our previous studies on hydrate selectivity, [START_REF] Patt | A grand canonical Monte Carlo study of the N 2 , CO, and mixed N 2 -CO clathrate hydrates[END_REF][START_REF] Petuya | Molecular Selectivity of CO-N 2 Mixed Hydrates: Raman Spectroscopy and GCMC Studies[END_REF][START_REF] Patt | Adsorption of CO and N 2 molecules at the surface of solid water. A grand canonical Monte Carlo study[END_REF] trapping isotherms have then been obtained by performing simulations, for given temperatures, at several increasing values of the chemical potentials corresponding to each guest species, namely CH 4 and C 2 H 6 . Thus, while the number of water molecules has been kept constant, the number of guest molecules has been allowed to fluctuate to reach the equilibrium of the chemical potential between the simulated hydrate phase and a virtual fluid phase. Notice that here, we have chosen to represent all the resulting trapping isotherms as the occupancy as a function of the pressure instead of the chemical potential (see below). As a reminder, the occupancy θ of an hydrate simply corresponds to the ratio of the average number of trapped molecules N guests divided by the number of cages N cages in the clathrate lattice, as :

θ = N guests N cages (1)
Notice that this equation can be also used to define a fractional occupancy per guest species (CH 4 and C 2 H 6 ) and/or per type of cage (small and large).

All the simulations have been performed using the Monte Carlo general purpose GIBBS software package. [START_REF] Ungerer | Applications of Molecular Simulation in the Oil and Gas Industry : Monte Carlo Methods[END_REF] A typical run has been split into equilibration and production stages of 8 × 10 7 and 2 × 10 7 MC steps, respectively. In the production stage, averages and molecular configurations have been saved every 5 × 10 3 steps, thus providing a set of 4 × 10 3 samples.

While all the molecules have been subjected to translation and rotation moves, only the guest molecules (methane and ethane) have been subjected to insertion/deletion and identity swap moves. The latter move, where a molecule of type A is replaced by a molecule of type B at the same position, allows to considerably speed up the convergence of GCMC simulations of mixtures, as long as molecules A and B are not too dissimilar. The probabilities of the translation and rotation moves have always been set to 20%. However, in the case of singlecomponent simulations, the insertion/deletion move has been given a probability of 60%, whereas in the case of mixture simulations, the insertion/deletion move had a probability of 55% and the identity swap move has been set to 5%.

Notice that, in the present work, as in adsorption studies on porous materials, the hydrate lattice has been treated as pre-existing solid substrate in which fluid adsorption could take place, a shared point of view with the van der Waals-Platteeuw theory. 1

Hydrate structures

As mentioned before, both sI and sII hydrate structures have been considered here. Regarding their molecular arrangement, the sI and sII unit cells contain 46 and 136 H 2 O molecules, respectively, and these molecules are assembled in 2 small (5 [START_REF] Mousis | Equilibrium composition between liquid and clathrate reservoirs on Titan[END_REF] ) and 6 large 5 [START_REF] Mousis | Equilibrium composition between liquid and clathrate reservoirs on Titan[END_REF] 6 2 cages for sI, and 16 small and 8 large cages (5 12 6 4 ) for sII. The atomic positions of both structures have been taken from the work of Takeuchi et al., [START_REF] Takeuchi | Water proton configurations in structures I, II, and H clathrate hydrate unit cells[END_REF] in which all possible proton configurations that obey the so-called "ice rules" 27 have been screened in order to find those displaying a zero net dipole moment and the lowest potential energies. The lattice constants of these resulting structures are 12.03 and 17.31 Å for sI and sII, respectively and they have been kept constant in the entire pressure and temperature ranges of the simulations. Such an assumption has allowed to save a significant computational cost and can be justified by the limited effects on trapping capacity caused by volume changes, [START_REF] Papadimitriou | The effect of lattice constant on the storage capacity of hydrogen hydrates: a Monte Carlo study[END_REF] which are themselves usually small. [START_REF] Costandy | Molecular dynamics simulations of pure methane and carbon dioxide hydrates: lattice constants and derivative properties[END_REF] Nevertheless, the hydrate structures have not been considered entirely rigid since H 2 O molecules have been allowed to translate and rotate. The simulation boxes have been built using 27 units cells (3 × 3 × 3) for sI and 8 unit cells (2 × 2 × 2) for sII, therefore containing 1242 and 1088 H 2 O molecules respectively. Thus, the sI simulated hydrates has formed a total of 54 small and 162 cages spanning in a cubic box of size 36.09 Å, while for the sII hydrates, 128 small and 64 large cages have composed the simulation box of size 34.61 Å.

A visualization of these structures made with the VESTA 3 software, 30 is given in Figure 1.

Molecular models and their interactions

The interaction potentials between water and guest molecules have been described classically, as the sum of the pairwise Lennard-Jones dispersion-repulsion and Coulomb electrostatic contributions of all pairs of interaction sites composing the H 2 O, CH 4 , and C 2 H 6 molecules.

The internal degrees of freedom have not been considered and, thus, all the molecules have been treated as rigid bodies. The water molecules have been represented by the TIP4P-Ew model, 31 already used in several hydrate studies, [START_REF] Tung | The Growth of Structure I Methane Hydrate from Molecular Dynamics Simulations[END_REF][START_REF] Ballenegger | Cage Occupancies in Nitrogen Clathrate Hydrates from Monte Carlo Simulations[END_REF] in which it has been shown to give satisfactory results compared to experimental data. In this representation, the oxygen atom is the only LJ interaction site and does not bear any electrostatic charge. The hydrogen atoms, on the other hand, only hold positive charges. The compensating negative charge is then placed on a virtual site M, which lies along the bisector of the H-O-H angle. Methane and ethane molecules have been represented by the TraPPE force field, in its united atom formulation, [START_REF] Martin | Transferable potentials for phase equilibria. 1. Unitedatom description of n-alkanes[END_REF] also used before to study hydrates. [START_REF] Ravipati | Calculation of chemical potentials and occupancies in clathrate hydrates through monte carlo molecular simulations[END_REF][START_REF] Ravipati | Calculation of three-phase methane-ethane binary clathrate hydrate phase equilibrium from Monte Carlo molecular simulations[END_REF] According to that force field, the methane and ethane molecules have been modeled by a single (CH 4 ) and two (CH 3 ) sites, respectively, which are in both cases, only LJ sites. All the interaction parameters are given in Table 1. The cross interaction parameters between unlike LJ sites have been calculated by using the standard Lorentz-Berthelot combining rules. The LJ contributions to the interaction potential have been truncated at a cutoff distance equal to half the simulation box length and long-range corrections have been taken into account. [START_REF] Frenkel | Understanding Molecular Simulation[END_REF] The calculation of the Coulomb interactions has been achieved by use of the Ewald summation technique. 37

Relation between pressure and chemical potential

As stated above, in GCMC simulations, the chemical potential is fixed, allowing the number of molecules to fluctuate. To express the results of such simulations in terms of pressure and not of chemical potential, it is possible to use an equation of state to relate these two quantities, as has been done in several previous studies. [START_REF] Papadimitriou | Evaluation of the Hydrogen-Storage Capacity of Pure H 2 and Binary H 2 -THF Hydrates with Monte Carlo Simulations[END_REF][START_REF] Papadimitriou | Hydrogen storage in sH hydrates: A Monte Carlo study[END_REF][START_REF] Glavatskiy | Toward a possibility to exchange CO 2 and CH 4 in sI clathrate hydrates[END_REF] However, such equations are not always suitable for simulations with empirical force fields for which they were not parametrised. Simulations of fluids in the isothermal-isobaric (N, p, T ) and grand canonical ensembles can also be performed, to interpolate the pressure of the former from the chemical potential of the latter, by using the densities obtained from both sets of simulations. This scheme has the advantage of providing a relation between the pressure and the chemical potential of a fluid phase which takes into account, in an adequate fashion, the force field that is used. However, it can require a huge number of simulations, and even more so when a mixture is considered.

Here, to determine a consistent and accurate relation between pressure and chemical potential with a reasonable computational cost, we have chosen a different approach. Thus, before the GCMC simulations concerning the hydrates have been performed, a series of Monte Carlo simulations in the (N, p, T ) ensemble have been launched, over the same ranges of pressure and temperature, in a cubic box containing at least 1000 CH 4 and/or C 2 H 6 molecules, for each single component and each mixture of interest. In these simulations, the molecules have been subjected not only to translation and rotation MC moves, but also to the identity swap move in the case of the mixtures, allowing to exchange the positions of two randomly selected molecules. These three moves have been performed with equal probabilities of 32% each. The system as a whole has also been subjected to the volume change move, with a 4% probability.

Once the systems have reached thermodynamic equilibrium, after at least 10 8 MC steps, the chemical potentials of CH 4 and/or C 2 H 6 constituting the fluid phase have been calculated

with the Widom insertion method, 41 during 2 × 10 8 additional MC steps. In this method, a test molecule has been randomly inserted into the system and the excess potential energy ∆U N +1 due to the addition has been derived. This operation is comparable to the insertion move from the grand canonical ensemble, but without actually accepting the insertion of the molecule. Thus, the excess chemical potentials, µ ex i , of each species, have been calculated by averaging the energetic contributions of the test molecules according to the following equation

µ ex i = µ i -µ 0 i = -k B T ln V exp - ∆U N i +1 k B T {N,p,T } (2)
with µ 0 i the reference chemical potential and V the fluctuating volume of the simulation box. Note that in the (N, p, T ) ensemble, the volume must be averaged together with the Boltzmann factor of the molecule addition. [START_REF] Frenkel | Understanding Molecular Simulation[END_REF] So, the obtained values of {µ ex i } have been used as inputs in the GCMC simulations, allowing us to make sure that the (virtual) fluid phase in contact with the hydrate structures has always been at the right pressure and composition.

Results and discussion

CH 4 and C 2 H 6 single guest clathrates

The trapping isotherms calculated for flexible single-guest clathrates, that give the average number of either CH 4 or C 2 H 6 molecules enclathrated as a function of the pressure are given in Figure 2, for the three temperatures and the two clathrate structures (sI and sII) considered here. Notice that the pressure axis is given on a logarithmic scale and that, on the y-axis, the number of trapped molecules has been divided by the total number of cages in the simulation box. Hence, an occupancy value of 1 corresponds to 64 and 192 enclathrated molecules in sI and sII clathrate structures, respectively.

First, it has to be mentioned that all these isotherms are, as expected, shifted toward larger pressures when the temperature is increased. Moreover, the comparison of these isotherms shows that larger pressure values are always required to start filling the clathrates with methane than with ethane molecules. In addition, while the simulated isotherms for methane exhibit a type I behavior (cf. IUPAC classification [START_REF] Sing | Reporting Physisorption Data for Gas/Solid Systems with Reference to the Determination of Surface Area and Porosity[END_REF][START_REF] Roquerol | Recommendations for the Characterization of Porous Solids[END_REF] ) irrespective of the temperature and the clathrate structure, there is a noticeable difference with the corresponding isotherms calculated for ethane. Indeed, for the C 2 H 6 molecule, all the trapping isotherms are characterized by an inflection corresponding to a partial occupancy of ∼ 0.75 and ∼ 0.33, for sI and sII structures, respectively. This particular form of isotherm is typical of an adsorption process governed by the presence of two different energetic sites. Taking into account that the unit cell of the clathrate structure is made of 1/4 of small and 3/4 of large cages, or of 2/3 of small and 1/3 of large cages, in the case of sI or sII structures, respectively, the shape of the trapping isotherms simulated for ethane strongly suggests that both the occupancy values of ∼ 0.75 and ∼ 0.33 are associated with the loading of the large cages only, the small cages requiring thus larger pressures for being occupied by the relatively large ethane molecule. By contrast, the Langmuir-like shape of the simulated trapping isotherm for methane indicates that this molecule can indifferently occupy both small and large cages of the clathrate structure.

It is noteworthy that C 2 H 6 appears to be more easily enclathrated than CH 4 irrespective of the temperature, because the isotherms simulated for ethane depart from zero at much lower pressure than those obtained for methane. However, methane clathrates seems to be more easily completely filled than ethane clathrates due to the apparent difficulty for the ethane molecules to fit in the small cages. This feature raises the question of a possible competition between these two molecules for occupying some cages of the clathrates, especially at large pressure values, when considering mixed clathrates of methane and ethane (see below).

To check these assumptions, we have also calculated the (partial) occupancy isotherms in small and large cages, separately. The results are given in Figure 3 for both methane and ethane single-guest clathrates, at the three temperatures considered here. As it can be seen, small and larges cages are concomitantly filled by the CH 4 molecule in structure sI, irrespective of the temperature, whereas the small cages appear to be slightly favored in the case of sII, especially at 70 and 91 K. By contrast, in the case of ethane, the complete filling of the large cages is observed at much lower pressure values than for the small cages, for both sI and sII clathrates and irrespective of the temperature. Moreover, the partial isotherms given in the right panels of Figure 3 clearly show that the onset of the small cage isotherm is observed well after the large cages are completly filled with the C 2 H 6 molecules.

In addition, whereas the partial isotherms of ethane in the small and large cages of sII are characterized by the same type I behavior, it appears clearly more difficult to insert the large C 2 H 6 molecules in the small cages of the structure sI, as indicated by the different shape of the corresponding isotherms. These simulated partial isotherms for ethane which clearly evidence the preferential filling of the large cages with respect to the small ones, are thus in accordance with the two-site behavior of the corresponding global isotherms.

For a better understanding of the respective behavior of methane and ethane, the isosteric heat of enclathration Q st has been calculated for the different systems and temperatures considered here. In a GCMC simulation, this quantity Q st can be easily calculated from both the fluctuations of the number N of molecules and of the internal energy in the clathrate U c as 44

Q st = RT - < U c N > -< U c >< N > < N 2 > -< N >< N > ( 3
)
where <> denotes an ensemble average in the grand canonical ensemble. R is the ideal gas constant and T is the temperature.

The heats of enclathration for the single-guest sI and sII methane and ethane clathrates are given in Figure 4, for the three temperatures considered here. In this figure, Q st is plotted as a function of the clathrate occupancy. Notice that because large fluctuations are usually obtained when calculating Q st for flexible clathrates, [START_REF] Patt | A grand canonical Monte Carlo study of the N 2 , CO, and mixed N 2 -CO clathrate hydrates[END_REF] the heats of enclathration have thus been represented here for rigid clathrates, only, for clarity.

For the single-guest methane clathrate of structure sI, the isosteric heat of enclathration remains more or less constant irrespective of both the clathrate occupancy and the temperature. When correlated with the corresponding isotherm behavior (see Figures 2 and3), this feature confirms that no significant difference is evidenced between small and large cages as far as their loading by the methane molecules in sI clathrate is considered. The situation appears slightly different for the trapping of methane molecules in sII clathrates because, in this case and for T=70 and 91 K, the Q st curve exhibits an inflexion point located at the value of the occupancy equal to about 2/3, i.e., the number ratio of the small cages in the sII clathrates. The corresponding values of Q st are slightly larger below this occupancy threshold than above, suggesting that there is a slight preferential trapping of the methane molecules in the small cages of the sII structure, as already infered from the analysis of the corresponding isotherm behavior. This can be related to the big size of the corresponding large cages in sII, that thus appear as the less favorable cages for the enclathration of CH 4 among all the types of cages present in sI and sII structures (we recall here that large cages of sII are in fact wider than large cages of sI). Notice that the difference between small and large cages for methane in sII clathrates is no longer observed at T=120 K, a temperature for which Q st remains almost constant irrespective of the clathrate occupancy.

By contrast, for the single-guest ethane clathrate, Q st shows an abrupt decrease of several kJ/mol at occupancy values of 3/4 and 1/3, for sI and sII structures, respectively, the heat of enclathration being larger below than above these thresholds. Because the occupancy values at which the decrease of Q st is observed correspond to the number ratios of large cages in the sI and sII clathrates, we can infer that large cages will be more easily stabilized by the ethane molecules than the small cages, in accordance with the behavior of the corresponding trapping isotherms, as analysed above. It is also interesting to mention that, when the small cages are occupied (i.e., at occupancy values above the threshold), Q st remains constant in sII structure whereas it continuously decreases in sI. Again this can be related to the behavior of the corresponding isotherms which show that it is clearly more difficult to insert the large C 2 H 6 molecules in the small cages of the structure sI than in those of the structure sII. Because the small cages are of same size in sI and sII structures, this difference certainly reflects the cooperative effect of the lateral interactions between the enclathrated C 2 H 6 molecules which have much more neighbors when they start to fill the small cages of sI than when they fill those of the sII structure.

In addition, overall, Figure 4 clearly shows that the isosteric heat of enclathration for ethane is always larger than for methane at low occupancy values (corresponding to low pressure values), suggesting that it would easier to trap C 2 H 6 than CH 4 , irrespective of the clathrate structure and of the temperature. However, at larger occupancy (i.e., larger pressure values), these differences vanish and the corresponding values of Q st are almost the same for methane and ethane molecules, which again rises the question of a possible competition between these two molecules for stabilizing the clathrates structures especially at high pressure, as already infered from the analysis of the trapping isotherms.

Mixed CH 4 -C 2 H 6 clathrate

The occupancy of clathrates in contact with a mixture of methane and ethane molecules of different compositions has also been investigated by using GCMC simulations. One should note that, especially at high pressure values, the mixed CH 4 -C 2 H 6 fluid could be not ideal, this nonideality also depending on the composition of the fluid. Thus, we recall that we have performed preliminary Monte Carlo simulations to determine the relation between chemical potential and pressure for each given composition of the methane-ethane mixture considered, ensuring that, in the following results, the mixture always keeps its fixed composition along the entire pressure range investigated here (see Methodology section above).

As an illustration of the results of these calculations, the simulated partial trapping isotherms (i.e., the isotherms giving separately the number of CH 4 and of C 2 H 6 molecules that are enclathrated) are shown in Figure 5 Indeed, in the case of the sI structure for instance, the ethane partial isotherm quickly departs from zero at very low pressure well before the methane isotherm starts, and then reaches a plateau after a sudden jump, corresponding to an occupancy of 0.75. This value being equal to the number ratio of large cages in this clathrate structure, the behavior of the ethane isotherm thus suggests that, first, all the large cages of the clathrate are occupied by the ethane molecules. Then, at larger pressures, the methane molecules start to fill the remaining (small) cages, as shown by the exponential shape of the methane partial isotherms which then reach also a plateau corresponding to a partial occupancy of 0.25. This value is equal to the number ratio of the small cages in sI structure, suggesting thus that the methane molecules fill entirely all the small cages at large pressure values. The situation appears at first sight similar when considering the sII structure, for which, again, the methane loading inside the small cages of the clathrate starts after the large cages are entirely filled by the ethane molecules, as indicated by the occupancy value of 0.33 corresponding to the plateau of the ethane partial isotherms (i.e., a value wich is equal to the number ratio of the large cages in the sII structures). Notice also that the loading of the methane molecules starts at pressure values which are slightly lower than those required to start loading of ethane in the small cages of the sII clathrate (as deduced from the analysis of the isotherms calculated for single-guest ethane clathrates). However, at the pressures for which the methane partial isotherms reach a plateau, a small inflexion is concomitantly obtained in the partial isotherms of ethane, which appears to vanish when the temperature increases. As a consequence, at large pressure values for which the sII clathrate appears saturated (i.e. for which no longer variations are obtained in the simulated isotherms), the ethane and methane occupancies are slightly higher than 0.33 and slightly lower than 0.66, respectively, indicating that contrarily to the case of sI clathrate, there is likely a competiton between C 2 H 6 and CH 4 for occupying some small cages of the sII clathrate lattice (at least at low temperature). These features can be related to the conclusions obtained for the single guest ethane clathrate which showed that ethane molecules are likely to occupy the small cages of sII more easily than those of sI.

To investigate more deeply this possible competition between methane and ethane molecules for the occupancy of the small cages, partial trapping isotherms of both species have also been simulated for various compositions of the mixed fluid in contact with sI and sII clathrates.

These compositions have been defined by the ratio between the number x of CH 4 and the number y = x of C 2 H 6 molecules in the fluid. Values of 1, 1.5, 3, and 9 have been chosen for

x and y, leading thus to six different situations, ranging from a fluid that is very enriched in methane with respect to ethane (composition 9:1) to the reverse case, i.e., a fluid very impoverished in methane with respect to ethane (composition 1:9). The simulated partial trapping isotherms are shown in Figures 6 and7, for the mixtures dominated by methane and ethane, respectively. When considering the sI structure, all the investigated situations give similar isotherms which do not differ from those simulated for the equimolar mixture (see Figure 5) and confirm that the trapping process always starts with the enclathration of the ethane molecules in the large cages of the sI clathrate, and is followed by the enclathration of the methane molecules, at larger pressures, in the small cages, when all the large cages are occupied by the ethane molecules. In addition, no competiton between methane and ethane is observed for the filling of these small cages, with perhaps the exception of fluid phases containing a very large amount of ethane molecules (ratio CH 4 :C 2 H 6 equal to 1:9). Indeed, in this situation only, a small perturbation of the ethane partial isotherm seems to appear when the loading of the sI clathrate with methane starts (right-top panel of Figure 7). For the sII clathrate, the ethane loading starts first, at much lower pressures than for methane, irrespective of the temperature and of the fluid composition. Moreover, the occupancy ratios corresponding to the plateau of the ethane and methane isotherms indicate that ethane and methane preferentially occupy the large and small cages, respectively, and that the filling of the small cages by methane occurs only after the large cages are filled by ethane. However, as the fluid is enriched in ethane with respect to methane (bottom panel of Figures 6 and7), an increasing perturbation is obtained in the ethane partial isotherms upon methane loading, which comes from the competition between these two molecules for occupying the small cages of the sII clathrate. This perturbation starts to be visible at the fluid mixing ratio CH 4 :C 2 H 6 equal to 1.5:1 and becomes very important when this ratio reaches the value of 1:9, i.e., when the fluid is highly enriched in ethane. Indeed, in this extreme case, the final occupancy of the clathrate corresponds to a nearly equimolar trapping of methane and ethane molecules (θ ∼ 0.5 for both molecules), indicating that the ethane molecules may also occupy up to 25 % of the small cages of the sII structure, at large pressure values. This competition between the two hydrocarbon molecules can be related to the isosteric heats of adsorption calculated in the small cages of the single-guest clathrates, which correspond to values equal to ∼ 22 and ∼ 20 kJ/mol for ethane and methane, respectively. Thus, in a pressure range where both molecules can be simultaneously incorporated into the sII clathrate, it is not really surprising to observe the preferential trapping of ethane molecules in the small cages that are not yet occupied by methane.

Conclusions

In this paper, we have investigated the trapping of CH 4 and C 2 H 6 molecules into sI and sII clathrates by performing Grand Canonical Monte Carlo simulations at 70, 91 and 120 K, i.e., temperatures that may be typical of Titan's surface. Various compositions of the fluid in contact with the clathrate phase have been considered, including pure methane, pure ethane and mixed fluids made of various methane:ethane ratios.

As far as single-guest clathrates are concerned, the trapping isotherms obtained from the present GCMC simulations clearly show that the methane molecules can occupy both small and large cages of the clathrate lattice, whereas the ethane molecules have a strong preference for the large cages, in accordance with the conclusions of, for instance, Raman and NMR spectroscopy experiments. [START_REF] Sum | Measurments of clathrate hydrates via Raman spectroscopy[END_REF][START_REF] Dec | Surface Transformation of Methane-Ethane sI and sII Clathrate Hydrates[END_REF] However, increasing the pressure may lead to the trapping of ethane in the small cages of the clathrates (as evidenced by the typical double-site shape of the corresponding isotherms), a feature that has also been experimentally observed using

Raman spectroscopy (although it should be mentioned that these experiments have been perfomed at much higher temperatures and, as a consequence, at much higher pressures than those considered here). [START_REF] Morita | Structure and stability of ethane hydrate crystal[END_REF][START_REF] Subramanian | Trends in vibrational frequency of guests trapped in clathrate hydrate cages[END_REF] This possible occupation of the small cages has also been more recently observed in coupled synchroton X-ray and neutron experiments, which have shown that while 100 % of the large cages of ethane sI clathrate are easily occupied, approximately 5 % of the small cages may also be filled by ethane under certain conditions. [START_REF] Murshed | Kinetics of methane-ethane gas replacement in clathrate-hydrates studied by time-resolved neutron diffraction and Raman spectroscopy[END_REF] In accordance with previous experimental 2, [START_REF] Murshed | Kinetic studies of methane-ethane mixed gas hydrates by neutron diffractions Raman spectroscopy[END_REF][START_REF] Ohno | Raman studies of methane-ethane hydrate metastability[END_REF][START_REF] Kida | Dissociation behavior of methane-ethane mixed gas hydrate coexisting structures I and II[END_REF][START_REF] Murshed | Kinetics of methane-ethane gas replacement in clathrate-hydrates studied by time-resolved neutron diffraction and Raman spectroscopy[END_REF][START_REF] Uchida | Spectroscopic observations and thermodynamic calculations on clathrate hydrates of mixed gas containing methane and ethane: determination of structure, composition and cage occupancy[END_REF] and theoretical 36 investigations conducted however at higher temperatures, the results of the GCMC simulations show that, when present together in the mixed fluid in contact with the clathrate phase, ethane and methane molecules preferentially occupy the large and small cages, respectively. However, the present results indicate that C 2 H 6 is always enclathrated at lower pressures than CH 4 .

In addition, the simulated isotherms also evidence the competition between methane and ethane for occupying the small cages at high pressures, especially inside sII clathrates. In this case, up to 25 % of the small cages may be filled by ethane molecules. This behavior could be related to the experimentally observed replacement of enclathrated methane by ethane, 49 a feature that has been invoked to explain the shape of Titan which could be modified by the presence of hydrocarbon clathrates at the poles of this satellite. [START_REF] Choukroun | Is Titan's shape caused by its meteorology and carbon cycle ?[END_REF] The above mentioned features could strongly influence the composition of a mixed methane:ethane fluid in contact with the clathrate phase, which might be thus first impoverished in ethane before methane starts to be trapped into the clathrate, in accordance with the conclusions of previous thermodynamic models. [START_REF] Mousis | The fate of ethane in Titan's hydrocarbon lakes and seas[END_REF] Thus, to better characterize this possible preferential trapping, we have calculated the selectivity α C2H6/CH4 of the sI and sII clathrates, defined as the ratio of the molar fractions of the C 2 H 6 and CH 4 molecules in the hydrate (x C2H6 , x CH4 ) and fluid (y C2H6 , y CH4 ) phases [START_REF] Petuya | Molecular Selectivity of CO-N 2 Mixed Hydrates: Raman Spectroscopy and GCMC Studies[END_REF][START_REF] Patt | Adsorption of CO and N 2 molecules at the surface of solid water. A grand canonical Monte Carlo study[END_REF]

α C2H6/CH4 = x C2H6 x CH4 y C2H6 y CH4 (4)
Worthwhile to point out, is that the higher the value of α C2H6/CH4 is, the more the system favors the incorporation of ethane at the expense of methane. Considering the pressure at the surface of Titan (1.5 bar), and the expected higher pressures beneath an hydrocarbon lake, it is reasonable to infer, from the isotherms simulated here, that the clathrate phases in contact with mixed methane:ethane fluids on Titan has reached the saturation. As a consequence, the selectivity values have been calculated at high pressures, only, i.e., when all the available cages of the clathrates are occupied The corresponding results are shown in Figure 8, as a function of the molecular fraction of ethane in the fluid phase, and for both sI and sII clathrates. As it could be inferred from the analysis of the isotherms, the calculated selectivity is similar for the three temperatures considered here, although a small temperature effect is visible for the sII clathrate, at very high concentrations of ethane in the fluid phase.

Nevertheless, it clearly appears from these calculations that, for both sI and sII structures, the clathrates remain selective with respect to ethane as far as the composition of the fluid phase does not exceed the proportions of the large cages in the clathrates. By contrast, above the threshold values corresponding to the number of large cages, methane molecules may be trapped in the small cages in such amounts that the clathrate phase becomes selective with respect to methane (i.e., α C2H6/CH4 < 1). These results indicate that in the region where α C2H6/CH4 > 1 , the fluid phase will be depleted in ethane with respect to methane and, as a consequence, the fluid phase, progressively, enriched in methane with respect to its initial concentration. On the contrary, when α C2H6/CH4 < 1, the clathrate phase will be enriched in methane with respect to the initial composition of the fluid phase and we can infer that the fluid phase will be, a contrario, more rich in ethane. These results perfectly agree with those of previous conclusions based on a thermodynamic model, which has however been used for sI clathrates, only. 9

However, it should be emphasized that the threshold for the change evidenced for the selectivity strongly depends on the clathrate structures. Although, available experimental results [START_REF] Takeya | Coexistence of structure I and II hydrates formed from a mixture of methane and ethane gases[END_REF][START_REF] Murshed | Kinetic studies of methane-ethane mixed gas hydrates by neutron diffractions Raman spectroscopy[END_REF][START_REF] Kida | Dissociation behavior of methane-ethane mixed gas hydrate coexisting structures I and II[END_REF] clearly show the possible coexistence of sI and sII structures when considering mixed methane-ethane clathrates, the corresponding experiments have been conducted at high temperatures. Extrapolating their conclusions to the temperature range observed at the surface of Titan is not obvious and thus, it is difficult to know wether the conditions of Titan can actually favor this coexistence or not. Anyway, the formation of the sII structure has always been related to an increased concentration of methane in the mixed fluid under consideration. Thus, if the composition of the lakes on Titan is enriched in ethane, as previously expected, [START_REF] Griffith | Evidence for polar ethane cloud on Titan[END_REF] this would rather favors the formation of sI clathrates and, as a consequence, the enclathration of ethane preferentially to methane (at least, as long as the initial concentration of ethane does not exceed 75 %). In this situation, the contact with the clathrate phase will lead to an impoverishment of the fluid phase in ethane, which could thus become more and more methane dominated, as recently inferred. 9 However, as soon as the concentration of methane in the remaining fluid will exceed a given threshold, formation of sII clathrate might be favored, which in turn, would impact on the selectivity of the clathrate phase. As a consequence, the present GCMC study clearly confirms that the composition of the lakes and seas on Titan could strongly depend on both their initial composition and on the clathrate phases they are in contact with.

Finally, it is worth noting that GCMC is a powerful and versatile tool for investigating preferential guest trapping in mixed clathrates, under (p,T) conditions that are not easily reachable in experimental investigations. This method thus appears especially suitable for investigating the composition of clathrates phases in various astrophysical environments, of course provided that accurate interaction potentials are available to describe the systems . [START_REF] Takeuchi | Water proton configurations in structures I, II, and H clathrate hydrate unit cells[END_REF] The cages formed by the water molecules are displayed explicitly. The blue polyhedra correspond to the small cages (5 [START_REF] Mousis | Equilibrium composition between liquid and clathrate reservoirs on Titan[END_REF] ), and the black and red polyhedra correspond to the large cages of the structure I (5 12 6 2 ) and II (5 12 6 4 ), respectively. This visualization was made using the VESTA 3 software.

Figures and tables

Structure I Structure

  for an equimolar composition of the mixture (i.e., 50 % CH 4 -50 % C 2 H 6 ) in contact with the clathrate phase of structure sI (top panel of the Figure) and sII (bottom panel of the Figure). For this equimolar mixture composition, the partial occupancy isotherms clearly show that C 2 H 6 always starts being enclathrated at lower pressures than CH 4 , irrespective of the clathrate structure and of the temperature.

Figure 1 :

 1 Figure 1: Visualization of the I and this work. The positions of the water have been taken from the work of Takeuchi et al..[START_REF] Takeuchi | Water proton configurations in structures I, II, and H clathrate hydrate unit cells[END_REF] The cages formed by the water molecules are displayed explicitly. The blue polyhedra correspond to the small cages (5[START_REF] Mousis | Equilibrium composition between liquid and clathrate reservoirs on Titan[END_REF] ), and the black and red polyhedra correspond to the large cages of the structure I (5[START_REF] Mousis | Equilibrium composition between liquid and clathrate reservoirs on Titan[END_REF] 6 2 ) and II (5 12 6 4 ), respectively. This visualization was made using the VESTA 3 software.[START_REF] Momma | VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data[END_REF]

30

 30 Figure 1: Visualization of the I and this work. The positions of the water have been taken from the work of Takeuchi et al..[START_REF] Takeuchi | Water proton configurations in structures I, II, and H clathrate hydrate unit cells[END_REF] The cages formed by the water molecules are displayed explicitly. The blue polyhedra correspond to the small cages (5[START_REF] Mousis | Equilibrium composition between liquid and clathrate reservoirs on Titan[END_REF] ), and the black and red polyhedra correspond to the large cages of the structure I (5[START_REF] Mousis | Equilibrium composition between liquid and clathrate reservoirs on Titan[END_REF] 6 2 ) and II (5 12 6 4 ), respectively. This visualization was made using the VESTA 3 software.[START_REF] Momma | VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data[END_REF]

Figure 4 :

 4 Figure4: Heats of enclathration for CH 4 (circles, left panels) and C 2 H 6 (squares, right panels) single-guest clathrates of structure sI (top) or sII (bottom) as a function of the total cage occupancy, at T=70 (red), 91 (blue) and 120 (grey) K. These heats of enclathration have been calculated using rigid hydrate structures to remove the noise due to water motions, for clarity. Note however, that the average values of Q st do not differ between the rigid and the flexible cases.

Figure 5 :

 5 Figure 5: Partial trapping isotherms of mixed CH 4 -C 2 H 6 clathrate hydrates as a function of the pressure, for structures I (top) and II (bottom), at 70 (red), 91 (blue) and 120 (gray) K as simulated for an equimolar composition (1:1) of the mixed fluid in contact with the clathrate phase. Squares and circles represent the results obtained for ethane and methane, respectively. Lines are guide for the eye, only.

Figure 6 :

 6 Figure 6: Partial trapping isotherms of mixed CH 4 -C 2 H 6 clathrate hydrates as a function of the pressure, for structures I (top panels) and II (bottom panels), at 70 (red), 91 (blue) and 120 (gray) K and various compositions of a CH 4 dominated fluid phase in contact with the clathrate phase. Squares and circles represent the results obtained for ethane and methane, respectively. Lines are guide for the eye, only.

Figure 7 : 33 Figure 8 :

 7338 Figure 7: Partial trapping isotherms of mixed CH 4 -C 2 H 6 clathrate hydrates as a function of the pressure, for structures I (top panels) and II (bottom panels), at 70 (red), 91 (blue) and 120 (gray) K and various compositions of a C 2 H 6 dominated fluid phase in contact with the clathrate phase. Squares and circles represent the results obtained for ethane and methane, respectively. Lines are guide for the eye, only.

Table 1 :

 1 Parameters of the interaction potentials used in the GCMC simulations, for the H 2 O, 31 CH 4 , 34 and C 2 H 6 molecules.[START_REF] Martin | Transferable potentials for phase equilibria. 1. Unitedatom description of n-alkanes[END_REF] The Lennard-Jones parameters σ (size parameter) is given in Å, like all the distances, whereas ε/k B (energy parameter) is given in K, with k B being the Boltzmann constant. The coulombic charge is denoted as q, and it is expressed in atomic units.

	Molecule	Site σ (Å) ε/k B (K) q (e) Geometry
	H 2 O (TIP4P-Ew)	O	3.164	81.90	r OH	0.9572
		H		+0.524	H-O-H 104.52°M
				-1.048	r OM	0.125
	CH 4 (TraPPE-UA) CH 4 3.73	148.0
	C 2 H 6 (TraPPE-UA) CH 3 3.75	98.0	r CH3-CH3	1.54

Acknowledgments

This project is supported by the French National Research Agency (ANR) in the framework of the MI2C project (ANR-15CE29-0016) and by the Region Bourgogne Franche-Comté through the grant PASCOA 2020-0052. Calculations have been performed thanks to the computational resources of UBFC (Dijon and Besançon-France).