Depth to Diameter Analysis on Small Simple Craters at the Lunar South Pole—Possible Implications for Ice Harboring

: In this paper, we present a study comparing the depth to diameter (d/D) ratio of small simple craters (200–1000 m) of an area between − 88.5° to − 90° latitude at the lunar south pole containing Permanent Shadowed Regions (PSRs) versus craters without PSRs. As PSRs can reach temperatures of 110 K and are capable of harboring volatiles, especially water ice, we analyzed the relationship of depth versus diameter ratios and its possible implications for harboring water ice. Variations in the d/D ratios can also be caused by other processes such as degradation, isostatic adjustment, or differences in surface properties. The conducted d/D ratio analysis suggests that a differentiation between craters containing PSRs versus craters without PSRs occurs. Thus, a possible direct relation between d/D ratio, PSRs, and water ice harboring might exist. Our results suggest that differences in the target’s surface properties may explain the obtained results. The resulting d/D ratios of craters with PSRs can help to select target areas for future In-Situ Resource Utilization (ISRU) missions.


Introduction
Volatile characterization on the lunar surface is ongoing research that started over 50 years ago. Watson et al. [1] suggested that ice might be trapped in Permanent Shadowed Regions (PSR) for short periods of time. Arnold [2] corroborated the idea of PSRs harboring volatiles, especially water ice, at the lunar south pole. Recent investigations also concur with this hypothesis (Nozette et al. [3], McClanahan et al. [4], and Susorney et al. [5]) however, others (Eke et al. [6], Haruyama et al. [7]) did not find any compelling results.
In polar regions, PSRs can be excellent cold trap candidates (Kokhanov et al. [8] and Rubanenko et al. [9]). Due to the low inclination of the spin axis of the Moon relative to the ecliptic (1.54 • ) and the extremely high topography in its south pole, PSRs can become excellent cold traps where ice might be present (Paige et al. [10]). Due to the lack of illumination, the characterization of water ice in these regions is mainly studied indirectly (e.g., relating PSR with low temperature) and by combining data from different instruments. A study by Mitrofanov et al. [11] registered a decrease of the neutron flux on the Lunar Exploration Neutron Detector (LEND) onboard the Lunar Reconnaissance Orbiter (LRO) in some PSR areas at the lunar south pole. This neutron flux decrease could be related to the volatile presence. In addition, Fisher et al. [12] found a positive correlation between reflectance measurements within PSRs, obtained by using the scattered sunlight from the inner crater walls, that the presence of ice might cause. Mitrofanov et al. [11] and Sanin et al. [13] registered a decrease of the neutron flux LEND, which might be related to possible traces of ice. The study Sanin et al. [14] showed that only three craters on the Moon, namely Shoemaker and Cabeus in the South and Rozhdestvensky in the North, had significant neutron suppression. Other PSRs were not showing relevant values, and some even showed an excess of neutron emission with respect to near-lit areas. Sanin et al. [14] concluded that except for the three craters mentioned above, no conclusive results could be extracted for the rest of PSRs. Furthermore, Fisher et al. [12] used reflectance measurements of the Lunar Orbiter Laser Altimeter (LOLA) tracks to find evidence for surface water ice. They found out that the reflectance increases with decreasing temperatures (around 110 K) which coincides with the presence of water ice. Finally, Hayne et al. [15] found a correlation between PSR crater temperatures and ultraviolet albedo spectra.
The Circular Polarization Ratio (CPR) is a measure of the quality of the reflected circularly-polarized radar signal. It is used to emphasize and ignore specific signals related to the texture of the reflecting material. Several studies on radar observation suggest that the CPR of the signal is directly related to the morphology of the radiating surface (Spudis et al. [16], Campbell et al. [17], and Thomson et al. [18]). Harmon et al. [19] identified that a high CPR (CPR > 1) in several locations of the south pole of Mercury coincided with areas where water ice might be present. This correlation was also found on Europa, Ganymede, and Callisto (Ostro et al. [20]). The problem of CPR is that a high CPR signal could also be produced by fresh regolith exposed to the surface, terrain ruggedness, and boulder distribution coinciding with the CPR scale (in the case of Mini-RF, 10 cm) (Baloga et al. [21]). Using CPR alone to characterize ice introduces significant uncertainties and needs to be supported by other data to eliminate possible false positives. On the Moon, several attempts to relate high CPR signatures with ice have been pursued, although no conclusive results have been achieved (Eke et al. [6] and Nozette et al. [3]). CPR was also found to have high signals in illuminated and shadowed areas, thus concluding that CPR depends on the morphology and nature of the terrain. Eke et al. [6] analyzed various polar craters and concluded that CPR alone is not sufficient to indicate the presence of ice. M-chi decompositions of radar sensors (Raney et al. [22]) is a recent technique used in surface analysis. By combining the Stokes products (Stokes and G. [23]) in an RGB composition, the backscattering behavior of the radar signal is shown in different colors. These colors reflect different characteristics of the surface (Raney et al. [24]). This RGB combination is helpful to identify the nature of high CPR signals inside craters, enhancing young fresh craters where ejecta and new material will return a high CPR value, thus allowing to discard the possibility of having ice.
The fact that PSRs are never directly illuminated makes it impossible to use an imaging spectrometer such as the Mineralogical Moon Mapper (M3) onboard Chandrayaan-1 (Pieters et al. [25]). However Li et al. [26] used scattered sunlight from crater walls within PSRs and analyzed spectra on those regions that showed consistency with the presence of ice. McCord et al. [27] detected absorption features outside PSRs from 2.8 to 3.0 µm, attributed to water. The bands near 2.8 µm in M3 are particularly sensitive to temperature, and images located at the lunar south pole need a specific temperature effect removal function (Combe et al. [28], Clark et al. [29]). Current temperature removal functions have been reported to be untrustworthy in images taken at the lunar south pole and might introduce uncertainty in the areas where water has absorption peaks in the spectra (Dr. Christian Wöhler, personal communication, May 2017). Furthermore, using the 3.0 µm band on M3 is inadvisable as it is too close to the detector boundary.
A crater's topography with steep walls and the low elevation of the Sun in south polar regions make craters ideal candidates to contain PSRs (Mazarico et al. [30]). Zuber et al. [31] extensively studied the Shackleton crater (where more than 70% of the area is a PSR) in order to find possible water ice by analyzing Lunar Orbiter Laser Altimeter (LOLA) observations. The study shows a possible water ice occurrence on the crater's floor, as such areas are brighter (at LOLA's wavelength of 1064 nm) than their surroundings. However, results are unclear as the brightness might also be caused by space weathering. Most of the studies conducted in this region focused on large craters such as Shackleton (Haruyama et al. [7], Thomson et al. [18], and Nozette et al. [3]) and a smaller number of studies concentrated on small simple craters (Kereszturi and Steinmann [32], Daubar et al. [33], and Stopar et al. [34]). The depth to diameter (d/D) ratio can be used to characterize crater age and evolutionary stage (Lagain et al. [35], Craddock et al. [36], and Robbins et al. [37]). Previous studies suggested that the d/D ratio for fresh primary craters is expected to be around 0.2 and decreases over time due to degradation (Elachi et al. [38] and Schultz et al. [39]). Recent studies also show d/D values from small simple craters to be slightly lower for fresh craters (∼0.1) than older ones. According to Stopar et al. [40], fresh simple craters are deeper than degraded craters regardless of their size. Changes in d/D are difficult to distinguish for minor and moderate crater degradations. Daubar et al. [33] suggest that small lunar craters are likely to be unrecognized secondary craters or some of those craters are primary craters that have undertaken degradation to similar shapes. In addition, small craters are less likely to show variation in surface heat flow after formation, making it easier to retain the possible volatiles carried by an impacting asteroid, as the crater's surface temperature might not reach temperatures capable of changing the volatile state of matter (e.g., water ice's sublimation temperature). If this hypothesis is true, small craters containing PSRs create an excellent environment for harboring water. Finally, according to Boyce et al. [41] slight variations in depths of small simple craters might be interpreted as volatile activity inside the crater. Stopar et al. [40], interprets the d/D ratio variations to be caused by other processes such as differences in target properties (e.g., layering, porosity, angle of repose, strength, fragmentation, cohesion, shear strength, and stratigraphy), resurfacing processes, or degradation.
This study focuses on the analysis of the depth to diameter characteristics of small simple craters that might contain water ice in a Region of Interest (RoI) at the lunar south pole between −88.5°to −90°latitude ( Figure 1). This area ( Figure 2) is a candidate landing site for future lunar exploration missions (De Rosa et al. [42], Kokhanov et al. [43], and Lemelin et al. [44]). The primary objective of our investigation is to analyze the d/D ratios small simple craters that fulfill good volatile harboring conditions on the Moon.

Data and Methods
The Lunar Reconnaissance Orbiter Camera (LROC) has a Wide Angle Camera (WAC) offering a resolution of 75 to 380 m and the Narrow Angle Camera (NAC) provides images with a full resolution of 0.5 m, both at an altitude of 50 km (Robinson et al. [45]). Previous studies (Boyce et al. [41] and Stopar et al. [40]) of d/D analysis focusing on small simple craters used the higher resolution available data. Several NAC and WAC images of the lunar south pole exist. However, due to the inclination of the rotation axis of the Moon and its extreme topography, most of the images only show large shadows. The presence of large shadows makes it almost impossible to create Digital Terrain Models (DTM) as stereophotogrammetry is not possible.
The Lunar Orbiter Laser Altimeter (LOLA) is a laser that creates a 65-m wide swath containing five laser dots with 10-m spacing along tracks (Robinson et al. [45]). DTMs using LOLA data are available on PDS. However, they contain artifacts (Gläser et al. [46] and Barker et al. [47]) that will influence the results of our study, especially the illumination.
In our study, we use a DTM provided by Gläser et al. [48] generated from Lunar Orbiter Laser Altimeter (LOLA) data ( Figure 3). The DTM was co-registered using optical data from LROC NAC, has a resolution of 20 m per pixel, and covers the area from −88.5°t o −90°S. We use the DTM to identify PSRs based on illumination conditions and digitize impact craters based on topography information.

Illumination
DTMs are a common data source to compute illumination on the lunar south pole (Mazarico et al. [30], Bussey et al. [49], Noda et al. [50], and Gläser et al. [48]). Illumination conditions based on elevation data can be computed through ray tracing or the horizons method. Both methods rely on computing the elevation of the Sun at a given position and comparing it with the elevation of the highest elevated terrain along the line from the observer to the Sun. If the Sun has a higher elevation, the pixel where the observer is located is illuminated. In the other case, the observer will be in shadow (see Gläser et al. [46]).
The ray-tracing method (Noda et al. [50] and Bussey et al. [49]) compares the elevation of the Sun versus the elevation of the highest point in the terrain at a given time along the direction of the Sun from a location on the terrain. This method is convenient when short periods of time are analyzed. The computation is pursued for each pixel at every position of the Sun. Therefore, the number of calculations increases for a long period of time. Hence, this method is only recommended for periods ranging from a few hours to one year if an hourly time step is used.
The horizons method (Mazarico et al. [30] and Gläser et al. [46]) introduces a preprocessing step where the highest elevated points around all pixels on the terrain are computed. Although very time-consuming, this initial step allows computing the illumination quickly when analyzing multiple azimuths and more extended periods. The result is about 720 files corresponding to the horizon around all the pixels every 0.5°. Choosing this step size of 0.5°e nsures that at least one azimuth direction lies within the Sun's angular diameter of 0.53°. In our study, we used the horizons method, as described in Gläser et al. [46], to compute the illumination as it was more convenient for our analysis. We chose to create one file per azimuth as this will reduce the input/output operations, which are highly time-consuming while computing the illumination.
We computed illumination over a 20-year period, which covers the precession cycle (18.6 years) of the Moon, as the orientation of its rotational axis is not fixed. This allows us to average the amount of light received over certain areas that receive high illumination during certain times and are less illuminated during other periods. We also modeled the Sun as a disk and not a simple dot, allowing us to assign a fraction of the visible disk, thus getting closer to an accurate representation. The illumination is computed by comparing the Sun's and Horizon's elevations over a certain period and assigning an illumination value to each pixel.
To minimize the computation time, Marco Figuera [51] and Marco Figuera et al. [52] created an illumination computation method adapted to work on GPUs using parallel computation. The process reduced the computation time significantly, allowing to analyze illumination for larger areas and along a more extended period of time. The resulting PSRs are shown in Figure 4.

Crater Selection
There are two methods to measure rim-to-rim diameter and rim-to-floor depth (Chappelow [53] and Stopar et al. [40]). One method uses DTMs, and the other the cast of shadows over craters. The shadow casting method is beneficial when craters are not nested, have clear and distinguishable rims, and the shadows lay within the crater (Chappelow [53] and Daubar et al. [33]). The second method for crater selection uses DTMs (Stopar et al. [40]). Using specific software, such as JMars (Christensen et al. [54]) and CraterTools (Kneissl et al. [55]), we can compute the rim-to-rim diameter and rim-to-floor depth. CraterTools is a toolbox running on ArcGIS and JMars is a java-based GIS client allowing the user to digitize craters in a straightforward way by selecting a fixed circle and hovering it over samesized craters. The results can be stored as a shapefile for later visualization in JMars or other GIS clients. The results can also be exported for further analysis in CraterStats (Michael and Neukum [56]). One of the main drawbacks of JMars is that the diameters of the circles have to be previously selected by hand and fitted to the craters. This makes the results to be approximated and not as precise as with CraterTools. In conclusion, the overall speed and accuracy drove the choice of CraterTools over JMars. As mentioned in Section 2, NAC images of the lunar south pole are mostly covered by shadows. Therefore, the DTM method using CraterTools was preferred over the shadow-casting method. To enhance the crater rims in the DTM to better digitize them, we computed a series of hillshaded reliefs locating the Sun at different angles and constant elevation. This was necessary as some of the craters at the lightly sloped areas faces are otherwise challenging to identify. The hillshades were computed using the GDAL library with a constant elevation of the Sun of 45°above the horizon and 45, 135, 225, and 315°azimuth, respectively. A hillshade at every 90°azimuth was considered sufficient to appreciate differences between the scenes. Once the hillshades were computed, we overlapped them with the DTM in ArcGIS and gave them a 50% transparency to help distinguish the craters.
To digitize the craters, we divided the region of interest, located in the south-westernmost part of the DTM, into smaller squared areas. The process starts by selecting one area, loading one hillshade, and digitizing all the visible craters, and repeating the process for the rest of the areas. In case a crater was not clearly visible, all the hillshades were loaded at the same time. If the crater edges were still not visible, the crater was not digitized.
In our study, we used CraterTools by drawing three points on the crater rim or by selecting two points along the crater's diameter. While the second method seems faster due to fewer points needed, better results are achieved when using the first. CraterTools stores the selected craters in a shapefile containing the diameter, longitude, and latitude of each crater.

Depth and Diameter Measurement
To derive the depth of the craters, we measured the highest elevated point on a crater's rim and the lowest elevated point on a crater's floor, using the raster calculator of ArcGIS. Craters on very steep areas (>25 • ) were discarded as the lowest point could not reflect the actual depth of the crater. Subtracting these two measurements, we computed the range of each crater and later included them into the craters table.
As mentioned in Robbins et al. [57], diameters of craters must be at least 10 pixels across as artifacts may occur within this distance. Taking the DTM's resolution of 20 m per pixel, the minimum diameter considered was 200 m. We digitized a total of 3317 craters and discarded 2564 after applying filters (slope <25º and diameter between 200 and 1000 m) leaving us with a total of 753 craters. Using ArcGIS's intersection tool, we selected the craters containing PSR regions and those without PSRs (Figure 4).

Uncertainties
The method used to select craters has a random error derived from measured crater diameters and depths. In order to evaluate the error, we conducted an experiment consisting of repetitively measuring diameters and depths. We divided the craters into groups of 100-m diameter steps, such that we had eight groups ranging from 200 m to 1000 m. We selected two craters randomly from each group and measured the diameter and depth three times.
The results and statistics of the experiment are shown in Table 1. The maximum standard deviation of diameters is 0.76 m for a crater diameter of 523.41 m with an associated error of 0.15%. However, the depth standard deviation of the crater is 0.01% for a depth of 112.21 m. For the depths, the maximum standard deviation is 0.08 m for a depth of 198.01 m, and a percentage of the standard deviation of 0.04%. The analysis shows good agreement with the values that were digitized.

Results
The results from the d/D measurement of craters containing PSRs vs. craters not containing PSRs are shown in Figure 5. Figure 5A shows the histogram grouped in d/D ratios of craters containing PSRs (blue) and not containing PSRs (orange). The results show that the average d/D ratio is smaller for craters with PSRs (d/D∼0.16) than without PSRs (d/D∼0.22). Figure 5B shows the slope versus the d/D ratio and a strong correlation between the wall slope and d/D ratio where the deepest craters have the steepest slopes. Craters without PSRs follow the expected trend of increasing the slope with an increasing d/D ratio. Craters with PSRs cluster at slopes smaller than 20°, denoting a possible increased shallowness of the crater profile. Figure 5C shows the depth versus diameter on a logarithmic scale. Although craters from both categories overlap, especially at a small diameter and depth, two distinguished classes can be identified. In general, we can see that craters in both groups tend to decrease their d/D ratio at around 400 m in diameter, whereas larger craters cluster around 0.2. This observation is consistent with the studies from Stopar et al. [40] and Fassett et al. [58].
According to Figure 5D, craters containing PSRs concentrate in the low slope region with an increasing diameter, whereas craters with no PSRs are primarily in the 200-300-m diameter range but with greater average slopes.
In general, our analysis shows that craters containing PSRs tend to be shallower than those without PSRs. Furthermore, craters with PSRs have lower slopes but higher variability than craters with no PSRs, which tend to follow a more linear trend ( Figure 5B).
We conducted a two-sided two-sample hypothesis testing to determine whether the difference between the two groups of craters is statistically significant. Several statistical tests can be performed for two-sample testing, for example, the Kolmogorov-Smirnov test, Student's t-test, or Mann-Whitney U test. Our study applies the Mann-Whitney U test since it is less sensitive to the different numbers of craters in the two populations (MacFarland et al. [59]). The null hypothesis (H 0 ) is that the d/D ratios in the two groups of craters are equal, whereas our alternate hypothesis (H a ) states a difference in d/D ratios between the two groups. We obtained a p-value smaller than 0.01, which is statistically significant (p < 0.05). We also ran the Mann-Whitney U test based on diameter ranges to see if the differences are also significant. We found craters with diameters smaller than 399 m and from 500 m to 599 m to be significant. However the remaining ranges have a p-value above 0.01. This is most likely caused by a very unbalanced population of craters containing PSRs vs. no PSRs in the range 400 m to 499 m. In the case of crater with a range of 700 m to 899 m, it might be due to a lack of no PSRs craters. Therefore, we can conclude that the d/D ratios of the craters, especially smaller than 400 m in diameter, containing PSRs are significantly different from the d/D ratios of craters without PSRs.

Discussion
We analyzed the d/D ratios of small simple craters at the lunar south pole. Our results show a close correlation between slopes and d/D ratios. In general, craters containing PSRs have average d/D ratios of about 0.16, whereas those without PSRs have average d/D ratios spanning from 0.20 to 0.23. Similarly, craters with and without PSRs form two distinguishable groups based on the slope gradients. Craters with PSRs are usually sparse and cluster at slopes lower than 20°, whereas craters without PSRs follow a linear tendency when increasing the slope. The same two groupings can be seen on the depth vs. diameter logarithmic plots. Craters without PSRs cluster around smaller diameters and are closer to d/D ratios of 0.2. On the contrary, craters with PSRs have a smaller d/D ratio of 0.1. Furthermore, the statistical test confirms the differences in d/D ratios between the two groups of craters. Our results support those obtained by Stopar et al. [40], Fassett et al. [58], and Kokhanov et al. [8], who used a different source DTM and included craters smaller than 200 m in diameter.
The difference in d/D ratios between craters containing PSRs and craters containing no PSRs can be caused by various processes related to diffusive crater degradation and target properties. For example, if PSRs contain a subsurface layer of ice, the surface dynamics should be different since the surface material would be harder than the regolith. Therefore, a subsurface layer of ice could cause a lower d/D ratio, as observed of craters containing PSRs. Furthermore, the flatter average slopes on craters containing PSRs coincide with McClanahan et al. [4], who observed that craters with an average wall slope of less than 20°are more likely to bound ice on the surface. Accordingly, if a layer of ice were present in these regions, we would expect the d/D ratio and the crater topography between the populations to be as distinguishable as we observed. However, different states of diffusive degradation due to crater age could show a similar effect (Fassett et al. [58]). If geologic units on craters containing PSRs would be significantly older, we would also expect lower d/D ratios and shallower crater slopes. However, since PSRs in the RoI are not located on unambiguously distinguishable geologic units, we can not study the effect of different degradation states on our interpretations using the given method. This would require further investigation and a detailed geologic analysis of the area.
As some craters containing PSRs might have temperatures not rising above the water ice sublimation temperature (110K), and they have, on average, a slope below the angle of repose, if any water ice exists, those craters fulfilling these constraints have the best morphological conditions to harbor water ice. However, according to the derived temperature maps from Paige et al. [10], the majority of the craters within the studied RoI are located in an area with temperatures, on average, that are above the water ice sublimation temperature. As the spatial resolution of Diviner ranges from 200 m to 1.3 km, small craters with PSRs might not be included in the map from Paige et al. [10]. Using the derived Diviner temperatures as a proxy in the crater diameter ranges used in this study is not sufficient enough to discard the craters containing PSRs as not suitable for harboring water ice. Therefore, the majority of the analyzed craters can be candidates for harboring water ice in its interior.
We also considered a possible bias of the results caused by secondary craters. However, it is challenging to unambiguously distinguish between primary and secondary craters in this region because secondaries do not appear in apparent clusters or rays. As we excluded an area encompassing two large impact structures from the analysis, a randomness analysis would most likely indicate a non-random crater distribution, even if no contamination was present. Since we assume that possible contamination by secondaries would affect craters containing PSRs and craters not containing PSRs to a similar degree and we observe d/D ratios between both populations are statistically distinguishable, we conclude that secondary contamination only marginally affects the results.
The analyzed area has been proposed as a candidate landing area for several missions, both from agencies and the private sector. The possibility of accessing a water ice source could lead, for example, to permanent lunar bases or refueling options. Moreover, as the craters containing PSRs are shallower with a slope that is generally below the maximum operational slope of a rover (<25 • ), this makes the analyzed craters perfect candidates for an In-Situ Resource Utilization (ISRU) mission.

•
We analyzed the d/D ratio of 753 craters with a diameter between 200 and 1000 m and a slope smaller than 25º at the lunar south pole and groupped the craters in two groups: Craters containing PSRs and without PSRs; • d/D ratios of craters containing PSRs vs. without PSRs, cluster in two distinct groups, which implies that the two groups might have some geomorphological differences; • The results from the Mann-Whitney U statistical test support the hypothesis that two distinct groups exist; • Possible secondary craters contamination will affect both crater population, therefore the contamination will only marginally affect the global results; • The difference in d/D ratios can be caused by a subsurface ice layer present within the PSRs. However, we cannot discard the influence of diffusive crater degradation or different geological units.