Vertical compositional variations of liquid hydrocarbons in Titan's alkanofers

According to clues left by the Cassini mission, Titan, one of the two Solar System bodies with a hydrologic cycle, may harbor liquid hydrocarbon-based analogs of our terrestrial aquifers, referred to as"alkanofers". On the Earth, petroleum and natural gas reservoirs show a vertical gradient in chemical composition, established over geological timescales. In this work, we aim to investigate the conditions under which Titan's processes could lead to similar situations. We built numerical models including barodiffusion and thermodiffusion (Soret's effect) in N_2+CH_4+C_2H_6 liquid mixtures, which are relevant for Titan's possible alkanofers. Our main assumption is the existence of reservoirs of liquids trapped in a porous matrix with low permeability. Due to the small size of the molecule, nitrogen seems to be more sensitive to gravity than ethane, even if the latter has a slightly larger mass. This behavior, noticed for an isothermal crust, is reinforced by the presence of a geothermal gradient. Vertical composition gradients, formed over timescales of between a fraction of a mega-year to several tens of mega-years, are not influenced by molecular diffusion coefficients. We find that ethane does not accumulate at the bottom of the alkanofers under diffusion, leaving the question of why ethane is not observed on Titan's surface unresolved. If the alkanofer liquid was in contact with water-ice, we checked that N_2 did not, in general, impede the clathration of C_2H_6, except in some layers. Interestingly, we found that noble gases could easily accumulate at the bottom of an alkanofer.


Introduction
Since the detection of its thick atmosphere by G. Kuiper (Kuiper 1944), Titan, the main satellite of Saturn, has been the subject of many studies (see, for instance, the monograph edited by Müller-Wodarg et al. 2014) and targeted by two major space missions: Voyager and Cassini, while an in situ exploration by the revolutionary rotorcraft Dragonfly (see, for example, Turtle et al. 2020) is in preparation. This unfailing interest led the planetary science community to numerous important discoveries. While expected for decades (Tyler et al. 1981;Owen 1982;Sagan & Dermott 1982;Flasar 1983;Eshleman et al. 1983), lakes and seas of liquid hydrocarbons were discovered by Cassini's instruments in Titan's polar regions (Stofan et al. 2007;Stephan et al. 2010). These observations have opened the door to a unique case of "exo-oceanography", connected to an "exo-hydrology" (or "alkanology") for which liquid methane is the main working fluid. In addition to seas and lakes, dry fluvial channels have been detected (Le Gall et al. 2010;Coutelier et al. 2021) and evidence of the presence of evaporite deposits (Barnes et al. 2011;Cordier et al. 2013;MacKenzie et al. 2014;Cordier et al. 2016) and underground alkanofers has been found (Corlies et al. 2017;Hayes et al. 2017). These latter geological structures were proposed to e-mail: daniel.cordier@univ-reims.fr be the Titanian analogs for terrestrial aquifers. While a porous water-ice matrix plays the role of Earth's porous rocks, the liquid contained in the pores is a mixture of methane and ethane, complemented by some amount of dissolved nitrogen (Cordier et al. 2017). The fast destruction rate of atmospheric methane predicted by pre-Cassini works (Flasar et al. 1981;Yung et al. 1984) lead scientists to hypothesize the existence of a subsurface source of methane. One potential endogenic source of methane could be the complete or partial dissociation of a clathrate hydrate reservoir (e.g., Lunine & Stevenson 1987;Tobie et al. 2006;Choukroun et al. 2010;Muñoz-Iglesias et al. 2018;Petuya et al. 2020) In addition, the possible presence of liquid reservoirs, more or less deeply buried under Titan's surface, has been proposed. For instance, Griffith et al. (2012) interpreted dark patches observed in equatorial regions as an emergence of liquid. A subsurface alkanofer is also mentioned by Sotin et al. (2012) in their discussion of Titan's organic cycle. Mousis et al. (2014) focused their study on the possible interplay between liquids and an icy matrix at least partially made of clathrates. In order to explain the cloud distribution observed by Cassini instruments, Turtle et al. (2018) proposed the existence of a widespread polar subsurface methane reservoir. In the same vein, MacKenzie et al. (2019) mentioned a drainage effect in a subsurface reservoir for interpretation of surface changes. However, the best indi-Article number, page 1 of 15 arXiv:2107.06348v1 [astro-ph.EP] 13 Jul 2021 rect evidence of alkanofer existence to date was found by Hayes et al. (2017). Indeed, the Cassini radar instrument has permitted altimetric measurements of polar regions. In this way, Hayes et al. (2017) found that several maria have their surfaces along the same gravity equipotential level, suggesting the existence of some subterranean connectivity. The study of possible Titan alkanofers is stimulating in many respects. First of all, the methane potentially stored in these reservoirs may participate significantly in the global carbon cycle of Titan (Horvath et al. 2016;Faulk et al. 2020). Secondly, the lack of ethane on the surface of Titan is a long-standing problem (Mousis et al. 2016;Gilliam & Lerman 2016); a response to this question could consist of a chemical vertical stratification of subsurface liquids. Since ethane (molecular weight: 30 g mol −1 ) is heavier than nitrogen (28 g mol −1 ) and methane (16 g mol −1 ), ethane-enriched layers could exist in Titan's deep alkanofers. Finally, the nature, the amplitude, and the temporality of exchanges between Titan's interior and its atmosphere is a "cold case", with an exobiological importance, for which any progress is welcome (Nixon et al. 2018;Kalousova & Sotin 2020). The variation of species in terrestrial hydrocarbon reservoirs is a well-established topic of great industrial interest (see reviews: Chilingar et al. 2005;Obidi 2014;Espósito et al. 2017). In most of the field measurements, a vertical variation of chemical composition is detected (see, for example, Metcalfe et al. 1988), with the lightest compounds floating above heavier ones. However, surprisingly, the reversed situation is also observed (Temeng et al. 1998). Here, we mainly address the question of a similar compositional grading in Titan's crustal hydrocarbon reservoirs. Our approach does not imply large-scale hydrodynamic models like in recent works (Horvath et al. 2016;Faulk et al. 2020), but rather local models focused on species diffusion in a porous icy matrix with a small hydraulic permeability that impedes convective transport. We propose a transposition of the physics of terrestrial oil reservoirs to the Titanian context. In this article, we first discuss the case of an isothermal system. We subsequently include the effect of the geothermal gradient, and we describe our in-depth study of the role of molecular diffusion coefficients. The last section is dedicated to a general discussion about some chemo-physical properties of alkanofers.

Case of an isothermal crust
Contrary to situations where a binary mixture is considered, the diffusion processes in a ternary mixture can no longer be described by a single equation (similar to Eq. 1 in Cordier et al. 2019). For a fluid containing N species, the vectorial diffusion flux − → J (kg m −2 s −1 ) has N − 1 components, and can be written (Ghorayeb & Firoozabadi 2000), for a 1D system related to a vertical z-axis, as with c being the total molar density (mol m −3 ), D M the (N − 1) × (N − 1) molecular diffusion tensor, D P the (N − 1) × 1 barodiffusion tensor, and D T the (N − 1) × 1 thermodiffusion tensor. Among these quantities, P is the pressure (in Pa) and T the temperature (in K). We also have x = (x 1 , x 2 , . . . , x N−1 ), with x i being the mole fraction of species i, taken under a condition of normalization. In the phenomenological Eq. 1, the three terms in the right hand side correspond to different microscopic transport processes. The first one represents the well-known molecular fickean diffusion, driven by compositional gradients. The sec-ond term corresponds to the diffusion induced by pressure variations (also known as barodiffusion), which may lead to gravity segregation. The last term represents thermal diffusion, alternatively called the "Soret effect" in the case of liquids (Soret 1879;Espósito et al. 2017). This effect is linked to the tendency for species, in a nonconvective mixture, to separate themselves under the influence of a temperature gradient. The generic forms of tensors D M , D P , and D T are available in the literature (Ghorayeb & Firoozabadi 2000); they specifically depend on the fickean diffusion coefficients, on the derivative of fugacities, and on the thermal diffusion coefficients. Here, the liquid under consideration involved three species: N 2 , C 2 H 6 and CH 4 . In such a case, at equilibrium (i.e., when J=0 everywhere), the vectorial Eq. 1 may be reformulated as a set of two (here N − 1 = 3 − 1 = 2) partial differential equations: with the D M i j s representing the elements of the diffusion tensor D M . Similarly, the D P i s and the D T i s are respective elements of D P and D T . The indices i and j are related to chemical species; we chose i = 1 for nitrogen, i = 2 for ethane, and i = 3 for methane. Since our system is idealized as a monodimensional reservoir, along the vertical direction, the chemical composition gradients are obtained by integrating the following equations: The abundance of species 3 (methane here) is implicitly computed according to the normalization condition of mole fractions. It is worth noting that the presence of methane is represented in the equation by physical quantities like molecular diffusion coefficient D i3 (see below) or the density of the liquid ρ liq, which depends on the chemical composition. All the detailed formulations of coefficients D M i j , D P i , and D T i are given in Appendix B. We emphasize that generalized diffusion coefficients, as they appear in Eqs 4 and 5, are noted with a "D" throughout the paper, while usual molecular diffusion coefficients are noted with a "D". The temperature gradient sensitivity of molecular fluxes, represented by the D T i s, may be explicitly written The terms a i3 , D i3, andM are defined in Appendix B, but their meaning is not required for the discussion that follows. The thermal diffusion ratio k T, i3 may be expressed as a function of the thermal diffusion coefficient α T, i3 . Unfortunately, the α T, i3 s relevant for our purpose are not available in the literature, and their estimation is not a straightforward task. This is why, in a first approach, we neglected the thermal diffusion contribution, as such the system is considered to be isothermal with a uniform temperature taken equal to the surface temperature T 0 = 90 K. In such a situation, the derivative ∂T/∂z in Eqs 4 and 5 vanishes, and the x N2 x C2H6 x CH4 Fig. 1. Vertical variations of mole fractions (denoted x i ) of the main components (i.e., N 2 , C 2 H 6 and CH 4 ) of Titan's alkanofer liquid. This figure depicts a scenario where the icy crust is assumed to be isothermal at the temperature T 0 = 90 K. The pressure varies between 1.5 bar at the surface (z = 0) and around 130 bar at the bottom of the simulated system (depth of z = 10 km). The three panels correspond to different surface compositions, expressed in mole fractions: (a) N 2 : 0.05, C 2 H 6 : 0.20, CH 4 : 0.75; (b) N 2 : 0.20, C 2 H 6 : 0.20, CH 4 : 0.60;, (c) N 2 : 0.20, C 2 H 6 : 0.30, CH 4 : 0.50. For all computations, the porosity of the water ice matrix was fixed to 5%, the diffusion theory formalism follows Ghorayeb & Firoozabadi (2000), the water ice equation of state was provided by Feistel & Wagner (2006), while the thermodynamic properties of liquids were computed using PC-SAFT. pressure derivative is provided by the hydrostatic equilibrium of the alkanofer: where ρ eff (kg m −3 ) is the effective density, which takes into account the matrix of water-ice, and g Tit (m s −2 ) is Titan's gravity. If Π is the average porosity of Titan's icy crust, then ρ eff = Π ρ liq + (1 − Π) ρ ice where ρ ice is the density of ice I h and ρ liq the density of the cryogenic liquid mixture. We neglected the depth dependency of Π, which can be found in other works (Kossacki & Lorenz 1996), since porosity has no direct influence on molecular diffusion (see also the results reported at the end of this section). Here, ρ liq is computed in the frame of the PC-SAFT 1 (Gross & Sadowski 2001) equation of state (EoS), successfully used in several past studies dealing with the Titan context (Tan et al. 2013;Luspay-Kuti et al. 2015;Cordier et al. 2016Cordier et al. , 2017. This equation of state is also employed to estimate the fugacity derivatives ∂ln f i /∂x j and the partial molar volumesV i , which appear in the expressions of D M i j s andV i s (see Appendix B). The water-ice density ρ ice is evaluated according to a dedicated I h ice equation of state (Feistel & Wagner 2006 Hildebrand (1971), and Vogel & Weiss (1981) methods (nicknamed "BHVW method" in the following) to estimate the viscosity of liquids. The whole procedure is summarized in Poling et al. (2007) (see their Eqs. 9-11.6, p. 9.72). We stress that the WC55 method has a validity restricted to diluted solutions, whereas our studied mixtures are not necessary diluted. We discuss the influence of fickean molecular diffusion coefficients on our results later in this paper. A few D i j estimations at P = 1.5 bar and T = 90 K, obtained with WC55 approach lead to D N 2 −CH 4 = 1.6454 × 10 −9 m 2 s −1 and D C 2 H 6 −CH 4 = 1.2439 × 10 −9 m 2 s −1 ; these values around 10 −9 m 2 s −1 appear consistent with what we can expect for liquids. The set of differential equations governing the alkanofer abundance profile, Eqs 4, 5, and 7, is recognized as an initial value problem, which is also called the "Cauchy problem". Thus, at the surface, meaning at z = 0, the pressure, temperature, and chemical composition have to be fixed and used as boundary conditions. This numerical problem is solved using a standard Runge-Kutta algorithm (e.g., Nougier 1987). For the surface pressure, we took Titan's ground pressure, measured in situ by the Huygens probe; it is close to 1.5 bar (Fulchignoni et al. 2005) and should not vary significantly over the satellite surface. The tem- Influence of PC-SAFT parameters m (number of segments) and σ (segment diameter), representing the size of the considered molecule. Our standard model, corresponding to Fig. 1 panel (b), is in blue, while red lines represent the model in which values of (m, σ) for N 2 and C 2 H 6 were exchanged. In order to avoid any misinterpretation, we emphasize that the red lines have no physical meaning, they correspond only to a toy model dedicated to testing the sensivity of parameters.
perature measured in a tropical region by Huygens is about 94 K, but the temperature in polar regions, where lakes and seas are located, should be a few degrees lower, that is around 90 K (Jennings et al. 2016). The chemical composition of seas remains relatively poorly constrained, this is the reason why we investigated several scenarios, keeping liquid methane as the solvent. In a first approach, with the aim of isolating the effect of gravity, and also because thermal diffusion coefficients are not well known, we chose to neglect the contribution of the geothermal gradient. This gradient is represented by the temperature derivative in Eqs. 4 and 5. As a consequence, the chemical composition gradients, represented by derivatives of compound mole fractions, are directly proportional to gravity. For instance, then, for a low gravity object like Titan, where g Tit = 1.352 m s −2 , we can expect a pretty smooth vertical stratification of chemical species. The widely accepted thickness of Titan's crust should be in the range 80 − 100 km (see e.g., Nimmo & Bills 2010;Tobie et al. 2012), but the surface liquids are likely to be present only within the first few kilometers of the subsurface layers. To recall  (2016) assumed a scale height of 5 km for the vertical hydraulic permeability variations law. On the Earth, fresh water aquifers seem to extend down to ∼ 2 km below the ground, and salty water is also found in oil wells at depths below ∼ 6 km (Baldwin & McGuinness 1976). Then, we decided to fix the lower boundary of our model at 10 km, a value that could reasonably represent the maximum expected thickness of a Titan alkanofer. However, a well-defined value for this limit is not strictly required for our discussion.
The results, obtained within the theoretical framework described above, are plotted in Fig. 1. A first glance at Fig. 1 reveals several striking features: (1) the progression of abundances with increasing depth is clearly linear; (2) the heaviest species tend to accumulate in the deepest layers as expected; (3) nitrogen (N 2 : 28.0134 g mol −1 ) seems to be more sensitive to gravitational effect than ethane (C 2 H 6 : 30.0690 g mol −1 ), although ethane is slightly heavier than nitrogen. This latter aspect can only be explained by nonideal effects within the liquid solution. The PC-SAFT parameters recalled in Table 1 show that molecules are individually associated with a specific set of parameters: the number of segments m, the segment diameter σ (Å), and the segment energy . According to the paradigm of PC-SAFT theory, individual molecules are idealized by a collection of hard spheres (called "segments"), and the parameters m and σ are the associated geometrical parameters. Of course, in realistic situations, individual molecules are not just a collection of identical hard spheres. Since the parameters m and σ are adjusted in order to fit experimental results, m are real numbers, not necessarily integers. Moreover, interaction parameters, denoted k i j , are introduced to account for inter-species interactions. Their values were taken from previous papers (Cordier et al. 2016;Tan et al. 2013): k C 2 H 6 −N 2 = 0.07, k C 2 H 6 −CH 4 = 0.00 and k N 2 −CH 4 = 0.03. We checked, by setting all these k i j s to zero, that the interaction parameters have no influence on the obtained abundances profiles. Similarly, by increasing the value of (N 2 ) to a value comparable to that of other species, we found the segment energy having no significant role. Finally, we exchanged the (m, σ) tuple between N 2 and C 2 H 6 . For a given molecule, the tuple (m, σ) represents its size. In Fig. 2, we compare a standard model with a model where we switched these (m, σ) values. This operation has a huge effect on the resulting profiles, with ethane becoming the dominant compound in the deepest layers of the reservoir. Therefore, the surprising behavior depicted in Fig. 1, leading to an alkanofer bottom more enriched in N 2 than in C 2 H 6 , is caused by nonideal effects due to the difference in size (see m and σ values in Table 1) of these molecules. Although the nitrogen molecule is slightly lighter than the ethane one; the smaller size of nitrogen molecule overcomes the effect of gravity.
Vertical variations of mole fractions (denoted x i ) of the main components (i.e., N 2 , C 2 H 6 and CH 4 ) of Titan's alkanofer liquid. This figure depicts a scenario where the icy crust undergoes a geophysical thermal gradient derived from a model based on Sohl et al. (2014). While at the surface the temperature is T = 90 K and the pressure has the ground value of 1.5 bar, at a depth of 10 km the pressure varies between 1.5 bar at the surface (z = 0) and around 130 bar at the bottom of the simulated system (depth of z = −10 km). The three panels correspond to different surface compositions, expressed in mole fractions: (a) N 2 : 0.05, C 2 H 6 : 0.20, CH 4 : 0.75; (b) N 2 : 0.20, C 2 H 6 : 0.20, CH 4 : 0.60; (c) N 2 : 0.20, C 2 H 6 : 0.30, CH 4 : 0.50. For all computations, the porosity of the water ice matrix has been fixed to 5%, the diffusion theory formalism follows Ghorayeb & Firoozabadi (2000), the water-ice EoS is provided by Feistel & Wagner (2006), while the thermodynamic properties of liquids are computed using PC-SAFT. For comparison, the results of Fig. 1 are added in thin gray lines.
To proceed to quantification, we introduced a vertical enrichment ratio for a given species i : with mole fractions x i subscripted with b or s corresponding to bottom or surface of the system. For the scenarios reported in Fig. 1, ∆ N 2 ranges between 34% and 45%, showing a clear enrichment. In contrast, ethane undergoes a lower enrichment with ∆ C 2 H 6 between 4.6% and 40%. For each scenario, the nitrogen enrichment is larger than ethane one.
We found that the surface pressure has only a negligible influence on the vertical composition profile. For instance, if we fix this pressure to 3 bar, a value that could be reached at Titan's sea bed (see Cordier et al. 2017), the profile remains unchanged. Equivalently, we obtain a globally similar picture by varying the temperature between 85 and 95 K or by changing the assumed uniform porosity from 5% to 20%. The latter parameter has only an indirect influence via the effective density ρ eff , which changes the value of the pressure through Eq. 7.

Contribution of thermal diffusion
Due to the vertical geothermal gradient, fluids confined in terrestrial oil or gas reservoirs undergo thermodiffusion, which can have, over geological timescales, consequences on the segregation of chemical species comparable to that of pressure and gravity (Galliero et al. 2017). In some cases, thermodiffusion shows even counterintuitive effects, leading heavy fluid mixtures floating on top of light fluids layers (Ghorayeb et al. 2003). Thus, it is important to tentatively estimate the possible amplitude of this effect in our context. Like many other moons, Titan has internal energy sources provided by radio elements and tidal effects (Tobie et al. 2006). This heat has to be dissipated and it is transported toward the surface through geological layers. According to numerical models (Sohl et al. 2014), the flux at the surface should be ∼ 1 mW m −2 , to be compared with the average geothermal flux at the surface of the Earth, which is more or less two orders of magnitude higher (Pollack et al. 1993). The thermal structure of Titan's crust can be reconstructed (see the description of our model in Appendix C), assuming an I h -ice viscosity (Sohl et al. 2014) high enough to impede solid-state convection (Nimmo & Bills 2010). The results show a linear evolution of the temperature within the layers of interest, corresponding to a gradient of ∼ 0.6 K km −1 , well below the terrestrial geothermal gradient that has a commonly accepted value around 30 K km −1 (Lowrie 2010). The existence of ice convection in Titan's outer crust is still debated (Nimmo & Bills 2010;Liu et al. 2016;Kalousova & Sotin 2020). However, a recent work (Kalousova & Sotin 2020) suggests the existence of some convection below ∼ 15 km, allowing for the possible icy alkanofer to remain in a static state.
In Equations 4 and 5, the geothermal gradient terms involve two coefficients: D T 1 and D T 2 , which contain thermal diffusion coefficients (see Eq. B.17 and B.18 in Appendix B). Here, we need two coefficients: for N 2 in CH 4 (denoted α N 2 −CH 4 ) and for C 2 H 6 in CH 4 (denoted α C 2 H 6 −CH 4 ), which are not available in the literature. In the general case, the thermal diffusion coefficient, α i j,T , of two species i and j is determined, in a complex manner, by the sizes and masses of molecules, the temperature and composition of the mixture, and the intermolecular interactions. This aspect is particularly relevant for nonideal fluids like those involved here. While performing accurate laboratory measurements of α i j,T s under microgravity conditions is difficult (Hu et al. 2014), theoretical estimations are possible. Consistently with the first model described in the previous section, for species (1) and (2) α 12,T can be derived from PC-SAFT via Haase's formula (Haase 1969;Pan et al. 2006): where the M i s and x i s are, respectively, the molecular weights and the mole fractions, and µ 1 is the chemical potential of component (1). The coefficient α 0 12,T represents the thermal diffusion coefficient for the corresponding ideal gas,h res 1 andh res 2 are the respective residual partial molar enthalpies of species (1) and (2). Even for most nonideal fluids α 0 i j,T α i j,T (Pan et al. 2006), the coefficients α 0 i j,T can be easily computed by applying the kinetic theory of gases (Chapman & Cowling 1970). On their side, the residual enthalpies are directly provided by PC-SAFT (see Eq. A.46 of Gross & Sadowski 2001): whereã res is the residual and reduced (i.e., molar) Helmholtz free energy, and Z is the compressibility factor (Gross & Sadowski 2001). The residual partial molar enthalpies are simply derived usinḡ In our framework, the term "residual" relates to a difference between the physical quantities for the actual fluid and the corresponding ideal gas. In the case of the chemical potential derivative, we also have with k B being the Boltzmann constant. The dimensionless ratio µ res /k B T is directly provided by PC-SAFT (see Eq. A.33 of Gross & Sadowski 2001). Of course, for any physical quantity Q, its residual counterpart Q res approaches zero when the system approaches an ideal gas state. This way, for an ideal gas ∂µ as expected. In Fig. 4, for illustration purposes, we report thermal diffusion factors k i j,T for the required pairs of species: N 2 -CH 4 and C 2 H 6 -CH 4 . Except for x 1 near zero or one, the sensitivity of nonideal systems (panel b) to temperature, represented 12,T is chosen, because it represents the sensitivity of the system to a gradient of temperature. In this panel, all quantities are related to the ideal gas state (marked by the "0" superscript). (b) Similar plot using the thermal diffusion coefficient provided by Haase's theory (see Eq. B.19). In both panels, the solid line is associated with the C 2 H 6 -CH 4 system, while the dashed line is related to N 2 -CH 4 . All computations represented in this figure assume a pressure of 1.5 bar and a temperature of 90 K. by k i j,T , is larger than the sensitivity of ideal systems (panel a). Moreover, according to these estimations, ethane seems to be more sensitive than nitrogen: |k Haase C 2 H 6 −CH 4 ,T | ≥ |k Haase N 2 −CH 4 ,T |. If we consider only thermodiffusion in a binary system, the flux j 1 of the species (1) along the vertical z-axis is given by (see Eq. 1 in Cordier et al. 2019) In the case of our icy crust model, we found |∂T/∂z| ∼ 0.6 K km −1 ; since k Haase C 2 H 6 −CH 4 ,T ≥ 0 and k Haase N 2 −CH 4 ,T ≤ 0, we have j C 2 H 6 ≥ 0 and j N 2 ≤ 0. As a consequence, the introduction of thermodiffusion should reinforce the abundance of nitrogen in the deepest parts of Titan's alkanofer. Of course, this approach is very simplified, if not naive, as it ignores the presence of the third species present in a ternary mixture, which we take into account in our full model. In Fig. 3, the profiles of molecular abundances obtained for the models including the geothermal gradient (colored lines) are compared to previous models corresponding to an isothermal icy matrix (gray lines). Clearly, the simple prediction made above seems to explain the results depicted in Figs. 3.b and 3.c, corresponding to scenarios where nitrogen has a relatively high surface-mole fraction (i.e., x N 2 ,s = 0.20). The first scenario, corresponding to panels (a) in the same figure, show a different behavior due to the low value of x N 2 ,s (0.05). In this case, ethane can undergo a substantial enrichment in the lowest layers of the alkanofer. This result may be regarded as rather obvious, the competition between nitrogen and ethane being minimized. Again, if the amount of dissolved nitrogen is not too small, the bottom of the alkanofer is enriched in nitrogen.
As part of their work, Kalousova & Sotin (2020) built a 2D model of Titan's icy crust. These authors assumed the existence of a layer of clathrates at the surface of the crust. Since clathrates have a lower thermal conductivity than regular water ice, the geothermal gradients they found are much higher than those considered up to this point. Depending on the value of the width h c of the clathrated layer, Kalousova & Sotin (2020) obtained |∂T/∂z| ranging between 2.8 K km −1 and 23.8 K km −1 . Included in our model, such high gradients reinforce the effect already seen in Fig 3 and discussed above, that is an enrichment in nitrogen and depletion in ethane. As an example, starting at the surface with a composition of 20% of N 2 , 20% of C 2 H 6 , and 60% of CH 4 , and assuming ∂T/∂z ∼ −23.8 K km −1 (corresponding to h c = 5 km in Kalousova & Sotin 2020), we obtained x N 2 ∼ 0.30 and x C 2 H 6 ∼ 0.0 at the depth of 1 km. Below this depth, our model is no longer valid since it is dedicated to ternary mixtures. This example raises the question of the partial vaporization of the mixture. In Fig. 5, we plot a phase diagram of the binary system N 2 +CH 4 for temperatures ranging from 110 K to 125 K. In the crust, at a depth of 1000 m, the pressure should be around 14 bar, and the temperature, according to Kalousova & Sotin's gradient, should be between 114 K and 120 K. As we can see in Fig. 5, the point (x N 2 ∼ 0.30, P = 14 bar), indicated by a star, is not very far from the liquidus at 120 K (triangle in Fig. 5). Even if our estimations do not favor the appearance of nitrogen enriched bubbles in the possible Titan alkanofer, the question has to be kept in mind for future studies. Indeed, a subsurface degassing process could be the origin of geological activities. We can imagine the formation of chemical composition gradients, driven over geological timescales by diffusion, feeding the formation of some pockets of nitrogen which could, at some point, be released abruptly bringing other materials to planetary surface. This kind of mechanism may explain morphological evidence for volcanic activity in Titan's north polar regions (Wood & Radebaugh 2020).

Role of molecular diffusion coefficients
According to the formalism developed by Ghorayeb & Firoozabadi (2000), which we adopted in the present work, the molecular "fickean" diffusion is represented by molecular diffusion coefficients, D i j s, which are key quantities. In the case of our ternary mixture, only two of them are explicitly included in the set of equations employed here: D 13 and D 23 , corresponding, respectively, to the diffusion of N 2 and C 2 H 6 in liquid methane. In models described in Sect. 2 and Sect. 3, the values of these coefficients are estimated using the approximation proposed by Wilke & Chang (1955), in this section we discuss the possible influence of these molecular diffusion coefficients.
A quick look at the expressions of the Onsager phenomenological coefficients (see Eqs. B.12, B.13, and B.14) leads to L 11 ∝ D 13 , L 22 ∝ D 23 and L 12 = L 21 ∝ D 13 , respectively. As a consequence, the ratio L 12 /L 11 is independent of all molecular diffusion coefficients, whereas L 21 /L 22 depends on D 13 /D 23 since L 21 /L 22 ∝ D 13 /D 23 . Therefore, the terms (see Eqs. B.15,B.16,B.17,and B.18). These dependencies have a quantitatively significant influence on derived chemical vertical profiles only if the values of other thermodynamical quantities do not inhibit this possible influence. This is why we performed sensitivity tests, multiplying D N 2 −CH 4 and D C 2 H 6 −CH 4 by an arbitrary large factor. Setting this factor to 10 3 , we found no influence on results concerning the "isothermal crust", whereas a relatively small influence arises when we take into account the geophysical thermal gradient. The variation of the vertical enrichment ratio (see Eq. 9) is about a few tenths of a percent for nitrogen. The abundance profile of ethane seems to be more sensitive, with ∆ C 2 H 6 going from −42% (for our "standard" model) to −39% (case D C 2 H 6 −CH 4 × 10 3 ) or −28% (case D N 2 −CH 4 × 10 3 ).
In order to investigate the actual influence of these molecular diffusion coefficients more precisely, we also evaluated the diffusion coefficients of N 2 and C 2 H 6 in CH 4 by molecular dynamics (MD) simulations. A similar technique has already been employed, in the context of Titan, by Kumar & Chevrier (2020). Here, we have employed the well-known open-source package GROMACS 2,3 2018 version (Abraham et al. 2015). Nitro-gen molecules have been described by the TraPPE 4 force field (Potoff & Siepmann 2001). Alkanes (i.e., methane and ethane) were modeled by adopting the TraPPE-UA 5 coarse-grained representation where pseudoatoms are introduced for methane and ethane methyl groups. (Martin & Siepmann 1998;Mansi et al. 2017). For all computations, the intermolecular interaction cutoff has been fixed to 1.4 nm, since it appears to be a well -suited value (Martin & Siepmann 1998;Shah et al. 2017). We split our MD calculations into three steps: (1) a 1-ns equilibration step in the NVT ensemble; (2) a 19-ns equilibration step in the N pT ensemble, where the temperature and pressure are fixed to the expected target values; (3) a 10 ns-accumulation step in the N pT ensemble. In our work, standard simulations were performed with a total number of molecules fixed to 5000 with 4000 CH 4 (80%) and 1000 N 2 or C 2 H 6 (20%), the infinite spatial extent of the liquid being reproduced by imposing periodic boundary conditions to the simulation box in the three space directions. In Table 2, we compare molecular diffusion coefficients derived 2). The differences between the coefficients obtained with the MD and WC methods remain below a factor of two, which is clearly not enough to significantly alter the alkanofer chemical vertical profiles. We can safely conclude that molecular diffusion coefficients have no influence on these profiles.

Liquid circulation between subsurface reservoirs
In this work, the most important assumption is the absence of macroscopic flows between underground reservoirs, and inside the considered reservoir. More precisely, the timescales associated with convection and flows through porous media have to be much larger than timescales related to diffusive processes, like those discussed in previous sections. The order of magnitude of the "diffusion timescale", here denoted τ diff , may be easily estimated using a relation similar to with L (m) being a characteristic length of the studied system, and D (m 2 s −1 ) a diffusion coefficient. According to values reported in Table 2, the diffusion coefficient is on the order of ∼ 10 −9 m 2 s −1 . For an alkanofer of depth L ∼ 100 m, these values lead to a timescale of τ diff ∼ 0.3 mega-year, and τ diff ∼ 30 mega-years for a 1 km deep reservoir. These 4 Transferable potentials for phase equilibria 5 Transferable potentials for phase equilibria -united atom timescales are comparable to those computed for terrestrial hydrocarbon reservoirs (Obidi 2014).
According to Darcy's law (Darcy 1856), the volume of liquid ∆V (m 3 ) exchanged between two reservoirs, with free surfaces altitudes difference of h, along a subsurface "channel" through a porous medium of length L and cross section S , is given by where K p,sc is the permeability (m 2 ) of the porous medium that connects two reservoirs or a reservoir and a lake, and the subscript "sc" stands for "subsurface channel". In Eq. 17, ρ is the density of the liquid (kg m −3 ), g is the gravity (m s −2 ), η is the dynamic viscosity (Pa s), ∆t is the time required to obtain free surfaces at the same equipotential. In other words, the transport of the volume ∆V of liquid has the duration ∆t. If we consider two reservoirs with identical permeability and geometries, one in a terrestrial context, the second on Titan, the relationship between Earth and Titan timescales is given by For the Earth, the working liquid is water, and in the case of Titan we considered a typical mixture N 2 +CH 4 +C 2 H 6 . We found a result that is in full accordance with an estimation made by Horvath et al. (2016). This longer timescale clearly favors the formation of liquid pockets, isolated from a large-scale subsurface environment. However, it is difficult to give a realistic estimate of horizontal seepage between two alkanofers, or between a lake and an alkanofer. In their work, Hayes et al. (2008) computed an interaction timescale τ 1/2 (see their Eq. 2) that depends on geometrical factors and on the terrains permeability if we put aside the role of surface evaporation. Nonetheless, we can reformulate Eq. 17 as where the ratio η/ρg represents the properties of the fluid, in the context of Titan (g), for this term we found roughly 10 −8 m s, since the dynamic viscosity η may be computed by using the BHVW method, already mentioned in Sect. 2. The factor L/S represents the geometry of the "subsurface channel", through which the liquid is supposed to circulate. If we take plausible values, for instance L ∼ 1 km and S ∼ 1 km 2 , then L/S ∼ 10 −3 m −1 . Moreover, we can consider that the composition of the targeted alkanofer should be significantly modified if ∆V is large enough; one can take, for instance, ∆V ∼ h 3 , letting the ratio ∆V/h be comparable to h 2 . On Titan, the topography is rather flat, h ∼ 1 km is a large value. Keeping this value, we obtain ∆V/h ∼ h 2 ∼ 10 6 m 2 , which leads to Thus, the effect of diffusion is not erased by inter-reservoir seepage if, roughly, K p,sc 10 −5 /∆t. For a diffusion timescale, for example ∆t ∼ 1 mega-year (i.e., 3 × 10 13 s), we have K p,sc 3 × 10 −19 m 2 . This result corresponds to impervious geological layers (see Table 5.5.1 in Bear 1972), which should surround the alkanofer, leaving the liquid trapped in a pocket for geological timescales. This possibility is compatible with topographic data analysis by Hayes et al. (2017), who observed lakes or dry lake beds at elevations that require either a localized aquifer or a low permeability. Nevertheless, geometrical parameters, presented in Eq. 20, are very uncertain and may conspire to yield a K p,sc value orders of magnitude higher or lower than what we find.

Convection in trapped liquids
When a fluid trapped in a porous matrix undergoes a temperature gradient, convection appears and the thermal convection flux overcomes the viscous forces; the ratio between these quantities is measured by the Rayleigh number (Ra) (Bories & Combarnous 1973): with α V being the volumetric thermal expansion coefficient (K −1 ), K p,m the permeability of the reservoir matrix (m 2 ), g the acceleration due to gravity (m s −2 ), (ρC p ) fluid the product of the fluid density (kg m −3 ) and the heat capacity (J K −1 kg −1 ), H the thickness of the reservoir (m), ν the kinematic viscosity of the fluid (Pa s kg −1 m 3 , i.e., m 2 s −1 ), and λ L the thermal conductivity (W m −1 K −1 ) of the medium. It is accepted that convection occurs when Ra 4π 2 (Nield & Bejan 2017). The volumetric coefficient of thermal expansion α V is given by with V m the molar volume of the fluid, and ρ the density. For a given composition, pressure, temperature, heat capacity C P, and density can be provided by PC-SAFT. The thermal conductivity λ L can be evaluated according to Latini et al. (1996) and Poling et al. (2007) for saturated hydrocarbons, and following Powers et al. (1954) for liquid nitrogen. Finally, the kinematic viscosity ν is easily derived from dynamic viscosity η = ρν, which is estimated according to the BHVW method mentioned in Sect. 2. As a robustness check for the validity of estimations performed accordingly to the above-mentioned methods, we also derived values of C P , α V, and η from MD simulations on the ternary mixture CH 4 (60%)-C 2 H 6 (20%)-N 2 (20%). The outputs of these compulations, made in the frame of the isothermal-isobaric (N pT ) ensemble, have been processed according to relevant fluctuation formulas, leading to results similar to those obtained with PC-SAFT and other cited methods. The latter have been adopted for routine calculations since they are much more computationally flexible and efficient. For an alkanofer 1 km deep and containing a typical mixture made up of 20% N 2 , 60% CH 4 , and 20% C 2 H 6 , under a thermal gradient of 0.6 K km −1 , we found Ra ∼ 3.5 (i.e., Ra/4π 2 < 1) (Nield & Bejan 2017) for a permeability of K p,m = 10 −12 m 2 , equivalent to that of very fine sand on Earth (see Table 5.5.1 in Bear 1972). Within the terrestrial context, the empirical relationship between grain diameter d (in micron) and K p,m are available in the literature (Bear 1972 convection no convection Fig. 6. Influence of the chemical composition of the fluid on its Rayleigh number. These computations were made with a fixed permeability K p,m = 10 −12 m 2 , P = 20 bar, T = 90 K, H = 1 km, and a thermal gradient of 0.6 K km −1 . The influence of ethane content (blue curve) has been explored by forcing the mole fraction of nitrogen to follow that of methane: x N 2 = 0.20 x CH 4 . The role of nitrogen abundance (red curve) has been studied by keeping ethane mole fraction constant: x C 2 H 6 = 0.01.
favorable factors for the absence of convective transport in a Titan alkanofer. The role of chemical composition does not appear explicitly, but physical quantities such as α V , ρ, C P , ν, and λ L depend on the composition. In Fig. 6, for typical values of these physical parameters, we varied the mole fractions of nitrogen and ethane independently. A high abundance of ethane would damp convective fluxes, while large nitrogen abundances would have the opposite effect. Finally, the low permeability requirement has implications concerning the formation of a "diffusive alkanofer". Liquids of course have to fill the reservoir either following a "diffusive regime" if the permeability is initially low, or following a "hydrodynamic regime" if the permeability is relatively large during the early ages of the reservoir. In the latter case, the permeability has to decrease with time due to geophysical processes like tectonic activity or sedimentary deposits; this way, the permability could reach values low enough to damp convection. In this work, the chemical composition of the solid matrix has been hypothesized to be I h -water ice. Of course, even if this kind of scenario is very plausible and supported by our knowledge of Titan's inner structure, the actual Titan reservoirs could be somewhat different in nature. Carbon -bearing species seem to be ubiquitous at the surface of the moon (Lorenz et al. 2008), Cassini radar measured a bulk dielectric constant = 2 (Elachi et al. 2005), inconsistent with water ice ( = 3.1) or ammonia ice ( = 4.5). The interpretation of radar observation is not unique and is still discussed (Hofgartner et al. 2020). Nonetheless, one can assert safely that the surface and the near subsurface of Titan are complex systems with horizontal and vertical heterogeneities. This is the reason why in some places, a porous deposited organic material a few hundred meters deep may host liquid hydrocarbons, leading to processes like those discussed in the present work. In our model of crust, the waterice thermal conductivity has a value around 5 W m −1 K −1 , and methane hydrates are much more insulating with values around 0.5 W m −1 K −1 , as emphasized by Kalousova & Sotin (2020). If we adopt bitumen as an analog of Titan's organic sediments, we find a thermal conductivity close to the clathrate one (Nazki et al. 2020), and, based on a study of Titan's crater relaxation, Abundances of clathrate hosts (in mole fraction) as a function of the alkanofer depth z: (a) for structure I clathrates; (b) for structure II clathrates. In both cases, the assumed liquid composition corresponds to the situation depicted in Fig. 3(b), for which the crust undergoes a geothermal gradient. This gradient, together with pressure variations, has been taken into account consistently. Schurmeier & Dombard (2018) found a thermal conductivity for organic-rich sand as low as 0.025 W m −1 K −1 . Hence, the discussion in Sect. 3 about Kalousova & Sotin's results, is also relevant in the case of an organics-dominated matrix.

Liquid-solid interaction
Interactions between trapped liquids and solid porous substrates may also be considered. In the case of an organic matrix, if the material resembles tholins produced in laboratories, the interaction should be weak since tholins seem to be very poorly soluble in cryogenic solvents (Carrasco et al. 2009). The situation is really different when the solid phase is dominated by water, which may form clathrate hydrates in the presence of molecules including methane, nitrogen, and ethane. This case has been considered many times in previous works (Lewis 1971;Lunine & Stevenson 1987;Loveday et al. 2001;Tobie et al. 2006;Choukroun & Sotin 2012;Kalousova & Sotin 2020;Vu et al. 2020a). The question that arises here is whether an enrichment of N 2 can inhibit, in some way, the clathration of CH 4 or C 2 H 6 . In search of an answer, we built an equilibrium model for clathrate, based on the van der Waals & Platteeuw (1959) approach. The clathration of species contained in the alkanofer liquid is measured by the fractional occupancy for a guest molecule i trapped in a clathrate crystal forming a cage of size S (i.e., small or large), with a structure of type T (type I orsII). The fractional occupancy is where the f i s are the fugacities of species in the liquid. These quantities have been evaluated according to the equation of state PC-SAFT already employed in this work. The factors C T S,i are the Langmuir constants of guest species i in a cage of size S and clathrate structure T. The expression of C T S,i is available in the literature (see, for instance, Thomas et al. 2008) and derived from Kihara potential (Sloan 1998). This theoretical approach enables the calculation of the mole fraction x Clath,T i of guest molecule i in a clathrate of structure type T: In Fig. 7, we represent these mole fractions for clathrate structures I and II, taking into account our "standard" alkanofer profile and including a geothermal gradient, as depicted in Fig. 3(b).
As we can see, the presence of N 2 in the clathrate cages remains mostly negligible. Except around z ∼ −3 km, the hydrocarbons CH 4 and C 2 H 6 are clearly the most abundant trapped molecules. For structure II, methane is the preponderant guest compound, although we observe an inversion of abundance with ethane at high depth, meaning deeper than ∼ 9 km. In the case of a structure I clathrate, ethane is the most abundant trapped species except at an interval between approximately 2000 and 4000 m. All these simulations assumed an equilibrium between the clathrate phase and the liquid, the latter being at the "diffusion equilibrium". Implicitly, the existence of a quantity of liquid large enough to fill the parent rock and clathrates, with respect to their porosities, is hypothesized. Of course, the straightforward approach described here would deserve a more detailed study. Finally, concerning the interaction solid-liquid, we cannot exclude the existence of more complex, and mainly unknown, physico-chemical processes leading to the formation of somewhat exotic species such as the co-crystals recently revealed by laboratory studies (Vu et al. , 2020aCable et al. 2014Cable et al. , 2019Cable et al. , 2020McConville et al. 2020;Czaplinski et al. 2020).

Possible implication for noble gases
In 2005, during the descent through Titan's atmosphere, the Gas Chromatograph Mass Spectrometer (GCMS) aboard the Huygens probe observed a strong depletion in noble gases (Niemann et al. 2005(Niemann et al. , 2010. For instance, the primordial argon isotope, 36 Ar, showed a ratio of 36 Ar/N 2 around 10 −7 , while Kr and Xe were not detected by the GCMS. Given the instrument detection threshold at 10 −8 , the atmospheric mole fractions of Kr and Xe should be below this upper limit. According to the literature, ratio in numbers Ar/N 2 , Kr/N 2 , and Xe/N 2 should be around 7 × 10 −2 , 5.6 × 10 −5 , and 3.9 × 10 −6 respectively, if solar values are adopted (Lodders 2003(Lodders , 2010(Lodders , 2019. Several scenarios have been investigated to explain this apparent lack of rare gases: some authors have considered the trapping in possible crustal clathrate layers. Osegovic & Max (2005), Thomas et al. (2007Thomas et al. ( , 2008, and Jacovi & Bar-Nun (2008) have tried to explain the depletion via a mechanism involving atmospheric haze, while Cordier et al. (2010) studied the dissolution in Titan's liquid phases. None of these works brought a fully satisfactory answer to the problem. Here, we can investigate whether noblegas atoms, which are small in size and heavier than N 2 or C 2 H 6 (Ar: 39.948 g mol −1 , Kr: 83.798 g mol −1 , Xe: 131.293 g mol −1 ), could accumulate in the deepest alkanofer layers. We performed preliminary computations by replacing N 2 from the ternary mixture and setting both surface mole fractions of the chosen noble gas "X" and ethane to 10 −3 , which is probably a huge value for a noble gas but a small one for C 2 H 6 . Adding the species X to our initial ternary mixture, N 2 +C 2 H 6 +CH 4 , would have required a large and in-depth modification of our model, then forming a quaternary mixture, this is the reason why we simply swapped N 2 and X. However, this simplified approach allows us to appreciate the global behavior of noble-gas species regarding their diffusion in a liquid solvent dominated by methane. As is evident from Table 3, all the derived vertical enrichment ratios (see definition given by Eq. 9) are significantly larger than 1. Noble Table 3. Vertical enrichment ratios ∆ X (X= Ar, Kr, Xe) for noble gases.

Ar
Kr Xe Model with an isothermal crust 1.64 11.5 119 Model including a geothermal gradient 1.95 8.98 33.0 gases tend to accumulate at the bottom of the alkanofer, leading to enrichment up to 2 orders of magnitude larger than the surface mole fraction fixed to 10 −3 in our test case. The tendency is clear, but a quantitative analysis would be required to solve the problem. The starting point of a follow-up study could be the development of an atmosphere-liquid phase equilibrium model based on the up-to-date PC-SAFT equation of state.

Conclusion
As stated in the introduction, the fate of ethane produced in the atmosphere poses a puzzling question. In this work, we addressed the possibility of an ethane storage in the deepest part of an alkanofer, due to the effect of molecular diffusion over geological timescales. This kind of scenario has been suggested by the existence of vertical chemical composition gradients through terrestrial hydrocarbon reservoirs. Our main hypothesis consists of the existence of an icy matrix, porous enough to enable molecular diffusion but with an hydraulic permeability below the threshold at which convective transport starts. We found that, even if ethane molecules are significantly heavier than their methane counterparts, C 2 H 6 does not accumulate at the bottom of the system. In addition, and somewhat surprisingly due to its relatively small molecular size, nitrogen seems to easily enrich deep layers. Many important alkanofer parameters are unknown; first of all, the composition and the total quantity of available liquids. The geographical and vertical extent of the reservoirs are not well constrained; the actual porosity, permeability, and connectivity of subsurface reservoirs are also poorly known. As a consequence, it is difficult, if not scientifically adventurous, to draw strong conclusions about the mentioned problem. However, arguments developed in this work clearly disfavor ethane accumulation in lower parts of subsurface reservoirs, under the action of molecular diffusion. Alternatively, if diffusion cannot account for the formation of significant chemical gradients, other transport processes may be more efficient. Several previous works have considered liquid subsurface flows through porous media (Hayes et al. 2008;Horvath et al. 2016). The scale taken into account may be regional (e.g., Horvath et al. 2016;Hayes et al. 2017) or global (Faulk et al. 2020); in these studies, the working fluid is chemically homogenous since it is liquid methane. The present investigations addresses the case of relatively small systems, meaning a few kilometers or a few tens of kilometers in size, hydrologically isolated to allow diffusion to be a dominant transport process. Future studies may focus on hydrodynamical mixing of local reservoirs with different chemical compositions.
The icy crust of Titan represents an interface between the internal global water ocean and the carbon-rich atmosphere where a plethora of organic molecules are produced. From an exobiological perspective, it is of interest to investigate the supply of molecules larger than C 2 -compounds to the underground aqueous environment. The deepest layers of alkanofers should have a low porosity caused by compaction, and then in these regions molecular diffusion should dominate the transport processes of these species. Even if the alkanofers have a maximum depth of a few kilometers, ice convection could subsequently bring this carbon-rich material down to the ocean (Kalousova & Sotin 2020) at a depth of about 80 − 100 km. Thus, the effect of molecular size on diffusion in liquid CH 4 +C 2 H 6 +N 2 mixtures should be explored.
To probe the possible polar Titan alkanofer would require techniques such as deep drilling, resistivity measurements (Griffiths & Barker 1993), and acoustic or electromagnetic sounding (Reynolds 2011). To date, one mission is planned to perform an in situ exploration of the surface of Titan: this is the octocopter Dragonfly 6 ). On board this probe, a scientific package called DraGMet 7 should contain instruments dedicated to geophysical measurements. However, to the best of our knowledge, the precise capabilities of these instruments are not publicly available. It would be surprising if they could probe the Titan subsurface layers down to a depth of several kilometers. Nonetheless, some particular circumstances could enable the search for clues concerning the behavior and the fate of liquids within underground layers. For example, Dragonfly will explore a region covered by large dunes and we can imagine a comparison of liquid composition analysis between samples collected at the base of a dune and at its top.

Data availability
We make a part of the data related to the present work publicly available on the Open Science platform Zenodo 8 , which belongs to the OpenAIRE project 9 . The dataset provided allows the reproduction of molecular dynamics simulations used in this work. Results from GROMACS molecular dynamics simulations performed on two binary mixtures, CH 4 (80%) + C 2 H 6 (20%) and CH 4 (80%) + N 2 (20%), and the ternary mixture CH 4 (60%) + C 2 H 6 (20%) + N 2 (20%), are available online 10 and it is also referenced with the DOI: 10.5281/zenodo.4975047. Each mixture is composed of 5000 molecules, and two values of temperature and pressure are considered, namely (T, p) = (90 K, 1.5 bar) and (95 K, 120 bar). This dataset collects information on trajectories (i.e., atomic positions and velocities, energies, and statistical data) and transport properties (diffusion and viscosity) for the six aforementioned cases (three mixtures under two (T, p) conditions). Our PC-SAFT implementation in FORTRAN 2008 has been made also publicly available on Github 11 and it has been archived on Zenodo 12 under the DOI: 10.5281/zenodo.

5085305.
Acknowledgements. The present research was supported by the Programme National de Planétologie (PNP) of CNRS-INSU co-funded by CNES, and also partially supported by the HPC center of Champagne-Ardenne (ROMEO) and the French HPC center CINES under project number A0090711976. D.C. thanks the University of Reims Champagne-Ardenne for granting him with the "dispositif Mobilité Courte Entrante et Sortante" for its stay at the Jet Propulsion Laboratory in 2019. Part of this work was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. This work has been done mainly using open-source softwares like Python, gfortran, kate, Jupyter Notebook, L A T E X under GNU/Debian Linux operating system, and also GROMACS under Red Hat Linux. The authors warmly acknowledge the whole Free Software community. . (B.14) The terms related to the barodiffusion are provided by Eqs. (46) and (47) of Ghorayeb & Firoozabadi (2000): The thermal diffusion ratio k T i, j (dimensionless) is a function of α T i, j , the thermal diffusion coefficient (dimensionless), according to the formula k T i, j = α T i, j x i x j . In all these equations, the density of the liquid ρ (kg m −3 ), used in c = ρ/M, the partial molar volumeV i (m 3 mol −1 ), the fugacities f i , and the thermal diffusion coefficient α T i, j are derived from the PC-SAFT equation of state. The equation that providesV i has been established in a previous work (see Eq. 19 in Cordier et al. 2019). The thermal diffusion coefficient α i j,T is derived from PC-SAFT thanks to the Haase formula (Haase 1969;Pan et al. 2006): where the M i s and x i s are, respectively, the molecular weights and the mole fractions, µ i is the chemical potential of component i, α 0 i j,T represents the thermal diffusion coefficient for the corresponding ideal gas of the couple of species (i, j), andh res i is the residual partial molar enthalpy of species i, also provided by PC-SAFT.