Issue 
A&A
Volume 653, September 2021



Article Number  A80  
Number of page(s)  14  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202140789  
Published online  10 September 2021 
Vertical compositional variations of liquid hydrocarbons in Titan’s alkanofers
^{1}
Université de Reims Champagne Ardenne, CNRS,
GSMA UMR CNRS 7331,
51097
Reims,
France
email: daniel.cordier@univreims.fr
^{2}
NASA Jet Propulsion Laboratory, California Institute of Technology,
Pasadena,
CA
91109,
USA
^{3}
Engineering Management of Additional Recovery, Mexican Petroleum Institute,
Eje Central Lázaro Cárdenas Norte 152,
07730
Mexico City,
Mexico
Received:
11
March
2021
Accepted:
22
June
2021
Context. According to clues left by the Cassini mission, Titan, one of the two Solar System bodies with a hydrologic cycle, may harbor liquid hydrocarbonbased analogs of our terrestrial aquifers, referred to as “alkanofers”.
Aims. 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.
Methods. We built numerical models including barodiffusion and thermodiffusion (Soret’s effect) in N_{2}+CH_{4}+C_{2}H_{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.
Results. 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 megayear to several tens of megayears, 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 waterice, we checked that N_{2} did not, in general, impede the clathration of C_{2}H_{6}, except in some layers. Interestingly, we found that noble gases could easily accumulate at the bottom of an alkanofer.
Key words: planets and satellites: composition / planets and satellites: surfaces / planets and satellites: interiors / equation of state / molecular processes
© D. Cordier et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 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üllerWodarg 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 “exooceanography”, connected to an “exohydrology” (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, 2016; MacKenzie et al. 2014) and underground alkanofers has been found (Corlies et al. 2017; Hayes et al. 2017). These latter geological structures were proposed to be the Titanian analogs for terrestrial aquifers. While a porous waterice 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 preCassini 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ñozIglesias 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 indirect 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 longstanding 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}), ethaneenriched 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 wellestablished 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 questionof a similar compositional grading in Titan’s crustal hydrocarbon reservoirs. Our approach does not imply largescale 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 indepth study of the role of molecular diffusion coefficients. The last section is dedicated to a general discussion about some chemophysical properties of alkanofers.
2 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 (kg m^{−2} s^{−1}) has N − 1 components, and can be written (Ghorayeb & Firoozabadi 2000), for a 1D system related to a vertical zaxis, as (1)
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 wellknown molecular fickean diffusion, driven by compositional gradients. The second 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 offugacities, 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: (2) (3)
with the s representing the elements of the diffusion tensor D^{M}. Similarly, the s and the 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: (4) (5)
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 (see below) or the density of the liquid ρ_{liq}, which depends on the chemical composition. All the detailed formulations of coefficients , and are given in Appendix B. throughout the paper, while usual molecular diffusion coefficients are noted with a “”. The temperature gradient sensitivity of molecular fluxes, represented by the s, may be explicitly written (6)
The terms a_{i3}, and 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 pressure derivative is provided by the hydrostatic equilibrium of the alkanofer: (7)
where ρ_{eff} (kg m^{−3}) is the effective density, which takes into account the matrix of waterice, 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 onmolecular diffusion (see also the results reported at the end of this section). Here, ρ_{liq} is computed in the frame of the PCSAFT^{1} (Gross & Sadowski 2001) equation of state (EoS), successfully used in several past studies dealing with the Titan context (Tan et al. 2013; LuspayKuti et al. 2015; Cordier et al. 2016, 2017). This equation of state is also employed to estimate the fugacity derivatives ∂ln f_{i} ∕∂x_{j} and the partial molar volumes , which appear in the expressions of s and s (see Appendix B). The waterice density ρ_{ice} is evaluated according to a dedicated I_{h} ice equation of state (Feistel & Wagner 2006).
The fickean molecular diffusion coefficients’ s are needed by s’ calculations. Given that we are dealing with liquids, we employed Wilke & Chang’s method (hereafter WC55 method; see Wilke & Chang 1955) combined with the Batschinski (1913), 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 estimations at P = 1.5 bar and T = 90 K, obtained with WC55 approach lead to m^{2} s^{−1} and 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 RungeKutta 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 temperature 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, (8)
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 orders of magnitude, according to bathymetric maps, derived from Cassini Synthetic Aperture Radar observations (Hayes 2016), the maximum depth of Titan’s seas should be around 200 m. In their subsurface liquid circulation model, Horvath et al. (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 welldefined 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 PCSAFT 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 PCSAFT 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_{ij}, are introduced to account for interspecies interactions. Their values were taken from previous papers (Cordier et al. 2016; Tan et al. 2013): , and . We checked, by setting all these k_{ij}s to zero, thatthe 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.
To proceed to quantification, we introduced a vertical enrichment ratio for a given species i: (9)
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, ranges between 34 and 45%, showing a clear enrichment. In contrast, ethane undergoes a lower enrichment with 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).
PCSAFT parameters.
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 pressurevaries 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 PCSAFT. 
3 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 solidstate 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 Eqs. (4) and (5), the geothermal gradient terms involve two coefficients: and , 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 ) and for C_{2} H_{6} in CH_{4} (denoted ), which are not available in the literature. In the general case, the thermal diffusion coefficient, α_{ij,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 α_{ij,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 PCSAFT via Haase’s formula (Haase 1969; Pan et al. 2006): (10)
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 represents the thermal diffusion coefficient for the corresponding ideal gas, and are the respective residual partial molar enthalpies of species (1) and (2). Even for most nonideal fluids (Pan et al. 2006), the coefficients can be easily computed by applying the kinetic theory of gases (Chapman & Cowling 1970). On their side, the residual enthalpies are directly provided by PCSAFT (see Eq. (A.46) of Gross & Sadowski 2001): (11)
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 using (12)
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 (13)
with k_{B} being the Boltzmann constant. The dimensionless ratio μ^{res}∕k_{B}T is directly provided by PCSAFT (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 and consequently (14)
as expected. In Fig. 3, for illustration purposes, we report thermal diffusion factors k_{ij,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 by k_{ij,T}, is larger than the sensitivity of ideal systems (panel a). Moreover, according to these estimations, ethane seems to be more sensitive than nitrogen: . If we consider only thermodiffusion in a binary system, the flux j_{1} of the species (1) along the vertical zaxis is given by (see Eq. (1) in Cordier et al. 2019) (15)
In the case of our icy crust model, we found ∂T∕∂z~ 0.6 K km^{−1}; since and , we have and . 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. 4, 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. 4b and c, corresponding to scenarios where nitrogen has a relatively high surfacemole fraction (i.e., ). The first scenario, corresponding to panels a in the same figure, show a different behavior due to the low value of (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, thecompetition 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 and 23.8 K km^{−1}. Included in our model, such high gradients reinforce the effect already seen in Fig. 4 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 and 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 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 and 120 K. As we can see in Fig. 5, the point , 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).
Fig. 2 Influence of PCSAFT 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. 
Fig. 3 (a) Variation of as a functionof the mole fraction x_{1} of the first compound. The quantity 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. 
Fig. 4 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’salkanofer 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 waterice EoS is provided by Feistel & Wagner (2006), while the thermodynamic properties of liquids are computed using PCSAFT. For comparison, the results of Fig. 1 are added in thin gray lines. 
Fig. 5 Phase diagrams of the N_{2}+CH_{4} system at 110, 115, 120, and 125 K. Circles represent laboratory measurements from Parrish & Hiza (1974), curves are computed according to the PCSAFT EoS (no experimental data are available at 125 K). Solid lines are the liquidus, while dashed lines represent the vaporus. The star represents the case where and P = 14 bar, which is tobe compared with the point marked by a triangle (see text). This figure is similar to Fig. 3b of Tan et al. (2013). 
4 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, 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: and , corresponding, respectively, to the diffusion of N_{2} and C_{2} H_{6} in liquid methane. In models described in Sects. 2 and 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.14)) leads to , and , respectively. As a consequence, the ratio L_{12}∕L_{11} is independent of all molecular diffusion coefficients, whereas L_{21}∕L_{22} depends on since . Therefore, the terms and are functions of the ratio and of , respectively.This way, the equations that structure the diffusive properties of our model, namely Eqs. (4) and (5), have their barodiffusion and thermodiffusion terms as functions of . Nonetheless, these dependencies are different since and have different behaviors regarding molecular diffusion coefficients (see Eqs. (B.15)–(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 and 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 going from − 42% (for our “standard” model) to − 39% (case ) or − 28% (case ).
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 wellknown opensource package GROMACS^{2}^{,}^{3} 2018 version (Abraham et al. 2015). Nitrogen molecules have been described by the TraPPE^{4} force field (Potoff & Siepmann 2001). Alkanes (i.e., methane and ethane) were modeled by adopting the TraPPEUA^{5} coarsegrained 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 1ns equilibration step in the NV T ensemble; (2) a 19ns equilibration step in the NpT ensemble, where the temperature and pressure are fixed to the expected target values; (3) a 10 nsaccumulation step in the NpT 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 from MD simulations, with estimations performed with Wilke & Chang’s (WC) method (see Sect. 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.
Diffusion coefficients derived from our MD simulations based on TraPPEUA force fields.
5 Discussion and conclusion
5.1 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 (16)
with L (m) being a characteristic length of the studied system, and (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 megayear, and τ_{diff} ~ 30 megayears for a 1 km deep reservoir. These 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 (17)
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 (18)
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 (19)
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 largescale 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 (20)
where the ratio η∕ρg represents the properties of the fluid, in the context of Titan (g), for this term wefound 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 (21)
Thus, the effect of diffusion is not erased by interreservoir seepage if, roughly, K_{p, sc} ≲ 10^{−5}∕Δt. For a diffusion timescale, for example Δt ~ 1 megayear (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.
5.2 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 () (Bories & Combarnous 1973): (22)
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}), 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 (Nield & Bejan 2017). The volumetric coefficient of thermal expansion α_{V} is given by (23)
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 PCSAFT. 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 abovementioned 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 isothermalisobaric (NpT) ensemble, have been processed according to relevant fluctuation formulas, leading to results similar to those obtained with PCSAFT 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 (i.e., (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), leading to d ≃ 40 μm. Of course, this kind of estimation must be regarded with caution, since empirical laws established for the Earth case may no longer be valid for Titan, even if it provides a plausible value. Equation (22) clearly shows that low permeability K_{p, m}, small depths H, and gentle thermal gradients ΔT∕H are 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, Schurmeier & Dombard (2018) found a thermal conductivity for organicrich 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 organicsdominated matrix.
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: . The role of nitrogen abundance (red curve) has been studied by keeping ethane mole fraction constant: . 
5.3 Liquidsolid 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 (24)
where the f_{i}s are the fugacities of species in the liquid. These quantities have been evaluated according to the equation of state PCSAFT already employed in this work. The factors are the Langmuir constants of guest species i in a cage of size S and clathrate structure T. The expression of 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 of guest molecule i in a clathrate of structure type T: (25)
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. 4b. 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 solidliquid, we cannot exclude the existence of more complex, and mainly unknown, physicochemical processes leading to the formation of somewhat exotic species such as the cocrystals recently revealed by laboratory studies (Vu et al. 2014, 2020a,b; Cable et al. 2014, 2019, 2020; McConville et al. 2020; Czaplinski et al. 2020).
Fig. 7 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. 4b, for which the crust undergoes a geothermal gradient. This gradient, together with pressure variations, has been taken into account consistently. 
Vertical enrichment ratios Δ_{X} (X = Ar, Kr, Xe) for noble gases.
5.4 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, 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 atmosphericmole 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, 2010, 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. (2007, 2008), and Jacovi & BarNun (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 indepth 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 noblegas 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 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 followup study could be the development of an atmosphereliquid phase equilibrium model based on the uptodate PCSAFT equation of state.
5.5 Conclusion
As stated inthe 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 clearlydisfavor 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 homogenoussince 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 carbonrich 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 carbonrich 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} (Lorenz et al. 2018). 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.
5.6 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 datasetprovided 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: https://doi.org/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 PCSAFT implementation in FORTRAN 2008 has been made also publicly available on Github^{11} and it has been archived on Zenodo^{12} under the DOI: https://doi.org/10.5281/zenodo.5085305.
Acknowledgements
The present research was supported by the Programme National de Planétologie (PNP) of CNRSINSU cofunded by CNES, and also partially supported by the HPC center of ChampagneArdenne (ROMEO) and the French HPC center CINES under project number A0090711976. D.C. thanks the University of Reims ChampagneArdenne 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 opensource softwares like Python, gfortran, kate, Jupyter Notebook, LATEX under GNU/Debian Linux operating system, and also GROMACS under Red Hat Linux. The authors warmly acknowledge the whole Free Software community.
Appendix A List of variables
P and T: local pressure (Pa) and temperature (K).
z: depth below Titan’s surface (m).
x_{i}: local mole fraction of compound i.
M_{i}: molecular mass of species i, mean molecular mass of considered material (kg mol^{−1}).
: partial molar volume of component i (m^{3} mol^{−1}).
f_{i}: fugacity of component i (Pa).
: fickean molecular diffusion coefficient (m^{2} s^{−1}) of compound i in species j.
D^{M} (kg mol^{−1} m^{2} s^{−1}), D^{P} (kg mol^{−1} m^{2} s^{−1} Pa^{−1}) and D^{T} (kg mol^{−1} m^{2} s^{−1} K^{−1}): tensors introduced by Ghorayeb & Firoozabadi (2000) are not to be confused with the ’s.
g, g_{Tit}, g_{Earth}: gravity (generic notation, Titan’s value, Earth’s value, m s^{−2}).
k_{ij,T}: thermal diffusion ratio (dimensionless).
α_{ij,T}: thermal diffusion coefficient, k_{ij,T} = x_{i}x_{j}α_{ij,T} (dimensionless).
ρ_{eff}: effective density (kg m^{−3}) of the alkanofer material, ρ_{eff} = Π ρ_{liq} + (1 − Π) ρ_{ice}, with Π being the porosity, ρ_{liq} the density of the liquid (kg m^{−3}), and ρ_{ice} the density of the icy matrix (kg m^{−3}).
μ_{i}: chemical potential of component i (J mol^{−1}), is the residual chemical potential of component i (J mol^{−1}) the term “residual” relates to a difference between the physical quantities for the actual fluid and the corresponding ideal gas.
h_{i}: residual molar enthalpy of component i (J mol^{−1}).
: residual partial molar enthalpy of component i (J mol^{−1}).
ã^{res}: residual and reduced (i.e., molar) Helmholtz free energy (J mol^{−1}).
Z: compressibility factor of the considered fluid (dimensionless).
K_{p,sc}: permeability (m 2) of the porous media that connects two reservoirs.
η: dynamic viscosity (Pa s); ν is the kinematic viscosity (η = ρν).
Δ_{i}: vertical enrichment ratios (dimensionless) for species i.
Appendix B Model of molecular, pressure, and thermal diffusion in a nonideal ternary mixture
Our model is essentially based on the formalism described in Ghorayeb & Firoozabadi (2000) for nonideal multicomponent mixtures. In this work, for the purpose of clarity we mainly kept the notation adopted in the mentioned reference. Thus, the derivative of the logarithm of the fugacity f_{i} of component i, with respect to the mole fraction x_{j} of the species j, is (B.1)
The fugacity f_{i} is evaluated according to the PCSAFT equation of state, already used in many Titanrelated papers (e.g., Tan et al. 2013; Cordier et al. 2017). Generally speaking, the fugacity may be written as f_{i} = Φ_{i}x_{i}P, where Φ_{i} is the fugacity coefficient directly provided by PCSAFT, and P stands for the pressure. The derivatives are obtained by finite differences. Employing a very usual notation, M_{i} represents the molecular mass of species i, while the average molecular mass is noted . In the following, we also use the reduced expression (B.2)
In a usual way, the fickean molecular diffusion coefficients of species i in component j are denoted (m^{2} s^{−1}). The gas constant is simply written as R. Four constants are introduced (see Eqs. (36)–(39), in Ghorayeb & Firoozabadi 2000): (B.3) (B.4) (B.5) (B.6)
The elements of the molecular diffusion matrix D^{M} are given by(see Eqs. (40)–(43), in Ghorayeb & Firoozabadi 2000) (B.8) (B.9) (B.10) (B.11)
The L_{ij}s are the socalled Onsager phenomenological coefficients given by (B.12) (B.13)
Using Eq. (49) of Ghorayeb & Firoozabadi (2000) for i = 1 and j = 2, we can derive (B.14)
The terms related to the barodiffusion are provided by Eqs. (46) and (47) of Ghorayeb & Firoozabadi (2000): (B.15) (B.16)
Finally, concerning the thermal diffusion, we have (see Eqs. (44) and ()45) in Ghorayeb & Firoozabadi 2000): (B.17) (B.18)
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 , the partial molar volume (m^{3} mol^{−1}), the fugacities f_{i}, and the thermal diffusion coefficient α_{T i,j} are derived from the PCSAFT equation of state. The equation that provides has been established in a previous work (see Eq. (19) in Cordier et al. 2019). The thermal diffusion coefficient α_{ij,T} is derived from PCSAFT thanks to the Haase formula (Haase 1969; Pan et al. 2006): (B.19)
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, represents the thermal diffusion coefficient for the corresponding ideal gas of the couple of species (i, j), and is the residual partial molar enthalpy of species i, also provided by PCSAFT.
Appendix C Structure of the crust
Following the Titan interior model published by Sohl et al. (2014), the properties of the icy crust may be retrieved by integrating (C.1) (C.2) (C.3)
where P, T, and r are, respectively, the pressure (Pa), the temperature (K), and the distance from Titan’s center of mass (m). The density of the medium (kg m^{−3}) is denoted ρ as usual, k_{c} represents the thermal conductivity (W m^{−1} K^{−1}) of the material, while is the local Nusselt number, which parameterizes the strength of the convective heat flux. is given by (see Eq. (7) in Sohl et al. 2014) (C.4)
where the “adiabatic gradient” is provided by (C.5)
and the “conductive gradient” may be defined as (C.6)
According to the mixing length theory (Kamata 2018), if (C.7)
otherwise, k_{v} = 0. In Eqs. (C.5) and (C.8), α, C_{P}, and η are, respectively, the thermal expansion coefficient (K^{−1}) of waterice I_{h}, the specific heat capacity (J kg^{−1} K^{−1}), and the viscosity (Pa s). The parameter l has been taken equal to the local pressure scale height H_{P} which is evaluated by (C.9)
The viscosity of the crust η is essentially unknown. Sohl et al. (2014) chose a value (10^{20} Pa s) large enough to impede solid convection, since arguments have been put forward concerning the absence of convection in Titan’s crust (Nimmo & Bills 2010). The system of ordinary differential equations formed by Eqs. (C.1)–C.3 can easily be solved by a RungeKutta method, taking the values at the surface as boundary conditions: P = P_{s} = 1.5 bar, q = q_{s} = 3.11 × 10^{−3} W m^{−2}, and T = T_{s} = 94 K (Sohl et al. 2014).
References
 Abraham, M. J., Murtola, T., Schulz, R., et al. 2015, SoftwareX, 12, 19 [Google Scholar]
 Baldwin, H. L., & McGuinness, C. L. 1976, A Primer on Ground Water, Technical Report, United States Departement of the Interior [Google Scholar]
 Barnes, J. W., Bow, J., Schwartz, J., et al. 2011, Icarus, 216, 136 [Google Scholar]
 Batschinski, A. J. 1913, Z. Physik. Chim., 84, 643 [Google Scholar]
 Bear, J. 1972, Dynamics of Fluids in Porous Media (New York: Dover Publications) [Google Scholar]
 Bories, S. A., & Combarnous, M. A. 1973, J. Fluid Mech., 57, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Cable, M. L., Vu, T. H., Hodyss, R., et al. 2014, Geophys. Res. Lett., 41, 5396 [NASA ADS] [CrossRef] [Google Scholar]
 Cable, M. L., Vu, T. H., Malaska, M. J., et al. 2019, ACS Earth Space Chem., 3, 2808 [NASA ADS] [CrossRef] [Google Scholar]
 Cable, M. L., Vu, T. H., Malaska, M. J., et al. 2020, ACS Earth Space Chem., 4, 1375 [NASA ADS] [CrossRef] [Google Scholar]
 Carrasco, N., SchmitzAfonso, I., Bonnet, J.Y., et al. 2009, J. Phys. Chem. A, 113, 11195 [CrossRef] [Google Scholar]
 Chapman, S., & Cowling, T. G. 1970, The Mathematetical Theory of NonUniform Gases (Cambridge University Press, Cambridge) [Google Scholar]
 Chilingar, G. V., Buryakovsky, L. A., Eremenko, N. A., & Gorfunkel, M. V. 2005, Geology and Geochemistry of Oil and Gas (Amsterdam: Elsevier) [Google Scholar]
 Choukroun, M., & Sotin, C. 2012, Geophys. Res. Lett., 39, 4201 [Google Scholar]
 Choukroun, M., Grasset, O., Tobie, G., & Sotin, C. 2010, Icarus, 205, 581 [NASA ADS] [CrossRef] [Google Scholar]
 Cordier, D., Mousis, O., Lunine, J. I., et al. 2010, ApJ, 721, L117 [NASA ADS] [CrossRef] [Google Scholar]
 Cordier, D., Barnes, J. W., & Ferreira, A. G. 2013, Icarus, 226, 1431 [NASA ADS] [CrossRef] [Google Scholar]
 Cordier, D., Cornet, T., Barnes, J. W., et al. 2016, Icarus, 270, 41 [Google Scholar]
 Cordier, D., GarcíaSánchez, F., JustoGarcía, D. N., & LigerBelair, G. 2017, Nat. Astron., 1, 0102 [NASA ADS] [CrossRef] [Google Scholar]
 Cordier, D., Bonhommeau, D. A., Port, S., et al. 2019, ApJ, 880, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Corlies, P., Hayes, A. G., Birch, S. P. D., et al. 2017, Geophys. Res. Lett., 44, 11 [Google Scholar]
 Coutelier, M., Cordier, D., Seignovert, B., et al. 2021, Icarus, 364, 114464 [Google Scholar]
 Czaplinski, E., Yu, X., Dzurilla, K., & Chevrier, V. 2020, Planet. Sci. J., 1, 76 [NASA ADS] [CrossRef] [Google Scholar]
 Darcy, R. P. G. 1856, Mémoires Présentés à l’Académie des Sciences de l’Institut de France, 14, 141 [Google Scholar]
 Elachi, C., Wall, S., Allison, M., et al. 2005, Science, 308, 970 [NASA ADS] [CrossRef] [Google Scholar]
 Eshleman, V. R., Lindal, G. F., & Tyler, G. L. 1983, Science, 221, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Espósito, R. O., Rodrigues Alijó, P. H., Scilipoti, J. A., & Tavares, F. W. 2017, Compositional Grading in Oil and Gas Reservoirs (Oxford, United Kingdom: Gulf Professional Publishing) [Google Scholar]
 Faulk, S. P., Lora, J. M., Mitchell, J. L., & Milly, P. C. D. 2020, Nat. Astron., 4, 390 [NASA ADS] [CrossRef] [Google Scholar]
 Feistel, R., & Wagner, Y. 2006, J. Phys. Chem. Ref. Data, 35 [Google Scholar]
 Flasar, F. M. 1983, Science, 221, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Flasar, F. M., Samuelson, R. E., & Conrath, B. J. 1981, Nature, 292, 693 [NASA ADS] [CrossRef] [Google Scholar]
 Fulchignoni, M., Ferri, F., Angrilli, F., et al. 2005, Nature, 438, 785 [NASA ADS] [CrossRef] [Google Scholar]
 Galliero, G., Bataller, H., Bazile, J.P., et al. 2017, npj Microgravity, 3, 20 [Google Scholar]
 Ghorayeb, K., & Firoozabadi, A. 2000, AIChE J., 46, 883 [Google Scholar]
 Ghorayeb, K., Firoozabadi, A., & Anraku, T. 2003, SPE J., 8, 114 [CrossRef] [Google Scholar]
 Gilliam, A. E., & Lerman, A. 2016, Icarus, 275, 252 [Google Scholar]
 Griffiths, D. H., & Barker, R. D. 1993, J. Appl. Geophy., 29, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Griffith, C. A., Lora, J. M., Turner, J., et al. 2012, Nature, 486, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Gross, J., & Sadowski, G. 2001, Ind. Eng. Chem. Res., 40, 1244 [CrossRef] [Google Scholar]
 Haase, R. 1969, Thermodynamics of Irreversible Processes (Reading, MA, USA: Addison Wesley), 201 [Google Scholar]
 Hayes, A. G. 2016, Ann. Rev. Earth Planet Sci., 44, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Hayes, A., Aharonson, O., Callahan, P., et al. 2008, Geophys. Res. Lett., 35, 9204 [NASA ADS] [CrossRef] [Google Scholar]
 Hayes, A. G., Birch, S. P. D., Dietrich, W. E., et al. 2017, Geophys. Res. Lett., 44, 11 [Google Scholar]
 Hildebrand, J. H. 1971, Science, 174, 490 [NASA ADS] [CrossRef] [Google Scholar]
 Hofgartner, J. D., Hayes, A. G., Campbell, D. B., et al. 2020, Nat. Commun., 11, 2829 [NASA ADS] [CrossRef] [Google Scholar]
 Horvath, D. G., AndrewsHanna, J. C., Newman, C. E., Mitchell, K. L., & Stiles, B. W. 2016, Icarus, 277, 103 [Google Scholar]
 Hu, W. R., Zhao, J. F., Long, M., et al. 2014, Microgravity Sci. Technol., 26, 159 [NASA ADS] [CrossRef] [Google Scholar]
 Jacovi, R., & BarNun, A. 2008, Icarus, 196, 302 [NASA ADS] [CrossRef] [Google Scholar]
 Jennings, D. E., Cottini, V., Nixon, C. A., et al. 2016, ApJ, 816, L17 [CrossRef] [Google Scholar]
 Kalousova, K.,& Sotin, C. 2020, Geophys. Res. Lett., 47, e2020GL087481 [NASA ADS] [CrossRef] [Google Scholar]
 Kamata, S. 2018, J. Geophys. Res. Planets, 123, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Kossacki, K. J., & Lorenz, R. D. 1996, PSS, 44, 1029 [NASA ADS] [Google Scholar]
 Kuiper, G. P. 1944, ApJ, 100, 378 [NASA ADS] [CrossRef] [Google Scholar]
 Kumar, P., & Chevrier, V. F. 2020, ACS Earth Space Chem., 4, 241 [NASA ADS] [CrossRef] [Google Scholar]
 Latini, G., Passerini, G., & Polonara, F. 1996, Int. J. Thermophys., 17, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Le Gall, A., Janssen, M. A., Paillou, P., et al. 2010, Icarus, 207, 948 [CrossRef] [Google Scholar]
 Lewis, J. S. 1971, Icarus, 15, 174 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, Z. Y.C., Radebaugh, J., Harris, R. A., et al. 2016, Icarus, 270, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Lodders, K. 2003, ApJ, 591, 1220 [Google Scholar]
 Lodders, K. 2010, ApSSP, 16, 379 [NASA ADS] [Google Scholar]
 Lodders, K. 2019, ArXiv eprints [arXiv:1912.00844] [Google Scholar]
 Lorenz, R. D., Mitchell, K. L., Kirk, R. L., et al. 2008, Geophys. Res. Lett., 35, L02406 [NASA ADS] [CrossRef] [Google Scholar]
 Lorenz, R. D., Turtle, E. P., Barnes, J. W., et al. 2018, Dragonfly: A Rotorcraft Lander Concept for Scientific Exploration at Titan, Technical Repoert 3, Johns Hopkins APL [Google Scholar]
 Loveday, J. S., Nelmes, R. J., Guthrie, M., et al. 2001, Nature, 410, 661 [CrossRef] [PubMed] [Google Scholar]
 Lowrie, W. 2010, Fundamentals of Geophysics (Cambridge UK: Cambridge University Press) [Google Scholar]
 Lunine, J. I., & Stevenson, D. J. 1987, Icarus, 70, 61 [NASA ADS] [CrossRef] [Google Scholar]
 LuspayKuti, A., Chevrier, V. F., Cordier, D., et al. 2015, EPSL, 410C, 75 [Google Scholar]
 MacKenzie, S. M., Barnes, J. W., Sotin, C., et al. 2014, Icarus, 243, 191 [Google Scholar]
 MacKenzie, S. M., Barnes, J. W., Hofgartner, J. D., et al. 2019, Nat. Astron., 3, 506 [NASA ADS] [CrossRef] [Google Scholar]
 Mansi, S. S., Siepmann, J. I., & Tsapatsis, M. 2017, AIChE J., 63, 5098 [Google Scholar]
 Martin, M. G., & Siepmann, J. I. 1998, J. Phys. Chem. B, 102, 2569 [CrossRef] [Google Scholar]
 McConville, C. A., Tao, Y., Evans, H. A., et al. 2020, Chem. Commun., 56, 13520 [Google Scholar]
 Metcalfe, R. S., Vogel, J. L., & Morris, R. W. 1988, SPE Reservoir Eng., 3, 1025 [Google Scholar]
 Mousis, O., Choukroun, M., Lunine, J. I., & Sotin, C. 2014, Icarus, 239, 39 [Google Scholar]
 Mousis, O., Lunine, J. I., Hayes, A. G., & Hofgartner, J. D. 2016, Icarus, 270, 37 [NASA ADS] [CrossRef] [Google Scholar]
 MuñozIglesias, V., Choukroun, M., Vu, T. H., et al. 2018, ACS Earth Space Chem., 2, 135 [Google Scholar]
 MüllerWodarg, I., Griffith, C. A., Lellouch, A., & Cravens, T. E. E. 2014, TITAN – Interior, Surface, Atmosphere, and Space Environment, eds. I. MüllerWodarg, C. A. Griffith, A. Lellouch, & T. E. Cravens (New York: Cambridge University Press) [Google Scholar]
 Nazki, M. A., Chopra, T., & Chandrappa, A. K. 2020, Constr. Build. Mater., 238, 117693 [Google Scholar]
 Nield, D. A., & Bejan, A. 2017, Convection in Porous Media, 5th edn. (Cham, Switzerland: Springer) [Google Scholar]
 Niemann, H. B., Atreya, S. K., Bauer, S. J., et al. 2005, Nature, 438, 779 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Niemann, H. B., Atreya, S. K., Demick, J. E., et al. 2010, J. Geophys. Res., 115, E12006 [Google Scholar]
 Nimmo, F., & Bills, B. G. 2010, Icarus, 208, 896 [NASA ADS] [CrossRef] [Google Scholar]
 Nixon, C. A., Lorenz, R. D., Achterberg, R. K., et al. 2018, Planet. Space Sci., 155, 50 [CrossRef] [Google Scholar]
 Nougier, J. P. 1987, Méthodes de Calcul Numérique (Paris: Masson) [Google Scholar]
 Obidi, O. C. 2014, PhD thesis, Imperial College London, UK [Google Scholar]
 Osegovic, J. P., & Max, M. D. 2005, J. Geophys. Res., 110, E08004 [NASA ADS] [Google Scholar]
 Owen, T. 1982, Sci. Am., 246, 98 [NASA ADS] [CrossRef] [Google Scholar]
 Pan, S., Jiang,C., Yan, Y., Kawaji, M., & Saghir, Z. 2006, J. NonEquilib. Thermodyn., 41, 47 [Google Scholar]
 Parrish, W. R., & Hiza, M. J. 1974, Adv. Cryog. Eng., 19, 300 [Google Scholar]
 Petuya, C., Choukroun, M., Vu, T. H., Sotin, C., & Davies, A. G. 2020, ACS Earth Space Chem., 4, 526 [NASA ADS] [CrossRef] [Google Scholar]
 Poling, B. E., Prausnitz, J. M., & O’Connell, J. 2007, The Properties of Gases and Liquids, 5th edn. (Englewood Cliffs: McGrawHill Professional) [Google Scholar]
 Pollack, H. N., Hurter, S. J., & Johnson, J. R. 1993, Rev. Geophys., 31, 267 [NASA ADS] [CrossRef] [Google Scholar]
 Potoff, J. J., & Siepmann, J. I. 2001, AIChE J., 47, 1676 [CrossRef] [Google Scholar]
 Powers, R. W., Mattox, R. W., & Johnston, H. L. 1954, J. Am. Chem. Soc., 76, 5968 [Google Scholar]
 Reynolds, J. M. 2011, An Introduction to Applied and Environmental Geophysics, 2nd edn. (Chichester, UK: John Wiley & Sons) [Google Scholar]
 Sagan, C., & Dermott, S. F. 1982, Nature, 300, 731 [NASA ADS] [CrossRef] [Google Scholar]
 Schurmeier, L. R., & Dombard, A. J. 2018, Icarus, 305, 314 [Google Scholar]
 Shah, M. S., Siepmann, J. I., & Tsapatsis, M. 2017, AIChE J., 63, 5098 [Google Scholar]
 Sloan, E. D. 1998, Clathrates Hydrates of Natural Gases (New York: Marcel Decker) [Google Scholar]
 Sohl, F., Solomonidou, A., Wagner, F. W., et al. 2014, J. Geophys. Res. Planets, 119, 1013 [NASA ADS] [CrossRef] [Google Scholar]
 Soret, C. 1879, Arch. Sci. Phys. Nat., 2, 48 [Google Scholar]
 Sotin, C., Lawrence, K. J., Reinhardt, B., et al. 2012, Icarus, 221, 768 [Google Scholar]
 Stephan, K., Jaumann, R., Brown, R. H., et al. 2010, Geophys. Res. Lett., 37, L07104 [Google Scholar]
 Stofan, E. R., Elachi, C., Lunine, J. I., et al. 2007, Nature, 445, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Tan, S. P., Kargel, J. S., & Marion, G. M. 2013, Icarus, 222, 53 [Google Scholar]
 Temeng, K. O., AlSadeg, M. J., & AlMulhim, W. A. 1998, Compositional Grading in the Ghawar Khuff Reservoirs (Society of Petroleum Engineers), 685 [Google Scholar]
 Thomas, C., Mousis, O., Ballenegger, V., & Picaud, S. 2007, A&A, 474, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Thomas, C., Picaud, S., Mousis, O., & Ballenegger, V. 2008, PSS, 56, 1607 [Google Scholar]
 Tobie, G., Lunine, J. I., & Sotin, C. 2006, Nature, 440, 61 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Tobie, G., Gautier, D., & Hersant, F. 2012, ApJ, 752, 125 [CrossRef] [Google Scholar]
 Turtle, E. P., Perry, J. E., Barbara, J. M., et al. 2018, Geophys. Res. Lett., 45, 5320 [NASA ADS] [CrossRef] [Google Scholar]
 Turtle, E. P., Trainer, M. G., Barnes, J. W., et al. 2020, in 51st Lunar and Planetary Science Conference (Houston: Lunar and Planetary Institute), Abstract #2288 [Google Scholar]
 Tyler, G. L., Eshleman, V. R., Anderson, J. D., et al. 1981, Science, 212, 201 [NASA ADS] [CrossRef] [Google Scholar]
 van der Waals, J. H., & Platteeuw, J. C. 1959, Clathrate Solutions. In: Advances in Chemical Physics, (New York: Interscience), 2, 1 [Google Scholar]
 Vogel, E., & Weiss, A. 1981, Ber. Bunsenges. Phys. Chem, 85, 539 [Google Scholar]
 Vu, T. H., Cable, M. L., Choukroun, M., Hodyss, R. & Beauchamp, P. 2014, J. Phys. Chem. A, 118, 4087 [NASA ADS] [CrossRef] [Google Scholar]
 Vu, T. H., Choukroun, M., Sotin, C., & MuñozIglesias, V. MaynardCasely, H. E. 2020a, Geophys. Res. Lett., 47, e2019GL086265 [Google Scholar]
 Vu, T. H., MaynardCasely, H. E., Cable, M. L., et al. 2020b, J. Appl. Crystallogr., 53, 1524 [Google Scholar]
 Wilke, C. R., & Chang, P. 1955, AIChE J., 1, 264 [Google Scholar]
 Wood, C. A., & Radebaugh, J. 2020, J. Geophys. Res., 125, e06036 [NASA ADS] [CrossRef] [Google Scholar]
 Yung, Y. L., Allen, M., & Pinto, J. P. 1984, ApJS, 55, 465 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
All Tables
Diffusion coefficients derived from our MD simulations based on TraPPEUA force fields.
All Figures
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 pressurevaries 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 PCSAFT. 

In the text 
Fig. 2 Influence of PCSAFT 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. 

In the text 
Fig. 3 (a) Variation of as a functionof the mole fraction x_{1} of the first compound. The quantity 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. 

In the text 
Fig. 4 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’salkanofer 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 waterice EoS is provided by Feistel & Wagner (2006), while the thermodynamic properties of liquids are computed using PCSAFT. For comparison, the results of Fig. 1 are added in thin gray lines. 

In the text 
Fig. 5 Phase diagrams of the N_{2}+CH_{4} system at 110, 115, 120, and 125 K. Circles represent laboratory measurements from Parrish & Hiza (1974), curves are computed according to the PCSAFT EoS (no experimental data are available at 125 K). Solid lines are the liquidus, while dashed lines represent the vaporus. The star represents the case where and P = 14 bar, which is tobe compared with the point marked by a triangle (see text). This figure is similar to Fig. 3b of Tan et al. (2013). 

In the text 
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: . The role of nitrogen abundance (red curve) has been studied by keeping ethane mole fraction constant: . 

In the text 
Fig. 7 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. 4b, for which the crust undergoes a geothermal gradient. This gradient, together with pressure variations, has been taken into account consistently. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.