Modeling the Distribution of Organic Carbon and Nitrogen in Impact Crater Melt on Titan

Titan is a chemically rich world that provides a natural laboratory for the study of the origin of life. Titan’s atmospherically derived C x H y N z molecules have been shown to form amino acids when mixed with liquid water, but the transition from prebiotic chemistry to the origin of life is not well understood. Investigating this prebiotic environment on Titan is one of the primary motivations behind NASA’s Dragonfly mission. One of its objectives is to visit the 80 km diameter Selk crater, where a melt sheet of liquid water would have formed during the impact cratering process. Organic molecules on Titan’s surface could have mixed with this water, forming molecules of prebiotic interest. Constraining how this material becomes trapped in the refreezing ice is necessary for Dragonfly to effectively target and interpret the samples it aims to acquire. In this work, we adapt the planetary ice model of Buffo et al. to Titan conditions to track how organic molecules will become trapped within the ice of the freezing melt sheet. We use HCN as a model impurity because of its abundance on Titan and its propensity to form amino acids in aqueous solutions. We show that without hydrolysis, HCN will be concentrated in the upper and middle portions of the resolidified melt sheet. In a closed system like Selk crater, the highest concentration of HCN appears 75% of the way into the frozen melt pond (relative to the surface), but HCN should be accessible at high concentrations nearer the surface as well.


Introduction
Titan is one of the most chemically rich places in the solar system. It houses a dense nitrogen-rich atmosphere with small amounts of methane (Lindal et al. 1983;Tomasko et al. 2005). This introduces a unique complexity to Titan's chemistry and geology (Hörst 2017;Malaska et al. 2020). Solar and cosmic radiation constantly dissociate the molecules in Titan's atmosphere, and, over time, the carbon, hydrogen, and nitrogen recombine to form a range of complex organic compounds referred to as tholins (Sagan & Khare 1979;Sagan & Thompson 1984;Krasnopolsky 2014). These tholins are the source of Titan's haze, layers of large molecules that scatter visible light and prevent it from reaching the surface (Tomasko et al. 2005;Hörst 2017). Any chemical arrangement that can form likely does, producing molecules with the general formula C x H y N z (Khare et al. 1986). Indeed, Titan possesses an abundance of organic molecules, the most prevalent of these being ethane (C 2 H 6 ), acetylene (C 2 H 2 ), propane (C 3 H 8 ), ethylene (C 2 H 4 ), and hydrogen cyanide (HCN; Hörst 2017). Furthermore, Titan's surface pressure is very similar to Earth's (1.5 bar versus 1 bar), and with temperatures of 94 K, methane is stable as both a gas and a liquid at the surface (Kouvaris & Flasar 1991). This produces a methane hydrological cycle that contributes to Titan's surface degradation (Neish et al. 2013(Neish et al. , 2016(Neish et al. , 2018Hedgepeth et al. 2020).
It is possible that Earth's early atmosphere was similar to Titan's present atmosphere (Trainer et al. 2006;He & Smith 2014). In 1953, the Miller-Urey experiment showed how conditions similar to Titan (an atmosphere rich in reducing molecules) could produce important amino acids relevant to life (Miller 1953). Later work demonstrated that similar amino acids could be produced in simulated Titan conditions (McDonald et al. 1994;Neish et al. 2010;Cleaves et al. 2014). Thus, understanding the processes occurring in Titan's atmosphere and on its surface may help us understand how life on Earth could have transitioned from abiotic chemistry to biology. Large quantities of organic compounds exist on Titan's surface in the form of sand dunes and methane lakes (Lorenz et al. 2008); the dunes alone cover ∼40% of the equatorial regions Lorenz & Zimbelman 2014). While the exact composition of the sand dunes is not well constrained, the consensus is that they are sourced from Titan's atmosphere and are organic in nature (Soderblom et al. 2007; Barnes et al. 2015;Janssen et al. 2016). These organics make Titan a natural laboratory for prebiotic chemistry. Neish et al. (2010) showed that Titan tholins produce amino acids when exposed to liquid water at room temperature. Titan has an abundance of water in the form of ice, but frigid surface temperatures (∼94 K) would impede any longterm life from taking hold (Clarke 2014). However, there exist a number of geophysical processes that may permit opportunities for Titan's tholins to mix with liquid water (Thompson & Sagan 1992;Neish et al. 2018).
Impact craters on Titan are potentially the most suitable geologic features to facilitate prebiotic chemical reactions of tholins in liquid water (Neish et al. 2018;Crósta et al. 2021). Impact craters produce massive amounts of energy and melt the bedrock ice at the site of the impact (Osinski et al. 2011). On Titan, this process would produce warm liquid water (>0°C), and depending on the size of the crater, could take hundreds to thousands of years to refreeze, far exceeding the timescales that can be examined in a laboratory (Thompson & Sagan 1992;Artemieva & Lunine 2005;O'Brien et al. 2005;Neish et al. 2018). Other mechanisms like cryovolcanos may contain water magma chambers that are sufficiently large to remain liquid over similar timescales, but the lava flows that reach the surface (where the organics are located) would be much shorter lived (Neish et al. 2006). Furthermore, cryovolcanic processes on Titan are likely to be rare, and the cold, ammonia-laced fluid erupted from depth is less conducive to the fast chemical reactions needed to produce biological molecules before the water source freezes (Thompson & Sagan 1992;Lorenz & Lunine 1996;Lopes et al. 2007;Neish et al. 2008Neish et al. , 2010Neish et al. , 2018. Therefore, impact craters are likely the most important astrobiological target on Titan, as advanced prebiotic chemistry could have occurred in their temporarily liquid melt sheets. Titan's biological potential is the primary motivation behind the newly selected Dragonfly mission . Dragonfly is a rotorcraft lander selected to be NASA's fourth New Frontiers mission. Dragonfly will directly sample and analyze the chemistry frozen in Titan's surface ice as well as its atmosphere, including at Selk crater (Lorenz et al. , 2021. Selk crater was chosen as a target because it is one of the larger craters on Titan, offering a maximal time span for chemistry to evolve in its putatively thicker melt deposit (Neish et al. 2018). It is unclear how organic impurities might become trapped in the residual melt and overlying ice as it freezes; however, longer timescales are likely the most conducive to the chemical evolution of life (Wald 1954). This suggests that the most astrobiologically interesting samples will be found in the deep interior of the melt layer, where the water freezes last.
There are 90 suspected impact craters identified on Titan, many of which exhibit high levels of degradation (Hedgepeth et al. 2020). Of particular importance for the present study is naturally forming methane rivers, which are known to incise into Titan craters (Neish et al. 2016(Neish et al. , 2018Hedgepeth et al. 2020). These rivers will incise into the frozen melt sheet, depositing ice (and the material trapped in it) in alluvial fans and exposing the interior of the melt deposits along the walls of the riverbed. The alluvial fans and exposed riverbed walls are two direct ways of accessing potentially prebiotic molecules with a landed asset on Titan (Neish et al. 2018). In addition to its size, Selk's depth to diameter ratio indicates that it is relatively fresh and is therefore less likely to be buried in a large amount of aeolian or fluvial deposited sediment (Neish et al. 2018;Hedgepeth et al. 2020). However, its relative depth suggests it is still more than 40% degraded relative to a fresh crater, and its rim heights are generally lower than expected (but within error; Hedgepeth et al. 2020). This is consistent with some fluvial erosion, but it should be noted, no fluvial features are discernible in RADAR or Imaging Science Subsystem data, given their limited resolutions. However, evidence of fluvial erosion has been proposed in the Visible and Infrared Mapping Spectrometer (VIMS) data (Soderblom et al. 2010).
In this work, we seek to determine where the products of aqueous prebiotic chemistry are most likely to be located within impact melt sheets on Titan. We use terrestrial sea ice as an endmember analog for the freezing of organic molecules within impact craters on Titan. Terrestrial studies have investigated sea ice biogeochemistry, but these results are inherently skewed toward biologic processes so are unlikely to be a good predictor for organics in Titan melt ponds (Collins et al. 2008;Santibáñez et al. 2019). Terrestrial sea ice naturally rejects most impurities as it freezes (Hobbs 1974), but the physics that control how much is rejected are dependent on the thermal gradient at the solidification interface, as well as the concentration, and chemistry of the brine mixture (Hunke et al. 2011;Buffo et al. 2018Buffo et al. , 2020. Therefore, it cannot be assumed that organic molecules on Titan will follow the same pattern as terrestrial impurities, because Titan has a different chemistry (organics versus salt), and the ambient temperature is much lower (94 K versus 273 K), leading to significantly different density and thermal gradient extremes. However, once a thin layer (∼10 m) of ice covering the melt sheet forms, the thermal gradients in impact melt sheets on Titan become comparable to those in terrestrial sea ice. For Dragonfly to adequately sample the frozen melt for prebiotic (or possibly biologic) products, the team must understand the geologic history of the water-organic mixture at Selk crater. Here, we constrain how much HCN will become trapped in the ice versus rejected into the remaining liquid reservoir under the conditions predicted for fresh impact crater melts on Titan. We use the planetary ice model SlushFund2.0-v2 (SF2; Buffo et al. 2020) and a heat transfer model (Chivers et al. 2021) to track the concentrations of molecular impurities, both within the ice and the residual melt, as the melt pond refreezes, providing improved estimates of organic material distribution in Titan environments for future planetary exploration.

Initial Conditions
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 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.
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 bowlshaped 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 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. 2008Neish et al. , 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%).

1D Ice Model
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(Worster , 1997Feltham 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, f 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 (f) (e.g., a hypothetical value y would be averaged, y y y 1 br ice . All variables and their values can be found in Table 1. 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 f). 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 (f) 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 f  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 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 j P 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(Buffo et al. , 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).

Application to Titan
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 (f 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).
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(Buffo et al. , 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 Figure 2. The SF2 model tracks the mushy water-ice interface (red box) within the finite region of space H (m) where each increment (dz) represents a volume average over some finite length (m). When the uppermost increment of the water-ice interface reaches a critical porosity (here f 0.05), the water-ice interface (H) shifts downward by one spatial increment (dz), and the concentration of the frozen increment (f 0.05) is assumed to be trapped in the ice and remains constant throughout the remainder of the simulation. This process continues throughout the depth of the pond h (m), beginning at some interface depth d (m). The model tracks the freeze front where buoyancy will drive impurities out of the ice. For salt, the model tracks the chamber roof's freezing front as imaged here, but for HCN, the model is tracking the floor's upward propagating freeze front. In our work, the model runs bottom-up because HCN is less dense than water and ice, resulting in upward gravity drainage (see the text for more details).
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).
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 Table 2 Constitutive Equations Developed by Buffo et al. (2020) with the Corresponding Constants Used to Fit to Our Results from the SF2 Model (<10 K m −1 ) 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%). (Bottom) We visualize four temporal steps of the ∼273 K water melt (blue) freezing into ice (white). As time (t n ) passes, the warmer ice insulates the melt, reducing the thermal gradient at the freeze front. Simultaneously, the concentration of HCN in the melt (C n ) is increasing.
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)).
We couple these results using the constitutive equations ( Table 2)   . Concentration of HCN (ppt) in the frozen melt pond, for initial concentrations of 1 ppt (0.1%, yellow) to 250 ppt (25%, dark blue), at different thermal gradients (a), and depths (b). The SF2 model results are shown with black points (they appear as disjointed lines because the spatial resolution is so fine). The analytical fits from Buffo et al. (2020) are applied to our data, where solid lines are fits formed using SF2 results, and dashed lines are formed by extrapolating from the lower concentration fits (see Table 2 and the text). 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 iceliquid 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.

1D Concentrations
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.
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.

2D Distributions
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 −1 thermal gradients. At thermal gradients below ∼5 K m −1 , 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 Figure 6. We use the Chivers et al. (2021) heat transfer model to calculate the final concentration distribution for an 80 km diameter crater (crater floor ∼55 km) using an initial concentration of 1 ppt (0.1%) and an assumed melt pond of 100 m depth. On the left is the concentration of HCN in each element, in the middle is the temperature of each element, and on the right is the phase of each element (0 is frozen, 1 is liquid). The model uses the thermal gradient on the interface boundary along with the constitutive equation of Table 2 to determine what the concentration will be upon freezing (see the text for further explanation). 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 xaxis (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.

Implications and Future Work
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. 2018Buffo et al. , 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. 2013Neish et al. , 2016Werynski 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(Neish et al. , 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.

Assumptions
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 −1 K −1 ) 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 firstorder 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.

Conclusions
We use the planetary ice model of Buffo et al. (2020) to constrain the spatial distribution of HCN (and by extension, other organic compounds) trapped within a freezing crater melt pond on Titan. We found that HCN will become entrained in the upper boundary of the melt sheet at the initial concentration of the melt reservoir, and HCN will be removed from the lower boundary due to density differences with water. As the system freezes, HCN will concentrate in the remaining fluid of the melt sheet, becoming most concentrated ∼75% into the melt sheet (relative to the surface). This portion of the melt sheet is likely to contain the most evolved aqueous chemistry, because it remains liquid for the longest time, but it is also the region where HCN becomes a less accurate representation of how impurities will act. Our results offer constraints for Dragonfly to use to predict how many organic molecules it will find frozen into the melt sheet at any depth. Furthermore, our results can be utilized to link near-surface sample observations to deeper reservoir properties to better understand the history of the melt sheet and the initial conditions when it formed. Finally, the freezing timescales we found will aid in our understanding of the extent to which prebiotic chemistry could have evolved within the freezing melt sheet. This is critical for understanding the astrobiological potential of Titan and provides necessary context to the samples Dragonfly will collect.