© D. Cordier et al. 2021

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ü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, 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 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 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 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⁻¹) is heavier than nitrogen (28 g mol⁻¹) and methane (16 g mol⁻¹), 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 questionof 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.

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⁻² s⁻¹) has N − 1 components, and can be written (Ghorayeb & Firoozabadi 2000), for a 1D system related to a vertical z-axis, as (1)

with c being the total molar density (mol m⁻³), 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₁, x₂, …, 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 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₂, C₂ H₆ and CH₄. 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₀ = 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⁻³) is the effective density, which takes into account the matrix of water-ice, and g_Tit (m s⁻²) 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 PC-SAFT¹ (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. 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 water-ice 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² s⁻¹ and m² s⁻¹; these values around 10⁻⁹ m² s⁻¹ 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 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⁻², 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 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₂: 28.0134 g mol⁻¹) seems to be more sensitive to gravitational effect than ethane (C₂ H₆: 30.0690 g mol⁻¹), 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_ij, are introduced to account for inter-species interactions. Their values were taken from previous papers (Cordier et al. 2016; Tan et al. 2013): , and . We checked, by setting all these k_ijs to zero, thatthe interaction parameters have no influence on the obtained abundances profiles. Similarly, by increasing the value of ϵ(N₂) 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₂ and C₂ H₆. 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₂ than in C₂ H₆, 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).

Fig. 1 Vertical variations of mole fractions (denoted xi) of the main components (i.e., N2, C2 H6 and CH4) of Titan’s alkanofer liquid. This figure depicts a scenario where the icy crust is assumed to be isothermal at the temperature T0 = 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) N2: 0.05, C2 H6: 0.20, CH4: 0.75; (b) N2: 0.20, C2 H6: 0.20, CH4: 0.60; (c) N2: 0.20, C2 H6: 0.30, CH4: 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.

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⁻², 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⁻¹, well below the terrestrial geothermal gradient that has a commonly accepted value around 30 K km⁻¹ (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₂ in CH₄ (denoted ) and for C₂ H₆ in CH₄ (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,Ts 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): (10)

where the M_is and x_i s are, respectively, the molecular weights and the mole fractions, and μ₁ 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 PC-SAFT (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_BT 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 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₂ -CH₄ and C₂ H₆-CH₄. Except for x₁ 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₁ of the species (1) along the vertical z-axis 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⁻¹; 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 surface-mole 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⁻¹. 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₂, 20% of C₂ H₆, and 60% of CH₄, and assuming ∂T∕∂z ~−23.8 K km⁻¹ (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₂ +CH₄ 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 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 N2 and C2 H6 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 x1 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 C2 H6-CH4 system, while the dashed line is related to N2-CH4. 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 xi) of the main components (i.e., N2, C2 H6 and CH4) 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) N2: 0.05, C2H6: 0.20, CH4: 0.75; (b) N2: 0.20, C2H6: 0.20, CH4: 0.60; (c) N2: 0.20, C2H6: 0.30, CH4: 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.

Fig. 5 Phase diagrams of the N2+CH4 system at 110, 115, 120, and 125 K. Circles represent laboratory measurements from Parrish & Hiza (1974), curves are computed according to the PC-SAFT 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₂ and C₂ H₆ 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₁₂∕L₁₁ is independent of all molecular diffusion coefficients, whereas L₂₁∕L₂₂ 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³, 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₂ and C₂ H₆ in CH₄ 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). Nitrogen molecules have been described by the TraPPE⁴ force field (Potoff & Siepmann 2001). Alkanes (i.e., methane and ethane) were modeled by adopting the TraPPE-UA⁵ 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 NV T ensemble; (2) a 19-ns equilibration step in the NpT ensemble, where the temperature and pressure are fixed to the expected target values; (3) a 10 ns-accumulation step in the NpT ensemble. In our work, standard simulations were performed with a total number of molecules fixed to 5000 with 4000 CH₄ (80%) and 1000 N₂ or C₂ H₆ (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.

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² s⁻¹) a diffusion coefficient. According to values reported in Table 2, the diffusion coefficient is on the order of ~ 10⁻⁹ m² s⁻¹. 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 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³) 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²) 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⁻³), g is the gravity (m s⁻²), η 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₂ +CH₄+C₂H₆. 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 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 (20)

where the ratio η∕ρg represents the properties of the fluid, in the context of Titan (g), for this term wefound roughly 10⁻⁸ 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², then L∕S ~ 10⁻³ m⁻¹. 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³, letting the ratio ΔV∕h be comparable to h². On Titan, the topography is rather flat, h ~ 1 km is a large value. Keeping this value, we obtain ΔV∕h ~ h² ~ 10⁶ m², which leads to (21)

Thus, the effect of diffusion is not erased by inter-reservoir seepage if, roughly, K_{p, sc} ≲ 10⁻⁵∕Δt. For a diffusion timescale, for example Δt ~ 1 mega-year (i.e., 3 × 10¹³ s), we have K_{p, sc} ≲ 3 × 10⁻¹⁹ m². 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⁻¹), K_{p, m} the permeability of the reservoir matrix (m²), g the acceleration due to gravity (m s⁻²), the product of the fluid density (kg m⁻³) and the heat capacity (J K⁻¹ kg⁻¹), H the thickness of the reservoir (m), ν the kinematic viscosity of the fluid (Pa s kg⁻¹ m³, i.e., m² s⁻¹), and λ_L the thermal conductivity (W m⁻¹ K⁻¹) 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 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₄(60%)-C₂H₆(20%)-N₂(20%). The outputs of these compulations, made in the frame of the isothermal-isobaric (NpT) 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₂, 60% CH₄, and 20% C₂ H₆, under a thermal gradient of 0.6 K km⁻¹, we found (i.e., (Nield & Bejan 2017) for a permeability of K_{p, m} = 10⁻¹² m², 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 water-ice thermal conductivity has a value around 5 W m⁻¹ K⁻¹, and methane hydrates are much more insulating with values around 0.5 W m⁻¹ K⁻¹, 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 organic-rich sand as low as 0.025 W m⁻¹ K⁻¹. Hence, the discussion in Sect. 3 about Kalousova & Sotin’s results, is also relevant in the case of an organics-dominated matrix.

Fig. 6 Influence of the chemical composition of the fluid on its Rayleigh number. These computations were made with a fixed permeability Kp, m = 10−12 m2, 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 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₂ can inhibit, in some way, the clathration of CH₄ or C₂ H₆. 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_is 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 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₂ in the clathrate cages remains mostly negligible. Except around z ~−3 km, the hydrocarbons CH₄ and C₂ H₆ 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. 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.

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, ³⁶Ar, showed a ratio of ³⁶Ar/N₂ around 10⁻⁷, while Kr and Xe were not detected by the GCMS. Given the instrument detection threshold at 10⁻⁸, the atmosphericmole fractions of Kr and Xe should be below this upper limit. According to the literature, ratio in numbers Ar/N₂, Kr/N₂, and Xe/N₂ should be around 7 × 10⁻², 5.6 × 10⁻⁵, and 3.9 × 10⁻⁶ 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 & 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 noble-gas atoms, which are small in size and heavier than N₂ or C₂ H₆ (Ar: 39.948 g mol⁻¹, Kr: 83.798 g mol⁻¹, Xe: 131.293 g mol⁻¹), could accumulate in the deepest alkanofer layers. We performed preliminary computations by replacing N₂ from the ternary mixture and setting both surface mole fractions of the chosen noble gas “X” and ethane to 10⁻³, which is probably a huge value for a noble gas but a small one for C₂ H₆. Adding the species X to our initial ternary mixture, N₂ +C₂H₆+CH₄, 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₂ 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 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⁻³ 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.

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₂ H₆ 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 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₂ -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₄+C₂H₆+N₂ 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⁶ (Lorenz et al. 2018). On board this probe, a scientific package called DraGMet⁷ 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⁸, which belongs to the OpenAIRE project⁹. 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₄ (80%) + C₂ H₆ (20%) and CH₄ (80%) + N₂ (20%), and the ternary mixture CH₄ (60%) + C₂ H₆(20%) + N₂ (20%), are available online¹⁰ 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 PC-SAFT implementation in FORTRAN 2008 has been made also publicly available on Github¹¹ and it has been archived on Zenodo¹² 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 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, 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⁻¹).

• : partial molar volume of component i (m³ mol⁻¹).

• f_i: fugacity of component i (Pa).

• : fickean molecular diffusion coefficient (m² s⁻¹) of compound i in species j.

• D^M (kg mol⁻¹ m² s⁻¹), D^P (kg mol⁻¹ m² s⁻¹ Pa⁻¹) and D^T (kg mol⁻¹ m² s⁻¹ K⁻¹): 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⁻²).

• k_ij,T: thermal diffusion ratio (dimensionless).

• α_ij,T: thermal diffusion coefficient, k_ij,T = x_ix_jα_ij,T (dimensionless).

• ρ_eff: effective density (kg m⁻³) of the alkanofer material, ρ_eff = Π ρ_liq + (1 − Π) ρ_ice, with Π being the porosity, ρ_liq the density of the liquid (kg m⁻³), and ρ_ice the density of the icy matrix (kg m⁻³).

• μ_i: chemical potential of component i (J mol⁻¹), is the residual chemical potential of component i (J mol⁻¹) 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⁻¹).

• : residual partial molar enthalpy of component i (J mol⁻¹).

• ã^res: residual and reduced (i.e., molar) Helmholtz free energy (J mol⁻¹).

• 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 PC-SAFT equation of state, already used in many Titan-related papers (e.g., Tan et al. 2013; Cordier et al. 2017). Generally speaking, the fugacity may be written as f_i = Φ_ix_iP, where Φ_i is the fugacity coefficient directly provided by PC-SAFT, 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² s⁻¹). 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)

In addition, we define (B.7)

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_ijs are the so-called 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⁻³), used in , the partial molar volume (m³ mol⁻¹), the fugacities f_i, and the thermal diffusion coefficient α_{T i,j} are derived from the PC-SAFT 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 PC-SAFT thanks to the Haase formula (Haase 1969; Pan et al. 2006): (B.19)

where the M_is 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 PC-SAFT.

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⁻³) is denoted ρ as usual, k_c represents the thermal conductivity (W m⁻¹ K⁻¹) 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)

then (C.8)

otherwise, k_v = 0. In Eqs. (C.5) and (C.8), α, C_P, and η are, respectively, the thermal expansion coefficient (K⁻¹) of water-ice I_h, the specific heat capacity (J kg⁻¹ K⁻¹), 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²⁰ 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 Runge-Kutta method, taking the values at the surface as boundary conditions: P = P_s = 1.5 bar, q = q_s = 3.11 × 10⁻³ W m⁻², and T = T_s = 94 K (Sohl et al. 2014).

References

All Tables

All Figures

Fig. 1 Vertical variations of mole fractions (denoted xi) of the main components (i.e., N2, C2 H6 and CH4) of Titan’s alkanofer liquid. This figure depicts a scenario where the icy crust is assumed to be isothermal at the temperature T0 = 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) N2: 0.05, C2 H6: 0.20, CH4: 0.75; (b) N2: 0.20, C2 H6: 0.20, CH4: 0.60; (c) N2: 0.20, C2 H6: 0.30, CH4: 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.

In the text

Fig. 2 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 N2 and C2 H6 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 x1 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 C2 H6-CH4 system, while the dashed line is related to N2-CH4. 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 xi) of the main components (i.e., N2, C2 H6 and CH4) 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) N2: 0.05, C2H6: 0.20, CH4: 0.75; (b) N2: 0.20, C2H6: 0.20, CH4: 0.60; (c) N2: 0.20, C2H6: 0.30, CH4: 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.

In the text

Fig. 5 Phase diagrams of the N2+CH4 system at 110, 115, 120, and 125 K. Circles represent laboratory measurements from Parrish & Hiza (1974), curves are computed according to the PC-SAFT 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 Kp, m = 10−12 m2, 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