Investigating the Stability and Distribution of Surface Ice in Mercury's Northernmost Craters

During an S-band (12.6 cm, 2380 MHz) radar observing run, the Arecibo Observatory would transmit a modulated, circularly polarized beam of light and receive echoes in both the same circular (SC) and opposite circular (OC) polarization as transmitted. Resulting data products included SC and OC delay-Doppler images, which are maps of radar backscatter intensity. In the case of Mercury, Harmon et al. (2011) also produced total power (i.e., SC+OC) images via a weighted sum of each polarization to account for differences in noise levels. These images were then projected onto hermiocentric coordinates with a resolution of 1.5 km. Higher-resolution maps (∼300 m) were also produced by Harmon et al. (2011) but without details of the process used to create the higher-resolution maps. Here we use the higher-resolution radar maps for visualization but use the base radar resolution maps for data analysis. The final map products were averaged over a multiyear observing campaign (1999–2005) encompassing over 20 nights of observations, which included subradar longitudes spanning 11°W–349°W, and thus offered a full 360° view of the northern polar deposits (Harmon et al. 2001, 2011). This helps to overcome some radar shadowing effects that can occur due to local topography during a single apparition.

To investigate variations in radar backscatter over Chesterton, Tolkien, Tryggvadóttir, and Kandinsky, we applied a k-means clustering algorithm (Lloyd 1982) to the multiyear, total power maps from Harmon et al. (2011), following the methods of Rivera-Valentin et al. (2022). This algorithm partitions data into clusters by minimizing within-cluster variances and allows for a robust comparison of local-scale variations in radar backscatter. Here we chose to cluster total power into six groups. Unlike in Rivera-Valentin et al. (2022), where k-means was only applied to radar returns identified as pertaining to bright features, clustering was done over the entire scene for each studied crater. In Figure 3, we show radar backscatter class maps with the PSR boundary from Figure 1(b) overlaid. Here backscatter intensity increases from Class 1 to Class 6, with Class 0 representing backscatter that was indistinguishable from noise.

Zoom In Zoom Out Reset image size

For each crater, the majority of high radar backscatter (i.e., Classes ≥4) occurs within the PSR boundary. Small discrepancies where high radar backscatter occurs outside of the PSR boundary, most notably in Figure 3(d), may be due to uncertainties in alignment or to converting the delay-Doppler radar images to hermiocentric coordinates, which can result in a 1.5–4.5 km offset (i.e., 1–3 pixels with respect to the native radar resolution). Of note, though, is that a several-kilometer wide region in Prokofiev has been shown to also host a radar-bright deposit just outside of the PSR boundary, consistent with a thermal environment that can support subsurface water ice (Barker et al. 2022). As radar-bright regions extending just beyond the PSRs exist at other locations on Mercury, discrepancies such as those in Figure 3(d) are not unrealistic.

In regard to low backscatter, we found that regions near central peaks in Tolkien and Kandinsky, where ice is not stable at the surface according to thermal models, typically correspond to radar returns of Class 0–2. Other than these features, low backscatter is not readily identifiable inside of Chesterton, Tolkien, Tryggvadóttir, and Kandinsky. In fact, Class 0–2 radar backscatter only corresponds to ∼6% of the radar return from inside these craters. Additionally, potential small-scale deposits outside of the larger craters are also identifiable in Figure 3, most notably in Figure 3(d), by their high radar backscatter. These features are also associated with the locations of small-scale PSRs.

We cross-correlated maps from Figure 3 with our new high-resolution maximum temperature maps. For every map in Figure 3, we found the average and corresponding standard deviation of the modeled maximum temperature associated with each cluster of radar backscatter. Results are shown in Figure 4. We found that low radar backscatter (Classes ≤2) is associated with average maximum temperatures ${\overline{T}}_{\max }=330$ ± 100 K (with 1σ uncertainty) and overall correlates well with non-PSR-hosting terrains. High backscatter (Classes ≥4) is associated with ${\overline{T}}_{\max }=102$ ± 27 K and correlates with PSR-hosting terrains. The intermediate backscatter cluster (Class 3) has a range of maximum temperatures with an average ${\overline{T}}_{\max }=131$ ± 86 K. Thus, the cross-correlation between radar and temperature maps demonstrates that there is a clear pattern between backscatter and temperature.

Zoom In Zoom Out Reset image size

Ice stability depth predictions were also compared against the maps of k-means clusters of radar backscatter inside of Chesterton, Tolkien, Tryggvadóttir, and Kandinsky. Results are shown in Figure 5. Low radar backscatter (i.e., Classes 1 and 2) is typically associated with ice stability depths near ∼1 m. The radar penetration depth in dry, lunar-like regolith is estimated to be on the order of 10 times the radar wavelength. For this study, accounting for the high radar incidence angle, the sampled depth normal to the surface is ∼0.3 m (Rivera-Valentin et al. 2022). As such, the predicted mantling of ∼1 m over these Class 1 and 2 regions is nearly three times the radar penetration depth, which would result in significant attenuation and therefore lower radar backscatter, as observed. On the other hand, high radar backscatter (i.e., Classes ≥4) is typically associated with areas predicted to have conditions that would allow surface or very shallow (i.e., 0.05 m) subsurface ice, as seen on Figure 5. Mantling much shallower than the radar wavelength (12 cm) results in negligible attenuation. This would be consistent with findings by Rivera-Valentin et al. (2022), who showed that backscatters from radar-bright features within craters that can host surface ice are indistinguishable from those that can only host ice beneath a thin mantling. As shown by Rivera-Valentin et al. (2022), the variation in radar backscatter for Classes 5 and 6 is well modeled by differences in ice purity.

Zoom In Zoom Out Reset image size

Interestingly, Figure 5 shows Class 3 radar backscatter at areas predicted to have conditions that would allow for surficial ice deposits within the four northernmost craters of this study. However, Class 3 radar backscatter values are on average nearly twice as low as Class 6. In fact, regions identified here as Class 3 backscatter were not distinguishable from noise in the work by Rivera-Valentin et al. (2022), which only included data from six consecutive nights in 2019. At the same time, though, Class 3 radar backscatter values are on average three times higher than Class 1. As such, Class 3 likely represents echoes from terrains with significant attenuation and/or lower backscatter enhancement with respect to Classes ≥4 and less attenuation than those associated with Classes ≤2. This may be due to differences in regolith metallic content (e.g., Campbell et al. 1997; Heggy et al. 2020) or another physical property that results in enhanced radar return.

Figure 3 shows that Class 3 radar backscatter primarily occurs from the floors of the four large craters where one would expect minimal backscatter enhancement from topography that decreases the local radar incidence angle. In comparison, Class 4–6 radar backscatter occurs on both the floor and walls of the crater, further suggesting that the backscatter differences are not due to topographic and viewing effects. Indeed, both SC and OC radar reflections from dry regolith have a flat response with respect to high incidence angles, such as those occurring in the study area (Fa et al. 2011; Thompson et al. 2011; Rivera-Valentin et al. 2022). Therefore, topography and surface roughness alone are not likely to produce the observed elevated radar returns. Structural subsurface irregularities and volume scattering, though, may result in higher radar backscatter at high incidence angles.

This leaves two potential explanations for regions with Class 3 radar backscatter: either (1) they are devoid of water ice but have heightened radar return with respect to the background terrain due to subsurface irregularities, or (2) they host some ice, but it is buried beneath a thick mantling, is less abundant, or contains increased amounts of contaminants in order to result in lower radar return with respect to Classes ≥4. Class 3 backscatter typically occurs from terrains surrounding high radar backscatter regions, even from small-scale PSRs outside of the four large craters. Such an association would not necessarily occur under scenario 1 above. As such, we propose that the intermediate values of total power backscatter occur from terrains with some icy component, such as dirty ice, pore filling ice, or buried ice. Such structures, though, are not readily distinguishable by radar (Fa et al. 2011).

There is already evidence that water ice is not uniformly distributed among the available cold traps at either Mercury's north (Deutsch et al. 2016) or south (Chabot et al. 2018b) poles given the apparent absence of strong radar backscatter in many PSRs, though potential biases in the composite radar products can complicate such investigations (Gläser & Oberst 2022). Independent of radar backscatter, surface roughness derived from MLA measurements is also suggestive of variations in ice abundance between and within individual north polar PSRs (Deutsch et al. 2022). The variation in radar backscatter observed within the PSRs of the four large craters in this study (Figure 3) could be further evidence of uneven water ice distribution in Mercury's cold traps that is not correlated with the thermal environment. One theory that has been offered to explain the uneven ice distribution at Mercury's cold traps is delivery of water ice by a large cometary impact (Ernst et al. 2018). Local-scale heterogeneities may be due to variations in ice purity within radar-bright features (as found by Rivera-Valentin et al. 2022), variations in ice burial depths, or variations in the ice abundances throughout each of the polar deposits. Finally, if the intermediate Class 3 radar backscatter is caused by ice diffused through the regolith at depth, these deposits may be the most relevant analog for ice within some lunar PSRs (Mitchell et al. 2018).

MESSENGER's MDIS Wide Angle Camera (WAC) captured images that revealed the permanently shadowed surfaces within some of Mercury's north polar craters using its broadband filter (Chabot et al. 2014, 2016). These detailed images of the PSRs might allow us to identify boundaries on the volatiles and assess the level of agreement with our models (Figure 2). The WAC broadband filter was designed to obtain calibration images of stars, making it very sensitive to light, with a central wavelength of ∼700 nm and a bandpass of ∼600 nm (Denevi et al. 2018). Since the PSRs are only minimally illuminated by light scattered off of nearby sunlit surfaces, the latter were saturated in these images in order to capture and resolve the details of the surface of the PSRs (Figure 6). When pixels of a charged coupled device (CCD; such as the one used for the MDIS WAC detector) saturate, the charge can spill into adjacent pixels in a process known as blooming. Blooming defaced many of the MDIS WAC images that were attempting to capture the surfaces within the PSRs of the northernmost craters, leaving just a few useful images to examine for each crater. An example of blooming can be seen in Figure 6(2a). Importantly, in these images that were successful at revealing the surface within the PSR, the readout edge of the CCD camera was situated along the PSR, which prevented blooming across the crater floor (Chabot et al. 2014). As a result, though, the whole PSRs could not be captured in a single image.

Zoom In Zoom Out Reset image size

The nature of MESSENGER's orbit proved to be simultaneously favorable and detrimental to collecting images to reveal the PSR within Mercury's northernmost craters. MESSENGER's high-eccentricity orbit (Solomon & Anderson 2018) permitted a close-up view of the northern region not obtained for the southern region, which allowed MDIS to acquire high-resolution images of the surface. However, MESSENGER's off-polar inclination diminished the amount of time it spent over the highest-latitude craters, reducing the total number of images ultimately obtained compared to the lower-latitude northern regions. Additionally, near Mercury's pole, the Sun's incidence angle is very near 90°, resulting in a terrain that is heavily shadowed with limited sunlit surfaces to illuminate the PSRs in comparison to lower-latitude locations.

Details of the MDIS WAC broadband images that best revealed the PSR for the four northernmost craters are given in Table 1. Figure 6 shows four of the MDIS images used in this study that best reveal the PSR for each crater in our study, with additional images from Table 1 available in the Appendix. On Figure 6, the yellow lines outline the PSR in each crater, determined by the white region shown in Figure 1(b), while the purple lines outline the region where ice is predicted to be stable at the surface, as described by the white region shown in Figure 1(f).

To aid in examining the MDIS WAC broadband images for evidence of surface volatiles, simulated images were constructed, with the methods from Mazarico et al. (2018) and as used in Hamill et al. (2020). These simulated images were made using the high-resolution DEMs for each crater and the illumination conditions corresponding to the MDIS images in Table 1. Simulated images were produced to predict the radiance (W m⁻² sr⁻¹ μm⁻¹) at each imaging time. In these images, a uniform surface albedo is assumed for each crater, allowing us to see how brightness variations emerge from light scattering in the crater owing solely to illumination geometry and topography. This allows us to compare these models to the actual MDIS WAC images and identify inconsistencies in brightness that can only be explained by deviations in surface reflectance, such as those due to the presence of surface volatiles. Figure 6 shows the simulated images and their corresponding MDIS images.

We employed a quantitative approach in our analysis of the MDIS images by finding the ratio of the MDIS images to their corresponding simulated images. With the MDIS images corrected to digital numbers (DN) as described in Hamill et al. (2020) and the simulated images in units of radiance, we first scaled the pixel values of the simulated image to match the values of the MDIS image to achieve a meaningful ratio. We did this by spatially mapping pixels in the MDIS image to the pixels in the simulated image, then with the simulated pixel values on the x-axis and the MDIS pixel values on the y-axis, we found a linear equation for the line of best fit to describe the pixel scaling. Pixels in the MDIS image that were at or near saturation were preemptively removed from the equation calculation in an effort to exclude compromised data. We did this by first excluding saturated pixels and then excluding pixels that were brighter than the mean pixel value. We applied the equation for the line of best fit to the simulated image and then took the ratio between the original MDIS image and the rescaled simulated image. Figure 6 (column (c)) shows the resulting ratio images from this approach. With these ratio maps, we could evaluate how much brighter or darker regions of the MDIS image are compared to regions of the simulated image. In case of a perfect representation, values >1 in the ratio image could be attributed to the presence of surface ice, blooming, or artifacts of the MDIS WAC, while values <1 would hint at low-reflectance volatiles at the surface.

All of the craters in this study have PSRs that cover nearly the entire crater floor, and most are predicted to have ice at the surface in a region that largely overlaps with the PSR. This leaves only a small region for low-reflectance volatiles at the surface, making boundaries between low- and high-reflectance volatiles more challenging to identify. Kandinsky is a special case whose geographic location and topography may allow for a large region of low-reflectance, organic-rich volatile compounds, represented by coronene, to be adjacent to a nearly equally large region of high-reflectance water ice, as shown in Figure 2. However, in these ratio images, we were unable to find conclusive evidence for the presence of high-reflectance volatiles that border low-reflectance volatiles at the surface. We do not observe a significant contrast in brightness along the purple outline in the MDIS images, nor the ratio in Figure 6 for any of the craters, though most notably at Kandinsky, where we would expect this contrast to be most evident.

It is worth noting that in each of the ratio images the ratio tends to be greater along the crater walls than in the center of the images. This brightening effect is likely due to artifacts of the MDIS camera that have been identified as being due to scattered light within the imaging system (Denevi et al. 2018), rather than the presence of exposed, high-reflectance compounds. This phenomenon is more pronounced in the direction that corresponds to the orientation of the readout edge of the CCD (while the readout edge is on the left, we observe a stronger gradient toward the top of the image). This trend of a brightness gradient that is consistent with an MDIS artifact becomes more apparent during these extremely challenging low-light PSR imaging opportunities. So while this feature can be seen in all four images in Figure 6, in Tryggvadóttir (Figures 6(2a) and 11(1a)), it is uniquely abrupt at about midway along the crater and is clearly visible directly in the MDIS WAC broadband image. This feature does not appear to be a consequence of blooming since there are no saturated pixels to the left that could compromise the image. Multiple MDIS imaging orientations would have helped to evaluate if the brightness variation seen within the MDIS WAC broadband images of Figures 6(2a) and 11(1a) of the PSR of Tryggvadóttir is real or rather an anomaly caused by the MDIS camera. However, due to the limited MESSENGER data for Tryggvadóttir, such data are not available.

Our analysis of the other viable MDIS images we identified for each crater (listed in Table 1) showed similar results, which are presented in the Appendix (see Appendix A.2). Our analysis of the MDIS images does not rule out the presence of volatiles on the surface, but rather just shows that there are no obvious features in the available MDIS images that can solely be attributed to the presence of surface volatiles.

MESSENGER's MLA measured the surface albedo at 1064 nm in the northern hemisphere of Mercury (Neumann et al. 2013; Deutsch et al. 2018). At this wavelength, the average reflectance of Mercury's surface is 0.17 (Chabot et al. 2018a). Exposed surface ice on Mercury is observed to be nearly a factor of 2 brighter, with a surface reflectance r_s ≥ 0.3 (Deutsch et al. 2017), while low-reflectance volatile materials are observed to be nearly a factor of 2 darker (r_s < 0.1; Chabot et al. 2018a).

We used the MLA data (Deutsch et al. 2017) to analyze the brightness variations between the regions where we anticipate ice and low-reflectance volatile compounds, represented by coronene, to be stable at the surface according to the thermal model prediction shown in Figure 2. Barker et al. (2022) found that MLA values that were captured at pulse widths <15 ns, ranges <100 m and >800 m, and emission angles >20° were anomalous or unreliable. We filtered the MLA data according to the pulse width and range constraints identified by Barker et al. (2022). However, due to the lack of data at latitudes >86°N, we chose to remove MLA points at emission angles >35° rather than 20° to retain some MLA data in the northernmost craters of this study for analysis.

All the MLA data between 825 and 90°N latitude with the specified MLA constraints applied are plotted in Figure 7. As is evident from this figure, the data are increasingly sparse with increasing latitude. Additionally, suspect MLA data, especially those due to high emission angles, were abundant at high latitudes where the craters of this study reside, further reducing the already limited data available to reliably examine the craters.

Zoom In Zoom Out Reset image size

We plotted the MLA data within the PSR of each of the northernmost craters in Figure 8, where the pink outline represents the region where we predict ice to be stable at the surface. The overall lack of data in these four northernmost craters is apparent in Figure 8. The number of filtered MLA points in the PSR of each crater is given in Table 2. The sum total of MLA points in the four northernmost craters is 485. By comparison, the total number of points in Figure 7 is on the order of 10⁶.

Zoom In Zoom Out Reset image size

While we do not observe a correlation between pixel values in the background MDIS images and MLA reflectance in Figure 8, the MDIS images are not calibrated to reflectance values and, as described in Section 3.2, the extreme low-light conditions make it difficult to image the surface within the PSRs. In looking at the MLA data in Figure 8, we notice that most of the MLA points within the pink outline are very bright (r_s > 0.3), consistent with ice being present at the surface, in agreement with our prediction (Figure 2). The MLA data within the region where we predict water ice to be unstable within Kandinsky's PSR do visually appear to be lower reflectance on average (despite the apparent uniformity in the MDIS image), though the number of points in this region, and in comparable regions within the other three craters, is inadequate to determine whether low-reflectance compounds are exposed at the surface.

Figure 9 shows the probability density as a function of the MLA reflectance values for various regions. The purple curve shows the distribution of all of the refined MLA data available between 55° and 90°N latitude. This curve peaks sharply just before 0.2. This is consistent with the average reflectance of Mercury, which we find to be 0.19 with the constrained MLA data set. While the purple curve includes PSRs, the PSRs make up <1% of Mercury's surface between 55° and 90°N latitude, so this curve largely represents the reflectance of Mercury's regolith. By comparison, the yellow curve is the distribution of MLA points in the sunlit region at high latitudes (86°–90°N latitude), which also represents Mercury's regolith, since we anticipate the sunlit regions to be too warm to sustain ice or other volatiles long-term. The relatively enhanced surface reflectance (r_s ≥ 0.3) of the adjacent sunlit areas outside the craters, which can be observed in Figure 7, can potentially be ascribed to the presence of small patches of exposed water ice (Deutsch et al. 2017), but the fact that the purple and yellow curves differ substantially, even though both should represent Mercury's regolith, is indicating that care must be taken in interpreting MLA reflectance data acquired at these high latitudes on Mercury owing to the high emission angles prevalent there (Neumann et al. 2013; Deutsch et al. 2017).

Zoom In Zoom Out Reset image size

In Figure 9, the green curve is the distribution of MLA points inside the PSR of the northernmost craters. Though subsets of the same purple curve, the green curve entirely excludes sunlit regions, while the yellow curve entirely excludes PSRs in the same latitude range as the green curve. Therefore, comparing the green and yellow curves enables comparison between the MLA data in the PSR, where we expect volatile compounds to be present at the surface, and the MLA data in the surrounding sunlit region, where we do not expect volatile compounds to be at the surface (with the exception of the potential for small patches), with similarly high emission angles.

We find that the PSR of the craters as a whole was found to be brighter than the sunlit region at similar latitudes by a factor of 1.65; this contrast is consistent with an abundance of high-reflectance material at the surface of the PSR rather than artificial brightening caused by the high emission angles in these regions. This result is consistent with the conclusion reached by Deutsch et al. (2017) that demonstrated evidence for high-reflectance surfaces within these four large northernmost craters.

In Figure 9, we also show the results of separating the limited MLA points into regions where ice is predicted to be stable at the surface (red curve) and regions where coronene is predicted to be stable at the surface (blue curve). Unfortunately, the green, red, and blue curves largely overlap, indicating that these regions are statistically indistinguishable owing to the highly limited number of points in the regions where low-reflectance volatiles represented by coronene are predicted to be stable at the surface.

As seen in Figure 8, the points inside the PSR of each crater are primarily contained inside the region where ice is predicted to be stable at the surface, and the blue curve has large error bars and spread, which is consistent with the overall lack of data in the coronene region. In fact, there are only 76 MLA points available for the regions where low-reflectance volatiles like coronene are predicted to be present on the surface, compared to 409 MLA points for the regions where water ice is predicted to be present on the surface as outlined in Table 2. Thus, available MLA data are insufficient to fully test our prediction of precisely where ice or low-reflectance volatiles might be present at the surface within these craters.

Table 3 provides a listing of the MDIS images used to create the DEMs for each crater in this study.

In addition to the MDIS images shown in Figure 6, we show the other MDIS images that were assessed in this study. Some of the MDIS images presented in this section were more affected by blooming than the ones presenting in Figure 6. While blooming can deteriorate the image quality, we present the images that provide a sufficiently detailed view of the PSR. Only the fraction of images presented in this section, of the more than 50 images that were available for each crater with resolutions better than 200 m pixel⁻¹, were adequate enough for analysis. This is because imaging regions in permanent shadow in an orientation that best prevents blooming is challenging.

Figures 10–13 show the comparison between the MDIS and simulated images for Kandinsky, Tryggvaddóttir, Chesterton, and Tolkein, respectively. The details of the MDIS images can be found in Table 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 10 depicts the comparison of the remaining high-resolution MDIS images that sufficiently detail the PSR of Kandinsky to the corresponding simulated images. In these ratios (Figures 10(1c) and (2c)), we see similar results to Figure 6(1c), in that the brightness appears to increase along the crater walls, likely due to the proximity to the saturated sunlit surface, and that no contrast in brightness that would indicate variations in surface composition is observed, particularly along the purple outline, where we expect to see such variations.

Figure 11 depicts the comparison of the remaining high-resolution MDIS images that sufficiently detail the PSR of Tryggvadóttir to the corresponding simulated images. Much like in Figure 6(2c), Figure 11(1c) depicts an abrupt shift in brightness midway through the crater. While this same shift is not observed in Figure 11(2c), this MDIS image is severely compromised by blooming. Since there are only two images with very similar views of Tryggvadóttir that clearly show this shift in brightness, we cannot make definitive claims about the nature of this feature, and so we look to BepiColombo to investigate further.

Figure 12 depicts the comparison of the remaining high-resolution MDIS images that sufficiently detail the PSR of Chesterton to the corresponding simulated images. In these ratio images (Figures 6(3c) and 12(1–8c)), we observe minor differences. However, we were not able to identify features that could indicate the presence of water ice at the surface in any of them.

Finally, Figure 13 depicts the comparison of the remaining high-resolution MDIS images that sufficiently detail the PSR of Tolkien to the corresponding simulated images. Figures 12(1c) and (2c) show similar results to Figure 6(4c), with little evidence for variations that agree with the prediction (Figure 2).