© The Authors 2025

1 Introduction

A notable feature of Titan’s atmosphere is the absence of heavy noble gases except argon, as revealed by the Gas Chromatograph Mass Spectrometer (GCMS) aboard the Huygens probe during its descent to the surface in January 2005. The GCMS detected argon in two forms: primordial ³⁶Ar and radiogenic ⁴⁰Ar, a decay product of ⁴⁰K. However, other primordial noble gases, including ³⁸Ar, krypton (Kr), and xenon (Xe), were not detected, with GCMS placing upper limits on their atmospheric mole fractions of <10⁻⁸ (Niemann et al. 2005, 2010).

Titan is believed to have formed from a mixture of ices, silicates, and organic materials (Tobie et al. 2014; Reynard & Sotin 2023). A plausible origin for Titan’s volatile phase is that its building blocks had a comet-like volatile composition (Alibert & Mousis 2007; Mousis et al. 2009a). This is supported by the D/H ratio in Enceladus’s water vapor matching the value observed in comets (Waite et al. 2009), which suggests similar building blocks for Titan and Enceladus (Mousis et al. 2009a,b,c). Titan’s noble gas depletion could be explained by the trapping of Ar, Kr, and Xe in clathrate hydrates (hereafter clathrates) within its icy shell (Osegovic & Max 2005; Thomas et al. 2008; Mousis et al. 2011), with sequestration occurring at the surface and facilitated by cryovolcanism (Mousis et al. 2011).

However, early in Titan’s history, high surface temperatures likely supported a global ocean in contact with its primordial atmosphere (Lunine & Stevenson 1987; Lunine & Nolan 1992), enabling clathrate formation at the ocean–atmosphere interface. This study explores how Titan’s hydrosphere acted as a clathrate reservoir. The compositions of Titan’s ocean and atmosphere depended on the volatile inventory of its building blocks, with the CO₂/NH₃ mass ratio influencing atmospheric composition and thickness (Marounina et al. 2018; Amsler Moulanier et al. 2025). Assuming a comet-like volatile inventory, the early atmosphere likely contained CH₄ and CO₂ (Lunine et al. 2009; Tobie et al. 2014), with NH₃ as a key hydrospheric volatile (Lunine & Stevenson 1987; Grasset & Sotin 1996; Tobie et al. 2012). NH₃ likely contributed to Titan’s current N₂ atmosphere (Atreya et al. 1978; Jones & Lewis 1987; McKay et al. 1988; Lunine et al. 1989; Glein et al. 2009) and may have acted as an antifreeze (Lunine & Stevenson 1987; Grasset & Sotin 1996).

In this study, we first analyzed the evolution of Titan’s primordial hydrosphere composition, assuming its volatile inventory was primarily delivered by icy building blocks. Second, we applied a statistical thermodynamic model to investigate whether clathrate formation at the surface of Titan’s ocean explains the depletion of primordial noble gases and its impact as well as the impact of clathrate formation on atmospheric composition.

2 Methods

We aimed to calculate the minimum thickness of the mixed clathrate crust formed during Titan’s primordial phase to explain the current noble gas depletion. The model used in this study simulates the evolution of the composition of Titan’s primordial hydrosphere as temperature decreases, leading to the formation of an ice crust. If the stability conditions for clathrate formation are met, the model calculates how this process could affect the distribution of volatiles in Titan’s primordial atmosphere.

We assumed that shortly after accretion, Titan’s surface temperature was high enough to support a primordial ocean in direct contact with the atmosphere (Lunine & Stevenson 1987; Lunine & Nolan 1992). These two thermodynamic phases were assumed to be in equilibrium, and the partitioning of volatiles between them was determined accordingly. The model considers both the liquid-vapor equilibrium and the chemical equilibrium between CO₂ and NH₃ within the ocean. Starting from a given bulk composition, surface temperature, and an initial guess for the atmospheric composition, we calculated the equilibrium surface pressure of the global ocean and the distribution of volatiles between the primordial atmosphere and the ocean. Once equilibrium is reached, a mass balance check ensures consistency with the initial volatile mass in the bulk composition. If convergence is not achieved, the process is repeated, adjusting the atmospheric partial pressures until the volatile partitioning between the atmosphere and ocean satisfies bulk mass conservation. For a detailed formulation of the equations governing this equilibrium, we refer to Appendix A and Amsler Moulanier et al. (2025). The model does not account for water–rock interactions between the rocky mantle and the ocean. If the total atmospheric pressure at the surface exceeds the dissociation pressure of the clathrates (Eq. (1)), mixed clathrates can form at the ocean surface. The dissociation pressure of the mixed clathrate is given as $${\left(P_{diss}^{mix}\right)}^{-1}=\sum _{i=1}^N\frac{y_i}{P_{diss}^i},$$(1) where $P_{diss}^i$ is the dissociation pressure of the pure clathrate of gas i, y_i is the atmospheric mole fraction of gas i, and $P_{diss}^{mix}$ is the dissociation pressure of the mixed clathrate (Lipenkov & Istomin 2001). When the stability conditions for clathrates are met, their composition is determined using a statistical thermodynamic formation model derived from Mousis et al. (2013) and detailed in Appendix B.

For a given pressure, temperature and gas composition, our model calculates the fractions of volatiles trapped in mixed clathrates. Since in our model CH₄ and CO₂ are the most abundant gases in the atmosphere, they are the primary clathrate formers, and their predominant clathrate structure is expected to be structure I (Sloan & Koh 2008). Consequently, all calculations were performed for structure I clathrates (SI). At each time step, the thickness of the mixed clathrate layer is determined based on a convergence criterion. While clathrate formation is generally considered to be a rapid process (Englezos et al. 1987; Uchida et al. 2004), its kinetics remain poorly constrained. Consequently, interpretation of the time step in terms of years is challenging, and the results are instead expressed as a function of clathrate crust thickness. We stress that, as observed on Earth, the kinetics of clathrate formation occurs on timescales between 1000 and 100 000 years. However, while this timescale is relevant in Earth’s context, it still needs to be evaluated for Titan. Then, as the crust thickens, the diffusion of the gas becomes less efficient, slowing growth (Ruppel & Waite 2020).

An SI clathrate consists of eight cages per unit cell, with a volume of 1.73 × 10⁻²⁷ m³ (Sloan & Koh 2008). In our case, mixed clathrates form, with each cage accommodating one gas molecule. The amount of mixed clathrates required to trap noble gases is expressed as the corresponding trapped gas volume per kilometer of crust. This volume is then subtracted from the primordial atmosphere, which serves as the gas reservoir for clathrate formation. Finally, the density of the formed clathrates is computed using the method outlined in Appendix B.3. We emphasize that the feedback effects on ocean composition due to changes in atmospheric composition were not calculated in this study.

Fig. 1 Equilibrium partial pressures of CO2 (dark blue) and NH3 (green) in Titan’s primordial atmosphere as a function of the total NH3 dissolved in the ocean. Each set of green and blue curves, distinguished by different line styles, represents a specified initial CO2 atmospheric pressure prior to chemical equilibrium between the ocean and atmosphere. The initial theoretical PCO2 values of 1, 5, 10, 20, and 30 bar correspond to CO2 bulk masses of 7.8 × 1019 kg, 3.79 × 1020 kg, 7.38 × 1020 kg, 1.37 × 1020 kg, and 1.92 × 1030 kg, and a dissolved CO2 mole fraction of 1.36 × 10−3, 6.78 × 10−3, 1.36 × 10−2, 2.71 × 10−2, and 4.07 × 10−2, respectively.

3 Results

We investigated the composition of Titan’s primordial hydrosphere during two distinct phases: the open-ocean phase, when the ocean directly interacted with the atmosphere, and the closed-ocean phase, when an icy crust formed, isolating the ocean from atmospheric exchange.

3.1 Combined evolution of Titan’s ocean and atmosphere compositions during the open-ocean phase

The composition and evolution of Titan’s primordial hydrosphere were examined based on a specific bulk composition, with particular emphasis on the role of NH₃, whose abundance ranges from 1.5 to 15 wt% in the ocean (Tobie et al. 2012). The dissolution of ammonia and carbon dioxide in the global ocean can significantly influence atmospheric composition and thickness through the chemical balance between NH₃ and CO₂ speciation products. In particular, the ratio of dissolved CO₂ to NH₃ is a key factor in determining whether a CO₂-rich atmosphere can be sustained (Marounina et al. 2018; Amsler Moulanier et al. 2025).

Figure 1 presents the equilibrium partial pressures of CO₂ and NH₃ in Titan’s primordial atmosphere at 273.15 K, as a function of the total NH₃ fraction in the ocean and the theoretical initial CO₂ partial pressure. The results indicate that when the NH₃ concentration exceeds 2%, even an initially high atmospheric CO₂ abundance (>10 bar) cannot be sustained above Titan’s primordial ocean. This finding challenges the hypothesis that Titan once hosted a thick CH₄- and CO₂-rich primordial atmosphere (Lunine et al. 2009; Tobie et al. 2014) in direct contact with an NH₃-rich ocean (Lunine & Stevenson 1987; Tobie et al. 2012). It also raises questions about the feasibility of surface clathrate formation during the open-ocean phase, since this process requires the total atmospheric pressure to exceed the clathrate dissociation pressure. As shown by Marounina et al. (2018), we suggest that, in the temperature range considered, a thick CO₂-rich primordial atmosphere in equilibrium with a NH₃-rich ocean would not have been sustainable.

The composition of Titan’s primordial atmosphere, and its relative tendency to favor CO₂ or NH₃, is inferred from its primordial bulk composition. Here, Titan’s volatile inventory is assumed to have originated from volatile-rich icy planetesimals with a comet-like composition, with the volatile fraction in the ice phase analogous to that measured in comet 67P/Churyumov-Gerasimenko (67P/C-G) (Rubin et al. 2019). Given Titan’s bulk density, the total mass of ice in its composition is assumed to be ~50% (Lunine et al. 2009; Tobie et al. 2012). In this study, we assumed that about 30% of Titan’s mass contributed to the melting of ice that formed the primordial atmosphere and ocean (Lunine et al. 2009), corresponding to a melted ice fraction of 60%. By computing the normalized volatile mass distribution in ice, we could estimate the total bulk mass of each volatile delivered to Titan’s hydrosphere (see Table 1). Based on this volatile inventory, Titan’s primordial ocean would be approximately 500 km deep. If the melted ice fraction were to reach 80% (Monteux et al. 2014), the partial pressures above the ocean would increase. Specifically, higher CO₂ pressures could support the existence of liquid CO₂. However, a higher surface pressure could promote clathrate formation at a temperature slightly higher than 282 K if 80% of the ice melts.

Using our model (Appendix A), we computed the thermodynamic equilibrium between the primordial atmosphere ocean and how volatiles are partitioned between the two phases. Assuming comet-like building blocks, a significant amount of CO₂ is released from the ice (4.12 × 10²¹ kg). Our computations show that such a mass of CO₂ introduced into the hydrosphere would lead to a considerable pressure at Titan’s surface. To model a physically accurate system, we capped the initial partial pressure of CO₂ at its saturation pressure as computed with Antoine’s law. As mentioned previously, both CO₂ and NH₃ should be present within the primordial hydrosphere of Titan in significant proportion. Although NH₃ has a higher solubility than CO₂, Table 1 shows that at 273.15 K a thick CO₂ atmosphere forms at the ocean surface, while most of NH₃ is speciated and retained in the ocean as ions (${\mathrm{N}\mathrm{H}}_4^{+}$ and NH₂COO⁻) leading to a limited fraction of free NH₃ present within the ocean. This conclusion particularly contrasts with the assumption of a primordial NH₃-rich global ocean on Titan (Tobie et al. 2012). The formation of a CO₂-rich atmosphere is, however, consistent with the findings of Ishimaru et al. (2011), which supports the idea that a CO₂-rich oxidized primordial atmosphere could have allowed the formation of N₂ from NH₃ via impact shock chemistry.

Figure 2 illustrates the effect of temperature on the composition of Titan’s primordial atmosphere and ocean at equilibrium. As the temperature decreases from 300 K to 273.15 K, the atmosphere becomes less dense and the fraction of dissolved volatiles increases. The magnitude of this change varies with the gas under consideration. For example, in this temperature range, the partial pressure of Xe decreases by 24%, while the abundance of Ar in the atmosphere is only slightly affected (8%). In contrast, the evolution of the CO₂ and NH₃ abundances in the atmosphere and ocean follows a distinct trajectory, determined by their chemical equilibrium in the ocean.

Fig. 2 Evolution of the composition of Titan’s primordial atmosphere and ocean at equilibrium, as a function of surface temperature.

3.2 Effect of clathrate formation on the composition of the primordial atmosphere during Titan’s closed-ocean phase

We focused on two stages of Titan’s evolution, corresponding to a surface temperature of 280 K and 273.15 K. At those temperatures, the equilibrium total surface pressure falls within the mixed clathrate stability region, with boundaries calculated using Eq. (1). This boundary is reached when the surface temperature drops to 282 K, allowing clathrates to form on top of the ocean, still uncovered by the ice crust. In our model, we assumed that the clathration process halts when the surface pressure falls below the dissociation pressure of mixed clathrates. At 280 K and 273.15 K, we calculated the clathrate composition and assessed its impact on the distribution of volatiles in the atmosphere. Figure 3 illustrates the distribution of partial pressures as a function of mixed clathrate crust thickness. The detection limit of the Huygens GCMS instrument for volatile abundances (indicated in dashed black on the figure), set at 10⁻⁸ (corresponding to a maximum partial pressure of approximately 1.5 × 10⁻³ Pa; Niemann et al. 2005), serves as the threshold for determining when a gas can be considered depleted from the atmosphere. At 280 K, close to the mixed clathrate stability limit, a 1.23 km thick crust can form until the surface pressure becomes too low to allow clathrate stability. At such a low thickness, the noble gases are not trapped in significant proportions. Titan’s surface must cool to reduce the dissociation pressure of the clathrates, thereby trapping more gases. At 273 K, Xe is almost completely removed from the atmosphere, becoming sufficiently rare to be undetectable by the Huygens GCMS – when the clathrate crust reaches a thickness of 3.3 km. In contrast, trapping all of the Kr in the atmosphere requires a thicker crust, as this noble gas is less efficiently trapped than Xe. Hence, Kr is almost completely trapped in clathrates with a crust thickness of 10.31 km.

However, Ar is much less efficiently trapped in clathrates compared to other noble gases and cannot be sufficiently captured before reaching the maximum crust thickness of 11 km at 273 K. For Ar to be completely removed from the primordial atmosphere, Titan’s surface would need to cool further. Based on the total surface pressure before crust formation and the Ar abundance shown in Table 1, our calculations suggest that complete trapping of Ar in clathrates – whether mixed or pure – likely did not occur during Titan’s early history.

Additionally, a shift in the composition of the mixed clathrates occurs as the crust thickens and time progresses. This transition is highlighted in Fig. 3 by the dashed red line, marking a change in the dominant gas forming the clathrate. To the left of this line, CH₄ is the primary gas captured in the clathrate, despite being less abundant than CO₂ in the atmosphere, due to thermodynamic conditions favoring its trapping. To the right of the red line, however, CH₄ becomes too scarce in the atmosphere, and CO₂ takes over as the driving force for clathrate formation. This shift in the dominant clathrate-forming gas significantly alters the overall composition of the mixed clathrates, as shown in Fig. 3, and consequently affects the composition of the crust over time. As shown in Fig. 4, this change in composition leads to a variation in clathrate density. As the crust thickens and clathrates shift from a CH₄-rich to a CO₂-rich composition, the density of the clathrates increases, surpassing a value of 1. As a result, while CH₄-rich clathrates would float on the ocean, the more recent CO₂-rich upper layers become denser than the older, lower layers, exceeding the density of liquid water when ρ > 1. This density gradient could eventually lead to structural remodeling of the clathrate crust over time.

Fig. 3 Evolution of atmospheric abundances as a function of clathrate crust thickness at 280 and 273.15 K. The dashed red line indicates the threshold at which the clathrate composition changes from CH4-rich to CO2-rich. The dashed black line represents the detection limit of Huygens’ GCMS.

Fig. 4 Density of clathrates as a function of crust thickness, assuming formation at 273 K. The step for both temperatures highlights the transition from CH4-rich to denser CO2-rich clathrates.

3.3 Formation of a crustal CH₄-rich reservoir

Titan’s atmosphere features a ~5% CH₄ fraction measured at its surface by Huygens (Niemann et al. 2005). Since CH₄ is destroyed by UV photochemistry on 10–100 Myr timescales (Yung et al. 1984; Wilson & Atreya 2004), this suggests replenishment from an internal reservoir containing at least ~2.8 × 10¹⁷ kg of CH₄ to maintain this level over 4.5 billion years (Niemann et al. 2005).

Assuming a comet-like primordial bulk composition, we estimated that Titan’s primordial atmosphere at 300 K contained at least 9.4 × 10¹⁹ kg of CH₄, providing a sufficient initial reservoir. Once clathrates form, our calculations show that at 273.15 K, with a 11 km-thick crust, about 7.76 × 10¹⁹ kg of CH₄ would be retained. Thus, a clathrate crust just below a dozen kilometers thick could store enough CH₄ to serve as a long-term reservoir to maintain Titan’s atmospheric composition. Based on the current mass of CH₄ in Titan’s atmosphere, we estimate that the reservoir could replenish the atmosphere for at least 2.8 Gyr.

4 Discussion and conclusions

This study explores clathrate formation during Titan’s primordial phase as a key factor in its present-day noble gas depletion. If Titan’s volatiles originated from comet-like building blocks, a CO₂- and CH₄-rich atmosphere would have formed over a liquid water ocean depleted of free NH₃. In this ocean, the CO₂–NH₃ equilibrium would have converted most NH₃ into ${\mathrm{N}\mathrm{H}}_4^{+}$. Notably, if the NH₃ concentration in the ocean did not exceeded a few percent (which is compatible with existing constrains on comets and Enceladus data), the primordial atmosphere should be dominated by CO₂. As temperatures fell to 282 K and below, clathrates could have formed a crust on the surface of Titan’s ocean, incorporating Xe and Kr into a CH₄- and CO₂-rich structure and thereby accounting for their depletion in the atmosphere. In particular, if this clathrate crust reached a thickness of about 10 km at 273.15 K during Titan’s early evolution, it could have efficiently sequestered Xe and Kr from the atmosphere. However, this process could not fully remove Ar, differing from Mousis et al. (2011) due to variations in atmospheric composition. We propose that Titan accreted from planetesimals that partially devolatilized in Saturn’s sub-nebula (Mousis et al. 2009b), which would explain the absence of primordial Ar. This clathrate formation may have also retained CH₄, contributing to its long-term atmospheric replenishment.

We emphasize that the modeling carried out in this study does not account for the clathrate formation kinetics, which is a critical aspect in this process. Specifically, our results highlighting the capture of Kr and Xe during Titan’s primordial phase at 273.15 K are valid only if the clathrate crust formation occurred while the surface temperature remained at 273.15 K. Indeed, a comparison between the cooling rate of Titan and the growth rate of the clathrate crust would be necessary but remains challenging given the current state of the art. Moreover, as the crust thickens, the transport of fresh water to the surface and/or the transport of gas through the clathrate crust should be limited, affecting the capture of noble gases from the atmosphere.

The case of argon can be considered a critical aspect in our study. Argon mainly exists as ⁴⁰Ar in Titan’s atmosphere (Niemann et al. 2010). This isotope likely originated in the rocky mantle as ⁴⁰K, which subsequently decayed into ⁴⁰Ar and was transported into the ocean (Tobie et al. 2012). Its presence in the atmosphere suggests outgassing from the icy crust, where it may have interacted with clathrates (Tobie et al. 2014). However, given the low trapping efficiency of Ar in clathrate hydrates (Thomas et al. 2007, 2008), we propose that, in competition with other dissolved species such as CH₄ and CO₂, ⁴⁰Ar most likely escaped from the crust without being incorporated into clathrates.

Finally, the bulk composition scenario and the results presented in this study raise the following questions:

• If a thick, primordial CO₂-rich atmosphere formed from comet-like building blocks shortly after Titan’s accretion, what processes could account for the low present-day abundance of CO₂? Beyond clathrate sequestration, several plausible mechanisms include surface condensation as Titan cooled (Lunine et al. 2009), early atmospheric escape (Kuramoto & Matsui 1994), photolysis, and loss through impact-driven chemistry (Ishimaru et al. 2011);

• If N₂ formation from NH₃ via shock chemistry is possible in a CO₂-rich atmosphere (Ishimaru et al. 2011), what mechanism could sustain NH₃ in the atmosphere, considering its expected retention in the ocean under such conditions? Possible processes include continuous outgassing from the interior, photochemical recycling, or impact-driven volatilization (Sekine et al. 2011);

• If primordial clathrates containing CH₄, Ar, Kr, and Xe act as a reservoir for CH₄, how can CH₄ be released without simultaneously reintroducing noble gases – specifically, 3.74 × 10¹⁶ kg of Xe and 5.65 × 10¹⁶ kg of Kr-into the atmosphere? One possibility is that Xe and Kr remain sequestered in residual surface clathrates during CH₄ release, or are reenclathrated shortly thereafter. This implies the existence of an additional mechanism that prevents their effective outgassing (Thomas et al. 2008; Mousis et al. 2011). The existence of such a mechanism is supported by the strong thermodynamic preference of Xe and Kr for clathrate formation compared to other gases (Sloan & Koh 2008).

Acknowledgements

O.M. and A.B. acknowledge support from CNES. The project leading to this publication has received funding from the Excellence Initiative of Aix-Marseille Université – A*Midex, a French “Investissements d’Avenir programme” AMX-21-IET-018. This research holds as part of the project FACOM (ANR-22-CE49-0005-01 ACT) and has benefited from a funding provided by l’Agence Nationale de la Recherche (ANR) under the Generic Call for Proposals 2022.

Appendix A Description of the atmosphere-ocean equilibrium

The model presented in this section follows the approach outlined in Marounina et al. (2018), which models the H₂O–CO₂–NH₃ system based on previous works (Göppert & Maurer 1988; Pazuki et al. 2006; Darde et al. 2010). Our ocean equilibrium model, developed in Amsler Moulanier et al. (2025) and in this study, shares several similarities with Marounina et al. (2018), particularly the use of the Peng-Robinson-Gasem equation of state (EOS) for fugacity coefficients of gases and the extended UNIQUAC model for activity coefficients of dissolved species. However, several key differences distinguish our model:

• The standard-state fugacity of dissolved species is approximated using Henry’s law constant, with an added Poynting correction factor.

• The van der Waals one-fluid mixing rule is incorporated in the Peng-Robinson EOS to account for gas interactions.

• Constants such as Henry’s constants, K_i dissociation constants, and T_c, P_c, and ω_c values are sourced from different references.

• The activity coefficient for CH₄ is computed.

• Clathrate formation at the surface of the primordial ocean is included. A key novelty of our model is simulating the impact of clathrate formation on the evolution of the primordial atmosphere. Additionally, the composition of clathrates is computed through fractional occupancies of guest molecules, using a statistical thermodynamic model.

The components of our model are outlined below.

A.1 Liquid-vapor equilibrium

We computed the liquid-vapor equilibrium for each species starting from the equality between the vapor and liquid phases’ fugacities, $f_i^V$ and $f_i^L$, respectively: $$f_i^V(T,P,y_i)=f_i^L(T,P,x_i).$$(A.1)

Here, the fugacity of the species i in the vapor phase, $f_i^V$, is a function of temperature T, pressure P, and mole fraction y_i in the vapor phase. Likewise, the fugacity of the species i in the liquid phase, $f_i^L$, is a function of T, P, and mole fraction x_i in the liquid phase.

To account for the nonideal behavior in the gas phase and the nonideal interactions between the dissolved species in the liquid phase, we introduced respectively the fugacity coefficient (ϕ) and the activity coefficient (γ). Hence, the fugacity of each phase for an individual species i can be written as $$\begin{array}{rl}& f_i^V={\varphi }_i{y}_{i}P,\\ & f_i^L={\gamma }_i{x}_i{f}_i^0.\end{array}$$(A.2)

In the vapor phase, the fugacity $f_i^V$ of the specie i is a function of the fugacity coefficient ϕ_i, the mole fraction y_i, and the total pressure P in Pascal. In the liquid phase, the fugacity $f_i^L$ is a function of the activity coefficient γ_i, the mole fraction x_i, and the standard-state fugacity $f_i^0$.

Then, we were able to compute the partitioning of volatiles between the atmosphere and the ocean, taking the nonideal behavior of both phases into account, by substituting these expressions into Eq. (A.1). The standard-state fugacity $(f_i^0)$ of H₂O is set to the saturation pressure of water, $P_{H_{2}O}^s(T)$, which is calculated using Antoine’s equation (Stull 1947): $${\log }_{10}\left(P_{H_{2}O}^s\right)=A-\frac{B}{T+C},$$(A.3) with A = 4.6543, B = 1435.264, C = −64.848, T in K and P in bar.

Equation (A.2) can then be rewritten as $$\begin{array}{rl}& {\varphi }_{H_{2}O}P{y}_{H_{2}O}={\gamma }_{H_{2}O}{x}_{H_{2}O}{P}_{H_{2}O}^s(T)\\ & {\varphi }_{i}P{y}_i={\gamma }_i{x}_i{H}_{H_{2}O,i}(T)\exp \left(\frac{v_{i,H_{2}O}^{\mathrm{\infty }}(P-P_{H_{2}O}^s)}{RT}\right),\end{array}$$(A.4) where the standard-state fugacity $f_i^0$ of dissolved species i is approximated as the Henry’s law constant, corrected by a Poynting factor (Kawazuishi & Prausnitz 1987; Darde et al. 2010). The Henry’s law constants and their temperature dependences are reported in Table A.1 and were taken from Kawazuishi & Prausnitz (1987) for CO₂ and NH₃, and from Warneck & Williams (2012) for the other species.

The molar volume of the species i at infinite dilution, $v_{i,H_{2}O}^{\mathrm{\infty }}$, is extrapolated from Rumpf & Maurer (1993a,b) and Garcia (2001) for CO₂ and NH₃. For the other species, $v_{i,H_{2}O}^{\mathrm{\infty }}$ is set to 0. R is the ideal gas constant.

We computed the fugacity coefficient (ϕ_i) for each species using the Peng-Robinson EOS (Peng & Robinson 1976) and the van der Waals one-fluid mixing rule. The molar volume of the mixture v was derived from the Peng-Robinson EOS: $$P=\frac{RT}{v-b}-\frac{a\alpha }{v^2+2vb-b^2},$$(A.5) where $$\begin{array}{rl}& a=\frac{0.457235R^2{T}_c^2}{P_c},\\ & b=\frac{0.077796R{T}_c}{P_c},\\ & \alpha ={(1+\kappa (1-\sqrt{\frac{T}{T_c}}))}^2,\end{array}$$(A.6) and κ = 0.37464 + 1.54226ω − 0.2699ω².

Here, T_c is the critical temperature, P_c is the critical pressure, and ω is the acentric factor of each volatile. Their values are shown in Table A.2.

The mixture parameters A_mix and B_mix required for the Peng-Robinson EOS and the van der Waals one-fluid mixing rule were computed using the following equations: $$\begin{array}{rl}{A}_i=& a_i\times {\alpha }_i,\\ A_{ij}=& (1-k_{ij})\sqrt{A_i{A}_j},\\ A_{\text{mix\hspace{0.17em}}}=& \sum _i^N\sum _j^N{y}_i{y}_j{A}_{ij},\\ \text{\hspace{0.17em}and\hspace{0.17em}}{B}_{mix}=& \sum _i^N{y}_i{b}_i,\end{array}$$(A.7)

where N is the total number of species. α_i, a_i and b_i are computed for each species i using Eq. A.6. In order to better capture the nonideal behavior of the volatile mixtures in the atmosphere and ocean, binary interaction coefficients k_ij are used to further account for the nonideal interactions between the volatile species, following the approach described in Gasem et al. (2001). If they are not available in the literature, they are set to 0. These binary interaction coefficients are determined empirically and taken from the literature (Dhima et al. (1999) for H₂O-CO₂, and Vrabec et al. (2009) for the other couples). After calculating the molar volume v from the Peng-Robinson EOS (Eq. A.5), the fugacity coefficient ϕ_i of a species i can be derived from Eq. A.5: $$\begin{array}{rl}\ln ({\varphi }_i(T,v,y))& =\frac{b_i}{B_{mix}}(Z-1)-\ln \left(\frac{P}{RT}(v-B_{mix})\right)+\\ & \frac{A_{mix}}{2\sqrt{2}RT{B}_{mix}}\ln \left(\frac{v+B_{mix}2.414}{v-B_{mix}0.414}\right)\end{array}$$(A.8) with, $Z=\frac{Pv}{RT}$.

The activity coefficients of H₂O, CO₂ and NH₃ are computed with the extended UNIQUAC model (Thomsen & Rasmussen 1999). The activity coefficient of CH₄ is taken from Kvamme (2021). For the other neutral species, we set γ_i = 1.

A.2 Chemical equilibrium

The liquid-vapor equilibrium is coupled with the chemical equilibrium taking place within the aqueous H₂O-CO₂-NH₃ system: $$H_{2}O\rightleftharpoons H^{+}+{OH}^{-},$$(A.9) $${CO}_2+H_{2}O\rightleftharpoons {HCO}_3^{-}+H^{+},$$(A.10) $${NH}_3+H_{2}O\rightleftharpoons {NH}_4^{+}+{OH}^{-},$$(A.11) $${HCO}_3^{-}\rightleftharpoons {CO}_3^{2-}+H^{+},$$(A.12) $${HCO}_3^{-}+{NH}_3\rightleftharpoons {NH}_2{COO}^{-}+H_{2}O.$$(A.13)

The equations describing the chemical equilibrium within the H₂O-CO₂-NH₃ system are derived from the dissociation reactions of these species: $$K_{H_{2}O}=\frac{m_{{OH}^{-}}{m}_{H^{+}}}{x_{H_{2}O}}\frac{{\gamma }_{{OH}^{-}}{\gamma }_{H^{+}}}{{\gamma }_{H_{2}O}},$$(A.14) $$K_{{CO}_2}=\frac{m_{{HCO}_3^{-}}{m}_{H^{+}}}{x_{H_{2}O}{m}_{{CO}_2}}\frac{{\gamma }_{{HCO}_3^{-}}{\gamma }_{H^{+}}}{{\gamma }_{H_{2}O}{\gamma }_{{CO}_2}},$$(A.15) $$K_{{NH}_3}=\frac{m_{{NH}_4^{+}}{m}_{{OH}^{-}}}{x_{H_{2}O}{m}_{{NH}_3}}\frac{{\gamma }_{{NH}_4^{+}}{\gamma }_{{OH}^{-}}}{{\gamma }_{H_{2}O}{\gamma }_{{NH}_3}},$$(A.16) $$K_{{HCO}_3^{-}}=\frac{m_{{CO}_3^{2-}}{m}_{H^{+}}}{m_{{HCO}_3^{-}}}\frac{{\gamma }_{{CO}_3^{2-}}{\gamma }_{H^{+}}}{{\gamma }_{{HCO}_3^{-}}},$$(A.17) $$K_{{NH}_2{COO}^{-}}=\frac{m_{{NH}_2{COO}^{-}}{x}_{H_{2}O}}{m_{{NH}_3}{m}_{{HCO}_3^{-}}}\frac{{\gamma }_{{NH}_2{COO}^{-}}{\gamma }_{H_{2}O}}{{\gamma }_{{NH}_3}{\gamma }_{{HCO}_3^{-}}},$$(A.18) where K_i(T) is the dissociation constant of each reaction taken from Kawazuishi & Prausnitz (1987), m_i (in mol/kg) is the molality of the aqueous species, x_H2O the mole fraction of water, and γ_i the activity coefficient of the species i.

To compute the activity coefficients involved in the chemical equilibrium equations, we used the extended UNIQUAC model, first introduced in Thomsen & Rasmussen (1999) and the parameters from Darde et al. (2010). This model calculates the activity coefficient using combinatorial, residual, and electrostatic terms, with the electrostatic terms based on the extended Debye-Hückel law. It is well established for studying this specific chemical equilibrium (Thomsen & Rasmussen 1999; Thomsen 2005; Darde et al. 2010). Finally, we could compute the whole set of variables using the mass and charge balance equations derived as follows: $$m_{{CO}_{2tot}}=m_{{CO}_2}+m_{{HCO}_3^{-}}+m_{{CO}_3^{2-}}+m_{{NH}_2{COO}^{-}},$$(A.19) $$m_{{NH}_{3tot}}=m_{{NH}_3}+m_{{NH}_4^{+}}+m_{{NH}_2{COO}^{-}},$$(A.20) $$m_{H^{+}}+m_{{NH}_4^{+}}=m_{{OH}^{-}}+m_{{HCO}_3^{-}}+2m_{{CO}_3^{2-}}+m_{{NH}_2{COO}^{-}},$$(A.21) where m_CO2tot and m_NH3tot are the total molalities of CO₂ and NH₃ incorporated into the liquid phase.

Appendix B Clathrate formation model

B.1 The statistical thermodynamic model

The method described in this section, detailed in Mousis et al. (2013), uses classical statistical mechanics to describe the macroscopic thermodynamic properties of clathrates to the molecular structure and interaction energies. It is based on the original model of van der Waals & Platteeuw (1959), which assumes that trapping of guest molecules into cages corresponds to the three-dimensional generalization of ideal localized adsorption. It is based on four hypotheses:

• The contribution of the host molecules to the free energy is independent of the clathrate occupancy. This implies that the guest species do not distort the cages.

• The cages are singly occupied, and guest molecules rotate freely within the cage.

• Guest molecules do not interact with each other.

• Classical statistics is valid, meaning that quantum effects are negligible.

Hence, the fractional occupancy of a guest molecule i for a given type of cage q (small or large) in a clathrate structure can be written as $$y_{i,q}=\frac{C_{i,q}{f}_{hydro,i}}{1+\sum _J{C}_{J,q}{f}_{hydro,J}},$$(B.1) where the sum in the denominator includes all the species present in the gas. C_i,q is the Langmuir constant of the species i in the cage of type q, f_hydro,i is the fugacity of the guest species i computed at the total gas pressure (see Appendix B.2).

The Langmuir constant C_i,q depends on the strength of the interaction between each guest species and each type of cage. It can be calculated by integrating the molecular potential energy within the cavity as $$C_{i,q}=\frac{4\pi }{k_{B}T}{\int }_0^{R_c}exp(-\frac{w_{i,q}(r)}{k_{B}T})r^{2}dr,$$(B.2) where R_c represents the radius of the cavity assumed to be spherical, k_B its the Boltzmann constant, and w_i,q(r) is the spherically averaged potential (here Kiahara or Lennard-Jones potential) representing the interactions between the guest molecules i and the water molecules forming the surrounding cage q. This potential w_i,q(r) can be written for a spherical guest molecule as (McKoy & Sinanoğlu 1963) $$w(r)=2zϵ[\frac{{\sigma }^{12}}{R_c^{11}r}({\delta }^{10}(r)+\frac{a}{R_c}{\delta }^{11}(r))-\frac{{\sigma }^6}{R_c^{5}r}({\delta }^4(r)+\frac{a}{R_c}{\delta }^5(r))],$$(B.3) with the mathematical function δ^N(r) in the form $${\delta }^N(r)=\frac{1}{N}[{(1-\frac{r}{R_c}-\frac{a}{R_c})}^{-N}-{(1+\frac{r}{R_c}-\frac{a}{R_c})}^{-N}].$$(B.4)

In Eq. B.3, z is the coordination number of the cell. Parameters z and R_c depend on the cage type q (small or large) and on the clathrate structure (I or II). Their values are given in Table B.1. The intermolecular parameters a, σ and ϵ describing the guest molecule i-water (i-W) interactions in the form of a Kihara or Lennard-Jones potential are listed in Table B.2.

Finally, the mole fraction $x_i^{clath}$ of a guest molecule i in a clathrate of a given structure (I or II) can be computed as $$x_i^{clath}=\frac{b_s{y}_{i,s}+b_l{y}_{i,l}}{b_s\sum _J{y}_{J,s}+b_l\sum _J{y}_{J,l}},$$(B.5) where b_s and b_l are the number of small and large cages per unit cell, respectively, for the clathrate structure we considered. The mole fractions of the enclathrated species is normalized to 1.

B.2 Calculation of the fugacities

To determine f_hydro,i, we first calculated the specific volume (v) of the considered mixture via the Redlich-Kwong EOS (Redlich & Kwong 1949): $$P=\frac{RT}{v-b}-\frac{a}{\sqrt{T}v(v+b)}$$(B.6) with $$a=0.4278\frac{R^2{T}_c^{2.5}}{P_c},b=0.08664\frac{R{T}_c}{P_c}$$(B.7) where R is the perfect gas constant, T the ambient temperature, T_c and P_c the critical temperature and pressure of the substance i Tab. A.2. In this case, a and b account for a pure substance. For mixtures, we took the average properties of the mixture into account by using the mixing rules $$\begin{array}{ll}{a}_m=\sum _i\sum _j{y}_i{y}_j{a}_{ij}& \text{\hspace{0.17em}with\hspace{0.17em}}{a}_{ij}=\sqrt{a_i{a}_j},\\ b_m=\sum _i{y}_i{b}_i,& \end{array}$$(B.8) where y_i the molar fraction of the gas in the atmosphere, normalized for the enclathrated species to 1. Then, we defined the mixture pressure P_m dissolved in water, which corresponds to the sum of the individual dissolved pressure P_i. This pressure is different than the partial pressure of the gas in the atmosphere. It is expressed as a function of Henry’s law coefficient as $$P_i=H_{hydro,i}\times y_i,$$(B.9) where H_hydro,i is the Henry’s law constant calculated for specie i at the surface (hydrostatic) pressure P_tot and temperature T. It is calculated as $$\ln \left(\frac{H_{hydro,i}}{H_{1atm,i}}\right)=\frac{{\overline{V}}_i}{RT}(P_{tot}-1.013\times {10}^5),$$(B.10) where H_1atm,i is the Henry’s constant of specie i calculated at 1 atmosphere and 0°C and ${\overline{V}}_i$ is the partial molar volume of specie i. The fugacity of the mixture at 1 atm pressure (f_1atm,m) can be computed as (Redlich & Kwong 1949) $$\begin{array}{rl}& \ln f_{1atm,m}=\frac{b_m}{v_m-b_m}+\ln \frac{RT}{v_m-b_m}-\\ & \frac{a_m}{R{T}^{3/2}}(\frac{1}{v_m+b_m}+\frac{1}{b_m}\ln \frac{(v_m+b_m)}{v_m})\end{array}.$$(B.11)

The fugacity of the mixture calculated at the surface pressure f_hydro,m is related to f_1atm,m via (Miller 1974) $$\ln \frac{f_{hydro,m}}{f_{1atm,m}}=-n\ln a_w-\frac{P_{tot}-1.013\times {10}^5}{RT}(v{V}_{H_{2}O}-V_{c,m}),$$(B.12) where a_w is the activity of water relative to pure liquid water (~1 here), n is the moles of water per mole of species i in the clathrate (n ~ 6), V_H2O is the molar volume of liquid water (1.8 × 10⁻⁵ m³) and V_c,m is the volume of clathrate that contains 1 mol of substance m. Finally, the fugacity coefficient ϕ of the mixture can be computed from $$f_{hydro,i}=\varphi \times P_i.$$(B.13)

The values of $T_c,P_c,{\overline{V}}_i,H_{1atm,i}$ and f_1atm,i used for each species i are taken from Mousis et al. (2013) Table 1.

B.3 Calculation of clathrate density

From the calculation of the fractional occupancies (y_i,q) in both small and larges cages, it is possible to derive the clathrate density based upon a unit crystal (Sloan & Koh 2008): $$\rho =\frac{N_w{M}_{H_{2}O}+\sum _{q=1}^C\sum _{i=1}^N{y}_{i,q}{b}_q{M}_i}{N_{Ava}{V}_{cell}},$$(B.14) where N_W is the number of water molecules per unit cell, N_Ana is Avogadro’s number (6.023 × 10²³), M_i the molecular weight of molecule i, y_i,q is the fractional occupancy of molecule i in type q cavity (small or large), b_q is the number of type q cavities per water molecules in unit cell, V_cell is the volume of an unit cell, N is the number of cavity type per unit cell (in this case 2) and C the number of component trapped in hydrate phase. N_w, b_q and V_cell values are given in Table B.1 for each structure (I and II).

References

All Tables

All Figures

Fig. 1 Equilibrium partial pressures of CO2 (dark blue) and NH3 (green) in Titan’s primordial atmosphere as a function of the total NH3 dissolved in the ocean. Each set of green and blue curves, distinguished by different line styles, represents a specified initial CO2 atmospheric pressure prior to chemical equilibrium between the ocean and atmosphere. The initial theoretical PCO2 values of 1, 5, 10, 20, and 30 bar correspond to CO2 bulk masses of 7.8 × 1019 kg, 3.79 × 1020 kg, 7.38 × 1020 kg, 1.37 × 1020 kg, and 1.92 × 1030 kg, and a dissolved CO2 mole fraction of 1.36 × 10−3, 6.78 × 10−3, 1.36 × 10−2, 2.71 × 10−2, and 4.07 × 10−2, respectively.

In the text

	Fig. 2 Evolution of the composition of Titan’s primordial atmosphere and ocean at equilibrium, as a function of surface temperature.
In the text

	Fig. 3 Evolution of atmospheric abundances as a function of clathrate crust thickness at 280 and 273.15 K. The dashed red line indicates the threshold at which the clathrate composition changes from CH4-rich to CO2-rich. The dashed black line represents the detection limit of Huygens’ GCMS.
In the text

	Fig. 4 Density of clathrates as a function of crust thickness, assuming formation at 273 K. The step for both temperatures highlights the transition from CH4-rich to denser CO2-rich clathrates.
In the text