**Key Points:**

* We studied the evolution of impact melt pools in Titan's ice shell using numerical simulations of two-phase thermal convection

* While most melt pools never reach the ocean, we observe surface-to-ocean exchange in a small part of the investigated parameter space

* Our results suggest that Dragonfly may be able to sample melt deposits where organics reacted with liquid water to produce biomolecules

**Correspondence to:**

K. Kalousová,

klara.kalousova@mff.cuni.cz

**Citation:**

Kalousová, K., Wakita, S., Sotin, C., Neish, C. D., Soderblom, J. M., Souček, O., & Johnson, B. C. (2024). Evolution of impact melt pools on Titan. *Journal of Geophysical Research: Planets*, 129, e2023JE008107. https://doi.org/10.1029/2023JE008107

Received 15 SEP 2023

Accepted 29 JAN 2024

© 2024. The Authors.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_1.png" /></figure>

# Evolution of Impact Melt Pools on Titan

Klára Kalousová<fnref n="1" /> ID, Shigeru Wakita²,³ ID, Christophe Sotin⁴ ID, Catherine D. Neish<fnref n="5" />,⁶ ID, Jason M. Soderblom² ID, Ondřej Souček⁷ ID, and Brandon C. Johnson³,⁸ ID

**Abstract** Titan is an ocean world with a dense atmosphere, where photochemistry produces complex organic molecules that fall to the surface. An important astrobiological question is whether this material can mix with water and form molecules of biological interest. Large impacts heat the moon's subsurface and create liquid water melt pools. A recent study investigated impacts into Titan's clathrate-covered ice shell. Methane clathrates are stable at Titan's surface conditions and have low thermal conductivity, making them efficient insulators that can lead to steep thermal gradients and a thin stagnant lid. The authors showed that the clathrate layer thickness primarily influences the melt distribution, while its volume is governed by the impactor size. Here, we investigate the fate of melt formed during an impact into a clathrate-covered ice shell. Our results show two different behaviors: in cases when less melt is produced, the subsurface melt pool remains close to the surface and freezes on timescales ≲25 kyr; in cases when larger volumes of melt are produced, a downward-oriented transport of the molten material occurs. As it descends, part of the melt freezes but some may reach the ocean within a few kyr under certain conditions; vertical impacts, high surface porosity, low viscosity, and tidal heating all favor this surface-to-ocean exchange. While providing insights on parameters that allow a subsurface melt pool to remain liquid beneath a Selk-sized crater for a few kyr, this study suggests that Dragonfly may be able to sample melt deposits where organics reacted with liquid water to produce biomolecules.

**Plain Language Summary** Titan, Saturn's largest moon, harbors a subsurface ocean beneath its ice shell. The moon also has an atmosphere, which is rich in large organic molecules that settle onto its surface. When atmospheric methane reacts with surface water ice, it forms methane clathrate. A clathrate layer atop Titan's ice shell affects the formation of Titan's impact craters, as it is both stronger than ice and a better insulator. Here, we study the effect of this clathrate layer on the formation and subsequent fate of melt pools produced by impacts on Titan. We investigate whether the melt descends through the ice shell to reach the ocean or remains near the surface and freezes. Our results show only a limited range of scenarios where impact melt reaches the ocean. In the majority of models, impact melt freezes near the surface within short timescales, ranging from a few thousands to tens of thousands of years. This implies that surface organic molecules may have interacted with subsurface melt pools. These results provide a positive outlook for NASA's Dragonfly mission, which will explore Titan's Selk crater in search of organic materials that have reacted with water to potentially form molecules of biological interest.

## 1. Introduction

Titan, Saturn's largest moon, harbors a deep ocean of liquid water below its ice shell, which raises the question of its habitability (Hendrix et al., 2019). The presence of a deep ocean has been inferred from two observations. First, the measured obliquity of 0.3° (Meriggiola et al., 2016; Stiles et al., 2008) is about three times larger than the value expected for an entirely solid Titan (0.12°; Bills & Nimmo, 2008), which has been interpreted as being indicative of the presence of an internal ocean (Baland et al., 2011, 2014; Bills & Nimmo, 2008, 2011). Second, the Cassini spacecraft was able to determine the long wavelength variations of gravity. The component of the degree 2 gravity field, which is due to the tidal forces during Titan's eccentric orbit around Saturn, has been found to be quite large. This is best explained by a decoupling between the outer ice shell and the interior (Durante et al., 2019). However, the depth and thickness of this ocean are not well constrained.

<footnote n="1">¹Department of Geophysics, Charles University, Prague, Czech Republic, ²Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA, USA, ³Department of Earth, Atmospheric, and Planetary Sciences, Purdue University, West Lafayette, IN, USA, ⁴LPG, Nantes University, Nantes, France,</footnote>

<footnote n="5">⁵Department of Earth Sciences, The University of Western Ontario, London, ON, Canada, ⁶The Planetary Science Institute, Tucson, AZ, USA, ⁷Mathematical Institute, Charles University, Prague, Czech Republic, ⁸Department of Physics and Astronomy, Purdue University, West Lafayette, IN, USA</footnote>
In addition to water, an environment needs essential elements such as carbon and nitrogen to be habitable (Cockell et al., 2016). Titan has an abundance of these elements; it is the only moon in the Solar System known to have a dense atmosphere, which is composed mainly of nitrogen with a few percent of methane. Organic molecules form by photochemistry of these two compounds (e.g., Loison et al., 2019; Vuitton et al., 2019; Yung et al., 1984). The Cassini mission, as well as telescopic observations, have identified more than 20 molecules, primarily organic, in Titan's upper atmosphere (MacKenzie et al., 2021). A flux of oxygen and OH from Enceladus may also provide some oxygen (Hörst, 2017). Laboratory analogs of Titan's organic molecules (known as “tholins”) are composed of thousands of different organic molecules (Hörst, 2017; Maillard et al., 2018) including polyaromatic hydrocarbons. Some of these molecules represent the building blocks for life as we know it, while others could provide the energy required for simple organisms. As the 94 K surface temperature of Titan is much too low for life as we know it to develop, those organics would have to be transported to more favorable conditions, such as the ones existing in subsurface liquid water reservoirs.

One way to transfer Titan's carbon-rich surface material to its interior is during an impact. Compared to other icy moons of similar size like Ganymede and Callisto, there are few impact craters on Titan's surface. Hedgepeth et al. (2020) identified only 90 potential impact craters with diameters between 4 and 500 km on the Cassini Synthetic Aperture Radar images that cover about 70% of the surface. Therefore, Titan's surface is geologically young, between 200 and 1,000 Myr (Hedgepeth et al., 2020; Neish & Lorenz, 2012). The impactor that created Menrva, Titan's largest impact crater, could have breached Titan's ice shell (Crosta et al., 2021), bringing organic-rich surface material to the ocean. Smaller impactors are not sufficient to breach Titan's ice shell, but they impart enough energy to melt a fraction of the ice crust. Such impactors can form a subsurface pool where ~10% of the ejected organics will be lightly shocked and deposited back inside the crater (Artemieva & Lunine, 2003). Such environments are interesting astrobiological targets, as they present an environment where biological molecules such as amino acids can form (e.g., Neish et al., 2010). One such crater, Selk crater, has been chosen as the target of the fourth NASA New Frontiers mission, Dragonfly (Turtle et al., 2017).

In their recent paper, Carnahan et al. (2022) showed that, on Europa, impact melt pools could sink into the ice shell due to negative phase buoyancy (since water is denser than ice). In the present study, we investigate whether melt pools could similarly sink in Titan's ice shell. Due to efficient tidal heating, Europa is expected to have a thinner stagnant lid and warmer interior than Titan, when a pure water-ice shell is considered. However, methane clathrates are stable at Titan's surface conditions (Lewis, 1971) and may form a few km thick outer layer (e.g., Choukroun & Sotin, 2012; Lunine & Stevenson, 1987; Tobie et al., 2006). Owing to their very low thermal conductivity (about one order of magnitude smaller than that of ice at Titan's surface temperature), the presence of clathrates results in high temperature gradients (Kalousová & Sotin, 2020b). This steep thermal gradient, in turn, reduces the thickness of the stagnant lid (as compared to ice shells without clathrates), which could help the sinking of the subsurface melt pool since the melt is embedded in warmer, less viscous ice. This may facilitate the transfer of material, including surface organics, through Titan's ice shell.

In this study, we use initial conditions for crustal temperature, melt volume and distribution, and methane-clathrate distribution, from impact simulations that resulted in a Selk-sized crater (Wakita et al., 2023). We tailored a previously developed two-phase numerical model (Kalousová & Sotin, 2018, 2020a; Kalousová et al., 2018) to conditions relevant to Titan's ice I shell and followed the evolution of the impact-generated melt pool. This model is described further in Section 2. We describe the results of the numerical simulations in Section 3, and discuss the model assumptions and the implications of our work in Section 4. Finally, we conclude this study in Section 5.

## 2. Model

### 2.1. Impact Model

Impacts can produce pools of molten bedrock during the crater formation process. To model impact-induced melt pools, we perform impact simulations using the iSALE-2D shock physics code (Amsden et al., 1980; Collins et al., 2004; Wünnemann et al., 2006). Wakita et al. (2023) studied 3–4 km-diameter-impactors striking a 5–15 km thick methane clathrate layer with a temperature profile corresponding to 1 mm ice grains from Kalousová and Sotin (2020b). In this work, we investigate additional model runs using the same resolution (50 m), impact velocity (10.5 km/s), and material input parameters as Wakita et al. (2023). We add runs for 3–4 km-diameter-impactors into pure water ice (i.e., no methane-clathrate layer), and 4 km-diameter-impactors into a 10 km thick
methane-clathrate layer with a temperature profile assuming 0.5 and 2 mm ice grains (Kalousová & Sotin, 2020b). Since the near-surface porosity in the target increases the amount of melt produced by an impact (Artemieva & Lunine, 2005; Wünnemann et al., 2008), we also examine a 4 km-diameter-impactor into a 10 km thick methane-clathrate layer of 10% and 20% porosity with a temperature profile assuming 1 mm ice grains (Kalousová & Sotin, 2020b). Note, that we only consider the effects of this near-surface, gas-filled porosity on the amount of melt (see below). It is not to be mistaken for the (water-filled) porosity of the two-phase code (described in the next section)—that is why we avoid using the term “porosity” for the quantity $\phi$ and we call it water volume fraction or water content. Following Johnson et al. (2017) and assuming a reference viscosity of $10^{15}$ Pa s for water ice at 273 K, we find that porosity in the ice below a 5 km thick clathrate layer will be removed by viscous pore closure (e.g., Fowler, 1985) on a time scale of ~50 kyr. This time scale is much shorter in the case of a 10 or 15 km thick clathrate layer. Since any porosity induced by earlier impacts (Wiggins et al., 2022) or other processes would be removed on a relatively short timescale, we take the porosity of the water-ice layer to be zero at the time of the impact. Note that the methane-clathrate layer can remain porous over longer time scales because it is much stronger (Durham et al., 2003), colder, and at lower overburden pressure than the underlying water ice. To account for the melted material, we followed the same methods as in Wakita et al. (2023), that is, we regard material as melted when the peak pressure of tracer particles exceeds the threshold of incipient melting (5.34 GPa for 0%, 3.45 GPa for 10%, and 2.31 GPa for 20% porosity). Note that the porous target has a larger volume change and more work is done during shock compression resulting in additional heating as compared to non-porous material shocked to the same pressure (Melosh, 1989). As a result, a porous target requires a lower peak pressure to melt than a non-porous target.

## 2.2. Two-Phase Flow Model

### 2.2.1. Governing Equations

To model the evolution of the melt pools forming in Titan's impact craters, we treat the ice shell as a mixture of solid ice (subscript i) and liquid water (subscript w). We define the volume fraction of water in the whole mixture as $\phi$ and, consequently, the ice volume fraction can be written as $1 - \phi$. To describe the evolution of the mixture, we use the two-phase flow formalism (e.g., Bercovici et al., 2001; McKenzie, 1984; Souček et al., 2014; Šrámek et al., 2007). We further assume that the two phases are in thermal equilibrium and that water is locked in the ice matrix and advected with it (i.e., there is no differential motion between the two phases, we call this approximation an “impermeable limit”). To capture the evolution of the thermal state as well as the phase change between the solid ice and liquid water, we use the enthalpy method (e.g., Carnahan et al., 2022; Jordan & Hesse, 2015; Katz, 2008, and Appendix A). We solve the conservation equations of mass, linear momentum, and enthalpy of the whole mixture. Assuming that (a) the mixture as a whole is incompressible, (b) it deforms viscously, and (c) the inertial forces are negligible with respect to viscous forces (infinite Prandtl number), we can write the governing equations in the following form (see Kalousová et al., 2016, for Equations 1a–1b and Appendix A1 for Equation 1c):

$$ \nabla \cdot \mathbf{v} = 0, \qquad (1a) $$

$$ -\nabla \Pi + \nabla \cdot \boldsymbol{\sigma} - \phi \Delta \rho \mathbf{g} - (1 - \phi) \rho_i \alpha_i \Delta T \mathbf{g} = 0, \qquad (1b) $$

$$ \rho \frac{\partial h}{\partial t} + \rho \mathbf{v} \cdot \nabla h - \phi \Delta \rho \mathbf{v} \cdot \mathbf{g} = \nabla \cdot (k \nabla T) + Q. \qquad (1c) $$

Here $\mathbf{v}$ is the mixture velocity, $\Pi$ is the dynamic pressure, $\boldsymbol{\sigma} = \eta(\nabla\mathbf{v} + \nabla\mathbf{v}^T)$ is the deviatoric stress tensor, with $\eta$ the effective viscosity (see Section 2.2.2 for more details), $\Delta\rho = \rho_i - \rho_w$ is the density difference between the ice and water, $\mathbf{g}$ is the gravity acceleration, $\alpha_i$ is the ice thermal expansivity, $\Delta T$ is the temperature difference from the reference value $T_0$, $\rho = \phi\rho_w + (1 - \phi)\rho_i$ is the mixture density, $h$ is the mixture enthalpy, $k$ is the thermal conductivity (Section 2.2.2), and $Q$ is internal heating (Section 2.2.2).

We define the specific (and gravity reduced) enthalpy of the mixture, $h(T, p, \phi)$, as (see derivation in Appendix A2):
$$h(T, p, \phi) = \phi L_0 + (\phi c_w^P + (1-\phi)c_i^P)(T-T_0) + \left[ \phi \frac{1-\alpha_w T}{\rho_w} + (1-\phi)\frac{1-\alpha_i T}{\rho_i} - \frac{1}{\rho_i} \right] (p-p_0), \quad (2)$$

with T the mixture temperature, p the pressure, L₀ the latent heat at reference temperature T₀ and pressure p₀, and cᵢ, αᵢ, and ρᵢ the heat capacity, thermal expansivity, and density of a particular phase k (ice or water). We also define the maximum enthalpy of solid ice, as:

$$h_i^{\max}(p) = \max_{T=T_m, \phi=0} h(T, p, \phi) = c_i^P(T_m - T_0) - \frac{\alpha_i T_m}{\rho_i}(p - p_0). \quad (3)$$

Depending on the actual value of enthalpy h, the mixture can either be completely solid (h < hᵢ<sup>max</sup>) or partially molten (h ≥ hᵢ<sup>max</sup>). Note that we do not treat the case of completely molten material (φ = 1) warmer than the melting temperature (T > T<sub>m</sub>). The temperature and water content are evaluated from the enthalpy as follows (see Appendix A3 for derivation):

$$h < h_i^{\max}, \quad T = T_m + \frac{h - h_i^{\max}}{c_i^P - \frac{\alpha_i}{\rho_i}(p - p_0)}, \quad \phi = 0, \quad (4a)$$

$$h \ge h_i^{\max}, \quad T = T_m, \quad \phi = \frac{h - h_i^{\max}}{L_0 + (c_w^P - c_i^P)(T_m - T_0) + \left[ \frac{1 - \alpha_w T_m}{\rho_w} - \frac{1 - \alpha_i T_m}{\rho_i} \right]}. \quad (4b)$$

### 2.2.2. Material Parameters

Depending on the local (p, T, σ) conditions, ice deforms by various mechanisms (volume and boundary diffusion, dislocation, basal slip, and grain-boundary sliding) that are combined in the composite creep (Goldsby & Kohlstedt, 2001) as (e.g., Kalousová & Sotin, 2020b):

$$\frac{1}{\eta_i} = \frac{1}{\eta_{\text{diff}}^{\text{v}}} + \frac{1}{\eta_{\text{diff}}^{\text{b}}} + \frac{1}{\eta_{\text{disl}}} + \frac{1}{\eta_{\text{gbs}} + \eta_{\text{bs}}}. \quad (5)$$

The particular creep viscosities are computed as:

$$\eta_j = \frac{1}{2} \frac{d^m T^{l_j}}{A_j \sigma_{\text{II}}^{n_j-1}} \exp\left(\frac{E_j}{RT}\right), \quad (6)$$

where the exponents m<sub>j</sub>, n<sub>j</sub>, and l<sub>j</sub>, prefactor A<sub>j</sub>, and activation energy E<sub>j</sub> are specific for each creep mechanism, d is the grain size, σ<sub>II</sub> is the second invariant of the deviatoric stress tensor, and R the universal gas constant. Note that the effective viscosity strongly depends on grain size d and therefore by choosing the grain size, we fix the reference viscosity. When the melt is present, viscosity can be reduced by several orders of magnitude (De La Chapelle et al., 1999). To account for this effect, we follow the approach of Tobie et al. (2003) but only decrease the viscosity by ~2 orders of magnitude in maximum (i.e., η<sub>min</sub> ~ 10<sup>12</sup> Pa s, see also Carnahan et al., 2022). Similarly, for numerical reasons, we limit the maximum viscosity to ~10<sup>24</sup> Pa s. We thus prescribe the mixture viscosity as follows:

$$\eta = \min\{\max[\eta_i \exp(-\gamma_m \phi), \eta_{\min}], \eta_{\max}\}. \quad (7)$$

The mixture thermal conductivity k is a function of temperature T, phase (solid or liquid), and composition (ice or clathrate). We first determine the conductivity of the solid phase (ice + clathrate) as (Kalousová & Sotin, 2020b)

$$k_{\text{h,c}} = (1 - C) \left( \frac{k_i^1}{T} + k_i^2 \right) + Ck_c, \quad (8)$$
<table><caption>**Table 1**
Model Parameters</caption><thead><tr><th>Symbol</th><th>Physical property</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>α<sub>i</sub></td><td>Ice thermal expansivity</td><td>1.6 × 10<sup>-4</sup></td><td>1/K</td></tr><tr><td>α<sub>w</sub></td><td>Water thermal expansivity</td><td>1.6 × 10<sup>-4</sup></td><td>1/K</td></tr><tr><td>c<sup>p</sup><sub>i</sub></td><td>Ice isobaric heat capacity</td><td>2,100</td><td>J/kg/K</td></tr><tr><td>c<sup>p</sup><sub>w</sub></td><td>Water isobaric heat capacity</td><td>4,184</td><td>J/kg/K</td></tr><tr><td>ρ<sub>i</sub></td><td>Ice density</td><td>920</td><td>kg/m<sup>3</sup></td></tr><tr><td>ρ<sub>w</sub></td><td>Water density</td><td>1,000</td><td>kg/m<sup>3</sup></td></tr><tr><td>g</td><td>Gravity</td><td>1.35</td><td>m/s<sup>2</sup></td></tr><tr><td>k<sup>1</sup><sub>i</sub></td><td>Constant of ice thermal conductivity (Equation 8)<sup>a</sup></td><td>488.12</td><td>W/m</td></tr><tr><td>k<sup>2</sup><sub>i</sub></td><td>Constant of ice thermal conductivity (Equation 8)<sup>a</sup></td><td>0.4685</td><td>W/m/K</td></tr><tr><td>k<sub>c</sub></td><td>Methane clathrate thermal conductivity<sup>b</sup></td><td>0.5</td><td>W/m/K</td></tr><tr><td>k<sub>w</sub></td><td>Liquid water thermal conductivity</td><td>0.5</td><td>W/m/K</td></tr><tr><td>L<sub>0</sub></td><td>Latent heat of melting</td><td>334</td><td>kJ/kg</td></tr><tr><td>H</td><td>Ice shell thickness</td><td>60</td><td>km</td></tr><tr><td>h<sup>c</sup><sub>c</sub></td><td>Clathrate layer thickness</td><td>0, 5, 10, 15</td><td>km</td></tr><tr><td>Q<sup>c</sup><sub>m</sub></td><td>Maximum heating amplitude</td><td>0, 5 × 10<sup>-8</sup></td><td>W/m<sup>3</sup></td></tr><tr><td>ω</td><td>Orbital frequency</td><td>4.6 × 10<sup>-6</sup></td><td>1/s</td></tr><tr><td>μ</td><td>Shear modulus of ice</td><td>3</td><td>GPa</td></tr><tr><td>η<sub>M</sub></td><td>Maxwell peak</td><td>6.5 × 10<sup>14</sup></td><td>Pa s</td></tr><tr><td>η<sub>min</sub></td><td>Minimum viscosity</td><td>10<sup>12</sup></td><td>Pa s</td></tr><tr><td>η<sub>max</sub></td><td>Maximum viscosity</td><td>10<sup>24</sup></td><td>Pa s</td></tr><tr><td>T<sub>s</sub></td><td>Surface temperature</td><td>94</td><td>K</td></tr><tr><td>T<sub>b</sub></td><td>Ocean temperature</td><td>266.9</td><td>K</td></tr><tr><td>R</td><td>Universal gas constant</td><td>8.314</td><td>J/K/mol</td></tr><tr><td>d<sup>c</sup></td><td>Grain size</td><td>0.5, 1.0, 2.0</td><td>mm</td></tr><tr><td>a</td><td>Melting curve parameter (Equation 11)<sup>d</sup></td><td>-395.2</td><td>MPa</td></tr><tr><td>c</td><td>Melting curve parameter (Equation 11)<sup>d</sup></td><td>9</td><td>-</td></tr><tr><td>T<sub>0</sub></td><td>Reference temperature<sup>d</sup> (Equation 11)</td><td>273.16</td><td>K</td></tr><tr><td>p<sub>0</sub></td><td>Reference pressure (Equation 11)<sup>d</sup></td><td>0</td><td>MPa</td></tr><tr><td>γ<sub>m</sub></td><td>Viscosity decrease due to melt<sup>e</sup></td><td>45</td><td>-</td></tr><tr><td>ε</td><td>Titan's surface emissivity</td><td>0.78</td><td>-</td></tr><tr><td>σ</td><td>Stefan-Boltzman constant</td><td>5.670367 × 10<sup>-8</sup></td><td>kg/s<sup>3</sup>/K<sup>4</sup></td></tr><tr><td>d<sup>c</sup><sub>i</sub></td><td>Impactor diameter</td><td>3, 3.5, 4</td><td>km</td></tr><tr><td>ξ<sup>c</sup></td><td>Surface porosity</td><td>0, 10, 20</td><td>%</td></tr></tbody></table>

with C the composition function, which is equal to 0 in ice and to 1 in clathrate, $k_i^1$ and $k_i^2$ the ice conductivity parameters (Hobbs, 1974), and $k_c$ the clathrate conductivity (considered constant). In the second step, we combine the solid conductivity with the water conductivity $k_w$ as follows:

$$k = \phi k_w + (1 - \phi) k_{i,c} \qquad (9)$$

In some simulations, a parameterized volumetric tidal heating considering Maxwell viscoelasticity is prescribed following Tobie et al. (2003) as:

<footnote><sup>a</sup>Hobbs (1974). <sup>b</sup>Sloan et al. (2007). <sup>c</sup>Indicates the varied parameters. <sup>d</sup>Chizhov (1993). <sup>e</sup>Tobie et al. (2003).</footnote>
<table><caption>**Table 2**
Creep Parameters of Ice Based on Goldsby and Kohlstedt (2001)</caption><thead><tr><th></th><th>T (K)</th><th>A<sub>j</sub> (m<sup>m</sup>/Pa<sup>q</sup>/s)</th><th>n<sub>j</sub></th><th>m<sub>j</sub></th><th>l<sub>j</sub></th><th>E<sub>j</sub> (kJ/mol)</th></tr></thead><tbody><tr><td>η<sub>diff</sub><sup>v</sup></td><td></td><td>9.1 × 10<sup>-4</sup></td><td>1.0</td><td>2.0</td><td>1.0</td><td>59</td></tr><tr><td>η<sub>diff</sub><sup>b</sup></td><td></td><td>6.4 × 10<sup>-4</sup></td><td>1.0</td><td>2.0</td><td>1.0</td><td>49</td></tr><tr><td rowspan="2">η<sub>disl</sub></td><td>≤258</td><td>4.0 × 10<sup>-19</sup></td><td>4.0</td><td>0.0</td><td>0.0</td><td>60</td></tr><tr><td>>258</td><td>7.9 × 10<sup>5</sup></td><td>4.0</td><td>0.0</td><td>0.0</td><td>180</td></tr><tr><td>η<sub>bs</sub></td><td></td><td>2.2 × 10<sup>-7</sup></td><td>2.4</td><td>0.0</td><td>0.0</td><td>60</td></tr><tr><td rowspan="2">η<sub>gbs</sub></td><td>≤258</td><td>6.2 × 10<sup>-14</sup></td><td>1.8</td><td>1.4</td><td>0.0</td><td>49</td></tr><tr><td>>258</td><td>5.6 × 10<sup>15</sup></td><td>1.8</td><td>1.4</td><td>0.0</td><td>192</td></tr></tbody></table>

$$Q = \frac{2Q_m}{\frac{\eta}{\eta_M} + \frac{\eta_M}{\eta}}, \qquad (10)$$

with $Q_m$ the maximum heating amplitude and $\eta_M = \mu/\omega$ the viscosity that maximizes the heating (the Maxwell peak) with $\mu$ the shear modulus and $\omega$ the orbital frequency. The tidal heating thus depends on stress $\sigma$, grain size $d$, the mixture temperature $T$ and water content $\phi$.

Finally, the melting temperature is evaluated as follows (Simon & Glatzel, 1929):

$$T_m = T_0 \left( \frac{P - P_0}{a} + 1 \right)^{1/c}. \qquad (11)$$

All parameter values are given in Tables 1 and 2.

### 2.2.3. Initial and Boundary Conditions

We use the results obtained by the iSALE shock physics code as an input for our two-phase code. In particular, we use the temperature, melt, and methane-clathrate distribution, as well as the crater shape, at 3,000 s after the impact. However, we do not evolve the surface shape in time as crater relaxation would occur on longer time scales (few tens of Myr, e.g. Schurmeier & Dombard, 2018; Schurmeier et al., 2023) than investigated here (few tens of kyr).

<table><caption>**Table 3**
Simulations and Their Results</caption><thead><tr><th>#</th><th>d<sub>i</sub> (km)</th><th>h<sub>c</sub> (km)</th><th>d (mm)</th><th>η<sub>b</sub> (Pa s)</th><th>ξ (%)</th><th>Q<sub>m</sub> (W/m<sup>3</sup>)</th><th>reg</th><th>d<sub>c</sub> (km)</th><th>V<sub>0</sub> (km<sup>3</sup>)</th><th>t<sub>f</sub> (kyr)</th><th>t<sub>1/2</sub> (kyr)</th><th>d<sub>m</sub><sup>max</sup> (km)</th></tr></thead><tbody><tr><td>1</td><td>3</td><td>0</td><td>1.0</td><td>6.2e14</td><td>0</td><td>0</td><td>1</td><td>3.0</td><td>178</td><td>12.1</td><td>0.31</td><td>5.1</td></tr><tr><td>2</td><td>3</td><td>5</td><td>1.0</td><td>6.2e14</td><td>0</td><td>0</td><td>1</td><td>1.5</td><td>157</td><td>6.5</td><td>0.17</td><td>4.7</td></tr><tr><td>3</td><td>3</td><td>10</td><td>1.0</td><td>6.2e14</td><td>0</td><td>0</td><td>2</td><td>1.1</td><td>150</td><td>10.2</td><td>0.15</td><td>16.8</td></tr><tr><td>4</td><td>3</td><td>15</td><td>1.0</td><td>6.2e14</td><td>0</td><td>0</td><td>1</td><td>2.1</td><td>140</td><td>12.7</td><td>0.35</td><td>1.7</td></tr><tr><td>5</td><td>3.5</td><td>0</td><td>1.0</td><td>6.2e14</td><td>0</td><td>0</td><td>1</td><td>3.3</td><td>293</td><td>16.3</td><td>0.36</td><td>6.2</td></tr><tr><td>6</td><td>3.5</td><td>5</td><td>1.0</td><td>6.2e14</td><td>0</td><td>0</td><td>1</td><td>1.5</td><td>245</td><td>2.1</td><td>0.08</td><td>1.5</td></tr><tr><td>7</td><td>3.5</td><td>10</td><td>1.0</td><td>6.2e14</td><td>0</td><td>0</td><td>2</td><td>1.0</td><td>262</td><td>9.2</td><td>0.10</td><td>17.4</td></tr><tr><td>8</td><td>3.5</td><td>15</td><td>1.0</td><td>6.2e14</td><td>0</td><td>0</td><td>2</td><td>1.6</td><td>249</td><td>19.5</td><td>0.59</td><td>14.3</td></tr><tr><td>9</td><td>4</td><td>0</td><td>1.0</td><td>6.2e14</td><td>0</td><td>0</td><td>1</td><td>3.8</td><td>446</td><td>21.7</td><td>0.42</td><td>7.6</td></tr><tr><td>10</td><td>4</td><td>5</td><td>1.0</td><td>6.2e14</td><td>0</td><td>0</td><td>1</td><td>1.7</td><td>364</td><td>1.1</td><td>0.03</td><td>0.9</td></tr><tr><td>11</td><td>4</td><td>10</td><td>1.0</td><td>6.2e14</td><td>0</td><td>0</td><td>2</td><td>1.1</td><td>382</td><td>22.1</td><td>0.15</td><td>49.6</td></tr><tr><td>12</td><td>4</td><td>15</td><td>1.0</td><td>6.2e14</td><td>0</td><td>0</td><td>2</td><td>1.2</td><td>369</td><td>25.2</td><td>0.45</td><td>19.7</td></tr><tr><td>13</td><td>4</td><td>10</td><td>0.5</td><td>1.6e14</td><td>0</td><td>0</td><td>3</td><td>1.1</td><td>375</td><td>37.5</td><td>0.29</td><td><strong>59.1</strong></td></tr><tr><td>14</td><td>4</td><td>10</td><td>2.0</td><td>2.5e15</td><td>0</td><td>0</td><td>2</td><td>1.2</td><td>384</td><td>14.8</td><td>0.31</td><td>11.7</td></tr><tr><td>15</td><td>4</td><td>10</td><td>1.0</td><td>6.2e14</td><td>10</td><td>0</td><td>2</td><td>1.4</td><td>496</td><td>33.2</td><td>0.53</td><td>33.6</td></tr><tr><td>16</td><td>4</td><td>10</td><td>1.0</td><td>6.2e14</td><td>20</td><td>0</td><td>3</td><td>1.8</td><td>653</td><td>55.0</td><td>0.59</td><td><strong>58.2</strong></td></tr><tr><td>17</td><td>4</td><td>0</td><td>1.0</td><td>6.2e14</td><td>0</td><td>5e-8</td><td>1</td><td>3.8</td><td>446</td><td>21.7</td><td>0.42</td><td>7.6</td></tr><tr><td>18</td><td>4</td><td>10</td><td>1.0</td><td>6.2e14</td><td>0</td><td>5e-8</td><td>3</td><td>1.1</td><td>382</td><td>37.6</td><td>0.15</td><td><strong>59.0</strong></td></tr></tbody></table>

Note. Columns show the: simulation number (#), impactor size ($d_i$ in km), clathrate layer thickness ($h_c$ in km), grain size ($d$ in mm), corresponding viscosity at the interface with the ocean ($\eta_b$ in Pa s), surface porosity ($\xi$ in %), the amplitude of internal tidal heating ($Q_m$ in W/m$^3$), melt behavior regime (1–3, see text), crater depth ($d_c$ in km), initial melt volume ($V_0$ in km$^3$), freezing time of all melt ($t_f$ in kyr), freezing time of half of initial melt volume $V_0$ ($t_{1/2}$ in kyr), and maximum depth reached by melt ($d_m^{max}$ in km)—bold numbers indicate the cases where melt reaches the underlying ocean.
We prescribe free slip on all boundaries, zero heat flux at the left and right boundaries, fixed temperature $T_b$ at the interface with the ocean, and a radiation condition at the surface:

$$q = -k \nabla T \cdot n = \epsilon \sigma (T^4 - T_s^4), \qquad (12)$$

with $q$ the normal conductive heat flux, $n$ the outer unit normal, $\epsilon$ Titan's surface emissivity, $\sigma$ the Stefan-Boltzman constant, and $T_s$ Titan's surface temperature. All parameter values are given in Table 1.

### 2.2.4. Numerical Method

We solve the equations described in Section 2.2.1 complemented by the boundary conditions (see above) in a 2d axisymmetric geometry with the left boundary being the axis of symmetry where the impactor has hit Titan's surface. The computational domain is 120 km wide and 60 km high. We use the finite element method implemented in the open source FEniCS software (Alnæs et al., 2015; Logg et al., 2012). The computational mesh consists of triangular elements and is significantly refined at the top and bottom boundaries and at the symmetry axis (left boundary). The code has been carefully tested in a series of numerical experiments for two-phase flow (Schmeling et al., 2016) as well as compared with our previous two-phase flow studies (Kalousová & Sotin, 2018, 2020a; Kalousová et al., 2016, 2018).

## 3. Results

We have performed a total of 18 numerical simulations, varying the impact diameter ($d_i = 3, 3.5, 4$ km), the clathrate layer thickness ($h_c = 0, 5, 10, 15$ km), the grain size ($d = 0.5, 1.0, 2.0$ mm), the surface porosity ($\xi = 0, 10\%, 20\%$), and the tidal heating amplitude ($Q_m = 0$ and $5 \times 10^{-8}$ W/m$^3$). All simulations and their results are summarized in Table 3. Note that for the simulations with tidal heating (#17–18), we did not perform additional impact simulations but started with the $T, \phi, C$ fields from the corresponding simulations without tidal heating (#9 and 11, respectively)—therefore, the initial volume $V_0$, melt distribution and crater depth $d_c$ are the same for the two pairs of simulations (#9 and 17, #11 and 18).

In Figure 1, we plot a few snapshots from the time evolution of temperature (left panels) and water content (right panels) for simulations with the largest impact diameter $d_i = 4$ km and four values of clathrate layer thickness $h_c$ (0, 5, 10, 15 km, simulations #9–12). In the case without clathrates ($h_c = 0$ km, left column) (a) the top half of the ice shell is much colder (see Kalousová & Sotin, 2020b, for a detailed discussion) and (b) the crater is much deeper than in the cases with some clathrate layer (see also Table 3 and Figure 4a)—this is the result of combination of the thermal profile difference (cases with clathrates are warmer) and clathrate rheology (clathrates are stronger than water ice, Durham et al., 2003). While the volume of generated melt is mainly dependent on the impactor diameter $d_i$ (see Table 3 and Figure 4b), the melt distribution depends strongly on the clathrate layer thickness (see also Wakita et al., 2023). In the case without clathrates ($h_c = 0$ km, first column), the impact generated melt pool is about 5 km deep and has a well defined hemispherical shape. For the thinnest non-zero clathrate layer considered in this study ($h_c = 5$ km, second column), the melt pool is broad and very shallow ($\lesssim 1$ km). For $h_c = 10$ km (third column), there are two connected molten regions—a broad shallow pool and a region of very deep (up to ~20 km) melt along the symmetry axis (see Section 4 for a discussion). Finally, for the largest considered clathrate layer thickness ($h_c = 15$ km, fourth column), the melt distribution is similar to the $h_c = 10$ km case, but the deep melt region has a substantially smaller volume and the two melt regions are not connected. The variations in the melt distribution due to the clathrate layer thickness can be explained by the following observation. The impact simulations with a 10 km thick clathrate layer have a higher pre-impact temperature (see Kalousová & Sotin, 2020b) and lower yield strength over much of the target compared to simulations with a 5 or 15 km thick clathrate layer (see Figure 2 in Wakita et al., 2023). The lower strength results in formation of a larger central uplift which, upon collapse, can also penetrate more deeply into the weak underlying ice. The case with a 5 km thick clathrate layer is the strongest and results in only near surface melt while the case of the 15 km thick clathrate layer has an intermediate strength resulting in some near surface melt and some deeper melt originating from collapse of the central uplift. Note that in the case without clathrates, the melt pool is embedded within the stagnant lid and thus underlain by cold ($T \lesssim 180$ K) and highly viscous ($\eta \gtrsim 10^{19}$ Pa s) ice, whereas in all cases with a clathrate layer, the molten region is underlain by warmer (and thus less viscous)
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Time evolution (top to bottom) of temperature (left half of each plot) and water content (right half of each plot) in the middle of the domain (±20 km from the symmetry axis) for $d_i = 4$ km and increasing values of clathrate thickness ($h_c = 0, 5, 10, 15$ km, from left to right). Black contours in the temperature fields correspond to 150, 200, and 250 K. The partially molten domain (i.e., where $T = T_m$) is marked by the dark red color in the temperature fields and highlighted by green contours in both fields. The orange contours in the water content fields border the clathrate layer. Note that Titan's surface is at $z = 60$ km and the ocean-ice interface is at $z = 0$ km.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The same as in Figure 1 but for $h_c = 10$ km and an increasing value of grain size ($d = 0.5, 1.0, 2.0$ mm, from left to right).</figcaption></figure>

ice. Also note that while the thinner clathrate layers ($h_c = 5$ and $10$ km) were breached by the impact, the thickest clathrate layer ($h_c = 15$ km) stayed mostly intact (see the orange contours in Figure 1).

The melt distribution after the impact governs its evolution. This is shown in Figure 1: For the case without clathrates (first column), the melt pool remains stationary within the stagnant lid and slowly freezes (freezing time
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Results for $d_i = 4$ km (simulations #9–18). Top: Time evolution of melt volume $V$ in the ice shell. The full circles show the initial volume $V_0$, and the empty circles indicate the time when all melt is frozen (#9–12, 14, 15, 17). The stars show the volume at the time of melt arrival at the ocean interface (#13, 16, 18). The inset shows the time <0.65 kyr in more detail. The empty squares mark the time when half of the initial amount of melt has frozen ($t_i$, see also Table 3). Bottom: Time evolution of the maximum melt depth $a_m^{\text{max}}$. The empty circles show the depth at the time of freezing (#9–12, 14, 15, 17). The stars show the time of melt arrival at the ocean interface (#13, 16, 18). See Figure 2 for a visual representation of the melt descending through Titan's ice shell into the ocean (simulation #13).</figcaption></figure>

$t_f = 21.7$ kyr, see Table 3). Even though the surrounding ice is warmed up by the freezing water, the largest melt depth does not vary significantly (see also the bottom panel of Figure 3). For $h_c = 5$ km (second column), the freezing time is significantly shorter (~1 kyr) because the melt pool is shallow (~1 km). For $h_c = 10$ km (third column), the melt evolution is much more complex—the shallow melt freezes quickly, but the column of deep melt around the symmetry axis starts to descend as it is located in very warm ice ($T > 250$ K) with low viscosity. The descent lasts several thousands of years and during that time, freezing competes with the melt transport. Eventually, all melt freezes in the warm ice layer before reaching the underlying ocean ($t_f = 22.1$ kyr). The deepest the melt reaches is about 10 km above the ocean-ice interface. Finally, for $h_c = 15$ km (fourth column), the deep melt region freezes within ~2 kyr while the shallow melt pool freezes in about 25 kyr, that is, slightly slower than in the no-clathrate case even though the initial melt volume is smaller—this is because the melt is fully embedded within the insulating clathrate layer (see the orange contours in Figure 1).

Next, we investigated the effect of ice viscosity by varying the ice grain size $d$ (larger ice grain results in larger viscosity and vice versa) while keeping the clathrate layer constant at 10 km thick. We find that the melt distribution after the impact does not depend as strongly on grain size as it does on the clathrate layer thickness. Figure 2 shows a few snapshots from the time evolution of temperature (left panels) and water content (right panels) for simulations with $d_i = 4$ km, $h_c = 10$ km, and three values of grain size $d$ (0.5, 1.0, 2.0 mm, from left to right, simulations #13, 11, 14). Note that the clathrate layer is breached in all three cases during impact (see the orange contours in Figure 2). Initially, the three molten regions are all composed of a deep melt region along the symmetry axis as well as a shallow melt pool—this is the effect of the clathrate layer thickness, $h_c = 10$ km. However, the depth reached by the deep melt region just after impact ($t = 0$ kyr) differs for different grain sizes, between >20 km ($d = 0.5$ mm, smallest viscosity, left), ~20 km ($d = 1.0$ mm, reference viscosity, middle), and ~10 km ($d = 2.0$ mm, largest viscosity, right). The effective viscosity of ice also plays a key role in the evolution of water content—the smaller the ice viscosity (grain size), the deeper the melt can penetrate. In the case of the smallest grain size considered here, $d = 0.5$ mm, the viscosity is smallest (at the bottom boundary, $\eta_b \sim 1.6 \times 10^{14}$ Pa s, see Table 3) and thus the flow of the ice and melt mixture is fastest—the downwelling reaches the bottom boundary at $t \sim 1$ kyr. For the reference grain size, $d = 1.0$ mm ($\eta_b \sim 6.2 \times 10^{14}$ Pa s), melt descends into the warm ice, but the progress is not fast enough and all melt freezes before

reaching the ocean-ice interface. Finally, for the largest grain size considered in this study, $d = 2.0$ mm ($\eta_b \sim 2.5 \times 10^{15}$ Pa s), the ice viscosity is too large for any significant movement to occur on the timescales of a few kyr. During this time, the melt is located in the coldest parts of the ice shell and thus freezes—for this particular simulation, the freezing time is $t_f = 14.8$ kyr.

We also studied the effects of surface porosity (simulations #15–16) and internal heating from tides (simulations #17–18). As shown in previous works (Artemieva & Lunine, 2005; Wünnemann et al., 2008), surface porosity strongly affects the amount of generated melt. The top panel of Figure 3 shows the time evolution of total melt volume $V$ for all simulations with $d_i = 4$ km (#9–18) and demonstrates the profound effect of surface porosity $\xi$ on the initial melt volume $V_0$: for $\xi = 0\%$ (simulation #11, the blue circle hidden below the cluster of circles just below 400 km$^3$), the initial melt volume is $382$ km$^3$; for $\xi = 10\%$ (simulation #15, dashed light green), $V_0 = 496$ km$^3$; and for $\xi = 20\%$ (simulation #16, dot-dashed dark green), $V_0 = 653$ km$^3$ (see also Table 3). In the inset, we also illustrate that the time to freeze half of the initial melt volume is <1 kyr for all simulations, indicating that a significant amount of melt is distributed close to the surface where freezing is fastest. Note that in
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Crater depth $d_c$ (panel a), initial amount of melt $V_0$ (panel b), freezing time $t_f$ (panel c), and maximum melt depth $d_m^{max}$ (panel d) as a function of clathrate layer thickness $h_c$ for all simulations (colors). The symbols correspond to the melt behavior regime: square (1) melt present only within the stagnant lid and freezes in place, diamond (2) deep melt present may move downward, but freezes before reaching the ocean, and star (3) melt transported to the ocean. Note that we only assumed $h_c = 0, 5, 10,$ and $15$ km but in the figure, we offset some symbols laterally around the given clathrate layer thickness for better visibility. We indicate the melt behavior regime in panel (d).</figcaption></figure>

some cases with deep initial melt, the freezing times (empty circles in the top panel of Figure 3) are out of range of the figure (>35 kyr).

The bottom panel of Figure 3 shows that the surface-to-ocean melt exchange occurred for 3 simulations: #13 (smallest viscosity, $d = 0.5$ mm, dashed cyan), #16 (largest surface porosity, $\xi = 20\%$, dot-dashed dark green), and #18 (additional tidal heating, dashed orange). In two cases—$h_c = 10$ km with $\xi = 0\%$ (#11, blue) and $h_c = 10$ km with $\xi = 10\%$ (#15, dashed light green), the melt descends into the warm ice but freezes during the descent. Note that the quick decrease in the maximum melt depth (at $t \sim 21$ and 28 kyr, respectively) corresponds to the time when the deep melt has frozen, but shallow melt is still present. While one could expect the shallow melt to freeze faster, this seeming swap occurs because the deep melt is embedded within warm ice with thermal conductivity of ~2.3 W/m/K while the shallow melt is surrounded by clathrates (see the orange contour in Figure 1) that have a thermal conductivity that is more than 4 times lower. The same thing happens in the case of the thickest clathrate layer ($h_c = 15$ km, #12, full light green) at $t \sim 2$ kyr. Also note that the gradual decrease in the maximum melt depth after ~10 kyr (#11, blue) and ~15 kyr (#15, dashed light green) is not due to melt upwelling but rather due to the fact that the freezing proceeds from the bottom up (where the water amount is smallest). Although the internal heating has no effect in the no-clathrate case ($h_c = 0$ km, #9 and 17), it prolongs the freezing time and allows the melt to reach the ocean in case of $h_c = 10$ km (#11 and 18). This indicates that the effect of internal heating is only important if the melt pools are located within the warm ice.

Finally, in Figure 4, we present the summary of all our results in a comparative way. Panel a shows that the crater depth $d_c$ primarily depends on the presence of a clathrate layer, while the particular clathrate layer thickness $h_c$, the impactor size $d_i$ (within the narrow range of 3–4 km), and the surface porosity play secondary roles. The crater depth is lowered from ~3–4 km for pure ice to ~1–2 km in the presence of a clathrate layer. The depth of Selk crater is ~0.5 km (Hedgepeth et al., 2020); such a shallow depth suggests the presence of a clathrate layer
combined with potential infilling by sand after the impact (Wakita et al., 2023). In panel b, we plot the initial melt volume $V_0$. This depends mainly on the impactor diameter $d_i$ and the surface porosity $\xi$. The presence or absence of the clathrate layer only slightly modulates the total melt volume. Panel c summarizes the freezing time $t_f$, which depends primarily on the clathrate layer thickness $h_c$ due to its role in shaping the melt distribution (deep/shallow molten regions). Melt in cases with $h_c = 5$ km freezes fastest because it forms a very shallow and broad melt pool. Freezing times of melt in other cases ($h_c = 0, 10, 15$ km) are similar, although various mechanisms affect the freezing process—strong melt localization in a well defined reservoir ~5 km deep and ~8 km wide for $h_c = 0$ km, presence of deep melt for $h_c = 10$ km, and the insulating effect of an intact clathrate layer for $h_c = 15$ km. In three cases—small grain size (#13, cyan star), high surface porosity (#16, dark green star), and internal tidal heating (#18, orange star)—the freezing times are longest because melt is located in warm parts of the ice shell. Finally, panel d depicts the maximum depth reached by the melt: we find that for the cases without clathrates ($h_c = 0$ km) and with a thin clathrate layer ($h_c = 5$ km), melt stays within the top 10 km of the ice shell. This is also the case for the smallest impactor ($d_i = 3$ km) and thickest clathrate layer ($h_c = 15$ km). In most of the cases with $h_c = 10$ and 15 km, some deep melt is present—either only in the initial setting or, in some cases, the initial deep melt moves further down into the warmer parts of the ice shell. However, only the simulations with a rather particular setting show an exchange of melt between the surface and the ocean. This behavior occurs in cases with initial deep melt due to a moderately thick clathrate layer ($h_c = 10$ km) and a second factor that further promotes the downwelling: either small viscosity due to small grain size (#13, cyan), large melt volume due to large surface porosity (#16, dark green), or internal tidal heating (#18, orange). Based on the results displayed in panel d, we describe 3 melt behavior regimes (see also Table 3):

1. ☐ Melt is initially present only within the stagnant lid and freezes in place without any significant movement (regime 1, simulations #1–2, 4–6, 9–10, 17). This is the case for thin ($h_c = 5$ km) or absent clathrate layers, as well as for the smallest impactor ($d_i = 3$ km) hitting the thickest clathrate layer ($h_c = 15$ km).

2. ◊ Deep melt is initially present (water below the stagnant lid) and may move downwards into the warmer domain, but eventually freezes either in the stagnant lid or on the way to the ocean (regime 2, simulations #3, 7–8, 11–12, 14–15). This is the case for thicker clathrate layers ($h_c = 10–15$ km) with little ($\xi = 10\%$) to no surface porosity, intermediate to large grain size ($d = 1–2$ mm), and no internal heating ($Q_m = 0$ W/m$^3$).

3. ☆ Deep melt is initially present and is transported all the way into the ocean (surface-to-ocean exchange, regime 3, simulations #13, 16, 18). This behavior only occurs in special cases of small viscosity ($d = 0.5$ mm), large surface porosity ($\xi = 20\%$), or significant internal heating ($Q_m = 5 \times 10^{-8}$ W/m$^3$) with a 10-km thick clathrate layer.

## 4. Discussion

In the present study, we investigated the fate of impact melt pools on Titan using numerical models. We found that surface-to-ocean exchange can only be achieved in a rather particular setting (the presence of deep melt coupled with low viscosity, high surface porosity, or internal tidal heating). This seems to be in contrast to a recent paper by Carnahan et al. (2022) who studied the fate of impact melt pools on Europa and found that, in the case of craters with transient cavity depths >50% of the ice shell thickness, the melt is drained very efficiently into the ocean on timescales of a few kyr. In both cases, surface-to-ocean drainage may help to deliver materials needed for habitability (such as oxidants and organic molecules). While the two studies bear many similarities (investigating the fate of impact melt pools, assuming cylindrical geometry and impermeable limit, and using the enthalpy method), there are also a few differences that we would like to discuss here. First, the ice shells of Titan and Europa differ from one another in at least two characteristics: (a) Europa's ice shell is likely thinner than Titan's—therefore Carnahan et al. (2022) investigate thicknesses between 10 and 40 km but here we assume that Titan's ice shell is 60 km thick (see the paragraph below on ice shell thickness). Note that based on geophysical measurements, the thickness is inferred to be at least 40 km (Baland et al., 2014) and can be as large as 100 km (Mitri et al., 2014). This means that while we are considering impacts with comparable kinetic energy (20–1,334 EJ in Carnahan et al. (2022) vs. 725–1,718 EJ here), in our study, no impact generates a transient crater that would be >30 km deep (Wakita et al., 2023, indicated that for the largest impactor with $d_i = 4$ km, the transient crater depth is ~20 km). Thus, none of our craters satisfy Carnahan et al.'s condition for melt delivery into the ocean. (b) Europa's crust is mostly made of water ice while Titan's surface may be covered by a few to tens of kilometers thick layer of thermally insulating methane clathrates. Consequently, Carnahan et al. (2022) prescribe constant thermal conductivity, but we take into account the thermal dependence of ice conductivity and combine it with the
presence of the clathrates (see Equation 8). Second, we include thermal buoyancy in our model (last term on the left-hand side of Equation 1b), which is neglected in Carnahan et al. (2022). This is an important process to take into account since the positive thermal buoyancy (warm material tends to ascend) acts to reduce the negative melt buoyancy (water tends to descend). For example, when molten material (T ~ 271 K) descends through ice at ~257 K, the third and fourth terms on the left-hand side of Equation 1b become comparable for φ ~ 2.5%. Therefore, the thermal buoyancy effect is most important when the water content is on the order of a few percent. However, even for larger water contents, the downwelling becomes less efficient if thermal buoyancy acts against it. Third, we employ the stress-dependent composite creep that combines the deformation by five mechanisms (see Equation 5), but Carnahan et al. (2022) prescribe temperature-dependent diffusion creep. Fourth, we implement the surface topography generated by the impact. Although we do not evolve this topography with time (the crater relaxation times being several orders of magnitude longer than the timescales investigated here), our approach does not add an artificial insulating layer of ice by filling the crater as implemented by Carnahan et al. (2022). Note, that 1 km of ice would delay the freezing by a few tens of thousands of years. Finally, we prescribe a radiation boundary condition (Equation 12) that allows for the realistic evolution of near-surface temperatures, which affect the melt evolution. These differences lead to different conditions for surface-to-ocean exchange between this study focused on Titan and the Carnahan et al. (2022) study focused on Europa.

In the present study, we assume Titan's ice shell to be 60 km thick. This value was chosen for the impact modeling by Wakita et al. (2023), because the ice shell thicker than 40 km has no effect on the crater formed by $d_i \le 5$ km impactors (Wakita, Johnson, et al., 2022). At the same time, thicker shells would need larger meshes, and thus the chosen value represents a sweet spot between the computational demands and the physical setting of the impact model. Since the data from Wakita et al. (2023) are used as an input for the two-phase simulations of impact melt evolution, we keep the same thickness. If Titan's ice shell is thinner, the possibility that the near-surface melt from the impact pool reaches the subsurface ocean should be higher than reported here as the travel time would be shorter (assuming the velocity would remain about the same) and thus the competition between the travel time and freezing time would become more favorable for the surface-to-ocean exchange. On the other hand, in the case of a thicker shell, it would be less likely for any impact melt to reach the underlying ocean as the travel time would become even longer while the freezing time would remain comparable with the presented results. Note that the chosen thickness value is in line with estimates of Titan's ice shell thickness from Cassini (e.g., Baland et al., 2014; Mitri et al., 2014; Sotin et al., 2021).

During the impact, the destabilization of clathrates leads to the formation of water (melt) and the release of methane gas. Methane clathrates that would not be destabilized during the impact would remain stable (solid) in the presence of water (melt). The solid clathrates have density very similar to that of ice and thus they would behave as ice, in a sense of buoyancy. Therefore, depending on their post-impact location, these solid clathrates would either stay in the stagnant lid or collapse downwards by the Rayleigh-Taylor instability (together with ice and the molten water). However, this is not treated in the numerical setting of the present paper where the clathrates distribution (given by the impact code) is prescribed constant in time (see orange lines in Figures 1 and 2).

The use of axisymmetric models for the vertical impacts might exaggerate the post-impact depth of the impact melt. When the largest impactor ($d_i = 4$ km) strikes the thick and warm methane-clathrate layer (≥10 km), a large central uplift is produced during the crater collapse (Wakita et al., 2023). The collapse of this uplift advects melt to great depth close to the symmetry axis. A similar distribution of melted material along the symmetry axis can be found in other two-dimensional impact simulations (e.g., Crosta et al., 2021; Potter et al., 2012). Although the mass of material near the symmetry axis is small, more realistic simulation of moderately oblique impacts or fully 3D simulations of vertical impacts may not result in melt reaching the same depths as the simulations presented here. Indeed, oblique impacts result in small central uplifts (e.g., Collins et al., 2020). The effect of impact angle on the impact heating and the volume of melt was investigated in several papers (e.g., Artemieva & Lunine, 2003; Pierazzo & Melosh, 2000; Wakita et al., 2019; Wakita, Genda, et al., 2022). Pierazzo and Melosh (2000) found that the volume of impact melt normalized by volume of the transient crater is nearly independent of impact angle for moderately oblique impacts. Wakita et al. (2019) and Wakita, Genda, et al. (2022) however found that by including the effects of shear heating, the volume of impact heated material is relatively insensitive to impact angle for moderately oblique impacts if impactor size is held constant. Note that for the same impact speed and final crater size, an oblique impact requires a larger impactor than the corresponding vertical impact. Thus, if shear heating dominates melt production, a crater produced by an oblique impact will result in more melt than a
similarly sized crater produced by a vertical impact. Although oblique impact may result in more impact melt for a similarly sized crater, our axisymmetric models may exaggerate post-impact melt depths. Thus, we expect our simulation result in favorable conditions for melt migration compared to oblique impacts. Future work using three dimensional codes can assess the distribution of deep melt pools in vertical and oblique impacts.

The description of the convecting two-phase ice-water system in this study is considered in a simplified impermeable limit. This means that the melt is trapped within the ice matrix and no differential motion between ice and water is considered. Consequently, the mobilization of the partially molten regions results solely from the formation of Rayleigh-Taylor (R.-T.) instabilities as a result of the negative buoyancy of the melt. While this approach can be justified for small porosities, for which the effective permeability of the matrix is low or negligible, it is only a crude approximation for the partially molten regions in this study, where porosities exceed tens of percent. Mobilization of the fluid phase would enhance phase separation and would lead to water accumulation at the bottom of the permeable zone. This would increase the density contrast and thus decrease the time scale for the formation of R.-T. instabilities (Turcotte & Schubert, 2014, see also Appendix B). Simultaneously, enhanced phase mixing and liquid convection may enhance heat extraction from the partially molten region, in analogy with convection in mushy regions (see Boury et al., 2021). This would probably shorten the time scale for freezing the partially molten regions. To account for these competing effects, simultaneous deformation of the ice matrix and the differential motion of the fluid would have to be modeled as in Kalousová et al. (2018) and Kalousová and Sotin (2018), Kalousová and Sotin (2020a). The additional uncertainties connected with parameterizations of such a system, for example, the form of permeability-porosity law, and substantial additional computational costs led us, in this study, to the adopted simplifications.

This work puts an upper limit on the amount of organics that can be transported from Titan's surface to its interior water ocean, which has implications for its habitability. To maintain a habitable environment, significant amounts of organic material must be delivered to Titan's ocean. Neish et al. (2024) assumed that melt transport through Titan's ice shell mirrored that of Europa, and found that at most 10¹⁴ mol/Gyr of glycine could be delivered from Titan's surface to its interior. This gives an upper limit for biomass productivity of ~10³ kgC/yr for a glycine fermentation metabolism, which is several orders of magnitude less than the hypothetical biomass production supported by Enceladus's subsurface ocean (Affholder et al., 2022). This new work suggests that significantly fewer impact craters on Titan can transport organics to Titan's interior than assumed in Neish et al. (2024). Thus, Titan's subsurface ocean is even less likely to be habitable, unless useful organic molecules can be sourced from its interior (Miller et al., 2023).

However, this work does suggest that it is more likely that Dragonfly will explore an impact whose subsurface melt pool froze in place and did not migrate downward. One of the primary goals of the Dragonfly mission is to study the prebiotic chemistry of impact melt deposits, where liquid water has mixed with surface organics to produce molecules of biological interest (Barnes et al., 2021). In the parameter space that we have investigated, a large number of the scenarios resulted in a melt deposit at the Selk crater that remained close to the surface. Evolved organic molecules including molecules of biological interest, would have been produced by the reaction of the surface organics with water (Neish et al., 2018). These molecules would have been frozen very rapidly on geological timescales. If rivers or other natural canyons have cut into these melt deposits, they could be sampled by Dragonfly.

The migration of melt formed during impacts would also be an important process for understanding the dynamics of the icy shells of the large Galilean moons, Ganymede and Callisto. These moons will be studied by the ESA-led JUICE mission (Grasset et al., 2013). Ganymede, the largest moon in the Solar System, has a radius 60 km larger than that of Titan and it also harbors a deep subsurface ocean (Kivelson et al., 2002). In contrast to Titan, it does not have a dense atmosphere and its surface is heavily cratered. The surface shows two types of terrains: dark older terrains and bright younger terrains (which are still >3 Gyr old). On these terrains, a large number of impact craters are present. The migration of melt to the ocean would modify the composition of the ocean by bringing non-ice compounds from both the surface and the impactor. The present work is not directly applicable to Ganymede as we do not expect a clathrate layer to be present and a dedicated study should be performed. The Ganymede laser altimeter (GALA) will precisely determine Ganymede's topography (Steinbrügge et al., 2015). Together with the gravity field, these observations will allow for the detection of subsurface mass anomalies like the frozen melt pools where non-ice material with densities different than that of water ice may be present. In addition, the ice penetrating radar (RIME) will characterize compositional, thermal, and/or structural variations in
the outer icy shell. These data will provide information on the composition beneath impact craters (Bruzzone et al., 2015) and will allow for the search of present or past reservoirs of liquid water by detecting strong reflections from large bodies as well as scatter from smaller reservoirs. This information will shed light on the stability of impact-generated subsurface melt pools within icy shells, a process that is of interest for all icy ocean worlds.

## 5. Conclusions

We investigated the evolution of impact melt pools on Titan using a numerical model that solves the problem of two-phase thermal convection of solid ice and liquid water mixture. We included the insulating effect of a few to tens of kilometers thick clathrate layer and confirmed that it significantly affects the distribution of generated melt, which then governs the melt pool evolution. Our results show that in the case of a thin (5 km) or absent clathrate layer, the melt pool freezes within the stagnant lid in a short time (≤12 kyr for $d_i = 3$ km, ≤16 kyr for $d_i = 3.5$ km, and ≤22 kyr for $d_i = 4$ km, respectively). For thicker clathrate layers ($h_c = 10–15$ km), deep melt is generated, which can be transported through the ice shell. We note that the melt depths obtained here should be considered as upper limits since we assumed only vertical impacts. In most of the investigated cases, the melt freezes without reaching the ocean interface—either in the stagnant lid or in the warmer ice below. Surface-to-ocean exchange only occurs in three special cases: $h_c = 10$ km and $d_i = 4$ km, combined with either (a) small viscosity, (b) large surface porosity, or (c) tidal heating. Note that in these cases, the freezing times are also longest (≥37 kyr). We used the initial melt pools that potentially lay beneath the Selk-sized crater, which is the 8th largest crater on Titan. This suggests that only a few impact craters on Titan could lead to the transport of surface organics to Titan's deep ocean, which strongly limits its habitability potential (unless the organic molecules are released from Titan's silicate core). On the other hand, our results suggest that NASA's Dragonfly mission may be able to sample the melt deposits at the Selk crater, where organics might have reacted with impact induced liquid water to produce complex organics or even biomolecules.

## Appendix A: Enthalpy Method

Here, we derive the mixture enthalpy balance (Appendix A1), the bulk enthalpy of the mixture (Appendix A2), and the melting-freezing relations to evaluate temperature $T$ and water content $\phi$ from the value of enthalpy $h$ (Appendix A3).

### A1. Enthalpy Balance

We start from the energy balance for the whole mixture which reads (e.g., Rudge et al., 2011; Šrámek et al., 2007)

$$ 
\begin{aligned}
\frac{\partial}{\partial t}(\rho_w \phi \varepsilon_w + \rho_i (1-\phi)\varepsilon_i) + \nabla \cdot (\rho_w \phi \varepsilon_w \mathbf{v}_w + \rho_i (1-\phi)\varepsilon_i \mathbf{v}_i) \\
= -\nabla \cdot \mathbf{q} + Q + \nabla \cdot (\phi \mathbf{v}_w \cdot \boldsymbol{\tau}_w + (1-\phi)\mathbf{v}_i \cdot \boldsymbol{\tau}_i) + \rho_w \phi \mathbf{v}_w \cdot \mathbf{g} + \rho_i (1-\phi)\mathbf{v}_i \cdot \mathbf{g},
\end{aligned}
\quad (A1) $$

where $\varepsilon_w$, $\varepsilon_i$ is the internal energy of each phase, respectively, $\mathbf{v}_w$, $\mathbf{v}_i$ is the velocity of each phase, respectively, $\mathbf{q}$ is the conductive heat flux, $Q$ is the relevant heat sources (e.g., tides), and $\boldsymbol{\tau}_w$, $\boldsymbol{\tau}_i$ is the Cauchy stress tensor of each phase, respectively. Note that we neglect the interfacial processes (see Šrámek et al., 2007, for more details) and also the kinetic energies of both phases.

The partial enthalpies of the fluid and the matrix are defined through

$$ h_w = \varepsilon_w + \frac{p}{\rho_w}, \quad h_i = \varepsilon_i + \frac{p}{\rho_i} \quad (A2) $$

where $p$ denotes the pressure, which is equal for both phases (due to the negligence of the effects of interfacial tension) and also equal to the interfacial pressure (Rudge et al., 2011; Šrámek et al., 2007).

The bulk mixture enthalpy can be obtained as (Carnahan et al., 2022; Katz, 2008):
$$ \rho h = \phi \rho_w h_w + (1 - \phi) \rho_i h_i, \qquad (A3) $$

with the mixture density defined as

$$ \rho = \phi \rho_w + (1 - \phi) \rho_i. \qquad (A4) $$

Assuming further that $v_w = v_i = v$ (impermeable limit, e.g. Kalousová et al., 2016) and noting that $\nabla \cdot v = 0$ (Equation 1a), we can rewrite the energy balance (Equation A1) as

$$ \rho \left( \frac{\partial h}{\partial t} + v \cdot \nabla h \right) + h \left( \frac{\partial \rho}{\partial t} + v \cdot \nabla \rho \right) - \left( \frac{\partial p}{\partial t} + v \cdot \nabla p \right) = - \nabla \cdot \mathbf{q} + Q + \tau : \nabla \mathbf{v} + v \cdot (\nabla \cdot \tau + \rho g), \qquad (A5) $$

with $\tau = \phi \tau_w + (1 - \phi) \tau_i$ the mixture stress tensor. In the following, we neglect viscous dissipation (third term on the right-hand side). Moreover, the second term on the left-hand side and the last term on the right-hand side are equal to zero due to the conservation of mixture mass and momentum, respectively (cf. e.g. Rudge et al., 2011; Šrámek et al., 2007). Furthermore, we approximate $p$ by the hydrostatic pressure in the pure ice phase (satisfying $\nabla p = \rho_w g$) and neglect its time derivative.

For numerical reasons, it is useful to rewrite the above-derived equations in terms of gravity reduced enthalpy (e.g., Katz, 2008):

$$ \tilde{h} = h - \frac{p - p_0}{\rho_i}, \qquad (A6) $$

for which it clearly holds that

$$ \frac{\partial \tilde{h}}{\partial t} = \frac{\partial h}{\partial t}, \qquad \nabla \tilde{h} = \nabla h - \frac{1}{\rho_i} \nabla p = \nabla h - \mathbf{g}, \qquad (A7) $$

where we have already neglected the time derivative of hydrostatic pressure. Using Equations A6 and A7, we can finally rewrite Equation A5 as

$$ \rho \frac{\partial h}{\partial t} + \rho v \cdot \nabla h - \phi \Delta \rho v \cdot g = \nabla \cdot (k \nabla T) + Q, \qquad (A8) $$

where we used Fourier's law for the heat flux $q$, $\Delta\rho = \rho_i - \rho_w$ for the density difference between the ice and water, and we omitted the tildes for the gravity reduced enthalpy.

## A2. Enthalpy of the Mixture

Assuming constant material parameters ($\alpha_k$, $c_k^P$, $\rho_k$ with $k = i$ or $w$), we can express the partial enthalpies of ice and water as functions of temperature $T$ and pressure $p$:

$$ h_i(T, p) = h_i^0 + c_i^P (T - T_0) + \frac{1 - \alpha_i T}{\rho_i} (p - p_0), \qquad (A9a) $$

$$ h_w(T, p) = h_w^0 + c_w^P (T - T_0) + \frac{1 - \alpha_w T}{\rho_w} (p - p_0), \qquad (A9b) $$

with $h_k^0$ the value at reference temperature $T_0$ and pressure $p_0$. Equation A3 can be simplified as (noting that $L = h_w - h_i$)
$$
\begin{align*}
h &= \frac{1}{\rho}(\phi\rho_w h_w + (1-\phi)\rho_i h_i) \\
&= \frac{1}{\rho}(\phi\rho h_w + (1-\phi)\rho h_i + \phi(\rho_w - \rho)h_w + (1-\phi)(\rho_i - \rho)h_i) \tag{A10} \\
&= \phi h_w + (1-\phi)h_i + \phi(1-\phi)\frac{\rho_w - \rho_i}{\rho}L \\
&\doteq \phi h_w + (1-\phi)h_i,
\end{align*}
$$

where we neglect the last term since $\phi(1-\phi)^{\frac{\rho_w - \rho_i}{\rho}}L \le 0.02L$. Subtracting the contribution from the hydrostatic pressure (gravity reduced enthalpy, see Katz, 2008, for more details), and using the same notation for the new quantity from now on we can write

$$h(T, p, \phi) = \phi h_w + (1 - \phi)h_i - \frac{p - p_0}{\rho_i}. \qquad (A11)$$

Inserting Equations A9a and A9b into Equation A11 and setting $h_i^0 = 0$, we obtain:

$$h(T, p, \phi) = \phi L_0 + (\phi c_w^P + (1-\phi)c_i^P)(T-T_0) + \left[ \phi \frac{1-\alpha_w T}{\rho_w} + (1-\phi) \frac{1-\alpha_i T}{\rho_i} - \frac{1}{\rho_i} \right] (p-p_0), \quad (A12)$$

where we define $L_0 = h_w^0 - h_i^0$.

### A3. Melting Freezing Relations

To evaluate the temperature and water content from the enthalpy, we first define the maximum (gravity reduced) enthalpy of solid ice as:

$$h_i^{\max}(p) = \max_{T=T_m, \phi=0} h(T, p, \phi) = c_i^P(T_m - T_0) - \frac{\alpha_i T_m}{\rho_i}(p - p_0). \qquad (A13)$$

Depending on the actual value of enthalpy $h$, the mixture can either be completely solid ($h < h_i^{\max}$) or partially molten ($h \ge h_i^{\max}$). The particular equations for temperature and water content can be obtained as follows:

If $h < h_i^{\max}$, then the mixture is completely frozen and there is no liquid water present. Using Equations A12 and A13 and assuming that $\phi = 0$, we can write:

$$
\begin{align}
h - h_i^{\max} &= c_i^P(T - T_0) - \frac{\alpha_i T}{\rho_i}(p - p_0) - c_i^P(T_m - T_0) + \frac{\alpha_i T_m}{\rho_i}(p - p_0) \tag{A14} \\
&= \left[ c_i^P - \frac{\alpha_i}{\rho_i}(p - p_0) \right] (T - T_m), \nonumber
\end{align}
$$

from which we obtain

$$T = T_m + \frac{h - h_i^{\max}}{c_i^P - \frac{\alpha_i}{\rho_i}(p - p_0)}. \qquad (A15)$$

If, on the other hand, $h \ge h_i^{\max}$, then the mixture is partially molten, some water is present and the temperature equals the melting temperature. Using again Equations A12 and A13 and assuming now that $T = T_m$, we can write:
$$
\begin{align}
h - h_i^{\text{max}} &= \phi L_0 + (\phi c_w^P + (1-\phi)c_i^P)(T_m - T_0) + \left[ \phi \frac{1-\alpha_w T_m}{\rho_w} + (1-\phi)\frac{1-\alpha_i T_m}{\rho_i} - \frac{1}{\rho_i} \right] (p-p_0) \nonumber \\
&\quad - c_i^P (T_m - T_0) + \frac{\alpha_i T_m}{\rho_i} (p-p_0) \tag{A16} \\
&= \phi \left[ L_0 + (c_w^P - c_i^P)(T_m - T_0) + \left\{ \frac{1-\alpha_w T_m}{\rho_w} - \frac{1-\alpha_i T_m}{\rho_i} \right\} (p-p_0) \right], \nonumber
\end{align}
$$

from which we obtain

$$
\phi = \frac{h - h_i^{\text{max}}}{L_0 + (c_w^P - c_i^P)(T_m - T_0) + \left[ \frac{1 - \alpha_w T_m}{\rho_w} - \frac{1 - \alpha_i T_m}{\rho_i} \right]} \quad (A17)
$$

Appendix B: Destabilization Time Scales and Comparison With the Terrestrial Setting

The phenomenon of gravitational overturn of the partially molten layer which, in this study, represents the main
mechanism for liquid water drainage, is the classical Rayleigh-Taylor (R.-T.) instability. The characteristic time
scale τ<sub>a</sub> for the formation of an instability with a wavelength λ (in a stratified medium with the denser fluid above
the less dense one) can be estimated analytically using the linear stability analysis (Turcotte & Schubert, 2014,
limit for small wavelengths, p. 456, Equation 6.160)

$$
\tau_a = \frac{4\mu}{(\rho_1 - \rho_2)gb}\left(\frac{2\pi b}{\lambda}\right), \qquad (B1)
$$

where μ is the viscosity, ρ₁ and ρ₂ are the densities of the two layers, g is the gravity acceleration and b is the
characteristic thickness of the layers. Taking parameters relevant to our setting, that is, μ = 10¹⁴ Pa s (ice viscosity
at the melting point), b = 30 km, g = 1.35 m s⁻², and considering, as the most extreme case, a water layer above
the ice (i.e., taking ρ₁ = 1,000 kg/m³ and ρ₂ = 920 kg/m³), and finally, estimating λ = 5 km as the characteristic
size of the destabilized water-saturated region, we obtain the (lower bound estimate of) time scale of the formation
of R.-T. instabilities: τₐ = 147 yr. This is in satisfactory correspondence with the time scale of the destabilization
observed in our simulations—cf. Figure 3 (bottom panel), where the melt descend (in those cases where it does
occur) starts in less than 1 kyr.

It is tempting to relate our modeling approach to icy environments on the Earth with interstitial water contents in
similar thermal gradients (cold top, warm base), such as hydraulic lows in firn aquifers and sea ice. In such
systems, downward viscous flow driven by density difference due to enhanced water content (as described in the
present paper) is not observed. However, we believe that no direct correspondence can be made between the two
systems, since the conditions for terrestrial aquifer bodies within glaciers differ significantly from those of
partially molten regions within planetary ice shells in the following aspects. First, the density profiles in firn
aquifers on Earth often do not reach the density of underlying compacted ice (see for instance Montgomery
et al., 2020), therefore, despite the presence of high-density liquid in the matrix, the overall average density yields
a gravitationally stable configuration. Second, a potential fully molten liquid body, for which density contrast
should favor instability, would be probably drained much faster by hydrofracturing since this mechanism, to our
knowledge, represents the most effective drainage mechanism there (e.g., Krawczynski et al., 2009). This is not
true for planetary ice shells, where the overburden pressure typically maintains the ice and potential cracks in a
compressional regime thus preventing downward fracture propagation. Finally, for sea ice with a thickness of the
order of meters to hundreds of meters (for ice shelves), the estimate of the characteristic time scale for the most
unstable mode (relevant in the setting of a large water body on top of an ice layer) can be expressed as (Turcotte &
Schubert, 2014, Equation 6.162)

$$
\tau_a = \frac{13.04\mu}{(\rho_1 - \rho_2)gb} \qquad (B2)
$$
which yields (for terrestrial gravity, the same viscosity as before, and $b = 10^0–10^2$ m) time scales of hundreds to tens of thousands of years, making this mechanism again irrelevant compared to hydraulic fracturing. We thus believe that our modeling approach is correct and well justified for the present application.

## Data Availability Statement

All data used for producing figures in this study are available in Kalousová et al. (2023). Note that iSALE-2D is distributed on a case-by-case basis to academic users in the impact community. It requires a registration from the iSALE webpage (https://isale-code.github.io/) and usage of iSALE-2D and computational requirements are also shown in there.

## Acknowledgments

We would like to thank two anonymous reviewers for their comments that helped to improve the clarity of this manuscript. This work was supported by the Czech Science Foundation through project No. 22-20388S (KK, OS) and through NASA Cassini Data Analysis Program Grant 80NSSC20K0382 (SW, CDN, JMS, BCJ). This research was triggered by discussions during the meeting in Bern of the International Space Science Institute (ISSI) International Team “The habitability of Titan’s subsurface water ocean.” We gratefully acknowledge the developers of iSALE-2D, including Gareth Collins, Kai Wünnemann, Dirk Elbeshausen, Tom Davison, Boris Ivanov, and Jay Melosh (https://isale-code.github.io). This research was supported in part through computational resources provided by Information Technology at Purdue, West Lafayette, Indiana. All data used for producing figures in the paper can be found in Kalousová et al. (2023). Open access publishing facilitated by Univerzita Karlova, as part of the Wiley - CzechELib agreement.

## References

Affholder, A., Guyot, F., Sauterey, B., Ferrière, R., & Mazevet, S. (2022). Putative methanogenic biosphere in Enceladus’s deep ocean: Biomass, productivity, and implications for detection. *The Planetary Science Journal*, 3(12), 270. https://doi.org/10.3847/PSJ/aca275

Alnaes, M., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., et al. (2015). The FEniCS project version 1.5. *Archive of Numerical Software*, 3, 100. https://doi.org/10.11588/ans.2015.100.20553

Amsden, A. A., Ruppel, H. M., & Hirt, C. W. (1980). SALE: A simplified ALE computer program for fluid flow at all speeds. https://doi.org/10.2172/5176006

Artemieva, N., & Lunine, J. (2003). Cratering on Titan: Impact melt, ejecta, and the fate of surface organics. *Icarus*, 164(2), 471–480. https://doi.org/10.1016/S0019-1035(03)00148-9

Artemieva, N., & Lunine, J. I. (2005). Impact cratering on Titan II. Global melt, escaping ejecta, and aqueous alteration of surface organics. *Icarus*, 175(2), 522–533. https://doi.org/10.1016/j.icanus.2004.12.005

Baland, R.-M., Tobie, G., Lefèvre, A., & Van Hoolst, T. (2014). Titan’s internal structure inferred from its gravity field, shape, and rotation state. *Icarus*, 237, 29–41. https://doi.org/10.1016/j.icanus.2014.04.007

Baland, R.-M., Van Hoolst, T., Yseboodt, M., & Karatekin, Ö. (2011). Titan’s obliquity as evidence of a subsurface ocean? *Astronomy & Astrophysics*, 530, A141. https://doi.org/10.1051/0004-6361/201116578

Barnes, J. W., Turtle, E. P., Trainer, M. G., Lorenz, R. D., MacKenzie, S. M., Brinckerhoff, W. B., et al. (2021). Science goals and objectives for the Dragonfly Titan rotorcraft relocatable lander. *The Planetary Science Journal*, 2(4), 130. https://doi.org/10.3847/PSJ/abfdefc

Bercovici, D., Ricard, Y., & Schubert, G. (2001). A two-phase model for compaction and damage: 1. General theory. *Journal of Geophysical Research*, 106(B5), 8887–8906. https://doi.org/10.1029/2000JB0900430

Bills, B. G., & Nimmo, F. (2008). Forced obliquity and moments of inertia of Titan. *Icarus*, 196(1), 293–297. https://doi.org/10.1016/j.icanus.2008.03.002

Bills, B. G., & Nimmo, F. (2011). Rotational dynamics and internal structure of Titan. *Icarus*, 214(1), 351–355. https://doi.org/10.1016/j.icanus.2011.04.028

Boury, S., Meyer, C., Vasil, G., & Wells, A. (2021). Convection in a mushy layer along a vertical heated wall. *Journal of Fluid Mechanics*, 926, A33. https://doi.org/10.1017/jfm.2021.742

Bruzzone, L., Plaut, J. J., Alberti, G., Blankenship, D. D., Bovolo, F., Campbell, B. A., et al. (2015). JUpiter ICy moon explorer (JUICE): Advances in the design of the radar for icy Moons (RIME). In *2015 IEEE international geoscience and remote sensing symposium (IGARSS) (pp. 1257–1260)*. https://doi.org/10.1109/IGARSS.2015.7326002

Carnahan, E., Vance, S. D., Cox, R., & Hesse, M. A. (2022). Surface-to-ocean exchange by the sinking of impact generated melt chambers on Europa. *Geophysical Research Letters*, 49(24), e2022GL100287. https://doi.org/10.1029/2022GL100287

Chizhov, V. E. (1993). In *Thermodynamic properties and thermal equations of the state of high-pressure ice phases* (Vol. 34). https://doi.org/10.1007/BF00852521

Choukroun, M., & Sotin, C. (2012). Is titan’s shape caused by its meteorology and carbon cycle? *Geophysical Research Letters*, 39(4), L04201. https://doi.org/10.1029/2011GL050747

Cockell, C., Bush, T., Bryce, C., Direito, S., Fox-Powell, M., Harrison, J., et al. (2016). Habitability: A review. *Astrobiology*, 16(1), 89–117. https://doi.org/10.1089/ast.2015.1295

Collins, G. S., Melosh, H. J., & Ivanov, B. A. (2004). Modeling damage and deformation in impact simulations. *Meteoritics & Planetary Sciences*, 39(2), 217–231. https://doi.org/10.1111/j.1945-5100.2004.tb00337.x

Collins, G. S., Patel, N., Davison, T. M., Rae, A. S. P., Morgan, J. V., Gulick, S. P. S., et al. (2020). A steeply-inclined trajectory for the Chicxulub impact. *Nature Communications*, 11(1), 1480. https://doi.org/10.1038/s41467-020-15269-x

Crosta, A. P., Silber, E. A., Lopes, R. M. C., Johnson, B. C., Bjonnes, E., Malaska, M. J., et al. (2021). Modeling the formation of Menrva impact crater on Titan: Implications for habitability. *Icarus*, 370, 114679. https://doi.org/10.1016/j.icanus.2021.114679

De La Chapelle, S., Milsch, H., Castelnau, O., & Duval, P. (1999). Compressive creep of ice containing a liquid intergranular phase: Rate-controlling processes in the dislocation creep regime. *Geophysical Research Letters*, 26(2), 251–254. https://doi.org/10.1029/1998GL900289

Durante, D., Hemingway, D. J., Racioppa, P., Iess, L., & Stevenson, D. J. (2019). Titan’s gravity field and interior structure after Cassini. *Icarus*, 326, 123–132. https://doi.org/10.1016/j.icanus.2019.03.003

Durham, W. B., Kirby, S. H., Stern, L. A., & Zhang, W. (2003). The strength and rheology of methane clathrate hydrate. *Journal of Geophysical Research*, 108(B4), 2184. https://doi.org/10.1029/2002JB001872

Fowler, A. C. (1985). A mathematical model of magma transport in the asthenosphere. *Geophysical & Astrophysical Fluid Dynamics*, 33(1–4), 63–96. https://doi.org/10.1080/03091928508245423

Goldsby, D. L., & Kohlstedt, D. L. (2001). Superplastic deformation of ice: Experimental observations. *Journal of Geophysical Research*, 106(B6), 11017–11030. https://doi.org/10.1029/2000JB0900336

Grasset, O., Dougherty, M. K., Coustenis, A., Bunce, E. J., Erd, C., Titov, D., et al. (2013). JUpiter ICy moons Explorer (JUICE): An ESA mission to orbit Ganymede and to characterise the Jupiter system. *Planetary and Space Science*, 78, 1–21. https://doi.org/10.1016/j.pss.2012.12.002

Hedgepeth, J. E., Neish, C. D., Turtle, E. P., Stiles, B. W., Kirk, R., & Lorenz, R. D. (2020). Titan’s impact crater population after Cassini. *Icarus*, 344, 113664. https://doi.org/10.1016/j.icanus.2020.113664
Hendrix, A. R., Hurford, T. A., Barge, L. M., Bland, M. T., Bowman, J. S., Brinckerhoff, W., et al. (2019). The NASA roadmap to ocean worlds. *Astrobiology*, **19**(1), 1–27. https://doi.org/10.1089/ast.2018.1955

Hobbs, P. V. (1974). *Ice physics*. Clarendon Press.

Hörst, S. M. (2017). Titan's atmosphere and climate. *Journal of Geophysical Research: Planets*, **122**(3), 432–482. https://doi.org/10.1002/2016JE005240

Johnson, B. C., Sheppard, R. Y., Pascuzzo, A. C., Fisher, E. A., & Wiggins, S. E. (2017). Porosity and salt content determine if subduction can occur in Europa's ice shell. *Journal of Geophysical Research: Planets*, **122**(12), 2765–2778. https://doi.org/10.1002/2017JE005370

Jordan, J. S., & Hesse, M. A. (2015). Reactive transport in a partially molten system with binary solid solution. *Geochemistry, Geophysics, Geosystems*, **16**(12), 4153–4177. https://doi.org/10.1002/2015GC005956

Kalousová, K., & Sotin, C. (2018). Melting in high-pressure ice layers of large ocean worlds—Implications for volatiles transport. *Geophysical Research Letters*, **45**(16), 8096–8103. https://doi.org/10.1029/2018GL078889

Kalousová, K., & Sotin, C. (2020a). Dynamics of Titan's high-pressure ice layer. *Earth and Planetary Science Letters*, **545**, 116416. https://doi.org/10.1016/j.epsl.2020.116416

Kalousová, K., & Sotin, C. (2020b). The insulating effect of methane clathrate crust on Titan's thermal evolution. *Geophysical Research Letters*, **47**(13), e2020GL087481. https://doi.org/10.1029/2020GL087481

Kalousová, K., Sotin, C., Choblet, G., Tobie, G., & Grasset, O. (2018). Two-phase convection in Ganymede's high-pressure ice layer— Implications for its geological evolution. *Icarus*, **299**, 133–147. https://doi.org/10.1016/j.icanus.2017.07.018

Kalousová, K., Souček, O., Tobie, G., Choblet, G., & Čadek, O. (2016). Water generation and transport below Europa's strike-slip faults. *Journal of Geophysical Research: Planets*, **121**(12), 2444–2462. https://doi.org/10.1002/2016JE005188

Kalousová, K., Wakita, S., Sotin, C., Neish, C. D., Soderblom, J. M., Souček, O., & Johnson, B. C. (2023). Evolution of impact melt pools on Titan —Data files. *Zenodo*. https://doi.org/10.5281/zenodo.8329373

Katz, R. F. (2008). Magma dynamics with the enthalpy method: Benchmark solutions and magmatic focusing at mid-ocean ridges. *Journal of Petrology*, **49**(12), 2099–2121. https://doi.org/10.1093/petrology/egn058

Kivelson, M. G., Khurana, K. K., & Volwerk, M. (2002). The permanent and inductive magnetic moments of Ganymede. *Icarus*, **157**(2), 507–522. https://doi.org/10.1006/icar.2002.6834

Krawczynski, M. J., Behn, M. D., Das, S. B., & Joughin, I. (2009). Constraints on the lake volume required for hydro-fracture through ice sheets. *Geophysical Research Letters*, **36**(10), 1–5. https://doi.org/10.1029/2008GL036765

Lewis, J. S. (1971). Satellites of the outer planets: Their physical and chemical nature. *Icarus*, **15**(2), 174–185. https://doi.org/10.1016/0019-1035 (71)90072-8

Logg, A., Mardal, K.-A., & Wells, G. (Eds.) (2012). *Automated solution of differential equations by the finite element method* (Vol. 84). Springer-Verlag. https://doi.org/10.1007/978-3-642-23099-8

Loison, J.-C., Wakelam, V., Gratier, P., & Hickson, K. M. (2019). Chemical nitrogen fractionation in dense molecular clouds. *Monthly Notices of the Royal Astronomical Society*, **484**(2), 2747–2756. https://doi.org/10.1093/mnras/sty3293

Lunine, J. I., & Stevenson, D. J. (1987). Clathrate and ammonia hydrates at high pressure: Application to the origin of methane on titan. *Icarus*, **70**(1), 61–77. https://doi.org/10.1016/0019-1035(87)90075-3

MacKenzie, S. M., Birch, S. P. D., Hörst, S., Sotin, C., Barth, E., Lora, J. M., et al. (2021). Titan: Earth-like on the outside, ocean world on the inside. *The Planetary Science Journal*, **2**(3), 112. https://doi.org/10.3847/PSJ/abf7c9

Maillard, J., Carrasco, N., Schmitz-Afonso, I., Gautier, T., & Afonso, C. (2018). Comparison of soluble and insoluble organic matter in analogues of Titan's aerosols. *Earth and Planetary Science Letters*, **495**, 185–191. https://doi.org/10.1016/j.epsl.2018.05.014

McKenzie, D. (1984). The generation and compaction of partially molten rock. *Journal of Petrology*, **25**(3), 713–765. https://doi.org/10.1093/petrology/25.3.713

Melosh, H. J. (1989). *Impact cratering: A geologic process*. Oxford University Press.

Meriggiola, R., Iess, L., Stiles, B. W., Lunine, J. I., & Mitri, G. (2016). The rotational dynamics of Titan from Cassini RADAR images. *Icarus*, **275**, 183–192. https://doi.org/10.1016/j.icanus.2016.01.019

Miller, K. E., Daswani, M. M., Sotin, C., Malaska, M. J., Neish, C. D., Farnsworth, K., et al. (2023). Combined geochemical and geophysical constraints on Titan interior composition. In *AAS/Division for planetary sciences meeting abstracts* (Vol. 55).

Mitri, G., Meriggiola, R., Hayes, A., Lefevre, A., Tobie, G., Genova, A., et al. (2014). Shape, topography, gravity anomalies and tidal deformation of Titan. *Icarus*, **236**, 169–177. https://doi.org/10.1016/j.icanus.2014.03.018

Montgomery, L., Miège, C., Miller, J., Scambos, T. A., Wallin, B., Miller, O., et al. (2020). Hydrologic properties of a highly permeable firn aquifer in the Wilkins Ice Shelf, Antarctica. *Geophysical Research Letters*, **47**(22), e2020GL089552. https://doi.org/10.1029/2020GL089552

Neish, C. D., & Lorenz, R. D. (2012). Titan's global crater population: A new assessment. *Planetary and Space Science*, **60**(1), 26–33. https://doi.org/10.1016/j.pss.2011.02.016

Neish, C. D., Lorenz, R. D., Turtle, E. P., Barnes, J. W., Trainer, M. G., Stiles, B., et al. (2018). Strategies for detecting biological molecules on Titan. *Astrobiology*, **18**(5), 571–585. https://doi.org/10.1089/ast.2017.1758

Neish, C. D., Malaska, M. J., Sotin, C., Lopes, R. M. C., Nixon, C. A., Affholder, A., et al. (2024). Organic Input to Titan's Subsurface Ocean Through Impact Cratering. *Astrobiology*. https://doi.org/10.1089/ast.2023.0055

Neish, C. D., Somogyi, A., & Smith, M. A. (2010). Titan's primordial soup: Formation of amino acids via low-temperature hydrolysis of tholins. *Astrobiology*, **10**(3), 337–347. https://doi.org/10.1089/ast.2009.0402

Pierazzo, E., & Melosh, H. J. (2000). Melt production in oblique impacts. *Icarus*, **145**(1), 252–261. https://doi.org/10.1006/icar.1999.6332

Potter, R. W. K., Collins, G. S., Kiefer, W. S., McGovern, P. J., & Kring, D. A. (2012). Constraining the size of the South Pole-Aitken basin impact. *Icarus*, **220**(2), 730–743. https://doi.org/10.1016/j.icanus.2012.05.032

Rudge, J. F., Bercovici, D., & Spiegelman, M. (2011). Disequilibrium melting of a two phase multicomponent mantle. *Geophysical Journal International*, **184**(2), 699–718. https://doi.org/10.1111/j.1365-246X.2010.04870.x

Schmeling, H., Dannberg, J., Dohmen, J., Kalousová, K., Maurice, M., Noack, L., et al. (2016). A benchmark initiative on mantle convection with melting and melt segregation. In *EGU general assembly conference abstracts* (p. EPSC2016-8449).

Schurmeier, L. R., Brouwer, G. E., Fagents, S. A., Kay, J. P., Marusiak, A. G., & Vance, S. D. (2023). Crater relaxation caused by an insulating methane clathrate crust on Titan. In *Lunar and planetary science conference* (Vol. 54).2813.

Schurmeier, L. R., & Dombard, A. J. (2018). Crater relaxation on Titan aided by low thermal conductivity sand infill. *Icarus*, **305**, 314–323. https://doi.org/10.1016/j.icanus.2017.10.034

Simon, F., & Glatzel, G. (1929). Bemerkungen zur schmelzdruckkurve. *Zeitschrift für Anorganische und Allgemeine Chemie*, **178**(1), 309–316. https://doi.org/10.1002/zaac.19291780123

Sloan, E. D., Koh, C. A., & Koh, C. (2007). Clathrate hydrates of natural gases. CRC Press.
Sotin, C., Kalousová, K., & Tobie, G. (2021). Titan's interior structure and dynamics after the Cassini-Huygens mission. *Annual Review of Earth and Planetary Sciences*, 49(1), 579-607. https://doi.org/10.1146/annurev-earth-072920-052847

Souček, O., Kalousová, K., & Čadek, O. (2014). Water transport in planetary ice shells by two-phase flow—A parametric study. *Geophysical & Astrophysical Fluid Dynamics, 108*(6), 639–666. https://doi.org/10.1080/03091929.2014.969251

Šrámek, O., Ricard, Y., & Bercovici, D. (2007). Simultaneous melting and compaction in deformable two-phase media. *Geophysical Journal International*, 168(3), 964–982. https://doi.org/10.1111/j.1365-246X.2006.03269.x

Steinbrügge, G., Stark, A., Hussmann, H., Sohl, F., & Oberst, J. (2015). Measuring tidal deformations by laser altimetry. A performance model for the Ganymede Laser Altimeter. *Planetary and Space Science*, 117, 184–191. https://doi.org/10.1016/j.pss.2015.06.013

Stiles, B. W., Kirk, R. L., Lorenz, R. D., Hensley, S., Lee, E., Ostro, S. J., et al. (2008). Determination Titan's spin state from Cassini radar images. *The Astronomical Journal*, 135(5), 1669–1680. https://doi.org/10.1088/0004-6256/135/5/1669

Tobie, G., Choblet, G., & Sotin, C. (2003). Tidally heated convection: Constraints on Europa's ice shell thickness. *Journal of Geophysical Research*, 108(E11), 5124. https://doi.org/10.1029/2003JE002099

Tobie, G., Lunine, J. I., & Sotin, C. (2006). Episodic outgassing as the origin of atmospheric methane on titan. *Nature*, 440(7080), 61–64. https://doi.org/10.1038/nature04497

Turcotte, D. L., & Schubert, G. (2014). *Geodynamics* (3rd ed.). Cambridge University Press.

Turtle, E. P., Barnes, J. W., Trainer, M. G., Lorenz, R. D., MacKenzie, S. M., Hibbard, K. E., et al. (2017). Dragonfly: Exploring titan's prebiotic organic chemistry and habitability. In *48th annual lunar and planetary science conference* (p. 1958).

Vuitton, V., Yelle, R. V., Klippenstein, S. J., Hörst, S. M., & Lavvas, P. (2019). Simulating the density of organic species in the atmosphere of Titan with a coupled ion-neutral photochemical model. *Icarus*, 324, 120–197. https://doi.org/10.1016/j.icarus.2018.06.013

Wakita, S., Genda, H., Kurosawa, K., & Davison, T. M. (2019). Enhancement of impact heating in pressure-strengthened rocks in oblique impacts. *Geophysical Research Letters*, 46(23), 13678–13686. https://doi.org/10.1029/2019GL085174

Wakita, S., Genda, H., Kurosawa, K., Davison, T. M., & Johnson, B. C. (2022). Effect of impact velocity and angle on deformational heating and postimpact temperature. *Journal of Geophysical Research: Planets*, 127(8), e2022JE007266. https://doi.org/10.1029/2022JE007266

Wakita, S., Johnson, B. C., Soderblom, J. M., Shah, J., & Neish, C. D. (2022). Methane-saturated layers limit the observability of impact craters on titan. *The Planetary Science Journal*, 3(2), 50. https://doi.org/10.3847/PSJ/ac4e91

Wakita, S., Johnson, B. C., Soderblom, J. M., Shah, J., Neish, C. D., & Steckloff, J. K. (2023). Modeling the formation of Selk impact crater on titan: Implications for Dragonfly. *The Planetary Science Journal*, 4(3), 51. https://doi.org/10.3847/PSJ/acbe40

Wiggins, S. E., Johnson, B. C., Collins, G. S., Jay Melosh, H., & Marchi, S. (2022). Widespread impact-generated porosity in early planetary crusts. *Nature Communications*, 13(1), 4817. https://doi.org/10.1038/s41467-022-32445-3

Wünnemann, K., Collins, G. S., & Melosh, H. J. (2006). A strain-based porosity model for use in hydrocode simulations of impacts and implications for transient crater growth in porous targets. *Icarus*, 180(2), 514–527. https://doi.org/10.1016/j.icarus.2005.10.013

Wünnemann, K., Collins, G. S., & Osinski, G. R. (2008). Numerical modelling of impact melt production in porous rocks. *Earth and Planetary Science Letters*, 269(3), 530–539. https://doi.org/10.1016/j.epsl.2008.03.007

Yung, Y. L., Allen, M. B., & Pinto, J. P. (1984). Photochemistry of the atmosphere of titan: Comparison between model and observations. *The Astrophysical Journal*, 3, 465–506. https://doi.org/10.1086/190963