Modeling the Formation of Selk Impact Crater on Titan: Implications for Dragonfly

NASA's upcoming Dragonfly mission is scheduled to arrive at Titan in the 2030s to characterize its habitability and search for potential chemical biosignatures (Barnes et al. 2021). One of the main science goals of the mission is to examine how liquid water may have mixed with organics on Titan's surface, to produce molecules of prebiotic interest (e.g., Neish et al. 2010). To accomplish these goals, the Dragonfly spacecraft will target the region around the Selk impact crater for exploration (Lorenz et al. 2021). Crater-forming impacts shock heat the target material, generating melt pools of molten bedrock; like all large craters on Titan, Selk likely possessed a liquid water melt pool (Neish et al. 2018). Such environments are conducive to the formation of prebiotic molecules such as amino acids (Neish et al. 2010). Thus, investigating the production of impact melt pools is essential for understanding Titan's surface as a biological environment.

Selk crater is 84 ± 2 km in diameter and ${0.47}_{-0.105}^{+0.125}$ km in depth (Hedgepeth et al. 2020), with materials rich in water ice observed around the crater (Soderblom et al. 2010; Janssen et al. 2016), and organic-rich dunes beyond the crater rim (Soderblom et al. 2010; Lorenz et al. 2021). Although impact melt formation has been examined previously for Titan (Artemieva & Lunine 2003, 2005; Crósta et al. 2021), those works assumed that the target is water ice, ignoring the presence of methane clathrate in the substrate. Methane clathrate hydrate is stable at Titan's present-day surface temperature and pressure conditions and is expected to constitute a significant portion of Titans upper crust (Choukroun et al. 2010; Vu et al. 2020). Since methane clathrate is much stronger than water ice near the melting temperature of water ice (Durham et al. 2003), the presence of methane clathrate in the target can affect the formation and evolution of Titan's impact craters. Indeed, previous work showed that a methane clathrate target results in slightly smaller craters than a pure water-ice target (Wakita et al. 2022b). As such, target materials would change the crater size for a given impactor size.

The presence of methane clathrates, however, also has a strong effect on the target temperature at depth. Methane clathrate has a lower thermal conductivity (0.5 W m⁻¹ K⁻¹) than water ice (2.2 W m⁻¹ K⁻¹) at 263 K (Sloan et al. 2007), and the presence of methane clathrate in the crust results in higher lithospheric thermal gradients than would exist for crusts made purely of water ice with the same heat flow (Kalousová & Sotin 2020). The target material weakens at these higher temperatures and becomes strengthless as the temperature approaches its melting temperature. Thus, the vertical temperature profile of the target significantly influences crater formation and acts as a primary control on crater morphology (Bray et al. 2014; Silber & Johnson 2017; Bjonnes et al. 2022). The effect of methane clathrate's thermal properties is also missing in previous works modeling the production of impact melts on Titan.

Here we explore the formation of a Selk-sized crater on Titan. We simulate impacts with different initial conditions, considering a top surface layer of methane clathrate, which is accompanied by a higher thermal gradient and increased strength. We explore how craters on Titan vary depending on the thickness of methane clathrate and the associated temperature profile. We also examine the distribution of the impact-induced melt. By investigating the formation and evolution processes of Selk crater, with a particular focus on the formation and source of melt pools and the thermal effects of methane clathrates, we provide crucial information needed to understand one of the main targets of Dragonfly's scientific investigations.

We perform crater-forming impact simulations on Titan using the iSALE-2D shock physics code. The code has been developed to model planetary impacts and cratering (Wünnemann et al. 2006; Collins et al. 2016). It is based on the SALE hydrocode (Amsden et al. 1980) and has been improved by including various equations of state (EOS) and strength models (Melosh et al. 1992; Ivanov et al. 1997; Collins et al. 2004). We assume that a spherical icy impactor hits a flat target with an impact velocity of 10.5 km s⁻¹, the average velocity of an impact on Titan (Zahnle et al. 2003). As we perform two-dimensional simulations in an axisymmetric coordinate system, all simulations explored here assume vertical impacts. We consider spherical icy impactors (D_imp) with diameters that range from 3 to 4 km. iSALE-2D has a high-resolution zone and an extension zone. We use a resolution of 50 m for the high-resolution zone that includes the impactor and the target down to a depth of 15 km and out to a radius of 70 km. The cell size in the extension zone is incrementally increasing by 2% from the previous cell up to a maximum cell size of 1000 m (20 times larger than the cell size in the high-resolution zone). While we consider water ice as the composition of the impactor, we consider methane clathrate and water ice for the composition of the target. Because shock physics EOS or Hugoniot data for methane clathrate is currently unavailable, we use an analytical equation of state (ANEOS) for water ice for both compositions (see below), and the strength model of methane clathrate developed in previous work (Wakita et al. 2022b). The material input parameters used in this work are shown in Table 1.

We explore the sensitivity of Titan's impacts to (i) the thickness of a methane clathrate layer overlying a solid water-ice layer and (ii) the crustal temperature profile. Since the thermal conductivity of methane clathrate is much lower than that of water ice (Sloan et al. 2007), the temperature gradient within the methane clathrate conductive lid is higher than for the water-ice case with the same heat flow, as shown in Kalousová & Sotin (2020). To explore the effect of the temperature gradient, we study a 10 km thick methane clathrate layer with various temperature gradients ranging from 5 to 20 K km⁻¹. When we assume that a conductive ice lid lays on top of a convective ice layer, the thickness of the conductive lid depends both on the temperature gradient and on the temperature of the convecting ice (Silber & Johnson 2017). Assuming a temperature of 255 K for the convecting ice (Kalousová & Sotin 2020), we model the thickness of the conductive ice lid to be 32.4 km, 16.15 km, 10.75 km, and 8.06 km for temperature gradients of 5 K km⁻¹, 10 K km⁻¹, 15 K km⁻¹, and 20 K km⁻¹, respectively. We take the lower limit of 5 K km⁻¹, which is similar to the no-methane-clathrate case (3 K km⁻¹; Wakita et al. 2022b), and the upper value from the thin-methane-clathrate case (23 K km⁻¹; Kalousová & Sotin 2020). The water-ice layer beneath the 10 km thick methane clathrate layer is therefore in the conductive lid for all temperature gradients considered except for the 20 K km⁻¹ case. Using a constant temperature profile, however, can lead to errors in the heat flux. By definition, the heat flux must be the same on each side of the boundary between the methane clathrate layer and the water-ice layer. Because of the different thermal conductivity of the two materials, however, this requires different temperature gradients on each side of the boundary. Thus, our simple-temperature-profile model that assumes the same temperature gradient at the boundary leads to an inaccurate heat flux. As such, we cannot trust the heat flux between the methane clathrate layer and the water-ice layer for these simple cases. Nevertheless, our simplified models allow us to directly explore the effects of the temperature gradient on the crater morphology. We explore and discuss the effects of the temperature of the convective layer below (see Section 4, Discussion). Next, we consider temperature profiles for 5, 10, and 15 km thick methane clathrate layers, based on thermal modeling in Kalousová & Sotin (2020); these temperature profiles are likely to be more accurate because they enforce a continuous heat flux at the methane clathrate–water-ice boundary. Note that the temperature gradient in the conductive lid of methane clathrate is ∼23, ∼15, and ∼10 K km⁻¹ for 5, 10, and 15 km thick cases, respectively. The temperature profiles and yield strengths of these seven models are shown in Figures 1 and 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To identify any material that becomes impact melt (i.e., liquid water), we rely on the peak pressure of tracer particles. Lagrangian tracer particles in iSALE-2D are placed in each cell at the start of the simulation prior to impact (t = 0), and record the pressure at regular intervals throughout the impact process. Because they follow the parcel of material in which they are embedded, they are useful to estimate the status of melt. Since ANEOS is a semi-analytical model derived from the first principles of thermodynamics, the thermodynamical quantities of materials (e.g., temperature, pressure, and density) are obtained self-consistently from the Gibbs free energy (Thompson 1990; Melosh 2007). In iSALE, we can use the Tillotson EOS, which is a simple analytical model using a simplified phase diagram. The Tillotson EOS, however, cannot compute thermodynamical quantities consistently, especially at phase transitions. Because ANEOS can treat multiple phases, we use it in this work (see also Section 4). Figure 3 shows the Hugoniot curve of water ice using ANEOS. To determine the fate of material we compare the entropies of the solidus (2360 J kg⁻¹ K⁻¹) and liquidus (3580 J kg⁻¹ K⁻¹) to entropies and pressures along the Hugoniot curve. Incipient melting occurs when the material pressure exceeds the intersect of the solidus and Hugoniot curve (5.34 GPa). Complete melting occurs when its pressure exceeds the intersect of the liquidus and Hugoniot curve (8.69 GPa). We adopt these thresholds to determine incipient and complete melting of materials during the impact process. Note that these thresholds are similar to previous work that also uses ANEOS (Artemieva & Lunine 2003), but differ from other work that uses thermochemical tables (Pierazzo et al. 1997). While there is a more complex EOS for water ice (five-phase EOS), its estimated amount of complete melting is similar to that of ANEOS (Kraus et al. 2011). Thus, the amount of melt predicted using the less complex ANEOS is reasonable.

Zoom In Zoom Out Reset image size

In this work, we consider the incipient and complete melt of tracer particles according to their peak pressure. We calculate their volumes based on their initial spacing and locations (e.g., Johnson & Melosh 2014). While we will discuss the volumes of incipient and complete melt, we do not treat the volume in between them (i.e., partial melting tracer particles).

To account for melting of methane clathrate, we apply the same approach. Although the Hugoniot data of methane clathrate is unavailable, the melting curve of methane clathrate is given experimentally (Sloan et al. 2007). Methane clathrate has a similar phase curve to water ice in a high-pressure environment (≳1.5 GPa; Journaux et al. 2020). Since our threshold pressure of incipient melt (5.34 GPa) exceeds 1.5 GPa, we assume that methane clathrate follows the same criteria for dissociation as that of water ice. Moreover, methane clathrate may transform to water ice at lower pressure by releasing methane gas. The released methane gas would change the target material's volume, resulting in compression that would lead to an increase of internal energy and additional heating (Melosh 1989). The compaction of any preexisting porosity in the target would also lead to similar compression and heating. As such, our melt estimates should be treated as lower limits.

3.1. Effect of a Temperature Gradient in 10 km Thick Methane Clathrate Layer

We first consider the sensitivity of crater formation on the temperature gradients assuming a 10 km thick methane clathrate layer over water ice and a 4 km diameter icy impactor. Figure 4 illustrates the time series of crater formation with a temperature gradient of 15 K km⁻¹. Soon after impact a transient crater (19 km in depth and 14 km in radius; Figure 4(a)) forms, and then collapses producing a central uplift (Figure 4(b)). This central uplift then collapses, pushing warm material over the apparent crater rim (Figure 4(c)). Because the overflow of this warm material covers the rim, the rim is obscured when the crater is finally formed (Figures 4(d) and 5, and supplemental movie Figure A1). Because this overflow is likely exaggerated by the axisymmetry of our numerical simulations, the rim after the overflow is not representative of the crater rim. Similar behavior is observed in modeling crater formation on other icy satellites with higher temperature gradients such as Europa (Silber & Johnson 2017). When such overflow occurs in these models, we must examine the crater diameter before this overflow covers the rim (i.e., the time of rim formation; Figure 4(b)); following Silber & Johnson (2017), we take a rim location between the time that the transient crater collapses and the central uplift collapses. Figure 5 depicts the surface profile before the overflow (t = 400 s) and after the overflow (t = 2000 s). Although their rim locations are coincidentally near, the profile at the earlier time (t = 400 s) is more pronounced than the later time (t = 2000 s). As such, defining such a topographic height, which clearly exceeds the original surface as a rim, we determine the crater diameter and depth before the overflow occurs (Silber & Johnson 2017). The crater depth is taken as the vertical distance between the rim height at this intermediate time step and the lowest location in the crater floor at the end of the simulation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The temperature gradient governs the diameter and depth of craters. The top four rows in Table 2 show the diameters and depths of craters formed in a 10 km thick methane clathrate layer over water ice assuming different linear temperature gradients and a convective temperature of 255 K. In the case of the 5 K km⁻¹ thermal gradient, the transient crater is shallow (16 km; see Figure 6(a)) and remains entirely in a relatively high-yield-strength region of the target. As such, craters formed in a 5 K km⁻¹ target are relatively small and deep. For the higher thermal gradients (15 and 20 K km⁻¹ cases), the transient crater is deeper (∼20 km in both cases; see Figure 6) and reaches a lower-yield-strength region that allows for greater uplift of the water-ice layer. Thus, craters formed in 15–20 K km⁻¹ targets are wider and shallower and, because the yield strength at 20 km depth is similar for both scenarios (Figure 1), the depths are similar. The 10 K km⁻¹ case is intermediate between these scenarios: most of the transient crater resides in the colder, higher-yield-strength material, with only the base reaching the warmer and weaker material. As a result, craters formed in this target have diameters that are comparable to the 5 K km⁻¹ case and depths that are similar to the 15 and 20 K km⁻¹ cases (see Table 2 and Figure 6). It is worth noting that the crater depth can change by a few meters to a kilometer or more if we use a different convective temperature (see Section 4). Note that the behavior of the rim relates to the crater size; the rim of the 20 K km⁻¹ case behaves like the 15 K km⁻¹ case (see Figure 5).

Zoom In Zoom Out Reset image size

Table 2. Crater Diameters and Depths

Impactor	Crater	Crater	Methane Clathrate	Temperature	Convective
Diameter (km)	Diameter (km)	Depth (km)	Thickness (km)	Gradient (K km−1)	Temperature (K)
4	61	3.2	10	5	255
4	60	1.4	10	10	255
4	90	1.0	10	15	255
4	97	1.2	10	20	255
4	96	1.7	5	Kalousová & Sotin (2020) a
4	90	1.1	10	Kalousová & Sotin (2020) a
4	85	1.2	15	Kalousová & Sotin (2020) a
3.5	59	1.5	5	Kalousová & Sotin (2020) a
3.5	76	1.0	10	Kalousová & Sotin (2020) a
3.5	51	1.6	15	Kalousová & Sotin (2020) a
3	54	1.5	5	Kalousová & Sotin (2020) a
3	64	1.1	10	Kalousová & Sotin (2020) a
3	44	2.1	15	Kalousová & Sotin (2020) a
4	62	3.1	10	5	250
4	61	1.7	10	10	250
4	89	1.3	10	15	250
4	90	1.2	10	20	250
4	62	3.2	10	5	260
4	56	0.5	10	10	260
4	NA b	NA b	10	15	260
4	NA b	NA b	10	20	260

Notes. All craters are given by vertical impacts with 10.5 km s⁻¹.

^a We use the temperature profile in Kalousová & Sotin (2020). See also Figure 2. ^b Almost flat with little topography.

Download table as:  ASCIITypeset image

3.2. Effect of Methane Clathrate Layer Thickness

3.3. Effect of Impactor Size

Next, we explore how smaller impactors affect crater formation and crater morphology on Titan. We consider the same range of targets as in the previous section: a 5, 10, or 15 km thick methane clathrate layer and the corresponding temperature profile (Figure 2; Kalousová & Sotin 2020). Smaller impactors have less impact kinetic energy, and therefore open smaller transient craters; for both 3 and 3.5 km diameter impactors, the transient craters are less than 15 km in depth. We still observe uplifting of the water-ice layer in the case of 5 km and 10 km thick methane clathrate layers, but the uplift of water ice fails to occur if the methane clathrate layer is 15 km thick. This is because the transient crater forms entirely within the methane clathrate layer, which has a higher yield strength than the water-ice layer (see blue dashed line in Figure 2).

The occurrence of water-ice uplift directly affects crater depth. Figure 7 summarizes the diameter and depth of craters formed for the ranges of impactor sizes and methane clathrates layers we consider. Impact simulation studies typically assume uncertainties of two or three cells; we take three cells for the uncertainty on each side and depth of the crater (i.e., six cells). Although we plot these uncertainties as error bars, they are relatively small and lie almost within the symbols. Since uplift of the water-ice layer occurs in the case of both 5 and 10 km thick methane clathrate layers, these crater depths are similar, regardless of the impactor size (see red and green symbols in Figure 7). The depth of craters in a 15 km thick methane clathrate layer, however, is inversely proportional to the diameter of the impactor. Smaller impactors produce smaller and shallower transient craters and, in turn, also produce smaller and shallower final craters. As discussed above, we determine the rim location before the crater collapse (e.g., Figure 5); for consistency, we do this regardless of the occurrence and absence of overflow.

Zoom In Zoom Out Reset image size

3.4. Production of the Melt Pool

Although similar-sized craters can form under different structural conditions, the resulting distribution of heated materials varies significantly. Figure 8 illustrates the distributions of tracer particles colored by their peak pressure for a 4 km diameter impactor into a target with a 5, 10, or 15 km methane clathrate layer and the corresponding temperature profile of Kalousová & Sotin (2020). We color the tracer particles only if the peak pressure exceeds the incipient melting pressure (5.34 GPa; see Section 2). The black dotted lines indicate the material boundary; methane clathrate and water-ice distribute near the surface in cases of 5 and 10 km thick methane clathrate layers, but methane clathrate dominates the 15 km thick cases. For the case of a 5 km thick methane clathrate layer, the heated materials are located below the crater floor as a layered structure and distributed along the crater wall and rim (Figure 8(a)). On the other hand, the distribution of heated materials in the 10 km and 15 km thick layers are relatively complex (Figures 8(b), (c)). The folding of methane clathrate during the crater collapse creates torus-like structures of heated material (e.g., Figures 4 and A1). Note that there is a gap between the surface (black dashed lines) and the tracer particles (i.e., Figure 8(a)). When material collapses inward to the symmetry axis, a single tracer particle represents several cells of material. After this material collapses outward, those tracer particles corresponding to several cells deviate from the surface that is represented by the cell, resulting in potentially erroneous results at the boundary of the simulation; this may be the cause of the heated materials (colored tracer particles) being distributed along the axis in the 10 km thick case. Figure 9 illustrates the provenance plots of melt with the same impact conditions as in Figure 8. This clearly shows that the heated material originates only from the methane clathrate layer in the 10 and 15 km thick layers. On the other hand, the water-ice layer also contributes to the heated material for the 5 km thick layer case. Overall, the thickness of methane clathrate affects the distribution of heated material and its origin.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our results confirm that impact-induced heat can produce melt pools in Selk-sized craters on Titan and demonstrate that such melt occurs in the presence of a methane clathrate layer. Figure 10 represents the incipient and complete melt volume corresponding to the horizontal direction with the same impact conditions as in Figure 8. Note that we calculate the volume of each cell using the initial location of the tracer particles (see Section 2). For a given impactor diameter, the total incipient volume of melt is similar regardless of the initial target structure (thickness of methane clathrate): a 4 km diameter impactor produces ∼360–400 km³ of melt, a 3.5 km diameter impactor produces ∼240–260 km³, and a 3 km diameter impactor produces ∼150–160 km³, respectively. ⁶ The total volume of complete melt is also similar (190–220 km³, 130–140 km³, and 80–90 km³ for 4, 3.5, and 3 km diameter impactors, respectively). We note that higher impact velocities will produce larger total melt volumes (Pierazzo et al. 1997; Artemieva & Lunine 2005). A 3 km diameter impactor with a higher impact velocity would produce a similar amount of melt as a 4 km diameter impactor at 10.5 km s⁻¹. And, finally, as we ignore melted ejecta excavated from the crater, these total volumes include only material within the crater. As Figure 9 shows, the total volume of (incipient/complete) melt is almost completely independent of the methane clathrate thickness. Nevertheless, as mentioned above, the spatial distribution of melt is affected, with melt being concentrated within 20 km of the crater center (see Figures 8(b), (c) and Figures 10(b), (c)) if the clathrate layer is ≥10 km thick, and being broadly distributed over the crater floor, wall, and rim, if the clathrate later is 5 km thick (see Figures 8(a) and 10(a)). These differences in the distribution of the melt will dramatically affect freezing time, because the structure (e.g., thickness) governs the thermal evolution of the melted material. Additionally, the distribution of melt (i.e., liquid water) will influence its ability to mix with the surface material, including possible organic materials, with the wide distribution of the 5 km clathrate layer case (Figure 8(a)) having the greatest potential for such mixing. This implies that the chemistry of the melt pool in a 5 km thick case might be different from other cases. Therefore, even if the total melt volume is similar, the structure and chemistry of the melt pool may differ depending on the thickness of methane clathrate.

Zoom In Zoom Out Reset image size

Our results indicate that impactors a few kilometers in diameter striking a warm target with a surface layer of methane clathrate that is 5–15 km thick can produce Selk-like craters. Furthermore, we find that impacts into targets with higher temperature gradients produce craters that are shallower and wider than those formed in targets with lower temperature gradients (Table 2). A 4 km diameter impactor striking a warm surface with a 10 km or 15 km thick methane clathrate layer forms the most analogous Selk-sized crater (see square green and blue symbols in Figure 7), but a 10 km thick methane clathrate layer with D_imp = 3.5 km can also be a candidate (see green circle in Figure 7). These results suggest that a target with a 10–15 km methane clathrate layer produces similar-sized craters if the impactor is large enough to make a transient crater that is >15 km deep. This is because the yield strength below 15 km depth is similar for such targets (Figure 2), because of the reasonable temperature profiles considering the lower thermal conductivity of methane clathrate. On the other hand, the yield strength below 15 km depth in the case of the 5 km methane clathrate layer is higher than other cases (Figure 2), thus the resultant crater differs. This similarity in size only holds when an impactor forms a transient crater deeper than 15 km, which we only see in the 4 km diameter impactor cases.

Postimpact fluvial and aeolian erosion may explain the difference in depth between the model outputs and the measurements of Selk. For the ∼1 km deep craters produced in a warm methane clathrate layer, we would require ∼0.5 km of infilling to reach the reported depth of Selk crater of ${0.47}_{-0.105}^{+0.125}$ km (Hedgepeth et al. 2020). Although we can resolve the rim height with two to three cells in our model (e.g., t = 400 s in Figure 5), here we attempt to compare it with the observed rim height. For our best case for Selk, we here consider a 4 km diameter impactor hitting a target with a 10 km thick methane clathrate layer with the associated temperature profile from Kalousová & Sotin (2020; Figures 7, 11 and Table 2). Our rim height would be ∼0.15 km and lower than Selk's rim height (${0.280}_{-0.054}^{+0.050}$ km; Hedgepeth et al. 2020). While we may underestimate the rim height, the depth from the original surface of ∼0.95 km would be deeper than the observed terrain depth of ${0.190}_{-0.130}^{+0.155}$ km (Hedgepeth et al. 2020). In this case, the required amount of infilling would be ∼0.75 km, which is slightly more than our estimate from the crater depth. This level of infill is consistent with moderate amounts of fluvial erosion on Titan (Neish et al. 2016). However, the sharp crater rim of Selk observed in the SARTopo data (Hedgepeth et al. 2020; Lorenz 2021) suggests that infilling by sand might be a more reasonable explanation, given that fluvial erosion quickly flattens crater rims (Neish et al. 2016). A thick deposit of sand in the crater can also insulate the crater floor, such that moderate amounts of viscous relaxation can further reduce the depth (Schurmeier & Dombard 2018). Indeed, evidence of sand is visible in the Selk crater cavity in Cassini RADAR, Visual and Infrared Mapping Spectrometer (VIMS), and Imaging Science Subsystem (ISS) data (Lorenz et al. 2021). The interpretation that Selk might have been deeper when it formed (∼1 km) is further supported by radarclinometry results that find variations in the height of the rim of Selk (Bonnefoy et al. 2022). Further information from Dragonfly will help us to constrain the processes that have occurred at Selk.

Zoom In Zoom Out Reset image size

It is possible that the ice-shell structure at the time Selk formed might be different from our model, which would change the morphology of the crater. For example, if the heat flux from Titan's core was higher when Selk formed, the temperature gradient would be higher. Since the icy shell moves more readily if the temperature gradient is higher, a higher thermal gradient would make the crater shallower than the results presented here (e.g., Bjonnes et al. 2022). Alternatively, the methane clathrate layer may have contained impurities when Selk formed, such as volatiles or liquids, which would weaken it. For example, if the crust contained more liquid methane and ethane than could diffuse into the ice shell to form clathrates (Choukroun et al. 2010; Vu et al. 2020), the crust would form a fluid-saturated layer. Fluid-saturated layers behave like weak sediment, limiting the topography of craters, making them shallower than observed here (Wakita et al. 2022b).

Moreover, the choice of the temperature in the warm convective layer (i.e., convective temperature) significantly affects crater depth. While we took a convective temperature of 255 K for simple-temperature-profile cases (Section 3.1 and Figure 1), thermal modeling suggests plausible values that range from 254 to 258 K (Section 3.2 and Figure 2; Kalousová & Sotin 2020). Since the temperature governs the yield strength, higher convective temperatures lead to weaker ice layers. Our results of changing the convective temperature profile to 250 K and 260 K are shown in Table 2. Although there is almost no effect on depth for different convective temperatures in the 5 K km⁻¹ cases, the depth of craters changes in 10–20 K km⁻¹ cases (Table 2). When we use the lower convective temperature of 250 K, the craters are deeper by a few hundred meters. When we increase the convective temperature to 260 K, the craters are shallower by a few hundred meters to more than a kilometer, depending on the temperature gradient. Additionally, in case of temperature gradients of 15 and 20 K km⁻¹, the impacts produce little topography, and likely would be unobservable from orbital images (see bottom two rows in Table 2). The effect of the temperature gradient and the convective temperature needs to be more thoroughly explored for a variety of craters on Titan.

Although we assume a fixed impact velocity of 10.5 km s⁻¹ in this study, different impact velocities are also possible. For different impactor velocities, we may use the crater scaling laws (Melosh 1989): lower impact velocities at a given size of impactor forms a smaller crater. However, our choice of impact velocity may also affect the assumed acoustic fluidization parameters. Acoustic fluidization models the weakening mechanism of the target during crater-forming impacts (Melosh 1979). The simple block model is widely used to describe this, assuming a block size of oscillation (Ivanov et al. 1997). Silber et al. (2017) explored the effects of scaling block size; one is based on the impactor size (Wünnemann & Ivanov 2003) and the other is based on transient crater size (i.e., π scaling; Ivanov & Artemieva 2002). They found that craters can be the same size but differ in morphology depending on the choice of scaling. They proposed the scaling of block size by the impactor size is preferable to explain observed lunar craters, especially in a range of 15–20 km diameter (i.e., simple-to-complex craters). We consider the parameters for the acoustic fluidization by using the block size scaling based on the impactor size, following the previous work of icy craters on Ganymede (Bray et al. 2014). Although their impact velocity of 15 km s⁻¹ is higher than ours of 10.5 km s⁻¹, we assume that the parameters for acoustic fluidization are independent of impact velocity. We know, however, that the parameters for acoustic fluidization increase with the impact velocity when the parameters are based on transient craters (Silber et al. 2017). As such, it is possible that parameters based on the impactor size might also be a function of impact velocity. The sensitivity of the crater diameter to acoustic fluidization parameters is minor (Bray et al. 2014), but the crater depth is sensitive to these parameters. Using the same acoustic fluidization parameters in this study (i.e., independent of impact velocity), we have tested a 3.5 km diameter impactor with a velocity of 15 km s⁻¹ hitting a 10 km thick methane clathrate layer with the thermal profile as in Figure 2. As expected, the resulting crater is larger in diameter (94 km) than the crater formed by the same-sized impactor hitting at 10.5 km s⁻¹ (76 km) and in fact is close to the diameter of a crater formed by a 4 km diameter impactor traveling at 10.5 km s⁻¹ (90 km; see Table 2). Interestingly, a crater formed by this higher-velocity impactor is not substantially deeper (1.2 versus 1–1.2 km; Table 2) than a crater formed by similarly sized impactors impacting with a lower velocity. This is because the uplift of the water-ice layer causes shallow craters in all cases (see Section 3.2). Thus, unlike other bodies, impactor velocity may exhibit a greater control on the diameter of Titan's impact craters than on the depth of its craters.

There are other numerical considerations in our simulation methodology that may also influence the crater depth. In this work, we used the ANEOS for water ice, but there are other EOSs for water ice. We tested the Tillotson EOS for water ice (Tillotson 1962; Ivanov et al. 2002) and found that runs with the Tillotson EOS tend to form slightly shallower craters compared to the ANEOS. While ANEOS includes the phase transition (i.e., melting of water ice), the Tillotson EOS does not. As a result, Tillotson overestimates the amount of heated (and melted) material and makes the material move more readily. To account for the melt pool formed by the crater-forming impacts, our results using ANEOS are more accurate.

The production of a melt pool at Selk is extremely important for aspects of astrobiology on Titan. Organic molecules can react with the liquid water to form prebiotic molecules over timescales of years (Neish et al. 2010), but the longer the liquid water remains before freezing the more time the system has to produce more complex molecules. Thus, we briefly discuss the lifetimes of the melt pool here. A crater formed in a 15 km thick methane clathrate layer has a melt pool 8–10 km wide and 2–3 km deep (see Figures 8 and 10). It would take a torus shape, but we can estimate a freezing timescale from the depth of the deposit, assuming it is losing its heat through thermal conduction to the surface. For 1 km and 4 km thick deposits of liquid water on Titan, Neish et al. (2006) estimated freezing timescales of 5000 yr and 90,000 yr, respectively. Therefore, these deposits may last for tens of thousands of years before freezing completely. These timescales are much longer than those calculated by Hedgepeth et al. (2022) for a melt pool a few hundred meters thick, who found that it takes just a few hundred years to freeze such melt pools. These estimates for melt-pool depth came from previous modeling work, which assumed impact into cold water ice. Our results show that the presence of both methane clathrate and a warmer temperature profile result in much deeper melt pools with significantly extended lifetimes. This would increase the likelihood of the formation of complex organics on Titan. These molecules would then concentrate within the melt sheet at different depths (Hedgepeth et al. 2022). If these melt sheets are later eroded by fluvial processes, exposing these layers, Dragonfly may be able to detect the prebiotic species formed there (Neish et al. 2018). However, subsequent simulations are necessary to better constrain the lifetime and concentration of complex organic (carbon-hydrogen-nitrogen) molecules in the torus-shaped melt pool.

Note that since our simulations are axisymmetric, we are unable to reproduce some features of Selk crater. For example, Selk exhibits a polygonal-shaped rim, which is most pronounced at its northeast (i.e., the "corner"; Soderblom et al. 2010). Such polygonal-shaped structures can be seen on other planets and moons (e.g., the Meteor Crater, Arizona; Shoemaker 1963). One hypothesis for their formation involves the preexistence of faults that are exploited during the crater formation process (Melosh 1989). Because we perform two-dimensional impacts in an axisymmetric coordinates systems, all resultant craters become circular shapes. Three-dimensional simulations capable of producing more complex-shaped craters would be necessary to model the polygonal shape of Selk crater. Note that oblique impacts in the three-dimensional simulations may enlarge the melt volume within the methane clathrate layer (Wakita et al. 2019, 2022a). Nevertheless, their change in the topography might be modest (Davison & Collins 2022). As we focus on the crater size (diameter and depth), the fate of impact ejecta is also beyond the scope of our work. Future work on impact ejecta may elucidate the origin of the several hundred kilometer long landform east-southeast of Selk (i.e., the "bench"; Soderblom et al. 2010) that has been hypothesized to be a fluidized-ejecta flow, similar to the impact melt flows that are found around craters on Venus (Phillips et al. 1991; Schaber et al. 1992) and the Moon (Neish et al. 2014).

We performed impact simulations where an icy impactor hits a methane clathrate crust over water ice on Titan. The lower thermal conductivity of methane clathrate leads to a higher temperature profile in the methane clathrate ice shell than in a pure water-ice target. If a 5–15 km thick warm methane clathrate layer exists on Titan's surface, a 4 km diameter impactor would produce a Selk-sized crater. While the resultant craters are shallower than the lower-temperature-gradient case, they are still slightly deeper than observed for Selk. Postimpact erosional processes such as fluvial erosion and aeolian infill likely explain the discrepancy in the crater depth; measurements by the DragonCam instrument on Dragonfly will help us to determine if such postimpact modification can explain this discrepancy or if additional factors must be considered in the crater formation process (Barnes et al. 2021).

We also confirmed the production of large amounts of impact melt beneath Selk-sized craters, particularly if there is an insulating methane clathrate layer that results in a warmer crustal temperature profile. While the total volume of the melt pool is similar for a given impactor size (360–400 km³ for 4 km diameter impactors) for all methane clathrate layer thicknesses we consider, the methane clathrate layer thickness affects the distribution of the melt. The torus-like melt pool produced in a 10–15 km thick methane clathrate layer would take a longer time to freeze—perhaps as long as tens of thousands of years—than a shallower melt pool produced in a target with a 5 km thick methane clathrate surface layer. If a thick layer of methane clathrate existed on Titan's surface when Selk was formed, there was ample time for complex prebiotic chemistry to occur in the melt pool, increasing the likelihood of producing molecules of biological interest. Measurements by the Dragonfly Mass Spectrometer instrument on Dragonfly will be able to confirm these predictions (Grubisic et al. 2021).

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/). We also thank Tom Davison, the developer of the pySALEPlot tool, which we used to create some plots in this work. This research was supported in part through computational resources provided by Information Technology at Purdue, West Lafayette, Indiana. This work was supported by Cassini Data Analysis Program grant No. 80NSSC20K0382. We thank the two reviewers for their helpful feedback. Our model data are given by using iSALE-2D, and our input and output files are available on DataVerse at doi:10.7910/DVN/XEGCCO.

Animated Figure A1 shows impact simulations of a 4 km diameter impactor hitting on 10 km thick methane clathrate with a temperature gradient of 15 K km⁻¹, which covers time sequences of the impacts that are missing in Figure 4.