During the impact cratering process, extreme levels of energy compress and excavate the target material (Melosh 1989; Osinski & Pierazzo 2013). Hundreds of gigapascals of pressure are converted to heat and kinetic energy, vaporizing and melting the bedrock ice (Melosh 1989). Artemieva & Lunine (2005) estimated that a 10 km crater would have a 100–200 m thick melt pond and that larger craters will produce melt at ∼2%–5% of the total volume of a crater. This means Selk crater (D ≅ 84 km) would have a melt pond ∼80 m in depth (Bray et al. 2012; Hedgepeth et al. 2020). O'Brien et al. (2005) were more liberal with their estimation of melt production, scaling up the ∼3% melt volume estimate of Artemieva & Lunine (2005) to ∼7%, which may produce melt ponds 100–200 m in depth for a crater of Selk's size (Figure 1). This may be a reasonable estimate considering lunar models predict that larger impact craters produce larger melt volume fractions, with nearly an order of magnitude higher volume fraction for 100 km craters than 10 km craters (Cintala & Grieve 1998). Furthermore, Elder et al. (2012) estimated melt volumes of ∼10% for a crater of Selk's size on Ganymede, supporting the higher estimate from O'Brien et al. (2005). This is an ideal situation, though; in reality, the icy substrate will become fractured by impact, allowing for drainage of the denser water into the subsurface ice (Elder et al. 2012). Elder et al. (2012) estimated that crater melt ponds on icy satellites will be significantly shallower than the estimates of O'Brien et al. (2005), with only larger craters (>50 km) being able to retain substantial melt.

 Zoom In Zoom Out Reset image size

The thickness of an impact crater melt pond will depend on the crater shape, which changes depending on its size. Small, simple craters form bowl-like shapes, but larger, complex craters will experience rebound in the center and collapse (Osinski et al. 2011). The largest craters collapse so severely that they become basins (Melosh 1989). On Titan, the transition diameter from simple to complex craters is smaller than most of the observed craters (∼5 km), assuming Titan craters form with the same crater morphology as those on Ganymede (Schenk 2002). The only basin crater observed on Titan is Menrva at ∼400 km in diameter. Unlike simple bowl-shaped craters, complex craters typically have flat floors, especially on icy moons (Schenk 2002; Neish et al. 2013). Therefore, the melt ponds can be considered effectively cylindrical in shape (but it may be more of a torus if a central uplift is present).

In terms of organic material, Artemieva & Lunine (2005) estimated the loss of ∼90% of the original surface organics in the impact cratering event, with the remaining 10% mixing with the water and sustaining minimal damage from the impact. However, the damage by the impact may have more severe effects on some molecules than others (Pierazzo & Chyba 1999). For our application, we examine the organic molecule HCN, motivated by its relative abundance on Titan (Krasnopolsky 2014) and its well-studied chemistry (Coates & Hartshorne 1931; Haynes 2011). However, HCN will vaporize at ∼25.5°C, and regions of the melt may range from −70°C to 225°C (Artemieva & Lunine 2005). Whether the volume of the melt pond is above its vaporization temperature, and for how long, will determine whether HCN is retained. HCN has been estimated to compose 4%–10% of the material condensed to Titan's surface (Vuitton et al. 2014). Nearly half of the equatorial region is covered with, on average, 30 m of sand dunes, which are suspected to be composed of organics (Lorenz et al. 2008). If we assume this is a good estimate of the dune height in the region of Selk crater when it formed, we can estimate the total volume of dune material that mixed into the Selk melt pond as the volume of a cylinder,

where r_crater = 40 km, and h_dunes = 30 m. Then the final concentration of dunes within the melt sheet can be estimated as,

where only 10% of the dunes are expected to make it into the melt unshocked (Artemieva & Lunine 2005), and V_melt varies based on crater size. For an 80 km crater with 7% melt production (∼100 m thick), the concentration of dunes in the melt pond is ∼7%. If we assume a conservative concentration (4%) of HCN within the dune material, the concentration of HCN in the melt would be ∼0.25%, or ∼2.5% for a much thinner melt sheet (10 m). As the melt sheet freezes and impurities are rejected back into it, concentrations can reach as high as 25% depending on the initial conditions used. We therefore use the SF2 model to study a range of HCN concentrations from 0.1%–25% by weight (1–250 parts per thousand, or ppt) for melt sheet thicknesses of 10–250 m.

However, HCN will not be the only organic molecule present in the melt pond. HCN, as with many other compounds, will quickly undergo hydrolysis into its oxygenated daughter products (e.g., amino acids and nucleobases) at a rate that varies by temperature and pH (Ferris et al. 1978; McDonald et al. 1994; Neish et al. 2010; Hörst et al. 2012; Ruiz-Bermejo et al. 2021). HCN's half-life (the time it takes for one-half of it to evolve into its daughter products) is reduced from weeks to minutes by increasing the temperature from 60°C to 120°C, but if the pH drops below 7, the half-life can go to as high as 100 yr (Miyakawa et al. 2002). Extrapolating the results of Miyakawa et al. (2002) to temperatures that we would expect in Titan's melt sheet, we would expect half-lives on the order of years at high pHs (≥8) and up to centuries for low pHs (≤7); these rates are consistent with other low-temperature hydrolysis measurements for HCN (Sanchez et al. 1966; Neish et al. 2008, 2009). Therefore, HCN's daughter products would likely dominate the chemistry deeper into the melt sheet and change the chemical properties of the melt in the process.

We therefore acknowledge that tracing organic compounds in the melt sheet is not as simple as modeling just HCN. Other compounds are present, and the total amount of HCN will be reduced with time as it hydrolyzes into different compounds. Therefore, our results are not entirely representative of what Dragonfly is going to observe. Instead, we use HCN as a tracer for the other organic molecules in the melt sheet. Our work is intended as a first-order approximation of the freezing process, guiding us to the best place(s) within the melt pond to investigate with Dragonfly. We have modeled at initial concentrations of HCN that are as low as the lowest predicted initial concentrations (∼0.1%), to as high a concentration that we would expect in the thickest, deepest part of the melt sheet where the rejected material will concentrate with time (10%–25%).

We use the 1D multiphase reactive transport model (SF2) of Buffo et al. (2020) to determine the freezing profile of Titan impact melt deposits. SF2 is adapted from the sea ice model of Buffo et al. (2018) to accommodate diverse ocean chemistries and planetary thermal environments. Reactive transport modeling is a powerful tool used to simulate physical, chemical, and biological processes at a range of spatial and temporal scales (Steefel et al. 2005). SF2 solves the heat and mass conservation equations governing a reactive porous media using the finite difference method. It employs the enthalpy method to simulate the phase evolution of the system.

In prior work, Buffo et al. (2018) developed a sea ice model that incorporated advection-reaction-diffusion processes using the theory of mushy layers for ice-ocean systems (Worster 1992, 1997; Feltham et al. 2006; Hunke et al. 2011) coupled with parameterized 1D density-driven convection. The model takes a continuum approach and spatially averages state variables across discretization cells and employs the enthalpy method of Huber et al. (2008) to iteratively calculate the liquid fraction of each cell. The mushy layer theory (Feltham et al. 2006) guarantees the conservation of heat and mass. "Mushy" refers to a system in which multiple components or phases coexist. It allows for different regimes (ice, brine, ice+brine) to be described with the same set of equations

Here, ρ is density, T is temperature, t is time, z is the vertical coordinate, k is thermal conductivity, L is the latent heat of fusion for the water to ice phase transformation, ϕ is the liquid fraction, S is HCN concentration, and D is salt diffusivity. Subscripts "ice" and "br" refer to characteristics of the ice and brine components of the cell, respectively, and variables with an over bar are volumetrically averaged quantities. That is, they are averaged using the respective amounts of ice and brine in the mixture, using the current liquid fraction (ϕ) (e.g., a hypothetical value y would be averaged, ). All variables and their values can be found in Table 1.

Table 1. All Variables Used in the Text, along with Their Definition, Values, and Units

Symbol	Definition	Value	Units
α	1D advection coefficient	1.56 × 10−1	⋯
	Vertical brine transported	Calculated	kg
β	Density coefficient	3.139 × 10−4	kg ppt−1
cbr	Brine heat capacity	Calculated	⋯
cice	Water ice heat capacity	2000	J kg−1 K−1
dt	Time discretization	1–10,000	s
dz	Spatial discretization	0.01–0.025	m
D	HCN diffusivity	Calculated a	m2 s−1
H	Enthalpy	Calculated	J kg−1
Hc	Enthalpy of solid cell	Calculated	J kg−1
g	Acceleration due to gravity	1.352	m s−2
h	Distance to interface	Varies	m
H	Melt pond thickness	Varies	m
kHCN	HCN thermal conductivity	0.2256	W m−1 K−1
kbr	Brine thermal conductivity	Calculated	W m−1 K−1
kice	Water ice heat conductivity	2.0	W m−1 K−1
κ	Thermal diffusivity	Varies a	m2 s−1
L	Latent heat of fusion (ice)	334,774	J kg−1
μ	Kinematic viscosity	1.55 × 10−3 a	m2 s−1
ϕ	Liquid fraction	Calculated	⋯
ϕc	Critical porosity	0.05 a	⋯
ρbr	Brine density	Varies	kg m−3
ρice	Water ice density	917	kg m−3
Π	Permeability	Calculated a	m2
Ra	Rayleigh number	Calculated	⋯
Rac	Critical Rayleigh number	1.01 × 10−2	⋯
S	Concentration	Calculated	ppt
T	Temperature	Calculated	K
Tm	Melting temperature	94	K
t	Time	⋯	s
w	Brine velocity	0	m s−1
	HCN liquid density	687 b	kg m−3
	HCN melting temperature	−13.3 b	°C
	HCN boiling temperature	25.6 b	°C
cHCN	HCN heat capacity	2.612 × 103 b	J kg−1 K−1
kHCN	HCN heat conductivity	2.28 × 10−9 b	W m−1 K−1

 

Notes.

^a Assumes the value is the same, or calibrated for, sea ice (Buffo et al. 2018, 2020). ^b Properties of HCN at 1 atm and ∼20°C are from limited experimental data (Coates & Hartshorne 1931; MatWeb 2018; Rumble 2020).

Download table as:  ASCIITypeset image

In their basic forms, Equation (3) is an expression of heat conservation, and Equation (4) is an expression of HCN conservation. Equation (4) tracks the bulk change in the diffusion of the solute (where D_ice = 0), and the concentration of the solute by a phase change. The work of Buffo et al. (2018) was adapted to model the Europan ice and ocean boundary (Buffo et al. 2020). The key modification Buffo et al. (2020) made to their 2018 model was to ignore rejection of salt by volumetric changes that occur during freezing in the mushy layer. In terrestrial sea ice, volumetric changes were shown to play a negligible role in the rejection of salt (Buffo et al. 2020); the primary factor controlling salt removal is gravity drainage (Hunke et al. 2011).

Thus, we move forward with two partial differential equations that contain three unknowns (T, S, and ϕ). We therefore need another constraint to solve this system of equations. Buffo et al. (2018) used the enthalpy method (e.g., see Huber et al. 2008) to calculate the liquid fraction (ϕ) of each cell

Here, H is enthalpy, T_m is the melting temperature, c is specific heat capacity, and H_S is the enthalpy if the cell was completely frozen (enthalpy of solid ice). The model iterates Equations (5) and (6) with Equations (3) and (4) until the T, S, and ϕ converge before moving onto the next time step of the simulation.

As water freezes, the solute will concentrate, and a density difference will be created in the form of a concentration gradient. Buffo et al. (2018) described this using the approach of Griewank & Notz's (2013) 1D model of gravity drainage by finding the amount of vertical brine transported by convection (br_j )

Here, α is a constant of proportionality optimized by Griewank & Notz (2013), Ra_j is the Rayleigh number of the jth layer, Ra_c is the critical Rayleigh number, dz and dt are the spatial and temporal discretization sizes, respectively, g is the acceleration due to gravity, ρ_sw is the density of meltwater, β is a density coefficient describing the relationship between density and concentration, ΔS_j is the difference in concentration of the brine from ambient seawater, h_j is the height of the jth layer above the basal surface of the ice, κ is the thermal diffusivity of seawater, μ is the kinematic viscosity of seawater, and  is the minimum permeability of any layer between the jth layer and the basal ice surface. It is found using the permeability function given by Griewank & Notz (2013).

In essence, Equation (7) is a measure of how much an element exceeds the critical Rayleigh number, which triggers the brine transport. The localized differences in the Rayleigh number fuel convection, which is the dominant process in the desalinization of sea ice (Buffo et al. 2018, 2020). The brine transport changes the concentration (S) and temperature (T) that feed into the conservation equations of the mushy layer model. Note that we assume that constants optimized for sea ice by Griewank & Notz (2013) are adequate for our study, and that the thermal diffusivity and viscosity are approximately the same as those for seawater (Table 1). As water freezes, impurities are concentrated within the residual liquid fraction (pore space) of the two-phase interfacial layer. In terrestrial sea ice, this translates to pockets and veins of concentrated brines (Buffo et al. 2018).

The 1D SF2 model described above uses active interface tracking of the temperature and chemistry of the mushy layer in finite increments (Figure 2). It focuses on the active two-phase region at the water–ice interface (H = 1 m). As each increment (e.g., dz = 1 cm) in the active layer freezes (reaches a critical porosity that is estimated from terrestrial ice studies; Buffo et al. 2018), it is removed from the domain and cataloged while a new increment with ambient liquid water reservoir properties is added to the bottom of the domain. Once an increment of the mushy layer reaches critical porosity (ϕ ≤ 5%), the ice is assumed to be frozen enough such that all of the liquid brine, trapped in pockets and veins, is closed off from the underlying liquid. The impurities become fixed in the ice, and the 1 m mushy layer moves down to the next increment of the melt pond. While Titan's conditions reach more extreme temperatures than the terrestrial sea ice the model was initially designed for, the temperatures within the mushy layer are likely comparable with terrestrial ice (∼0°C).

 Zoom In Zoom Out Reset image size

In sea ice, the removal of salt from the forming ice layer is gravity driven and due to buoyancy differences between the denser interstitial brine and the underlying ocean (Buffo et al. 2018, 2020). Unlike salt, HCN is less dense than water and ice (Table 2; Haynes 2011). This relationship will invert the system such that impurities are buoyantly removed at the base of the melt pond, while at the top, gravity will not drive the removal of impurities (Chivers et al. 2021). Similar scenarios have been studied in terrestrial magma chambers, where low-density melt is concentrated in the porous roof magma–rock interface during solidification (Worster 1990). We propose that HCN in a Titan melt pond would follow the same pattern as less-dense impurities in terrestrial lavas because of its buoyancy within water and ice, with HCN becoming entrapped in the ice at the top freeze front. HCN's density relationship with water and ice suggests that as a melt pond freezes from all directions, the HCN dynamics will operate differently at the bottom and upper interfaces. A key feature of our impact melt system is that it is a closed finite body of water, similar to the proposed pockets of water thought to cause chaos terrain on Europa (Schmidt et al. 2011; Buffo et al. 2020; Chivers et al. 2021). SF2 assumes the mushy layer is atop an infinite ocean and that the rejected impurities will never significantly affect the bulk concentration of the ocean; however, the bulk concentration of a melt pond will become progressively more concentrated as it freezes. SF2 provides an estimate of impurities frozen into the ice for a range of temperature gradients (or depths) at a given concentration of HCN in the melt. Therefore, we produced multiple impurity entrainment versus thermal gradient relationships at a range of initial HCN concentrations in the melt to encapsulate the evolving conditions (i.e., thermal gradients and melt concentrations) of the melt pond as it freezes (Figure 3).

 Zoom In Zoom Out Reset image size

Table 2. Constitutive Equations Developed by Buffo et al. (2020) with the Corresponding Constants Used to Fit to Our Results from the SF2 Model

Constitutive Equations	a	b	c	d
	a1 = 0.1674	a1 = −27.98	a1 = −14.26
	a5 = 00.4654	a5 = −35.87	a5 = −6.127
	a10 = 0.6802	a10 = −43.35	a10 = −3.960
	a15 = 0.9876	a15 = −43.64	a15 = −2.535
	a20 = 1.120	a20 = −49.62	a20 = −2.242
	a25 = 1.328	a25 = −51.15	a25 = −1.710
	a50 = 2.346	a50 = −65.06	a50 = −0.7333
	a75 = 2.849 a	a75 = −76.75 a	a75 = 0.5926 a
	a100 = 3.348 a	a100 = −87.61 a	a100 = 0.7731 a
	a150 = 4.166 a	a150 = −107.1 a	a150 = 0.9729 a
	a200 = 4.838 a	a200 = −124.6 a	a200 = 1.086 a
	a250 = 5.417 a	a250 = −140.8 a	a250 = 1.343 a
	a1 = 0.3317	a1 = 0.0625
	a5 = 0.03900	a5 = 0.3587
	a10 = −0.1483	a10 = −0.1483
	a15 = −0.1847	a15 = 0.6002
	a20 = −0.2948	a20 = 0.6427
	a25 = −0.2805	a25 = 0.6387
	a50 = −0.2533 a	a50 = 0.6272 a
	a75 = −0.2533 a	a75 = 0.6272 a
	a100 = −0.2533 a	a100 = 0.6272 a
	a150 = −0.2533 a	a150 = 0.6272 a
	a200 = −0.2533 a	a200 = 0.6272 a
	a250 = −0.2533 a	a250 = 0.6272 a
	a1 = 1.982	a1 = −59.30	a1 = −2.381	a1 = 0.2337
	a5 = 5.942	a5 = −4.028	a5 = 1.885	a5 = 27.05
	a10 = 10.79	a10 = −7.847	a10 = 1.761	a10 = 38.29
	a15 = 15.77	a15 = −11.79	a15 = 1.766	a15 = 49.80
	a20 = 20.63	a20 = −15.61	a20 = 1.866	a20 = 58.65
	a25 = 26.60	a25 = −20.47	a25 = 1.936	a25 = 70.89
	a50 = 60.70	a50 = −49.17	a50 = 2.145	a50 = 103.4
	a75 = 76.24	a75 = −71.31	a75 = −0.3654	a75 = 65.06
	a100 = 167.2	a100 = −147.40	a100 = 1.566	a100 = 94.20
	a150 = 223.6 a	a150 = −193.8 a	a150 = 1.784 a	a150 = 96.89 a
	a200 = 301.4 a	a200 = −262.2 a	a200 = 1.771 a	a200 = 97.64 a
	a250 = 379.1 a	a250 = −300.7 a	a250 = 1.748 a	a250 = 98.09 a

 

Notes. These equations are fit to the data for 1–250 ppt.

^a Indicates values determined using extrapolation and averaging.

Download table as:  ASCIITypeset image

The 2D profiles of the impurity entrainment relationship are derived from the results of the 2D thermophysical model. Results are produced in increments of ≥5 m sections of the melt pond over a broad range of depths and thermal gradients sufficient to infer a complete relationship between the impurity entrainment and the depths or thermal gradients anticipated in Titan impact crater melt sheets. We solve the 2D heat transfer model by coupling it with the SF2 results. We fit impurity entrainment results to complete profiles using a set of constitutive equations produced by Buffo et al. (2020) to construct quasi-analytical concentration dependencies that can be utilized without explicit simulations (Table 2). We then reproduced the results for a range of initial concentrations of HCN in the melt. This accounts for the changing fluid concentration as the melt pond freezes, which is necessary to produce an impurity entrainment relationship for Titan crater melt ponds as opposed to an infinite ocean. The 2D heat transfer model uses the impurity entrapment–thermal gradient profiles, interpolating between them using a linear least-squares method to create a reference to track how much impurity is entrained into the ice in a closed system like a melt pond. This reference is crucial in tracking how much impurity is entrained in the ice because the thermal gradient at the freeze front and the impurity concentration in the melt change as the melt pond freezes, and both variables are inextricably linked to how much impurity is entrapped in the ice (Figure 3).

As previously outlined, we studied a range of concentrations from 0.1%–25%. However, the model was unable to converge at higher concentrations (>7.5%) and at low thermal gradients (<10 K m⁻¹) for moderate concentrations (5%–7.5%). Without these results, we could not construct full profiles using the constitutive equations (Table 2); these are necessary in tracking a 2D closed system where a range of concentrations will exist. Therefore, we extrapolated from the results at lower concentrations (≤7.5%).

There are three areas where extrapolation is performed: the linear thermal gradient domain, the shallow thermal gradient domain, and the entire depth domain (Figure 4). First, the linear domain is clearly distinguishable in the results for the lower initial concentrations (Figure 4(a)). Because the modeled concentrations converge at these lower thermal gradients, we assume the larger initial concentration profiles converge onto this line too. Therefore, their constants in the constitutive equation are taken to be the average of the constants derived from the lower concentrations. This extrapolation was done for all concentrations ≥5% (Figure 4(a)). Second, the shallow domain was best constrained using a linear extrapolation, done for all concentrations ≥15%. The depth relationship is defined by a single constitutive equation, and it was best constrained using a power-law extrapolation, where the extrapolation is performed for initial concentrations ≥7.5% (Figure 4(b)).

 Zoom In Zoom Out Reset image size

We couple these results using the constitutive equations (Table 2) with a 2D heat transfer model in a stand-alone code built by Chivers et al. (2021) to simulate a simple closed system melt sheet. The Chivers et al. (2021) model is a finite difference model representing a rectangular (cylindrical) melt sheet that tracks the transfer of heat, phase change, and chemistry as the melt sheet freezes. As each element freezes, the model references the SF2 results that correspond to the liquid melt concentration and the thermal gradient of the cell at the time of freezing. The Chivers et al. (2021) model uses this to account for how much HCN is frozen into the ice at the bottom ice–liquid interface, and the rejected HCN increases the liquid melt concentration. As each cell freezes at the upper ice–liquid interface, the model assumes that the amount of HCN that will be trapped in the ice is the same as the liquid melt concentration at the point of freezing. Unlike the SF2 model, the Chivers model uses a temperature-dependent ice conductivity to account for the large-scale variations of ice temperatures that exist across the entire model domain, which extends deep beneath the melt sheet into the ice shell.

We model a closed 2D melt pond of an 80 km diameter crater (∼55 km floor; Bray et al. 2012) mixed with HCN using the same initial conditions as in the SF2 model (Figure 5). We ran the model for melt sheets of 10–100 m thickness and an initial concentration of HCN ranging from 0.1% to 7.5% by weight HCN, where the highest concentrations in the melt will reach 25% as the melt freezes. The thicker 250 m melt sheet was modeled at just 0.1% because larger initial values caused the melt concentration to exceed 25% as the melt lens shrinks; this is the maximum concentration that we studied with the SF2 model. The heat transfer model begins with 1–10 m of ice frozen at the roof of the melt pond. We do this to avoid the extreme conditions encountered in this region of the melt sheet, which are likely to be chaotic and difficult to model. This initiates the model with a thin layer of ice insulating the melt sheet, assuming that once this ice layer exists, the model will proceed to freeze via conduction.

 Zoom In Zoom Out Reset image size

We have created a complete profile of how much HCN will freeze into a solidified water ice melt pond for an array of initial concentrations between 1 ppt (0.1%) and 250 ppt (25%) (Figure 4). We fit our results to the constitutive equations of Buffo et al. (2020; Table 2). These constitutive equations provide a complete profile of bulk HCN in the ice across the melt pond depth for variable temperature gradients (at the instant of freezing). Note that our fits are not perfect; they become less precise at higher thermal gradients and concentrations. The observed imprecision in the fits is not related to the model's physical assumptions, but rather its ability to be fit to the constitutive equations at the extremes found on other worlds. Furthermore, these extremes are occurring at the boundaries of our model (i.e., the highest concentrations occur near the end of the melt pond's evolution, and the highest thermal gradients occur at the very beginning of the melt pond's evolution). The intermediate stages are likely to be the most reliable because the thermal gradients are far less extreme and more akin to those in terrestrial environments. Fortunately, this is the region of greatest interest to Dragonfly because any biomolecules trapped in the upper portions of the melt sheet would have had very little time to evolve into something chemically and astrobiologically interesting. Our fits are also less precise for the domain where the amount of HCN that is trapped in the ice forms a linear relationship with the thermal gradient (Table 2). The transition into the linear domain appears to occur at progressively higher thermal gradients for higher initial concentrations. The results at higher concentrations (>7.5%) are largely extrapolated, as detailed in Section 2.3, but these profiles are mostly relevant in the deep interior of the melt sheet as it freezes. The point at which the melt sheet reaches said concentrations will depend on the initial concentration and the melt thickness (Figures 5 and 6). The primary purpose of the extrapolated profiles is to provide an estimate for these large initial concentrations, which only exist in a small region of the melt sheet, but the results in regions of lower concentration are more reliable.

 Zoom In Zoom Out Reset image size

Our results from the SF2 model exhibit clear trends. As expected, higher bulk concentrations will lead to a higher concentration of HCN trapped in the ice. The magnitude of this difference is largest at the highest thermal gradients. We observe that higher HCN concentrations in the liquid melt result in a higher concentration in the ice, but all of the concentration profiles soon begin to converge to ∼1 ppt and below at the lower thermal gradients (more ice insulating the melt). With HCN as the impurity, the SF2 model tracks the HCN entrainment of the freeze front at the base of the crater, and, as stated in Section 2.3, buoyancy is assumed to hold the HCN in place at the top interface of the melt pond, resulting in ice compositions that reflect that of the underlying melt.

The freeze fronts create insulating ice layers that grow, slowing the cooling process. However, the upper ice layer will not insulate as effectively as the ice at the base of the melt pond, causing the upper boundary of the melt pond to freeze faster than the lower boundary. Furthermore, the base crustal ice is likely to be warmer than the background temperature (94 K) because of the residual heat left over from the impact, but we assumed the basal ice was the same temperature as the background for this initial study. Consequently, the melt pond's final, most concentrated fluid solidifies just below the geometric center of the melt pond (as previous studies have suggested, e.g., Thompson & Sagan 1992; Chivers et al. 2021), offset by ∼25% of the melt sheet thickness in every run (Figures 5; 6). While this region exhibits the highest concentration of HCN, ∼70% of the total HCN content is still stored in the top half of the melt sheet. The position of this peak will shift depending on how much the basal ice is heated from the impact process, but it will also shift as the freezing temperature changes (i.e., gets lower at higher concentrations). The effect of melting temperature is observed at the highest initial concentrations of HCN; melt sheets of 50 and 100 m thickness peak slightly deeper for initial concentrations of 5% and 7.5%. In these cases, the melt sheet is sufficiently concentrated such that the lower melting temperature prolongs the freezing process.

HCN preferentially freezes into the ice in the upper boundary, and the concentration of HCN increases with depth as HCN is removed from the lower boundary. Conversely, the bottom boundary sees a small spike in HCN trapped in the ice before its concentration drops significantly. That is because the initial conditions assume the underlying ice is 94 K, but realistically, the thermal gradient would be much lower because of residual heat from the impact. Therefore, our results are a lower bound for the lifetime of the melt pond because warmer ice will insulate the melt pond prolonging the time it takes to freeze. Furthermore, higher thermal gradients will lead to more HCN freezing into the melt pond's basal ice, so the upper profile of HCN concentrated in the ice is a lower bound for the amount of HCN Dragonfly is likely to find in the ice. However, the increase in HCN in the basal ice that results from higher thermal gradients is negligible relative to the concentrations observed near the surface.

The lower boundary faces two conflicting environmental pressures: (1) a drop in the thermal gradient will decrease how much HCN brine freezes into the ice, but (2) the increasing bulk concentration will increase how much HCN brine freezes into the ice. Our SF2 results (Figure 4) suggest that an increase in the initial concentration of HCN would have a significant effect on how much HCN freezes into the ice for the first several meters of the melt sheet or up to ∼5 K m⁻¹ thermal gradients. At thermal gradients below ∼5 K m⁻¹, the initial concentration appears to matter much less, as the concentration in the ice across different initial conditions begins to converge. This trend is observed in the results of the Chivers et al. (2021) model as well. The higher initial concentrations of HCN produce a slightly higher concentration of HCN trapped in the ice at the bottom of the melt pond (where the SF2 model is applied). But the more the melt sheet freezes, the smaller the thermal gradients become such that less HCN is trapped in the ice. In the innermost region of the melt sheet, we see that a higher initial concentration of HCN in the melt produces profiles of HCN concentration in the ice that overlap (and drop below) profiles of lower initial HCN concentration in the melt (Figure 6), but this is not representative of reality. It is unlikely that higher impurity concentrations in the melt would entrain less impurity in the ice.

These unrealistic results are derived from how we construct and apply the trend lines for each concentration profile. Each profile converges upon the same linear concentration line (Figure 4(a)) at low thermal gradients. The higher the initial concentration, the sooner the profile reaches that linear region. The linear fits we use are only applicable for a specific region (up to ∼15 m) beyond which they breakdown. These fits allow us to track the concentration until it reaches the critical minimum concentration, and while we see the concentration profiles dip below our lowest initial concentration line (1 ppt), it is safe to assume the fit is only applicable until it reaches the lowest concentration line where it would converge with the profile.

For example, the lowest initial melt concentration profile that we study is 1 ppt (0.1%) concentration. In our thickest melt sheet (100 m), this concentration profile decreases at a much slower rate than the higher concentrations, reaching, at its lowest point, about 30% of the initial concentration. In this region of very low thermal gradients, the highest initial concentrations begin to produce concentrations in the ice as low as melt with an initial concentration of 1 ppt; this is the point at which the higher concentration linear fits cease to be applicable. This suggests that at very low concentrations, very little material will be rejected, even with the lowest of thermal gradients. Logically, this is to be expected, as lower concentrations lead to smaller density differences and thus smaller forces driving the movement of HCN within the forming ice lattice. Therefore, HCN may freeze at one-third its initial concentration (0.25%–2.5%; Section 2.1); if HCN undergoes hydrolysis into more complex biomolecules, these daughter products may freeze near their initial concentration given that they would be a fraction as abundant as HCN (at least for melt sheets <50 m). However, it should be noted that volumetric effects become a bigger factor as concentrations have an increasingly negligible effect on the density of the melt.

In a bulk scenario, where the total concentration of organics is much larger (∼7%; Section 2.1), more material will be rejected if the organics are assumed to act as a whole (rather than segregate). However, the lowest concentrations in the ice still seem to be ∼0.03% (0.3 ppt), well within the parts per trillion that Dragonfly will need for its investigation.

The results in Figure 6 are for the Chivers 2D heat transfer model, which accounts for lateral heat transfer. However, the x-axis (radius of the melt pond, ∼55 km) is orders of magnitude larger than the y-axis (depth of the melt pond, 50–250 m; see Chivers et al. 2021 for more on heat loss in the x-direction). Therefore, the horizontal freezing is not overly apparent in Figure 6, but if you look at the ice fraction column, you can see the corners begin to curve in. Our measured freezing timescale of tens to hundreds of years is consistent with previous timescale estimates for the freezing of a melt sheet at the thicknesses we modeled (Figures 5 and 7; Artemieva & Lunine 2005; O'Brien et al. 2005). However, our largest freezing timescales are much smaller than the most liberal estimates of O'Brien et al. (2005); they consider larger depths (up to d/D = 0.05) and volume fractions (up to 20%) than we deem reasonable. Their estimates would lead to depths multiple kilometers deep for the largest craters, but there is evidence that the depths of large craters on icy worlds reach a limit where relaxation produces shallower craters (Schenk 2002). In addition, they neglect the significant drainage into the underlying fractured bedrock ice that will likely occur, thinning the melt sheets (Elder et al. 2012). With a conservative depth to diameter ratio (d/D = 0.01, as opposed to the slightly higher estimate we used in Section 2.1 from Bray et al. 2012) and volume fraction (f = 2%–5%), a crater the size of Selk would produce a melt pond 60 m in depth (excluding any drainage), resulting in a freezing timescales of ∼10 yr. If we are more liberal with the volume fraction (f ≅ 15%), a crater the size of Selk would produce a melt pond ∼350 m in depth and freezing timescales of up to ∼600 yr. However, we focused on melt sheets no thicker than 250 m, as thicker melt ponds are unlikely to persist given the effects of drainage (Figure 7; Elder et al. 2012). The timescales of each melt sheet thickness become progressively longer but retain the same pattern of freezing when presented with the same initial conditions.

 Zoom In Zoom Out Reset image size

Our results suggest that Dragonfly should be able to access HCN at the top of the melt sheet in concentrations similar to those present at the point of its formation (assuming the melt pond is thoroughly mixed). Of course, these results do not account for the hydrolysis of HCN, but rather use it as a representative tracer for its byproducts. The byproducts of HCN aqueous chemistry may freeze into the ice in a different manner; our results suggest density would dictate the path they follow. If a particular byproduct of HCN has an inverse density relationship, the ratio of HCN to its byproduct may be lower than what is expected through aqueous reactions because the denser byproduct will freeze into the ice at lower concentrations than HCN. Based on the half-lives discussed in Section 2.1, it is unlikely that HCN's daughter products will dominate a melt sheet less than ∼50 m thick, and potentially as thick as 100 m depending on the pH, because the half-life of HCN in the best case is on the order of hundreds of years (Figure 7; Miyakawa et al. 2002). HCN may continue to segregate as we have predicted, while the daughter products act independently in the mushy layer freezing process. In this case, the trends we present are still accurate, but the magnitude of the concentration in the ice would be lower than we predict. Alternatively, such a mixture of compounds may not segregate at all, but instead be governed by its bulk properties, with much broader implications for the distribution of impurities within the melt sheet.

As HCN hydrolyzes and becomes less concentrated, its chemical properties would become decreasingly relevant to the physics of the SF2 model. Once the melt sheet has remained liquid beyond HCN's half-life, the chemical properties of its daughter products will dominate the physics. Therefore, our results are the least reliable beyond this point (perhaps as quickly as a few years and almost certainly after centuries). The most drastic difference would be if its daughter products (e.g., glycine) are denser than HCN and water, thus inverting the forces that drive how much, and where, impurity freezes into the ice. In such a scenario, the current pattern of organic distribution with a large segregation of organics in the top half and orders of magnitude less concentration of organics in the bottom half would still be present. But once HCN's daughter product becomes dominant (closer to the melt sheet's center), the same pattern would be reproduced but mirrored (much like salt; Chivers et al. 2021). Eventually, the two freeze fronts would meet, and the concentration would spike at the center. Realistically, this transition could happen as shallow as 10 m below the surface (with the fastest hydrolysis) or to as deep as 100 m into the surface (for slower rates). This is within the region Dragonfly may be able to access. Therefore, the team should be prepared for this potential reality.

An additional complication may arise as the HCN hydrolyzes into potentially more dense daughter products. Going from HCN's buoyant density to a denser daughter product might introduce a period where the bulk melt density is near that of water. In this case, buoyancy-driven rejection of impurities may become less important than volumetric phase-change-driven expulsion. These are ignored in our work because it has been shown that volumetric changes are negligible to those driven by density differences (Griewank & Notz 2013; Buffo et al. 2018, 2020). But if buoyancy differences do not exist, volumetric changes may become more important in the diffusion of impurities out of the ice.

There are further limitations to this work that must be considered as well. It is difficult to make precise approximations of the near-surface concentration of chemical impurities in the ice because of the chaotic nature of the freezing process in this region. Naturally formed methane rivers on Titan could incise deeper into the melt sheet, where the freezing conditions were less extreme, and the SF2 model would provide a more accurate numerical solution. Selk crater was chosen, in part, because of its relative freshness (Neish et al. 2018), but there is evidence to suggest erosion plays an ongoing role in the degradation of Titan's impact craters (Soderblom et al. 2010; Neish et al. 2013, 2016; Werynski et al. 2019; Hedgepeth et al. 2020; Solomonidou et al. 2020). Erosion will help to expose subsurface material, creating easily accessible samples for Dragonfly to target (Neish et al. 2018). However, sediment covering the melt sheet may still impede Dragonfly's access to these samples (Neish et al. 2013). If rivers do not incise in the frozen melt pond, or incise very little, the aqueous chemistry found near the surface of the melt sheet could only have evolved on very short timescales (minutes to hours), limiting its potential to evolve into more complex biochemistry (Neish et al. 2006, 2018).

Overall, we predict that the best samples for Dragonfly to target would be near the center depth of the melt sheet. Ignoring other chemistry, the region of maximum HCN concentration does not appear to shift based on different initial HCN concentrations or melt sheet thicknesses (although the absolute depth of this region will vary with the initial thickness of the melt sheet, which is unknown). The easiest location to access samples at would be in any available river plains, where portions of the melt sheet may have eroded out onto the valley floor (Neish et al. 2018). Our results suggest that even in low HCN concentrations (0.1%), there is likely to be plenty of HCN for Dragonfly to sample. However, the dynamics that govern the emplacement will likely complicate as HCN undergoes hydrolysis during freezing, potentially changing what may be observed. However, our results suggest even in the lowest of concentrations, there is likely to be plenty of HCN in the ice to sample with Dragonfly's instruments. Whether this is true for its daughter products requires further study.

This work offers a suggestion of how much organic material Dragonfly will likely find in Selk crater, but there are several assumptions that are associated with this work. The most important assumption is that the planetary ice model SF2 is applicable for alternate chemistries and environments other than those of terrestrial sea ice. Overall, we deem this to be a good assumption because the fundamental physics governing ice-brine environments will not change across the solar system. However, the extreme conditions found at Titan do make it more difficult to find a numerical solution, but only for a thin layer at the surface. A necessary part of adapting this model for Titan is the assumption that some of the unique parameters taken from sea ice research are the same for Titan (Table 1; Golden et al. 2007; Maus et al. 2020). Even in the terrestrial sea ice community, some of these parameters (e.g., permeability) are still heavily debated (Golden et al. 2007).

Furthermore, we assume our system is analogous to terrestrial sea ice, but as previously stated, the temperatures drop below those found in terrestrial sea ice. This would affect the SF2 model in its initial stages, but that is only because we assume the basal ice begins at background temperatures (i.e., 94 K). Once the temperatures adjust from the extreme initial conditions, the temperatures within the mushy layer will not vary nearly as significantly, and the assumption of a constant conductivity of ice will have a much smaller impact. The 2D heat transfer model accounts for temperature-dependent conductivities (Chivers et al. 2021), but our initial parameters still assume the basal ice is at the background temperature (94 K). Figures 6 and 7 show that the top boundary freezes at a rate almost twice that of the base of the melt sheet. This suggests that the timescales may extend another decade or more from what we predict, depending on the level of heating the ice experiences from impact. It also suggests the impurities may concentrate farther down into the melt sheet. This is because the upper freeze front will continue at the same rate that we observed, but the warmer basal ice will slow the progression of the lower freeze front. The impact is also likely to lead to initial temperatures well above the freezing point of water. This would produce higher thermal gradients (and therefore, higher impurity entrainment) and longer timescales over which the freezing will occur. Overall, we anticipate that the temperature assumptions we made in this work represent a lower limit on the actual freezing timescales encountered on Titan.

Even though the conditions in a Titan melt sheet are not extreme compared to those of terrestrial sea ice for most of the model, the SF2 model still faced limits on how well it could model the scenarios we gave it. This does not seem to be a physical problem, however. The limitations we observed appear to be specific to HCN's chemical properties, and as the initial concentration is increased, the problems are exacerbated. The combination of HCN's lower melting temperature (∼−13°C) with its specific heat (2612 J kg⁻¹ K⁻¹) at the extremely low thermal gradients within the melt sheet makes it progressively more difficult for the model to converge. This limited our SF2 model runs to relatively high thermal gradients and relatively low concentrations. It should be noted that most of the concentrations of interest (0.1%–7%) produced results that converged even at low thermal gradients. The extrapolation to higher concentrations was needed for the Chivers et al. (2021) model to simulate the melt when it reaches these concentrations. While we do not place strong weight onto the extrapolations, they do serve the purpose of producing a more cohesive picture. Further, while the regions where the melt concentrations rely on extrapolated SF2 results may be the least accurate, these are also the regions that are of least interest. That is because they are less accessible to Dragonfly, and if they are accessed, we can at least say that there will be ample organic samples for Dragonfly to analyze. Outside of this region, the extrapolation is less important.

Our work also assumes a well-mixed pond. This may be true during its energetic origin, but the initial energy of the system will dissipate with time. It is possible that large-scale convection may continue to mix the melt pond, as observed in terrestrial magma chambers (Sparks et al. 1984; Turner & Campbell 1986), and the model of Chivers et al. (2021) assumes this type of convection will be enough to keep the melt sheet well-mixed.

Another assumption is that HCN will be present and is not lost during the initial stages of impact cratering, where temperatures may exceed the vaporization point of water (373 K). This is above the vaporization point of HCN (299 K). The SF2 model cannot account for gases in the melt, and this problem is not unique to HCN. This points to the larger problem of the modification and destruction of organics that may occur during impact (Artemieva & Lunine 2005). The destruction or loss of HCN will largely affect the initial concentration of HCN in the melt pond, and the high pressures and temperatures experienced during impact may alter the chemistry as well (Pierazzo & Chyba 1999; Blank et al. 2001). Therefore, the estimates we give for HCN are likely an upper bound, and HCN's half-life may be significantly smaller for the much warmer temperature conditions. Therefore, the daughter products may become a bigger factor in the melt sheet sooner than we predicted.

We have assumed that HCN will act as a proxy for the bulk mixture of material, and we have discussed significantly how the presence of its daughter products (and other chemistry) would complicate this process. This relates not just to the traditional mushy layer problem being solved in this work but to other chemical processes we do not consider. For example, HCN will begin to dissociate when in high pH conditions, which is also when rates of hydrolysis are known to be especially high (Miyakawa et al. 2002; Gail et al. 2011). The final pH of the melt solution will depend on the other compounds present in the melt, and it will determine the stability and rate of hydrolysis of HCN, as well as for the other molecules that are present. Therefore, when we discuss the distribution of HCN, it is more accurately thought of as a first-order approximation of the distribution of organic molecules within the melt sheet. Future studies are needed to determine the full breadth of factors that will influence the final distribution of biomolecules in Titan's frozen melt ponds.