Stratification Dynamics of Titan's Lakes via Methane Evaporation

To investigate the thermophysical evolution of a hydrocarbon mixture on Titan, we developed the TITANPOOL numerical model in Matlab using the Reference Fluid Thermodynamic and Transport Properties Database (REFPROP). TITANPOOL is purpose-built to study the thermodynamic, chemical, and physical evolution of a nonideal, ternary mixture of methane, ethane, and nitrogen under the conditions present at the surface of Titan. Since REFPROP cannot directly compute the equilibrium composition and material properties of a mixture, TITANPOOL contains a series of subroutines that call REFPROP to obtain the properties of a specified mixture at a specified temperature and pressure. TITANPOOL uses these functional calls to linearly interpolate the equilibrium liquid composition and its properties. REFPROP natively computes the liquid density of a ternary mixture, along with compositions in vapor–liquid equilibrium (VLE).

REFPROP is produced, updated, maintained, and benchmarked by the National Institute of Standards and Technology (NIST; Lemmon et al. 2010) and has already been used to study the behavior of Titan's lakes (Hartwig et al. 2018). REFPROP is an industry standard program that first computes the Helmholtz free energy of pure and multicomponent material mixtures to compute thermodynamic and material properties (such as density and equilibrium composition). For hydrocarbon and nitrogen mixtures, REFPROP uses the GERG-2008 equation of state (Kunz et al. 2007; Kunz & Wagner 2012), which is commonly used in the oil and gas industry. The GERG-2008 equation of state most accurately models the material properties of hydrocarbon and nitrogen mixtures in its "Normal Range of Validity" (90K ≤ T ≤ 450 K, p ≤ 35 MPa), for which the availability of extensive experimental data has reduced uncertainties in liquid density computations to less than 0.5% and uncertainties in the pTxy relationships (pressure, temperature, and gas and liquid phase compositions) to less than ∼5% (Kunz & Wagner 2012). Outside of this range, GERG-2008 necessarily must extrapolate from the available experimental data.

In this "Extended Range of Validity" (60 K ≤ T ≤ 700 K and p ≤ 70 MPa), little experimental data exist to accurately assess model uncertainties. Kunz & Wagner (2012) estimate that the uncertainties in liquid component density and pTxy relationships are expected to be no more than a few percent. Recent experimental work of methane–ethane–nitrogen mixtures found that GERG-2008 indeed has errors in the densities of such mixtures of no more than ∼7% for conditions present at the surface of Titan (∼1.5 bar, ∼90–94 K; Richardson et al. 2018). Comparing GERG-2008 with our experimental work (∼1.5 bar and 83–87 K), we confirm that errors in pressure in this "Extended Range of Validity" are similarly no more than 7% (see Section 2.3). However, recent experimental work on methane–ethane–nitrogen mixtures found that GERG-2008's errors became quite significant (up to ∼18% in density) at pressures much higher than those at the surface of Titan or the depths of its lakes (Richardson et al. 2018). Nevertheless, since we are only considering shallow lakes at Titan's surface, we restrict our concern to GERG-2008's accuracy at Titan's surface pressure (∼1.5 bar), which is sufficient for our purposes.

We interfaced REFPROP with Matlab to construct our model of the thermodynamic properties of a lake of methane, ethane, and nitrogen on Titan. We assume that the entire lake is at a uniform temperature (thermal equilibrium—a reasonable assumption for ponds that we anticipate to be only a few meters deep) and that there are no sensible heat flows into or out of the lake. Due to the low volatility of ethane relative to methane under Titan surface conditions (Stevenson & Potter 1986; Tokano 2009), we also assume that no ethane evaporates and that the entire mixture is in VLE at 1.5 bar with the gas boundary layer above the surface of the lake. We use REFPROP to compute the equilibrium mole fraction of nitrogen dissolved in the liquid as a function of lake temperature and methane hydrocarbon fraction (f_methane) as a function of temperature at Titan surface conditions:

$$\begin{eqnarray*}&&{f}_{\mathrm{methane}}\equiv \displaystyle \frac{{x}_{\mathrm{methane}}}{{x}_{\mathrm{methane}}+{x}_{\mathrm{ethane}}}.\end{eqnarray*}$$

We then use REFPROP to compute the liquid density of the lake, also as a function of the temperature and methane hydrocarbon fraction.

This model neglects the effects of hydrostatic pressure and temperature gradients/heating on stratification. Nitrogen becomes increasing soluble in methane–ethane mixtures at increasing hydrostatic pressures and, thus, would drive up liquid densities significantly in the depths of deep liquid bodies (e.g., Kraken Mare); this would further stabilize lake stratification. However, these effects are negligible over the depths of the ephemeral ponds and shallow lakes that we consider here. Similarly, solar heating would tend to warm the surface layers of liquid bodies, causing epilimnia to lose dissolved nitrogen and become even more buoyant, allowing stratification to occur at even warmer temperatures. Therefore, our simple model may err conservatively on the conditions required for stratification, which places strong constraints on the dynamic behavior of Titan's lakes.

Nitrogen solubility in liquid alkanes depends strongly on alkane complexity, being highest in the simplest alkane, which is methane (Battino et al. 1984). It also shows strong temperature dependence, with cooler temperatures leading to increased nitrogen solubility (Battino et al. 1984; Malaska et al. 2017; Hartwig et al. 2018). In alkane mixtures, nonideal interactions between alkane species strongly affect nitrogen solubility (Figure 1). Because nitrogen is considerably denser than either methane or ethane, its abundance strongly influences the density of mixtures and can cause methane-rich mixtures to become denser than ethane-rich mixtures under certain conditions (Figure 1).

Zoom In Zoom Out Reset image size

Titan's shallow, ephemeral (Barnes et al. 2013; Soderblom et al. 2016), and polar "phantom" (MacKenzie et al. 2019) lakes are likely cooler than the ∼89–96 K surface (Cottini et al. 2012; Jennings et al. 2016). Methane evaporation allows falling rain to reach the surface ∼4 K cooler (Graves et al. 2008) than the surrounding atmosphere. Recent mesoscale modeling of the interactions between air and lakes shows that methane evaporation may further cool lake temperatures by up to 5 (or more), depending on the initial conditions (Rafkin & Soto 2020). This level of cooling occurs when the atmosphere has less than 20% relative methane humidity and the lake is shallow (∼1 m deep; Rafkin & Soto 2020). The methane humidity across Titan is generally greater than 50%, and even higher near the poles (Lora & Ádámkovics 2016). Nevertheless, a higher relative humidity would correlate with a 37.5% reduction in evaporative cooling, thus lowering the lake temperature by up to 3–4 K. Thus, the combination of these effects allows the temperatures of these shallow lakes to plausibly be as low as 83K.

We find that the density of methane–ethane–nitrogen mixtures at typical Titan surface temperatures (>87 K) decreases monotonically with increasing methane hydrocarbon fraction. However, significant deviations from this monotonic behavior appear at temperatures below 87 K, producing a local density maximum and minimum at methane alkane mole fractions (${f}_{\mathrm{methane}}=\tfrac{{x}_{\mathrm{methane}}}{{x}_{\mathrm{methane}}+{x}_{\mathrm{ethane}}}$) of 70–90 mole% and 40–70 mole%, respectively.  Below 84 K, the methane-rich local maximum composition (f_methane = 70–90 mole%) is denser than a methane-free, ethane–nitrogen binary mixture (f_methane = 0%). This surprising density behavior has important implications for the evaporation-dependent evolution of Titan's lakes, the stratification of Titan's lakes, and the limnological layer from which ethane ice might precipitate.

This composition and temperature dependence allows Titan's lakes to stratify. Above 87 K, evaporative methane loss enriches the lake's surface layer in ethane, producing a denser liquid regardless of composition (Figure 1). Thus, the lake's surface layer will sink into, and mix with the hypolimnion, preventing stratification. Similarly, the surface layer of cool (<87 K) methane-poor lakes (methane alkane fraction less than the composition of the local density minimum, Figure 1, left) will become denser from evaporation, mixing the lake and preventing stratification.

However, methane-rich lakes (relative to the local minimum's composition) below 87 K evolve differently (Figure 2). Preferential methane evaporation from the surface layer decreases the epilimnion's density (if also poorer in methane than the local density maximum), producing a buoyant epilimnion that fails to mix with the hypolimnion (Figure 2(a)). Further evaporation (Figure 2(b)) evolves the epilimnion's composition toward the local density minimum (Figure 2(c)), producing a relatively methane-poor epilimnion and methane-rich hypolimnion, separated by a strong chemocline. Continued evaporation (Figure 2(d)) increases the epilimnion's density above the local density minimum. Eventually, evaporation will drive the epilimnion's density above that of the hypolimnion, causing the lake to overturn (epilimnion sinks and mixes with the hypolimnion, Figure 2(e)). Evaporation from the well-mixed lake would then form a new buoyant epilimnion, and this stratification cycle would repeat.

Zoom In Zoom Out Reset image size

Each overturn event mixes the ethane-rich epilimnion with the methane-rich hypolimnion; the new bulk lake composition is a linear combination of the epilimnion and hypolimnion compositions. The addition of the epilimnion's ethane-rich material to the hypolimnion may result in a bulk lake that is supersaturated in nitrogen (see Figure 3). This supersaturated nitrogen would likely exsolve, forming bubbles; such bubbles have been detected experimentally under similar conditions (Farnsworth et al. 2019; Richardson et al. 2019). On Titan's larger seas, overturn-driven bubble formation may produce the observed "Magic Islands" (Hofgartner et al. 2016). However, it is unclear if these large lakes cool sufficiently to facilitate stratification-driven bubbling.

Zoom In Zoom Out Reset image size

Further evaporative methane loss depletes the bulk lake of methane (composition moves leftward in Figures 1–2). Eventually, the bulk lake composition will evolve past the local density minimum, inhibiting further stratification.

Below 84 K, methane evaporation drives the composition of the epilimnion toward the local density minimum (Figure 2(g)). Although further evaporation (Figure 2(h)) increases the epilimnion density, the epilimnion never sinks and never mixes with the methane-rich hypolimnion because the methane-free, ethane–nitrogen binary mixture (Figure 2(f)) is still less dense than the methane-rich local density maximum (Figure 2(i)). The lake becomes meromictic. These dynamic paths are illustrated in Figure 4.

Zoom In Zoom Out Reset image size

Although the shape of the methane–ethane–nitrogen liquidus curve is poorly understood under Titan conditions, methane-rich mixtures can remain liquid below ∼80 K (Hanley et al. 2017), while ethane-rich mixtures begin to freeze around ∼89 K (Farnsworth et al. 2017, 2020; Richardson et al. 2019). Thus, evaporative methane loss will eventually evolve the epilimnion across the liquidus curve, causing ethane ice to precipitate out of solution and, thus, stabilizing the epilimnion's composition.

We also carried out laboratory studies to investigate the interactions between methane, ethane, and nitrogen, conducting experiments at the Northern Arizona University Astrophysical Materials Laboratory. A detailed description of this facility is available in Tegler et al. (2012) and Grundy et al. (2011). We used a similar setup to Roe & Grundy (2012) and Protopapa et al. (2015). We prepared our samples by first mixing methane and ethane gases at the desired ratio in a 2 liter mixing chamber at room temperature. We then cooled that sample using a helium compressor to cool the cell to the desired temperature. Once the cell was at the desired temperature, we obtained a background spectrum of the empty cell with our Raman spectrometer (Kaiser Optical Systems model Rxn1 using a 785 nm laser). We then opened the line between the cell and the mixing chamber, allowing the sample to condense into the cell. We allowed the sample to reach equilibrium. We then slowly added nitrogen gas until we reached the desired pressure and waited until the temperature and pressure stabilized to ensure the sample reached equilibrium. We measured the composition of layers in the liquid with Raman spectroscopy using a model in which the molar fraction of a species in the liquid is linearly proportional to the band area in the spectrum. Photographs and videos were taken throughout the experiment.

For this specific experiment, we used a sample of 50% methane and 50% ethane at a fixed temperature, stepping from 84 to 87 K in 1 K increments. We added a total of 53.3 kPa (400 Torr) of hydrocarbons: 26.7 kPa (200 Torr) each of methane and ethane. At 83 K, nitrogen was added to reach a pressure at or slightly above 1.47 bar (1100 Torr), until equilibrium was achieved. We then increased the temperature in 1 K increments, consequentially increasing pressure. At each temperature, we agitated the cell to ensure the sample was homogeneous. We then reduced the pressure to 145 kPa (1090 Torr; Titan's surface pressure), which resulted in nitrogen and methane exsolving from the sample's surface. The sample stratified into two separate layers separated by a visible chemocline (Figure 5). We performed Raman spectroscopy at 1 mm intervals in each layer. Measurements in a layer were of nearly uniform composition. As a result, we computed the average and standard deviation of the composition in each layer.

Zoom In Zoom Out Reset image size

As the sample equilibrated, the separation between the layers weakened, demonstrating that nitrogen and methane exsolution was required for layer formation. We then agitated the sample to check for liquid–liquid–vapor equilibrium (e.g., Hanley et al. 2016; Cordier et al. 2017). Upon agitation, the layers mixed, demonstrating that the layering was the result of exsolution-driven stratification rather than the formation of two immiscible liquid phases (as would be present in a liquid–liquid–vapor equilibrium).

We carried out an additional experiment to determine if layering can form via temperature-driven exsolution. In this experiment, we mixed a liquid sample of 64% methane and 36% ethane at 90K. Then, we added nitrogen until we reached a pressure of 1090 torr (i.e., Titan's surface pressure). We reduced the temperature to 82 K and then raised the temperature to 83 K, keeping the pressure at 1090 torr, which resulted in nitrogen and methane exsolving from the sample surface. The liquid again separated into two layers, and we took spectra of each layer.

In Table 1, we tabulate our nitrogen and alkane compositional results. The first two rows are from the temperature-driven exsolution experiments, and the following eight rows are from the pressure-driven exsolution experiments. For each experiment, we provide X_N2, X_CH4/X_C2H6, and X_CH4/(X_CH4+X_C2H6) for the epilimnion (the upper layer) and hypolimnion (the lower layer). Notice that in all five experiments X_N2, X_CH4/X_C2H6, and X_CH4/(X_CH4+X_C2H6) are smaller in the epilimnion (the upper layer) than the hypolimnion (the lower layer). In other words, the upper layer is richer in ethane and the lower layer is richer in methane and nitrogen as predicted by our numerical model.

Table 1.  Data from Experimental Effort

T	P	XCH4/(XCH4+XC2H6)	Layer	XN2	XCH4/(XCH4+XC2H6)	${P}_{\mathrm{REFPROP}}$ (pred.)	Model Error
(K)	(Torr)	(Initial)		(Layers)	(Layers)	(Torr)	(%)
83	1090	0.64	U	0.455 ± 0.021	0.606 ± 0.013	1072	1.65
83	1090	0.64	L	0.560 ± 0.006	0.635 ± 0.007	...	...
84	1093	0.50	U	0.197 ± 0.006	0.487 ± 0.022	1110	1.02
84	1093	0.50	L	0.251 ± 0.003	0.503 ± 0.004	...	...
85	1090	0.50	U	0.161 ± 0.005	0.477 ± 0.020	1123	3.03
85	1090	0.50	L	0.191 ± 0.004	0.501 ± 0.003	...	...
86	1090	0.50	U	0.138 ± 0.004	0.477 ± 0.025	1123	3.03
86	1090	0.50	L	0.162 ± 0.005	0.498 ± 0.003	...	...
87	1090	0.50	U	0.124 ± 0.010	0.477 ± 0.026	1140	6.48
87	1090	0.50	L	0.135 ± 0.003	0.498 ± 0.002	...	...

Note. This table reports our experimental data for five experimental runs. The first two rows (83 K) present our temperature-driven exsolution experiment results, while the next eight rows (84–87 K) present our pressure-driven exsolution experiment results. The lower layer (hypolimnion) is enriched in methane for all five experiments. Uncertainties are the standard deviations of measurements in a layer. An ethane-enriched upper layer and a methane–nitrogen-enriched lower layer confirms our numerical model in this temperature range of interest. We also compute the expected pressure from the temperature and composition of the sample's upper layer to validate our TITANPOOL code and check that errors are small (of the order of a few percent).

Download table as:  ASCIITypeset image

Lastly, we used these experimental results to check the accuracy of the GERG-2008 equation of state in its "Extended Range of Validity" (see Section 2.1). We focused on the composition of the epilimnia (the upper layer), since the upper layer is equilibrated with the vapor phase (indeed, this is why the upper layer grows from the surface boundary). We input the temperature and composition of the upper layer into REFPROP (which implements GERG-2008) and output the expected pressure. We then compare this expected pressure to the measured pressure and find that GERG-2008 produces errors no larger than 6.48% for our experimental results (see Table 1). This is consistent with GERG-2008's expected error (Kunz & Wagner 2012) and previous experimental results (Richardson et al. 2018), which confirm that GERG-2008 and REFPROP are sufficiently accurate for our purposes.