The Effects of Methane Clathrates on the Thermal and Seismic Profile of Titan's Icy Lithosphere

We investigate the effects of methane clathrates on Titan’s thermal and seismic structure. The Dragonfly mission is planned to arrive at Titan in 2033 with a payload that includes a seismic package. The seismic instruments are tasked with recording seismic events and recovering the internal structure. Here, we explore whether differences in thermal and seismic profiles between a pure water ice shell and an ice shell with a clathrate lid could be detectable with seismic instrumentation. Due to their lower thermal conductivity, clathrates reduce the conductive lid thickness thus altering the thermal profile. The differences between seismic velocities of clathrates and pure water ice, coupled with changes in the thermal profile, indicate the clathrate lid will create lower seismic velocities, particularly for the upper 10 km of the surface ice shell. The differences in P and S velocity at the surface are 2.9% and 4.5%, respectively, and reach up to 8.4% (for both P and S) at a depth of 9.6 km. Due to changes in thermal profile, the seismic attenuation of the ice shell will change such that clathrates will suppress surface wave amplitudes relative to the pure water ice model. The clathrate lid will further create minor changes (≤2.0%) in the surface wave dispersion curves. Dragonfly, or other future seismic investigations, could provide evidence for or against the presence of a clathrate lid by constraining the thermal and seismic profile of Titan’s ice shell, by measuring the relative amplitudes of the surface to body waves, or by constraining the surface wave dispersion with high accuracy and precision.


Introduction
Titan, Saturn's largest moon, is of great interest for several reasons. Like other icy ocean worlds, such as Jupiter's Europa or Saturn's Enceladus, Titan has an icy shell overlying a subsurface ocean ). Titan's subsurface ocean may provide the requirements for life to currently exist. Beneath the subsurface ocean, could lie an additional layer of high-pressure ices such as Ice V and Ice VI (Journaux et al. 2020). Titan is unique compared to other icy ocean worlds, due to its thick methane-rich atmosphere (Niemann et al. 2005). Temperature and pressure conditions at the surface allow for liquid hydrocarbons to stably exist (Stofan et al. 2007). In addition to its presence in Titan's atmosphere, and on its surface, methane may also be present within Titan's icy lithosphere in clathrate structures (Choukroun et al. 2010).
Methane clathrates are nonstoichiometric compounds with a methane molecules hosted inside a cage-like structure of hydrogenbounded water molecules. For the remainder of the paper, we will use the term clathrates to refer to methane clathrates. Other forms of clathrates possibly found inside Titan in smaller amounts, such as ethane and carbon dioxide clathrates, have similar bulk and seismic properties, so we use measurements made on methane clathrates to calculate the clathrate layers' seismic properties. Due to their large stability range, it has been hypothesized that clathrates should remain stable throughout Titan's hydrosphere and cryosphere, and inside the liquid water ocean (Lunine & Stevenson 1985, 1987Choukroun et al. 2010). For this reason, CH 4 clathrates could play a major role in Titan's methane cycle (Atreya et al. 2006;Tobie et al. 2009;Mousis et al. 2011;Choukroun & Sotin 2012). Clathrates forming in Titan's interior and migrating upward would subsequently dissociate when reaching the surface, releasing methane into the atmosphere. In addition to contributing to Titan's methane cycle, clathrates may also influence the dynamics of the surface ice shell. Indeed, clathrate thermal and bulk transport properties vary significantly from those of pure water ice. Notably, the thermal conductivity of clathrates is ≈0.5 W m −1 K −1 compared to water ices 2.5 W m −1 K −1 , while the viscosity is nearly 20 times higher for clathrates than pure water ice Ih (Helgerud et al. 2003(Helgerud et al. , 2009Hand et al. 2006;Waite et al. 2007). Due to the important differences in transport properties compared to ice Ih, clathrates will likely alter the conductive ice shell thickness (Kamata et al. 2019;, thermal profile, and heat flux through the ice shell. Convection in the ice shell will affect how materials are exchanged between the surface, ice shell, and subsurface ocean, and could have implications for the habitability of the ocean. A thinner, warmer ice shell would allow for material to be more easily exchanged between the ice and the ocean. Furthermore, previous studies (e.g., Hemingway et al. 2013;Čadek et al. 2021;Carnahan et al. 2022) have invoked or shown that clathrates can help explain Titan's topography and density in the ice shell.
Future missions, such as NASA's Dragonfly mission (Barnes et al. 2021), will investigate Titan's surface and interior to assess the potential habitable environments. Dragonfly's payload includes a seismic package that has the goal of constraining Titan's seismicity and internal structure (Lorenz et al. 2019). Here, we model the effects methane clathrates have on Titan's thermal and seismic velocity profile to assess if a mission, such as Dragonfly, could help constrain the presence of methane clathrates within Titan's ice shell through seismology.

Methods
Previous studies investigated how a clathrate lid or the presence of clathrates may alter the thermal profiles of icy ocean worlds (Kamata et al. 2019;. These 2D numerical studies solve for thermal convection with set depths for total ice thickness. Here, we construct radial models of Titan's ice using the PlanetProfile  open-source software package, integrating the numerical results of  to include a 10 km clathrate lid at the top of a 100 km total thickness outer shell. We create a set of geophysically consistent models based on input parameters including moment of inertia, mass, total radius, and temperature at the ice-ocean boundary. We update PlanetProfile to allow for calculations using clathrate properties instead of pure water ice. This additional capability allows for computing the bulk and seismic properties of the resulting icy shells, and direct comparisons between the icy shell compositions and structures. For the purposes of this comparative study, we maintain the same ocean and silicate compositions, surface conditions, and ice-ocean temperatures for both models we describe below. We set the surface temperature to 94 K and the surface pressure to 0 MPa. The ocean composition is set at 10% MgSO 4 . The silicate interior is roughly 48.82% SiO 2 , 28.3%c MgO, 11% FeO, 8% Al 2 O 3 , and 3% CaO. PlanetProfile uses Perple X (Connolly 2009) to determine the equation of state for the silicate interior. We assume no metallic core. The modeled mineral assemblages do not include lower density minerals such as phlogopite that might explain the low density of Titan's interior (Néri et al. 2020), so the total mass of each Titan model exceeds the actual mass but would only affect calculations at the very center of Titan. This discrepancy does not affect the physical parameters of the ice, which assumes depth-dependent gravity based on the actual mass of Titan. We set the ice-ocean temperature to 260 K, which creates ice shells of roughly 100 km (Figure 1). To determine ice shell thickness, we calculated the pressure at which the ocean would freeze, given the ice-ocean temperature and ocean composition. The architecture of the PlanetProfile code builds the planet from the surface down; thus the ice shell thickness is determined by the depth at which the pressure from the overlying ice reaches the freezing pressure of the ocean. Since clathrates are slightly denser than water ice Ih, the clathrate lid reduces the ice shell thickness, but by less than 1 km. Once the thickness and boundary temperatures are set, we use the methods of Deschamps & Sotin (2001) to determine if convection occurs, and if it does, what are the resulting conductive and convection thicknesses and temperature profile. This approach is straightforward for pure water ice. We modify the approach for the clathrate lid as follows.
When the ice shell includes a clathrate lid, we use the values of heat flux (Q) and temperature of the convective layer (T conv ) from  to guide our approach. We benchmarked our pure water model against their pure water model to ensure they were in close agreement, justifying our use of heat flux and convective temperatures for the clathratelid model. For simplicity we assume a 2D geometry, rather than spherical when calculating heat flux. We determine the overall thickness of the conductive lid using Fourier's law (Equation (1)), where k is the thermal conductivity, and dT dZ is the thermal gradient between the surface and the convective temperature; Because clathrates have a thermal conductivity of ≈0.5 W m −1 K −1 and pure water ice has a thermal conductivity value closer to 2.5 W m −1 K −1 , we split Equation (1) into two calculations. The first calculation (Equation (2)) determines the temperature at the base of the clathrate lid (T clath ) using the heat flux, thermal conductivity of clathrates (k clath ), surface temperature (T surf ), and the thickness of the clathrate lid (h clath ); We then use T clath along with the convective temperature (T conv ), and thermal conductivity of ice (Andersson & Inaba 2005) to determine the thickness of the conductive ice  (h ice ). The total conductive lid thickness (e TBL ) is the sum of h ice and h clath . At the base of the ice shell there is a second conductive boundary layer between the convecting ice and the liquid ocean. We use Equation (3) to determine the thickness of this layer by finding the temperature difference between T conv and the temperature at the ice-ocean interface.
Once the pressure and temperature profile of the icy shell is determined, the bulk properties are calculated. We use the SeaFreeze library (Journaux et al. 2020) to calculate the pressure and temperature dependent properties of pure water ice phases (density, heat capacities, thermal expansion, seismic velocities). We use Helgerud et al. (2009) to determine bulk and seismic properties of clathrates based on pressure and temperature. With this approach we produce thermal and seismic velocity profiles for the two models, one with a pure water ice shell, and a second with a 10 km thick clathrate lid overlying a 90 km thick pure water ice shell. We select 10 km for the clathrate lid, as it represents a large lid thickness in which clathrates are still completely contained within the conductive portion of the ice shell, which simplifies our calculations and models.  showed the conductive lid thickness and heat flux are highly dependent on grain size for clathrate lids smaller than about 10 km; thus by using a 10 km ice shell we do not have to make assumptions about grain size. For thin ice shells, the grain size determines convection rigor, which can further reduce the stagnant lid thickness. Thicker ice shells are less dependent on convection rigor, grain size, and viscosity but more sensitive to the effects of thermal insulation. Lastly, and most importantly, a 10 km lid has been shown to produce the greatest change in thermal profile due to a thinning of the stagnant lid, and increase in convective temperature; thus we can measure the maximum effect on seismic velocities and resulting ground motions. We only alter the composition of the surface ice shells, meaning we assume any high-pressure ices beneath the ocean are pure water ices (e.g., Ice V, Ice VI).
Once we produce PlanetProfile seismic velocity profiles, the files are used as inputs for AxiSEM and InstaSEIS (Nissen-Meyer et al. 2014;van Driel et al. 2015), open-source codes that create seismic waveforms. To generate the waveforms, we set a dominant period of 2 s, a moment magnitude (M w ) of 3.1, source depth of 3 km, and a time length of 3600 s, which was sufficient to investigate seismic phases that pass through the ice shell, ocean, and any reflections off the interfaces between the surface ice, ocean, high-pressure ices, and silicate interior. The chosen magnitude represents a plausible event size (Hurford et al. 2020), but is completely arbitrary as we are comparing seismic events without any background noise and without instrument responses added. Receivers are placed across Titan's surface spanning 180°and spaced every 1°, allowing us to investigate seismic events from nearly any distance with high spatial resolution. No noise is added to the waveforms, to allow more direct comparison of the arrival times and waveforms. We use the open-source code TauP (Crotwell et al. 1999) to calculate the predicted arrival times of seismic phases.

Results
We investigate how a model with a clathrate lid compares to a pure water ice shell model. We find that clathrates will significantly reduce the thickness of the conductive lid from ≈30.7 to ≈10.5 km. Due to the reduction in the conductive lid thickness, the thermal profiles for pure water ice compared to a clathrate lid are significant (Figure 2(a)).
Seismic velocities within the icy shell are largely controlled by temperature, such that seismic velocities decrease with increasing temperature at at rate of ≈2.5 km s −1 K −1 for compressive waves and 1.5 km −1 s −1 K −1 for shear waves. Thus, changes in the thermal profile result in changes in the seismic profile (Figure 2(b)). At surface conditions where both models have the same pressure and temperature, clathrate seismic velocities are 2.9% lower for shear (V S ) and 4.6% lower for compressive (V P ) velocities compared to pure water ice (Figure 3). Due to the higher temperature gradient, the seismic velocity gradient is also larger for a clathrate lid than for pure water ice. This change causes the differences in velocities between the two models to increase with depth, maximizing at a depth of ≈10km, and then to decrease and eventually reach nearly zero once both models have pure water compositions, and both are convecting at similar (within 3°) temperatures. The clathrate lid also creates a small discontinuity at the clathrate-water ice interface, where the seismic velocity shows a slight increase before continuing with a negative velocity gradient. At its maximum, the difference in velocities between the models reached ≈8.4% for both P and S waves.
Temperature not only affects seismic velocities, but also seismic quality factors (Q s ): the inverse measurement of attenuation (Figure 2(c)). Following the approach of Cammarano et al. (2006), attenuation in the ice shell depends on temperature, Figure 4. Vertical component of an M w 3.1 event at epicentral distances of (a) 3°, (b) 50°, and at (c) 160°. Pure water (blue) is compared to a clathrate lid (red). To show both body and surface waves, the gray shaded areas have different y-axis limits. In the first plot, the clathrate-lid model has higher amplitudes than the pure water ice model. In the bottom two plots, the pure water ice model has greater amplitudes for surface waves but comparable amplitudes for body waves. such that warm ice will attenuate more seismic energy than cold brittle ice will. For the pure water ice model, Q s is relatively high for the upper 30 km. The high Q s allows surface waves to retain high amplitudes even at great distances (Figure 4). For the clathrate-lid model, Q s decreases rapidly within the upper 10 km. This reduction in Q s , greatly reduces surface wave amplitudes as distance between the event and receiver increases. The suppression of surface waves means a pure water ice model will produce average ground acceleration values up to 75 times greater than the clathrate-lid model ( Figure 5). Because clathrates act to suppress surface waves, while pure water ice maintains high amplitude surface waves, the amplitudes of body waves relative to surface waves could be used to identify a thin conductive lid, thus pointing to the presence of clathrates within the lid.
Changes in seismic velocity structure create changes in observable distances and arrival times for key seismic phases. For example, at a distance of 20°, the clathrate lid delays the arrival of P waves reflecting off the ice-ocean interface by 2 s; at 30°the S reflection is likewise delayed by ≈1.5 s. It is worth noting that differences on the order of a few seconds are typical uncertainties seen by the InSight Mars Quake Service (MQS; Giardini et al. 2020) so differences of a few seconds may be difficult to detect. See Section 4 for further details.
To further investigate the effects on surface waves, we invoke the open-source code, Mineos (Masters et al. 2011), to investigate surface wave dispersion. We supply the Mineos code with the two interior structure models from PlanetProfile. We look at both spheroidal and torodial modes with angular orders up to 3000 and frequencies up 200 Hz. We show in Figure 6 that there are slight differences in the dispersion characteristics between pure ice and clathrate-lid models. However, the differences between clathrate model values and pure ice model values are within 2% of each other, indicating that surface wave dispersion could be another tool to determine whether or not clathrates are present in Titan's surface ice shell, but would require very high precision and accuracy in measurements.

Discussion
Using PlanetProfile with updated parameters from  we show that clathrates will reduce the thickness of the conductive lid, thus altering the thermal profile of Titan. Because seismic velocities and quality factors depend on temperature, a clathrate lid overlying pure water ice will produce different seismic responses than a model composed purely of water ice. At most, the seismic velocities in the clathrate lid will be reduced by 8.4% compared to pure water. On Earth, the average seismic velocity is well known (Dziewonski & Anderson 1981;Kennett & Engdahl 1991) and tomographic studies have been able to discern variations down to a few percentages (5; Anderson & Dziewonski 1984;Grand et al. 1997). However, previous planetary missions like InSight and the Apollo seismic experiments have not benefited from numerous and dense seismic arrays nor frequent teleseismic events that occur on Earth. Due to the limited numbers of stations and scarcity of large events, the inferred internal structures for the Moon (Garcia et al. 2019) and Mars Knapmeyer-Endrun et al. 2021;Stähler et al. 2021) have uncertainties in velocities that exceed 10%, and corresponding errors in boundary locations on the order of several kilometers. To distinguish between the models presented here, a seismic investigation would need to recover an internal structure with uncertainties less than 8%, ideally less than 5%. It may be possible to reach those levels of uncertainty, but it would require Figure 5. Rms comparison of pure water ice shell to clathrate-lid waveforms. The rms ratio for the radial (blue) and vertical (red) components are plotted vs. distance. Values greater than 1.0 (dashed gray) indicate the clathrate-lid model has greater amplitudes than the pure water ice model. recording numerous high-quality seismic events for which the source location could be derived, and the seismic phases need to be well constrained. The InSight MQS team typically reports uncertainties in seismic phase arrivals ranging from ≈1 s up to tens of seconds for some reflection phases. The locations and depths also have uncertainties on the order of several degrees in epicentral distance and several kilometers in depth, even for the highest quality events.
In addition to arrival times of body waves, the surface waves will also be affected by the presence of a clathrate lid. The increased attenuation of the comparatively warmer ice will reduce the amplitudes of surface waves traveling through Titan's clathrate lid compared to those traveling through colder, brittler, pure water ice. Because body waves do not exclusively travel through the surface ice layers, the amplitudes of body waves for both models are comparable. The relative amplitudes of surface waves and body waves may help indicate if a clathrate lid is present, particularly at longer distances. In the shallow parts of the ice shell, the temperature-based model of seismic attenuation from Cammarano et al. (2006) predicts high Q values, but as discussed in Panning et al. (2018, attenuation in very low temperature ice is not well constrained. If the cold ice has a much lower quality factor, the difference between the two models would be reduced. Surface waves can also be suppressed if the source event has a deep source. If the quality of waveforms is sufficiently high, the arrival times and amplitudes of seismic phases may help indicate if clathrates are present. Figures 4 and 5 indicate there are subtle differences between the body waves, but significant differences in surface wave. Surface wave dispersion could also be implemented, but the differences between the two models are subtle. However, these waveforms and Mineos models include no added background noise, precisely known location, and timing of the event, and assume Titan is radially consistent with no porosity or lateral heterogeneity. The atmosphere, Titan's lakes, and surface roughness will create seismic background noise and create scattering effects Lorenz et al. 2021). Topography along the ice-atmosphere and ice-ocean boundary can also alter the relative arrival times of seismic phases. The thermal profile of Titan might also be altered by organic sands, which can lower the thermal conductivity more than clathrates can (Schurmeier & Dombard 2018). The sands may mimic the effects of clathrates, thus making it hard to distinguish between the materials. Lastly, physical and/or chemical heterogeneity within Titan's ice shell will also create variations in travel times and waveform amplitudes. The anomalies may be caused by changes in porosity, changes in concentration or depth of the clathrate lid, organic sand infilling craters, or entrained liquid water or methane. For these reasons, the differences in the waveforms of the two models may be trivial and nonunique once background noise and the above factors are considered.
For a seismometer to constrain the presence of clathrates within Titan's ice shell, numerous high-quality events need to be recorded. The arrival times and waveforms can be used to recover the seismic velocity profile of Titan, which can inform on the thicknesses of the conductive lid. If the seismic velocity structure can be recovered with high precision, it might be possible to determine whether or not clathrates are present. A better Figure 6. We compare the dispersion curves for group (solid) and phase (dashed) velocity using the pure ice model (blue) and the clathrate-lid model (red). The differences are subtle and minor. understanding of Titan's current internal structure, thermal profile, and composition will help inform how Titan might have evolved to its current state and place constraints on Titan's history.

Conclusions
Methane clathrates will insulate the icy shell of Titan. Compared to a shell composed of pure water ice, an ice shell with a clathrate lid will have a reduced conductive lid thickness and smaller seismic velocities. The difference in seismic velocities could reach up to ≈8.4%. Changes in the seismic velocity will result in changes in arrival times and amplitudes of seismic phases. However, these changes are subtle and would require high-quality seismic events. Background noise, heterogeneity, and topography may make it difficult to pick arrival times and will introduce uncertainty in the seismic velocity profile. If the seismic profile is known with 5% uncertainty, or if the conductive lid thickness and total ice thickness is well constrained, seismology can help constrain whether clathrates are present in Titan's surface ice shell.
A part of the research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). Copyright 2022. We thank two anonymous reviewers for their constructive remarks that helped improve this manuscript.