Evolution of Impact Melt Pools on Titan

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.


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(Baland et al., , 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.
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 organicrich 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 methaneclathrate 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.

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-diameterimpactors 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 methaneclathrate 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 ϕ 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.

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 ϕ and, consequently, the ice volume fraction can be written as 1 ϕ.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): Here v is the mixture velocity, Π is the dynamic pressure, σ = η(∇v + ∇v T ) is the deviatoric stress tensor, with η the effective viscosity (see Section 2.2.2 for more details), Δρ = ρ i ρ w is the density difference between the ice and water, g is the gravity acceleration, α i is the ice thermal expansivity, ΔT is the temperature difference from the reference value T 0 , ρ = ϕρ w + (1 ϕ)ρ 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, ϕ), as (see derivation in Appendix A2): with T the mixture temperature, p the pressure, L 0 the latent heat at reference temperature T 0 and pressure p 0 , and c p k , α k , and ρ k 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: Depending on the actual value of enthalpy h, the mixture can either be completely solid h < h max i ) or partially molten h ≥ h max i ).Note that we do not treat the case of completely molten material (ϕ = 1) warmer than the melting temperature (T > T m ).The temperature and water content are evaluated from the enthalpy as follows (see Appendix A3 for derivation):

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): The particular creep viscosities are computed as: where the exponents m j , n j , and l j , prefactor A j , and activation energy E j are specific for each creep mechanism, d is the grain size, σ II 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., η min ∼ 10 12 Pa s, see also Carnahan et al., 2022).
Similarly, for numerical reasons, we limit the maximum viscosity to ∼10 24 Pa s.We thus prescribe the mixture viscosity as follows: 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) with C the composition function, which is equal to 0 in ice and to 1 in clathrate, k 1 i and k 2 i 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: In some simulations, a parameterized volumetric tidal heating considering Maxwell viscoelasticity is prescribed following Tobie et al. (2003) as: with Q m the maximum heating amplitude and η M = μ/ω the viscosity that maximizes the heating (the Maxwell peak) with μ the shear modulus and ω the orbital frequency.The tidal heating thus depends on stress σ, grain size d, the mixture temperature T and water content ϕ.
Finally, the melting temperature is evaluated as follows (Simon & Glatzel, 1929): All parameter values are given in Tables 1 and 2.

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 2
Creep Parameters of Ice Based on Goldsby and Kohlstedt ( 2001) Table 3 Simulations and Their Results 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: with q the normal conductive heat flux, n the outer unit normal, ϵ Titan's surface emissivity, σ the Stefan-Boltzman constant, and T s Titan's surface temperature.All parameter values are given in Table 1.

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 (Alnaes 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., 2016Kalousová et al., , 2018)).

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 (ξ = 0, 10%, 20%), and the tidal heating amplitude (Q m = 0 and 5 × 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, ϕ, 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 (≲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 ≲ 180 K) and highly viscous (η ≳ 10 19 Pa s) ice, whereas in all cases with a clathrate layer, the molten region is underlain by warmer (and thus less viscous) 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 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 oceanice 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, η b ∼ 1.6 × 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 ∼ 1 kyr.For the reference grain size, d = 1.0 mm (η b ∼ 6.2 × 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 (η b ∼ 2.5 × 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 ξ on the initial melt volume V 0 : for ξ = 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 ξ = 10% (simulation #15, dashed light green), V 0 = 496 km 3 ; and for ξ = 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  (#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 (t1 2 , see also Table 3).Bottom: Time evolution of the maximum melt depth d max m .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).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, ξ = 20%, dot-dashed dark green), and #18 (additional tidal heating, dashed orange).In two cases-h c = 10 km with ξ = 0% (#11, blue) and h c = 10 km with ξ = 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 ∼ 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 ∼ 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 ξ.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 (ξ = 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 (ξ = 20%), or significant internal heating (Q m = 5 × 10 8 W/m 3 ) with a 10-km thick clathrate layer.

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'stherefore 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-toocean 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 ≤ 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 Journal of Geophysical Research: Planets 10.1029/2023JE008107 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 14 mol/Gyr of glycine could be delivered from Titan's surface to its interior.This gives an upper limit for biomass productivity of ∼10 3 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.

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-toocean 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.
with the mixture density defined as Assuming further that v w = v i = v (impermeable limit, e.g.Kalousová et al., 2016) and noting that ∇ ⋅ v = 0 (Equation 1a), we can rewrite the energy balance (Equation A1) as with τ = ϕτ w + (1 ϕ)τ 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 ∇p = ρ i 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): for which it clearly holds that where we have already neglected the time derivative of hydrostatic pressure.Using Equations A6 and A7, we can finally rewrite Equation A5 as where we used Fourier's law for the heat flux q, Δρ = ρ i ρ 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 (α k , c p k , ρ k with k = i or w), we can express the partial enthalpies of ice and water as functions of temperature T and pressure p: with h 0 k the value at reference temperature T 0 and pressure p 0 .Equation A3 can be simplified as (noting that where we neglect the last term since ϕ(1 ϕ) ρ w ρ i ρ L ≤ 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 Inserting Equations A9a and A9b into Equation A11 and setting h 0 i = 0, we obtain: where we define L 0 = h 0 w h 0 i .

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(T,p,ϕ) = c p i (T m T 0 ) Depending on the actual value of enthalpy h, the mixture can either be completely solid h < h max i ) or partially molten h ≥ h max i ) .The particular equations for temperature and water content can be obtained as follows: If h < h max i , then the mixture is completely frozen and there is no liquid water present.Using Equations A12 and A13 and assuming that ϕ = 0, we can write: from which we obtain If, on the other hand, h ≥ h max i , 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: 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 τ a 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) where μ is the viscosity, ρ 1 and ρ 2 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 14 Pa s (ice viscosity at the melting point), b = 30 km, g = 1.35 m s 2 , and considering, as the most extreme case, a water layer above the ice (i.e., taking ρ 1 = 1,000 kg/m 3 and ρ 2 = 920 kg/m 3 ), 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: τ a = 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) 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.

Figure 1 .
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.

Figure 2 .
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).

Figure 3 .
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 (t1 2 , see also Table3).Bottom: Time evolution of the maximum melt depth d max m .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 Figure2for a visual representation of the melt descending through Titan's ice shell into the ocean (simulation #13).

Figure 4 .
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 max m (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).
Time Scales and Comparison With the Terrestrial Setting